From Salvatore:
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Jul 2013 09:00:04 +0000 (09:00 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Jul 2013 09:00:04 +0000 (09:00 +0000)
- AliAnalysisTaskSAQA: new QA histograms for centrality, transition to THnSparse
- AliAnalysisTaskSAJF: transition to THnSparse
- AliAnalysisTaskEmcalJet: GetSortedArray returns number of accepted items
- AliJetModelBaseTask: flow modulation in random phi
- AliJetResponseMaker: transition to THnSparse
- AliJetEmbeddingFromAODTask: select events using given track pt range, attach embedding file name to current ESD/AOD event list
- AliHadCorrTask: fix matching rate histogram for tracks (only tracks in fiducial emcal acceptance go in it)
- AliEmcalJetTask: check if MC energy fraction of clusters is set before using it
- AliNamedString: new class to have a TString object findable in a TList collection
- AliEmcalPicoTrackMaker: increase class version to silent errors from the streamer

21 files changed:
PWG/CMakelibPWGTools.pkg
PWG/EMCAL/AliAnalysisTaskEmcal.cxx
PWG/EMCAL/AliEmcalPicoTrackMaker.h
PWG/PWGToolsLinkDef.h
PWG/Tools/AliNamedString.cxx [new file with mode: 0644]
PWG/Tools/AliNamedString.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJet.cxx
PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJet.h
PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
PWGJE/EMCALJetTasks/AliHadCorrTask.cxx
PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx
PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h
PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
PWGJE/EMCALJetTasks/AliJetModelBaseTask.h
PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
PWGJE/EMCALJetTasks/AliJetResponseMaker.h
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.h
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.h
PWGJE/EMCALJetTasks/macros/AddTaskSAQA.C

index a60107c..4112f69 100644 (file)
@@ -35,6 +35,7 @@ set ( SRCS
     Tools/AliFigure.cxx
     Tools/AliHelperPID.cxx
     Tools/AliNamedArrayI.cxx
+    Tools/AliNamedString.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
index 06d84cd..d797628 100644 (file)
@@ -390,12 +390,16 @@ Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track) const
     }
   }
 
-  if (track->Pt() < fTrackPtCut)
+  if (track->Pt() < fTrackPtCut) {
+    AliDebug(2,Form("Track not accepted because of minimum pt cut (%f < %f).", track->Pt(), fTrackPtCut));
     return kFALSE;
+  }
 
   if (track->Eta() < fTrackMinEta || track->Eta() > fTrackMaxEta || 
-      track->Phi() < fTrackMinPhi || track->Phi() > fTrackMaxPhi)
+      track->Phi() < fTrackMinPhi || track->Phi() > fTrackMaxPhi) {
+    AliDebug(2,"Track not accepted because of acceptance cut.");
     return kFALSE;
+  }
   
   return kTRUE;
 }
@@ -843,6 +847,7 @@ Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
       }
     } else {
       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
+      fCent = 99;
       fCentBin = 3;
     }
     AliEventplane *aliEP = InputEvent()->GetEventplane();
index 11271d3..3dab820 100644 (file)
@@ -55,6 +55,6 @@ class AliEmcalPicoTrackMaker : public AliAnalysisTaskSE {
   AliEmcalPicoTrackMaker(const AliEmcalPicoTrackMaker&);            // not implemented
   AliEmcalPicoTrackMaker &operator=(const AliEmcalPicoTrackMaker&); // not implemented
 
-  ClassDef(AliEmcalPicoTrackMaker, 5); // Task to make PicoTracks in AOD/ESD events
+  ClassDef(AliEmcalPicoTrackMaker, 6); // Task to make PicoTracks in AOD/ESD events
 };
 #endif
index 8bece86..e2868e6 100644 (file)
@@ -9,6 +9,7 @@
 #pragma link C++ class AliHelperPID+;
 #pragma link C++ class AliLatexTable+;
 #pragma link C++ class AliNamedArrayI+;
+#pragma link C++ class AliNamedString+;
 #pragma link C++ class AliPWGFunc+;
 #pragma link C++ class AliPWGHistoTools+;
 #pragma link C++ class AliTHn+;
diff --git a/PWG/Tools/AliNamedString.cxx b/PWG/Tools/AliNamedString.cxx
new file mode 100644 (file)
index 0000000..f20498a
--- /dev/null
@@ -0,0 +1,27 @@
+// $Id$
+//
+// Named string object.
+//
+// Author: S.Aiola
+
+#include "AliNamedString.h"
+
+ClassImp(AliNamedString)
+
+//________________________________________________________________________
+AliNamedString::AliNamedString() : 
+  TObjString(),
+  fName()
+{
+  // Dummy constructor.
+
+}
+
+//________________________________________________________________________
+AliNamedString::AliNamedString(const char *name, const char *string) :
+  TObjString(string),
+  fName(name)
+{
+  // Standard constructor.
+
+}
diff --git a/PWG/Tools/AliNamedString.h b/PWG/Tools/AliNamedString.h
new file mode 100644 (file)
index 0000000..fb70313
--- /dev/null
@@ -0,0 +1,27 @@
+#ifndef ALINAMEDSTRING_H
+#define ALINAMEDSTRING_H
+
+// $Id$
+
+#include <TObjString.h>
+
+class AliNamedString : public TObjString {
+ public: 
+  AliNamedString();
+  AliNamedString(const char *name, const char *string="");
+
+  void Clear(Option_t* /*option*/="")           { SetString("")      ; }
+
+  const char* GetName()               const { return fName       ; }
+  void        SetName(const char* n)        { fName = n          ; }
+
+ protected:
+  TString fName; // name of the string object
+  
+ private:
+  AliNamedString(const AliNamedString&);             // not implemented
+  AliNamedString& operator=(const AliNamedString&);  // not implemented
+  
+  ClassDef(AliNamedString, 1); // Named string object
+};
+#endif
index 8c672fc..db2b53d 100644 (file)
@@ -237,7 +237,7 @@ void AliAnalysisTaskEmcalJet::ExecOnce()
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskEmcalJet::GetSortedArray(Int_t indexes[], TClonesArray *array, Double_t rho) const
+Int_t AliAnalysisTaskEmcalJet::GetSortedArray(Int_t indexes[], TClonesArray *array, Double_t rho) const
 {
   // Get the leading jets.
 
@@ -248,6 +248,8 @@ Bool_t AliAnalysisTaskEmcalJet::GetSortedArray(Int_t indexes[], TClonesArray *ar
 
   const Int_t n = array->GetEntriesFast();
 
+  Int_t nacc = 0;
+
   if (n < 1)
     return kFALSE;
 
@@ -267,7 +269,8 @@ Bool_t AliAnalysisTaskEmcalJet::GetSortedArray(Int_t indexes[], TClonesArray *ar
       if (!AcceptJet(jet))
        continue;
       
-      pt[i] = jet->Pt() - rho * jet->Area();
+      pt[nacc] = jet->Pt() - rho * jet->Area();
+      nacc++;
     }
   }
 
@@ -287,7 +290,8 @@ Bool_t AliAnalysisTaskEmcalJet::GetSortedArray(Int_t indexes[], TClonesArray *ar
       if (!AcceptTrack(track))
        continue;
       
-      pt[i] = track->Pt();
+      pt[nacc] = track->Pt();
+      nacc++;
     }
   }
 
@@ -310,16 +314,15 @@ Bool_t AliAnalysisTaskEmcalJet::GetSortedArray(Int_t indexes[], TClonesArray *ar
       TLorentzVector nPart;
       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
       
-      pt[i] = nPart.Pt();
+      pt[nacc] = nPart.Pt();
+      nacc++;
     }
   }
 
-  TMath::Sort(n, pt, indexes);
-
-  if (pt[indexes[0]] == -FLT_MAX) 
-    return 0;
+  if (nacc > 0)
+    TMath::Sort(nacc, pt, indexes);
 
-  return kTRUE;
+  return nacc;
 }
 
 //________________________________________________________________________
index f150eb6..56007d3 100644 (file)
@@ -43,7 +43,7 @@ class AliAnalysisTaskEmcalJet : public AliAnalysisTaskEmcal {
   Double_t                    GetLeadingHadronPt(AliEmcalJet* jet)                                     const;
   void                        ExecOnce()                                                                    ;
   AliRhoParameter            *GetRhoFromEvent(const char *name)                                             ;
-  Bool_t                      GetSortedArray(Int_t indexes[], TClonesArray *array, Double_t rho=0)     const;
+  Int_t                       GetSortedArray(Int_t indexes[], TClonesArray *array, Double_t rho=0)     const;
   Bool_t                      IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted = kFALSE)       const;
   Bool_t                      IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted = kFALSE)      const;
   Bool_t                      RetrieveEventObjects()                                                        ;
index 4ff5eef..364cabd 100644 (file)
@@ -455,7 +455,7 @@ void AliEmcalJetTask::FindJets()
           maxNe = cPt;
 
         if (c->GetLabel() > fMinMCLabel) // MC particle
-          mcpt += cPt * c->GetMCEnergyFraction();
+          mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
 
         if (cPhi<0) 
           cPhi += TMath::TwoPi();
index 598de8d..ccaeda4 100644 (file)
@@ -533,10 +533,10 @@ void AliHadCorrTask::DoTrackLoop()
     if (!AcceptTrack(track))
       continue;
     
-    if (track->GetTrackEtaOnEMCal() < fGeom->GetArm1EtaMin() + fEtaMatch ||
-       track->GetTrackEtaOnEMCal() > fGeom->GetArm1EtaMax() - fEtaMatch || 
-       track->GetTrackPhiOnEMCal() < fGeom->GetArm1PhiMin() * TMath::DegToRad() + fPhiMatch || 
-       track->GetTrackPhiOnEMCal() > fGeom->GetArm1PhiMax() * TMath::DegToRad() - fPhiMatch)
+    if (track->GetTrackEtaOnEMCal() < fGeom->GetArm1EtaMin() - fEtaMatch ||
+       track->GetTrackEtaOnEMCal() > fGeom->GetArm1EtaMax() + fEtaMatch || 
+       track->GetTrackPhiOnEMCal() < fGeom->GetArm1PhiMin() * TMath::DegToRad() - fPhiMatch || 
+       track->GetTrackPhiOnEMCal() > fGeom->GetArm1PhiMax() * TMath::DegToRad() + fPhiMatch)
       continue;
     
     Int_t Nclus = emctrack->GetNumberOfMatchedObj();
index 273d095..8fa8a62 100644 (file)
@@ -39,6 +39,7 @@
 #include "AliFJWrapper.h"
 #include "AliLog.h"
 #include "AliInputEventHandler.h"
+#include "AliNamedString.h"
 
 ClassImp(AliJetEmbeddingFromAODTask)
 
@@ -69,6 +70,9 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
   fJetType(0),
   fJetAlgo(1),
   fJetParticleLevel(kTRUE),
+  fParticleMinPt(0),
+  fParticleMaxPt(0),
+  fParticleSelection(0),
   fIncludeNoITS(kFALSE),
   fCutMaxFractionSharedTPCClusters(0.4),
   fUseNegativeLabels(kTRUE),
@@ -89,6 +93,7 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
   fAODCaloCells(0),
   fAODMCParticles(0),
   fCurrentAODEntry(-1),
+  fAODFilePath(0),
   fHistFileMatching(0),
   fHistAODFileError(0),
   fHistNotEmbedded(0),
@@ -134,6 +139,9 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t
   fJetType(0),
   fJetAlgo(1),
   fJetParticleLevel(kTRUE),
+  fParticleMinPt(0),
+  fParticleMaxPt(0),
+  fParticleSelection(0),
   fIncludeNoITS(kFALSE),
   fCutMaxFractionSharedTPCClusters(0.4),
   fUseNegativeLabels(kTRUE),
@@ -154,6 +162,7 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t
   fAODCaloCells(0),  
   fAODMCParticles(0),
   fCurrentAODEntry(-1),
+  fAODFilePath(0),
   fHistFileMatching(0),
   fHistAODFileError(0),
   fHistNotEmbedded(0),
@@ -250,6 +259,13 @@ Bool_t AliJetEmbeddingFromAODTask::ExecOnce()
   else
     fEsdTreeMode = kTRUE;
 
+  fAODFilePath = static_cast<AliNamedString*>(InputEvent()->FindListObject("AODEmbeddingFile"));
+  if (!fAODFilePath) {
+    fAODFilePath = new AliNamedString("AODEmbeddingFile", "");
+    AliDebug(3,"Adding AOD embedding file path object to the event list...");
+    InputEvent()->AddObject(fAODFilePath);
+  }
+
   return AliJetModelBaseTask::ExecOnce();
 }
 
@@ -394,7 +410,8 @@ Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
 
   } while (!IsAODEventSelected());
 
-  fHistRejectedEvents->Fill(attempts);
+  if (fHistRejectedEvents)
+    fHistRejectedEvents->Fill(attempts);
 
   if (!fCurrentAODTree)
     return kFALSE;
@@ -411,7 +428,7 @@ Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
     AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
 
     // Trigger selection
-    if (fTriggerMask != AliVEvent::kAny) {
+    if (fTriggerMask != 0) {
       UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
       
       if ((offlineTrigger & fTriggerMask) == 0) {
@@ -452,6 +469,12 @@ Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
       
   }
 
+  // Particle selection
+  if ((fParticleSelection == 1 && FindParticleInRange(fAODTracks)==kFALSE) ||
+      (fParticleSelection == 2 && FindParticleInRange(fAODClusters)==kFALSE) ||
+      (fParticleSelection == 3 && FindParticleInRange(fAODMCParticles)==kFALSE))
+    return kFALSE;
+
   // Jet selection
   if (fJetMinPt > 0) {
     TLorentzVector jet;
@@ -482,6 +505,39 @@ Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
 }
 
 //________________________________________________________________________
+Bool_t AliJetEmbeddingFromAODTask::FindParticleInRange(TClonesArray *array)
+{
+  if (!array) return kFALSE;
+
+  if (array->GetClass()->InheritsFrom("AliVParticle")) {
+    const Int_t nentries = array->GetEntriesFast();
+    for (Int_t i = 0; i < nentries; i++) {
+      AliVParticle *part = static_cast<AliVParticle*>(array->At(i));
+      if (!part) continue;
+      if (part->Pt() > fParticleMinPt && part->Pt() < fParticleMaxPt) return kTRUE;
+    }
+  }
+  else if (array->GetClass()->InheritsFrom("AliVCluster")) {
+    const Int_t nentries = array->GetEntriesFast();
+    for (Int_t i = 0; i < nentries; i++) {
+      AliVCluster *clus = static_cast<AliVCluster*>(array->At(i));
+      if (!clus) continue;
+      
+      TLorentzVector vect;
+      Double_t vert[3] = {0,0,0};
+      clus->GetMomentum(vect,vert);
+
+      if (vect.Pt() > fParticleMinPt && vect.Pt() < fParticleMaxPt) return kTRUE;
+    }
+  }
+  else {
+    AliWarning(Form("Unable to do event selection based on particle pT: %s class type not recognized.",array->GetClass()->GetName()));
+  }
+
+  return kFALSE;
+}
+
+//________________________________________________________________________
 void AliJetEmbeddingFromAODTask::Run() 
 {
   if (!GetNextEntry()) {
@@ -496,6 +552,8 @@ void AliJetEmbeddingFromAODTask::Run()
   if (fHistEmbeddingQA)
     fHistEmbeddingQA->Fill("OK", 1);
 
+  fAODFilePath->SetString(fCurrentAODFile->GetName());
+
   if (fOutMCParticles) {
 
     if (fCopyArray && fMCParticles)
index 977b2e5..a37c798 100644 (file)
@@ -12,6 +12,7 @@ class AliVHeader;
 class TH2;
 class TH1;
 class TLorentzVector;
+class AliNamedString;
 
 #include "AliJetModelBaseTask.h"
 
@@ -50,6 +51,7 @@ class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
   void           SetJetAlgo(Byte_t t)                              { fJetAlgo            = t     ; }
   void           SetZVertexCut(Double_t z)                         { fZVertexCut         = z     ; }
   void           SetMaxVertexDist(Double_t d)                      { fMaxVertexDist      = d     ; }
+  void           SetParticlePtRange(Double_t min, Double_t max, Byte_t t=1) { fParticleMinPt = min; fParticleMaxPt = max; fParticleSelection = t; }
 
  protected:
   Bool_t          ExecOnce()            ;// intialize task
@@ -59,6 +61,7 @@ class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
   virtual Bool_t  GetNextEntry()        ;// get next entry in current tree
   virtual Bool_t  IsAODEventSelected()  ;// AOD event trigger/centrality selection
   TLorentzVector  GetLeadingJet(TClonesArray *tracks, TClonesArray *clusters=0);  // get the leading jet
+  Bool_t          FindParticleInRange(TClonesArray *array);// Find particle in array within range (fParticleMinPt, fParticleMaxPt)
 
   TObjArray     *fFileList            ;//  List of AOD files 
   Bool_t         fRandomAccess        ;//  Random access to file number and event
@@ -84,6 +87,9 @@ class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
   Byte_t         fJetType             ;//  Jet type (0=full, 1=charged, 2=neutral)
   Byte_t         fJetAlgo             ;//  Jet algorithm (0=kT, 1=anti-kT)
   Bool_t         fJetParticleLevel    ;//  Trigger, look at particle level jets
+  Double_t       fParticleMinPt       ;//  Select events with a particle pt between fParticleMinPt and fParticleMaxPt (see fParticleSelection)
+  Double_t       fParticleMaxPt       ;//  Select events with a particle pt between fParticleMinPt and fParticleMaxPt (see fParticleSelection)
+  Byte_t         fParticleSelection   ;//  Particles used to select events (def=0=none, 1=tracks, 2=clusters, 3=MC particles)
   Int_t          fAODfilterBits[2]    ;//  AOD track filter bit map
   Bool_t         fIncludeNoITS        ;//  True = includes tracks with failed ITS refit
   Double_t       fCutMaxFractionSharedTPCClusters;  // max fraction of shared TPC clusters
@@ -105,6 +111,7 @@ class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
   AliVCaloCells *fAODCaloCells        ;//! AOD cell collection
   TClonesArray  *fAODMCParticles      ;//! AOD MC particles collection
   Int_t          fCurrentAODEntry     ;//! Current entry in the AOD tree
+  AliNamedString *fAODFilePath        ;//! Current AOD file path being embedded
   TH2           *fHistFileMatching    ;//! Current file ID vs. AOD file ID (to be embedded)
   TH1           *fHistAODFileError    ;//! AOD file ID (to be embedded) error
   TH1           *fHistNotEmbedded     ;//! File ID not embedded
@@ -115,6 +122,6 @@ class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
   AliJetEmbeddingFromAODTask(const AliJetEmbeddingFromAODTask&);            // not implemented
   AliJetEmbeddingFromAODTask &operator=(const AliJetEmbeddingFromAODTask&); // not implemented
 
-  ClassDef(AliJetEmbeddingFromAODTask, 10) // Jet embedding from AOD task
+  ClassDef(AliJetEmbeddingFromAODTask, 11) // Jet embedding from AOD task
 };
 #endif
index 04a9048..936e919 100644 (file)
@@ -11,6 +11,7 @@
 #include <TLorentzVector.h>
 #include <TRandom3.h>
 #include <TList.h>
+#include <TF2.h>
 
 #include "AliVEvent.h"
 #include "AliAODCaloCluster.h"
@@ -55,8 +56,10 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fNTracks(0),
   fMarkMC(99999),
   fPtSpectrum(0),
+  fPtPhiEvPlDistribution(0),
   fDensitySpectrum(0),
   fQAhistos(kFALSE),
+  fPsi(0),
   fIsInit(0),
   fGeom(0),
   fClusters(0),
@@ -107,8 +110,10 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fNTracks(1),
   fMarkMC(99999),
   fPtSpectrum(0),
+  fPtPhiEvPlDistribution(0),
   fDensitySpectrum(0),
   fQAhistos(drawqa),
+  fPsi(0),
   fIsInit(0),
   fGeom(0),
   fClusters(0),
@@ -186,9 +191,9 @@ void AliJetModelBaseTask::UserExec(Option_t *)
   }
 
   if (fDensitySpectrum) {
-    fNTracks = fDensitySpectrum->GetRandom();
-    fNCells = fDensitySpectrum->GetRandom();
-    fNClusters = fDensitySpectrum->GetRandom();
+    fNTracks = TMath::Nint(fDensitySpectrum->GetRandom());
+    fNCells = TMath::Nint(fDensitySpectrum->GetRandom());
+    fNClusters = TMath::Nint(fDensitySpectrum->GetRandom());
   }
   
   // Clear map
@@ -205,6 +210,9 @@ void AliJetModelBaseTask::UserExec(Option_t *)
     }
   }
 
+  if (fPtPhiEvPlDistribution)
+    fPsi = gRandom->Rndm() * TMath::Pi();
+
   Run();
 
   if (fCaloCells && !fCopyArray) {
@@ -662,20 +670,25 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t labe
 AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Double_t ptemc, Bool_t ise, Int_t label, Short_t charge)
 {
   // Add a track to the event.
-
-  const Int_t nTracks = fOutTracks->GetEntriesFast();
   
-  if (pt < 0) 
-    pt = GetRandomPt();
-  if (eta < -100) 
-    eta = GetRandomEta();
-  if (phi < 0) 
-    phi = GetRandomPhi();
+  if (pt < 0 && eta < -100 && phi < 0) {
+    GetRandomParticle(pt,eta,phi);
+  }
+  else {
+    if (pt < 0) 
+      pt = GetRandomPt();
+    if (eta < -100) 
+      eta = GetRandomEta();
+    if (phi < 0) 
+      phi = GetRandomPhi(pt);
+  }
 
   if (label >= 0)
     label += fMarkMC+fMCLabelShift;
   else if (label < 0)
     label -= fMarkMC+fMCLabelShift;
+  
+  const Int_t nTracks = fOutTracks->GetEntriesFast();
 
   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, 
                                                                  eta, 
@@ -893,7 +906,9 @@ Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal)
     if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
   }
 
-  return gRandom->Rndm() * (phimax - phimin) + phimin;
+  Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
+
+  return result;
 }
 
 //________________________________________________________________________
@@ -908,6 +923,46 @@ Double_t AliJetModelBaseTask::GetRandomPt()
 }
 
 //________________________________________________________________________
+void AliJetModelBaseTask::GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal)
+{
+  // Get a random particle.
+  
+  eta = GetRandomEta(emcal);
+
+  if (fPtPhiEvPlDistribution) {
+    Double_t phimax = fPhiMax;
+    Double_t phimin = fPhiMin;
+    
+    if (emcal) {
+      const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
+      const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
+      
+      if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
+      if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
+      if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
+      if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
+    }
+    
+    if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) {
+      AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName()));
+      pt = GetRandomPt();
+      phi = GetRandomPhi(emcal);
+    }
+    else {
+      do {
+       fPtPhiEvPlDistribution->GetRandom2(pt,phi);
+       phi += fPsi;
+       if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2;
+      } while (phi > phimax || phi < phimin);
+    }
+  }
+  else {
+    pt = GetRandomPt();
+    phi = GetRandomPhi(emcal);
+  }
+}
+
+//________________________________________________________________________
 void AliJetModelBaseTask::Run() 
 {
   // Run.
index c47a292..d92e5f7 100644 (file)
@@ -4,13 +4,13 @@
 // $Id$
 
 class TClonesArray;
-class TH1I;
 class AliEMCALGeometry;
 class AliVCluster;
 class AliPicoTrack;
 class AliVCaloCells;
 class AliAODMCParticle;
 class AliNamedArrayI;
+class TF2;
 
 #include <TH1F.h>
 #include <TF1.h>
@@ -28,10 +28,11 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
 
   void                   SetEtaRange(Float_t min, Float_t max) { fEtaMin       = min;  fEtaMax = max; }
   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                   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",1000,f->GetXmin(),f->GetXmax()); 
                                                                  fPtSpectrum->Add(f); }
+  void                   SetPtPhiEvPlDistribution(TF2 *f)      { fPtPhiEvPlDistribution   = f;    }
   void                   SetDensitySpectrum(TH1 *f)            { fDensitySpectrum = f;    }
   void                   SetDensitySpectrum(TF1 *f)            { fDensitySpectrum = new TH1F("densitypectrum","densitypectrum",1000,f->GetXmin(),f->GetXmax()); 
                                                                  fDensitySpectrum->Add(f); }
@@ -71,6 +72,7 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   Double_t               GetRandomEta(Bool_t emcal=kFALSE);                                     // generate a random eta value in the given range
   Double_t               GetRandomPhi(Bool_t emcal=kFALSE);                                     // generate a random phi value in the given range
   Double_t               GetRandomPt();                                                         // generate a random pt value in the given range
+  void                   GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal=kFALSE);  // generate a particle with random eta,phi,pt values
   virtual Bool_t         ExecOnce();                                                            // intialize task
   virtual void           Run();                                                                 // do jet model action
 
@@ -97,8 +99,10 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   Int_t                  fNTracks;                // how many tracks are being processed
   Int_t                  fMarkMC;                 // which MC label is to be used (default=100)
   TH1                   *fPtSpectrum;             // pt spectrum to extract random pt values
+  TF2                   *fPtPhiEvPlDistribution;  // pt vs. (phi-psi) distribution to extract random pt/phi values
   TH1                   *fDensitySpectrum;        // particle density spectrum to extract random density values
   Bool_t                 fQAhistos;               // draw QA histograms
+  Double_t               fPsi;                    //!simmetry plane for the elliptic flow
   Bool_t                 fIsInit;                 //!=true if initialized
   AliEMCALGeometry      *fGeom;                   //!pointer to EMCal geometry
   Double_t               fVertex[3];              //!event vertex
@@ -121,6 +125,6 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   AliJetModelBaseTask(const AliJetModelBaseTask&);            // not implemented
   AliJetModelBaseTask &operator=(const AliJetModelBaseTask&); // not implemented
 
-  ClassDef(AliJetModelBaseTask, 9) // Jet modelling task
+  ClassDef(AliJetModelBaseTask, 10) // Jet modelling task
 };
 #endif
index 24987bd..7df5658 100644 (file)
@@ -8,7 +8,7 @@
 
 #include <TClonesArray.h>
 #include <TH2F.h>
-#include <TProfile.h>
+#include <THnSparse.h>
 #include <TLorentzVector.h>
 
 #include "AliAnalysisManager.h"
@@ -45,39 +45,41 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fMatchingPar2(0),
   fUseCellsToMatch(kFALSE),
   fMinJetMCPt(1),
+  fHistoType(0),
+  fDeltaPtAxis(0),
+  fDeltaEtaDeltaPhiAxis(0),
+  fDoJet2Histogram(0),
   fTracks2(0),
   fCaloClusters2(0),
   fJets2(0),
   fRho2(0),
   fRho2Val(0),
   fTracks2Map(0),
-  fMCEnergy1vsEnergy2(0),
+  fHistLeadingJets1PtArea(0),
+  fHistLeadingJets1CorrPtArea(0),
+  fHistLeadingJets2PtArea(0),
+  fHistLeadingJets2CorrPtArea(0),
+  fHistLeadingJets2PtAreaAcceptance(0),
+  fHistLeadingJets2CorrPtAreaAcceptance(0),
+  fHistJets1(0),
+  fHistJets2(0),
+  fHistJets2Acceptance(0),
+  fHistMatching(0),
   fHistJets1PhiEta(0),
   fHistJets1PtArea(0),
   fHistJets1CorrPtArea(0),
-  fHistLeadingJets1PtArea(0),
-  fHistLeadingJets1CorrPtArea(0),
   fHistJets1NEFvsPt(0),
   fHistJets1CEFvsCEFPt(0),
   fHistJets1ZvsPt(0),
   fHistJets2PhiEta(0),
   fHistJets2PtArea(0),
   fHistJets2CorrPtArea(0),
-  fHistLeadingJets2PtArea(0),
-  fHistLeadingJets2CorrPtArea(0),
   fHistJets2PhiEtaAcceptance(0),
   fHistJets2PtAreaAcceptance(0),
   fHistJets2CorrPtAreaAcceptance(0),
-  fHistLeadingJets2PtAreaAcceptance(0),
-  fHistLeadingJets2CorrPtAreaAcceptance(0),
   fHistJets2NEFvsPt(0),
   fHistJets2CEFvsCEFPt(0),
   fHistJets2ZvsPt(0),
-  fHistNonMatchedJets1PtArea(0),
-  fHistNonMatchedJets2PtArea(0),
-  fHistNonMatchedJets1CorrPtArea(0),
-  fHistNonMatchedJets2CorrPtArea(0),
-  fHistMissedJets2PtArea(0),
   fHistCommonEnergy1vsJet1Pt(0),
   fHistCommonEnergy2vsJet2Pt(0),
   fHistDistancevsJet1Pt(0),
@@ -85,9 +87,9 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fHistDistancevsCommonEnergy1(0),
   fHistDistancevsCommonEnergy2(0),
   fHistCommonEnergy1vsCommonEnergy2(0),
+  fHistDeltaEtaDeltaPhi(0),
   fHistJet2PtOverJet1PtvsJet2Pt(0),
   fHistJet1PtOverJet2PtvsJet1Pt(0),
-  fHistDeltaEtaPhi(0),
   fHistDeltaPtvsJet1Pt(0),
   fHistDeltaPtvsJet2Pt(0),
   fHistDeltaPtOverJet1PtvsJet1Pt(0),
@@ -99,10 +101,10 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fHistDeltaPtvsArea2(0),
   fHistDeltaPtvsDeltaArea(0),
   fHistJet1PtvsJet2Pt(0),
-  fHistDeltaCorrPtvsJet1CorrPt(0),
-  fHistDeltaCorrPtvsJet2CorrPt(0),
   fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
   fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
+  fHistDeltaCorrPtvsJet1CorrPt(0),
+  fHistDeltaCorrPtvsJet2CorrPt(0),
   fHistDeltaCorrPtvsDistance(0),
   fHistDeltaCorrPtvsCommonEnergy1(0),
   fHistDeltaCorrPtvsCommonEnergy2(0),
@@ -110,10 +112,10 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fHistDeltaCorrPtvsArea2(0),
   fHistDeltaCorrPtvsDeltaArea(0),
   fHistJet1CorrPtvsJet2CorrPt(0),
-  fHistDeltaMCPtvsJet1MCPt(0),
-  fHistDeltaMCPtvsJet2Pt(0),
   fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
   fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
+  fHistDeltaMCPtvsJet1MCPt(0),
+  fHistDeltaMCPtvsJet2Pt(0),
   fHistDeltaMCPtvsDistance(0),
   fHistDeltaMCPtvsCommonEnergy1(0),
   fHistDeltaMCPtvsCommonEnergy2(0),
@@ -151,39 +153,41 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fMatchingPar2(0),
   fUseCellsToMatch(kFALSE),
   fMinJetMCPt(1),
+  fHistoType(0),
+  fDeltaPtAxis(0),
+  fDeltaEtaDeltaPhiAxis(0),
+  fDoJet2Histogram(0),
   fTracks2(0),
   fCaloClusters2(0),
   fJets2(0),
   fRho2(0),
   fRho2Val(0),
   fTracks2Map(0),
-  fMCEnergy1vsEnergy2(0),
+  fHistLeadingJets1PtArea(0),
+  fHistLeadingJets1CorrPtArea(0),
+  fHistLeadingJets2PtArea(0),
+  fHistLeadingJets2CorrPtArea(0),
+  fHistLeadingJets2PtAreaAcceptance(0),
+  fHistLeadingJets2CorrPtAreaAcceptance(0),
+  fHistJets1(0),
+  fHistJets2(0),
+  fHistJets2Acceptance(0),
+  fHistMatching(0),
   fHistJets1PhiEta(0),
   fHistJets1PtArea(0),
   fHistJets1CorrPtArea(0),
-  fHistLeadingJets1PtArea(0),
-  fHistLeadingJets1CorrPtArea(0),
   fHistJets1NEFvsPt(0),
   fHistJets1CEFvsCEFPt(0),
   fHistJets1ZvsPt(0),
   fHistJets2PhiEta(0),
   fHistJets2PtArea(0),
   fHistJets2CorrPtArea(0),
-  fHistLeadingJets2PtArea(0),
-  fHistLeadingJets2CorrPtArea(0),
   fHistJets2PhiEtaAcceptance(0),
   fHistJets2PtAreaAcceptance(0),
   fHistJets2CorrPtAreaAcceptance(0),
-  fHistLeadingJets2PtAreaAcceptance(0),
-  fHistLeadingJets2CorrPtAreaAcceptance(0),
   fHistJets2NEFvsPt(0),
   fHistJets2CEFvsCEFPt(0),
   fHistJets2ZvsPt(0),
-  fHistNonMatchedJets1PtArea(0),
-  fHistNonMatchedJets2PtArea(0),
-  fHistNonMatchedJets1CorrPtArea(0),
-  fHistNonMatchedJets2CorrPtArea(0),
-  fHistMissedJets2PtArea(0),
   fHistCommonEnergy1vsJet1Pt(0),
   fHistCommonEnergy2vsJet2Pt(0),
   fHistDistancevsJet1Pt(0),
@@ -191,9 +195,9 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fHistDistancevsCommonEnergy1(0),
   fHistDistancevsCommonEnergy2(0),
   fHistCommonEnergy1vsCommonEnergy2(0),
+  fHistDeltaEtaDeltaPhi(0),
   fHistJet2PtOverJet1PtvsJet2Pt(0),
   fHistJet1PtOverJet2PtvsJet1Pt(0),
-  fHistDeltaEtaPhi(0),
   fHistDeltaPtvsJet1Pt(0),
   fHistDeltaPtvsJet2Pt(0),
   fHistDeltaPtOverJet1PtvsJet1Pt(0),
@@ -205,10 +209,10 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fHistDeltaPtvsArea2(0),
   fHistDeltaPtvsDeltaArea(0),
   fHistJet1PtvsJet2Pt(0),
-  fHistDeltaCorrPtvsJet1CorrPt(0),
-  fHistDeltaCorrPtvsJet2CorrPt(0),
   fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
   fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
+  fHistDeltaCorrPtvsJet1CorrPt(0),
+  fHistDeltaCorrPtvsJet2CorrPt(0),
   fHistDeltaCorrPtvsDistance(0),
   fHistDeltaCorrPtvsCommonEnergy1(0),
   fHistDeltaCorrPtvsCommonEnergy2(0),
@@ -216,10 +220,10 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fHistDeltaCorrPtvsArea2(0),
   fHistDeltaCorrPtvsDeltaArea(0),
   fHistJet1CorrPtvsJet2CorrPt(0),
-  fHistDeltaMCPtvsJet1MCPt(0),
-  fHistDeltaMCPtvsJet2Pt(0),
   fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
   fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
+  fHistDeltaMCPtvsJet1MCPt(0),
+  fHistDeltaMCPtvsJet2Pt(0),
   fHistDeltaMCPtvsDistance(0),
   fHistDeltaMCPtvsCommonEnergy1(0),
   fHistDeltaMCPtvsCommonEnergy2(0),
@@ -239,51 +243,26 @@ AliJetResponseMaker::~AliJetResponseMaker()
   // Destructor
 }
 
+
 //________________________________________________________________________
-void AliJetResponseMaker::UserCreateOutputObjects()
+void AliJetResponseMaker::AllocateTH2()
 {
-  // Create user objects.
-
-  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
-
-  // General histograms
-  if (fIsEmbedded) {
-    fMCEnergy1vsEnergy2 = new TH2F("fMCEnergy1vsEnergy2", "fMCEnergy1vsEnergy2", fNbins, fMinBinPt, fMaxBinPt*4, fNbins, fMinBinPt, fMaxBinPt*4);
-    fMCEnergy1vsEnergy2->GetXaxis()->SetTitle("#Sigmap_{T,1}^{MC}");
-    fMCEnergy1vsEnergy2->GetYaxis()->SetTitle("#Sigmap_{T,2}");
-    fMCEnergy1vsEnergy2->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fMCEnergy1vsEnergy2);
-  }
-
-  // Jets 1 histograms
-  fHistLeadingJets1PtArea = new TH2F("fHistLeadingJets1PtArea", "fHistLeadingJets1PtArea", 
-                                    50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistLeadingJets1PtArea->GetXaxis()->SetTitle("area");
-  fHistLeadingJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
-  fHistLeadingJets1PtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistLeadingJets1PtArea);
+  // Allocate TH2 histograms.
 
   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
   fOutput->Add(fHistJets1PhiEta);
   
-  fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1, fNbins, fMinBinPt, fMaxBinPt);
   fHistJets1PtArea->GetXaxis()->SetTitle("area");
   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
   fOutput->Add(fHistJets1PtArea);
-
   if (!fRhoName.IsNull()) {
-    fHistLeadingJets1CorrPtArea = new TH2F("fHistLeadingJets1CorrPtArea", "fHistLeadingJets1CorrPtArea", 
-                                          50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
-    fHistLeadingJets1CorrPtArea->GetXaxis()->SetTitle("area");
-    fHistLeadingJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
-    fHistLeadingJets1CorrPtArea->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistLeadingJets1CorrPtArea);
-
     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 
-                                   50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                   fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
@@ -310,19 +289,25 @@ void AliJetResponseMaker::UserCreateOutputObjects()
 
   // Jets 2 histograms
 
-  fHistLeadingJets2PtAreaAcceptance = new TH2F("fHistLeadingJets2PtAreaAcceptance", "fHistLeadingJets2PtAreaAcceptance", 
-                                              50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistLeadingJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
-  fHistLeadingJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
-  fHistLeadingJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistLeadingJets2PtAreaAcceptance);
+  fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
+  fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
+  fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistJets2PhiEta);
 
-  fHistLeadingJets2PtArea = new TH2F("fHistLeadingJets2PtArea", "fHistLeadingJets2PtArea", 
-                                    50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistLeadingJets2PtArea->GetXaxis()->SetTitle("area");
-  fHistLeadingJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
-  fHistLeadingJets2PtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistLeadingJets2PtArea);
+  fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets2PtArea->GetXaxis()->SetTitle("area");
+  fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
+  fHistJets2PtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets2PtArea);
+
+  if (!fRho2Name.IsNull()) {
+    fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 
+                                   fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
+    fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
+    fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistJets2CorrPtArea);
+  }
 
   fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
   fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
@@ -330,51 +315,19 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fOutput->Add(fHistJets2PhiEtaAcceptance);
 
   fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 
-                                       50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
+                                       fNbins/2, 0, 1, fNbins, fMinBinPt, fMaxBinPt);
   fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
   fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
   fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
   fOutput->Add(fHistJets2PtAreaAcceptance);
 
-  fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
-  fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
-  fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistJets2PhiEta);
-
-  fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistJets2PtArea->GetXaxis()->SetTitle("area");
-  fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
-  fHistJets2PtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistJets2PtArea);
-
   if (!fRho2Name.IsNull()) {
     fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 
-                                             50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                             fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
     fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
     fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
     fOutput->Add(fHistJets2CorrPtAreaAcceptance);
-
-    fHistLeadingJets2CorrPtAreaAcceptance = new TH2F("fHistLeadingJets2CorrPtAreaAcceptance", "fHistLeadingJets2CorrPtAreaAcceptance", 
-                                                    50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
-    fHistLeadingJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
-    fHistLeadingJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
-    fHistLeadingJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistLeadingJets2CorrPtAreaAcceptance);
-
-    fHistLeadingJets2CorrPtArea = new TH2F("fHistLeadingJets2CorrPtArea", "fHistLeadingJets2CorrPtArea", 
-                                          50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
-    fHistLeadingJets2CorrPtArea->GetXaxis()->SetTitle("area");
-    fHistLeadingJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
-    fHistLeadingJets2CorrPtArea->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistLeadingJets2CorrPtArea);
-
-    fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 
-                                   50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
-    fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
-    fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
-    fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistJets2CorrPtArea);
   }
 
   fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
@@ -397,45 +350,6 @@ void AliJetResponseMaker::UserCreateOutputObjects()
 
   // Matching histograms
 
-  fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 
-                                       50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
-  fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
-  fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistNonMatchedJets1PtArea);
-
-  fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 
-                                       50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
-  fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
-  fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistNonMatchedJets2PtArea);
-
-  if (!fRhoName.IsNull()) {  
-    fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 
-                                             50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
-    fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
-    fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
-    fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistNonMatchedJets1CorrPtArea);
-  }
-
-  if (!fRho2Name.IsNull()) {  
-    fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 
-                                             50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
-    fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
-    fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
-    fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistNonMatchedJets2CorrPtArea);
-  }
-
-  fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 
-                                   50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");  
-  fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
-  fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistMissedJets2PtArea);
-
   fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
   fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
   fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");  
@@ -478,6 +392,12 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
   fOutput->Add(fHistCommonEnergy1vsCommonEnergy2);
 
+  fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2);
+  fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)");
+  fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)");  
+  fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaEtaDeltaPhi);
+
   fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
   fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
   fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
@@ -490,12 +410,6 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
   fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
 
-  fHistDeltaEtaPhi = new TH2F("fHistDeltaEtaPhi", "fHistDeltaEtaPhi", 200, -0.995, 1.005, 200, -TMath::Pi()*99/200, TMath::Pi()*301/200);
-  fHistDeltaEtaPhi->GetXaxis()->SetTitle("#delta#eta");
-  fHistDeltaEtaPhi->GetYaxis()->SetTitle("#delta#phi");
-  fHistDeltaEtaPhi->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistDeltaEtaPhi);
-
   fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", 
                                  fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
   fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");  
@@ -546,21 +460,21 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fOutput->Add(fHistDeltaPtvsCommonEnergy2);
 
   fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1", 
-                                50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
   fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
   fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
   fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
   fOutput->Add(fHistDeltaPtvsArea1);
 
   fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2", 
-                                50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
   fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
   fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
   fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
   fOutput->Add(fHistDeltaPtvsArea2);
 
   fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea", 
-                                    100, -1.98, 2.02, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                    fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
   fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
   fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
   fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
@@ -633,21 +547,21 @@ void AliJetResponseMaker::UserCreateOutputObjects()
     fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2);
 
     fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1", 
-                                      50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                      fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
     fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
     fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
     fOutput->Add(fHistDeltaCorrPtvsArea1);
     
     fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2", 
-                                      50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                      fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
     fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
     fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
     fOutput->Add(fHistDeltaCorrPtvsArea2);
     
     fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea", 
-                                          100, -1.98, 2.02, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                          fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
     fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
     fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
@@ -671,7 +585,7 @@ void AliJetResponseMaker::UserCreateOutputObjects()
 
   if (fIsEmbedded) {
     fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt", 
-                                     fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                       fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");  
     fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
@@ -720,21 +634,21 @@ void AliJetResponseMaker::UserCreateOutputObjects()
     fOutput->Add(fHistDeltaMCPtvsCommonEnergy2);
 
     fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1", 
-                                    50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                    fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
     fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
     fOutput->Add(fHistDeltaMCPtvsArea1);
 
     fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2", 
-                                    50, 0, 2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                    fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
     fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
     fOutput->Add(fHistDeltaMCPtvsArea2);
 
     fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea", 
-                                        100, -1.98, 2.02, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+                                        fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
     fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
     fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
@@ -747,11 +661,459 @@ void AliJetResponseMaker::UserCreateOutputObjects()
     fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
     fOutput->Add(fHistJet1MCPtvsJet2Pt);
   }
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::AllocateTHnSparse()
+{
+  // Allocate THnSparse histograms.
+
+  TString title[20]= {""};
+  Int_t nbins[20]  = {0};
+  Double_t min[20] = {0.};
+  Double_t max[20] = {0.};
+  Int_t dim = 0;
+
+  title[dim] = "#phi";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
+  dim++;
+
+  title[dim] = "#eta";
+  nbins[dim] = fNbins/4;
+  min[dim] = -1;
+  max[dim] = 1;
+  dim++;
+
+  title[dim] = "p_{T}";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "A_{jet}";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1;
+  dim++;
+
+  title[dim] = "NEF";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "Z";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  Int_t dim1 = dim, dim2 = dim;
+
+  if (!fRhoName.IsNull()) {
+    title[dim1] = "p_{T}^{corr}";
+    nbins[dim1] = fNbins*2;
+    min[dim1] = -250;
+    max[dim1] = 250;
+    dim1++;
+  }
+
+  if (fIsEmbedded) {
+    title[dim1] = "p_{T}^{MC}";
+    nbins[dim1] = fNbins;
+    min[dim1] = 0;
+    max[dim1] = 250;
+    dim1++;
+  }
+
+  fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
+  for (Int_t i = 0; i < dim1; i++) 
+    fHistJets1->GetAxis(i)->SetTitle(title[i]);
+  fOutput->Add(fHistJets1);
+
+  if (!fRho2Name.IsNull()) {
+    title[dim2] = "p_{T}^{corr}";
+    nbins[dim2] = fNbins*2;
+    min[dim2] = -250;
+    max[dim2] = 250;
+    dim2++;
+  }
+
+  if (fDoJet2Histogram) {
+    fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
+    for (Int_t i = 0; i < dim2; i++) 
+      fHistJets2->GetAxis(i)->SetTitle(title[i]);
+    fOutput->Add(fHistJets2);
+  }
+
+  fHistJets2Acceptance = new THnSparseD("fHistJets2Acceptance","fHistJets2Acceptance",dim2,nbins,min,max);
+  for (Int_t i = 0; i < dim2; i++) 
+    fHistJets2Acceptance->GetAxis(i)->SetTitle(title[i]);
+  fOutput->Add(fHistJets2Acceptance);
+  
+  // Matching
+
+  dim = 0;
+
+  title[dim] = "p_{T,1}";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "p_{T,2}";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "A_{jet,1}";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1;
+  dim++;
+
+  title[dim] = "A_{jet,2}";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1;
+  dim++;
+
+  title[dim] = "distance";
+  nbins[dim] = fNbins/2;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "CE1";
+  nbins[dim] = fNbins/2;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "CE2";
+  nbins[dim] = fNbins/2;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  if (fDeltaPtAxis) {
+    title[dim] = "#deltaA_{jet}";
+    nbins[dim] = fNbins/2;
+    min[dim] = -1;
+    max[dim] = 1;
+    dim++;
+
+    title[dim] = "#deltap_{T}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (!fRhoName.IsNull()) {
+    title[dim] = "p_{T,1}^{corr}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (!fRho2Name.IsNull()) {
+    title[dim] = "p_{T,2}^{corr}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (fDeltaPtAxis && (!fRhoName.IsNull() || !fRho2Name.IsNull())) {
+    title[dim] = "#deltap_{T}^{corr}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (fDeltaEtaDeltaPhiAxis) {
+    title[dim] = "#delta#eta";
+    nbins[dim] = fNbins/2;
+    min[dim] = -1;
+    max[dim] = 1;
+    dim++;
+
+    title[dim] = "#delta#phi";
+    nbins[dim] = fNbins/2;
+    min[dim] = -TMath::Pi()/2;
+    max[dim] = TMath::Pi()*3/2;
+    dim++;
+  }
+  if (fIsEmbedded) {
+    title[dim] = "p_{T,1}^{MC}";
+    nbins[dim] = fNbins;
+    min[dim] = 0;
+    max[dim] = 250;
+    dim++;
+
+    if (fDeltaPtAxis) {
+      title[dim] = "#deltap_{T}^{MC}";
+      nbins[dim] = fNbins*2;
+      min[dim] = -250;
+      max[dim] = 250;
+      dim++;
+    }
+  }
+      
+  fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
+    
+  for (Int_t i = 0; i < dim; i++)
+    fHistMatching->GetAxis(i)->SetTitle(title[i]);
+   
+  fOutput->Add(fHistMatching);
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::UserCreateOutputObjects()
+{
+  // Create user objects.
+
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+  // Jets 1 histograms
+  fHistLeadingJets1PtArea = new TH2F("fHistLeadingJets1PtArea", "fHistLeadingJets1PtArea", 
+                                    fNbins/2, 0, 1, fNbins, fMinBinPt, fMaxBinPt);
+  fHistLeadingJets1PtArea->GetXaxis()->SetTitle("area");
+  fHistLeadingJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
+  fHistLeadingJets1PtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistLeadingJets1PtArea);
+
+  if (!fRhoName.IsNull()) {
+    fHistLeadingJets1CorrPtArea = new TH2F("fHistLeadingJets1CorrPtArea", "fHistLeadingJets1CorrPtArea", 
+                                          fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistLeadingJets1CorrPtArea->GetXaxis()->SetTitle("area");
+    fHistLeadingJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
+    fHistLeadingJets1CorrPtArea->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistLeadingJets1CorrPtArea);
+  }
+
+  fHistLeadingJets2PtAreaAcceptance = new TH2F("fHistLeadingJets2PtAreaAcceptance", "fHistLeadingJets2PtAreaAcceptance", 
+                                              fNbins/2, 1, 2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistLeadingJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
+  fHistLeadingJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
+  fHistLeadingJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistLeadingJets2PtAreaAcceptance);
+
+  fHistLeadingJets2PtArea = new TH2F("fHistLeadingJets2PtArea", "fHistLeadingJets2PtArea", 
+                                    fNbins/2, 0, 1, fNbins, fMinBinPt, fMaxBinPt);
+  fHistLeadingJets2PtArea->GetXaxis()->SetTitle("area");
+  fHistLeadingJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
+  fHistLeadingJets2PtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistLeadingJets2PtArea);
+
+  if (!fRho2Name.IsNull()) {
+    fHistLeadingJets2CorrPtAreaAcceptance = new TH2F("fHistLeadingJets2CorrPtAreaAcceptance", "fHistLeadingJets2CorrPtAreaAcceptance", 
+                                                    fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistLeadingJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
+    fHistLeadingJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
+    fHistLeadingJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistLeadingJets2CorrPtAreaAcceptance);
+
+    fHistLeadingJets2CorrPtArea = new TH2F("fHistLeadingJets2CorrPtArea", "fHistLeadingJets2CorrPtArea", 
+                                          fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistLeadingJets2CorrPtArea->GetXaxis()->SetTitle("area");
+    fHistLeadingJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
+    fHistLeadingJets2CorrPtArea->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistLeadingJets2CorrPtArea);
+  }
+
+  if (fHistoType==0)
+    AllocateTH2();
+  else 
+    AllocateTHnSparse();
 
   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 }
 
 //________________________________________________________________________
+void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t Z, Double_t CorrPt, Double_t MCPt, Int_t Set)
+{
+  if (fHistoType==1) {
+    THnSparse *histo = 0;
+    if (Set==1)
+      histo = fHistJets1;
+    else if (Set==2)
+      histo = fHistJets2;
+    else if (Set==3)
+      histo = fHistJets2Acceptance;
+
+    if (!histo)
+      return;
+
+    Double_t contents[20]={0};
+
+    for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
+      TString title(histo->GetAxis(i)->GetTitle());
+      if (title=="#phi")
+       contents[i] = Phi;
+      else if (title=="#eta")
+       contents[i] = Eta;
+      else if (title=="p_{T}")
+       contents[i] = Pt;
+      else if (title=="A_{jet}")
+       contents[i] = A;
+      else if (title=="NEF")
+       contents[i] = NEF;
+      else if (title=="Z")
+       contents[i] = Z;
+      else if (title=="p_{T}^{corr}")
+       contents[i] = CorrPt;
+      else if (title=="p_{T}^{MC}")
+       contents[i] = MCPt;
+      else 
+       AliWarning(Form("Unable to fill dimension %s!",title.Data()));
+    }
+
+    histo->Fill(contents);
+  }
+  else {
+    if (Set == 1) {
+      fHistJets1PtArea->Fill(A, Pt);
+      fHistJets1PhiEta->Fill(Eta, Phi);
+      fHistJets1ZvsPt->Fill(Z, Pt);
+      fHistJets1NEFvsPt->Fill(NEF, Pt);
+      fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
+      if (!fRhoName.IsNull())
+       fHistJets1CorrPtArea->Fill(A, CorrPt);
+    }
+    else if (Set == 2) {
+      fHistJets2PtArea->Fill(A, Pt);
+      fHistJets2PhiEta->Fill(Eta, Phi);
+      if (!fRho2Name.IsNull())
+       fHistJets2CorrPtArea->Fill(A, CorrPt);
+    }
+    else if (Set == 3) {
+      fHistJets2PtAreaAcceptance->Fill(A, Pt);
+      fHistJets2PhiEtaAcceptance->Fill(Eta, Phi);
+      fHistJets2ZvsPt->Fill(Z, Pt);
+      fHistJets2NEFvsPt->Fill(NEF, Pt);
+      fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
+      if (!fRho2Name.IsNull())
+       fHistJets2CorrPtAreaAcceptance->Fill(A, CorrPt);
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, Double_t MCPt1)
+{
+  if (fHistoType==1) {
+    Double_t contents[20]={0};
+
+    for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
+      TString title(fHistMatching->GetAxis(i)->GetTitle());
+      if (title=="p_{T,1}")
+       contents[i] = Pt1;
+      else if (title=="p_{T,2}")
+       contents[i] = Pt2;
+      else if (title=="A_{jet,1}")
+       contents[i] = A1;
+      else if (title=="A_{jet,2}")
+       contents[i] = A2;
+      else if (title=="distance")
+       contents[i] = d;
+      else if (title=="CE1")
+       contents[i] = CE1;
+      else if (title=="CE2")
+       contents[i] = CE2;
+      else if (title=="#deltaA_{jet}")
+       contents[i] = A1-A2;
+      else if (title=="#deltap_{T}")
+       contents[i] = Pt1-Pt2;
+      else if (title=="#delta#eta")
+       contents[i] = Eta1-Eta2;
+      else if (title=="#delta#phi")
+       contents[i] = Phi1-Phi2;
+      else if (title=="p_{T,1}^{corr}")
+       contents[i] = CorrPt1;
+      else if (title=="p_{T,2}^{corr}")
+       contents[i] = CorrPt2;
+      else if (title=="#deltap_{T}^{corr}")
+       contents[i] = CorrPt1-CorrPt2;
+      else if (title=="p_{T,1}^{MC}")
+       contents[i] = MCPt1;
+      else if (title=="#deltap_{T}^{MC}")
+       contents[i] = MCPt1-Pt2;
+      else 
+       AliWarning(Form("Unable to fill dimension %s!",title.Data()));
+    }
+
+    fHistMatching->Fill(contents);
+  }
+  else {
+    fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1);
+    fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2);
+
+    fHistDistancevsJet1Pt->Fill(d, Pt1);
+    fHistDistancevsJet2Pt->Fill(d, Pt2);
+
+    fHistDistancevsCommonEnergy1->Fill(d, CE1);
+    fHistDistancevsCommonEnergy2->Fill(d, CE2);
+    fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
+
+    fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2);
+
+    fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1);
+    fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2);
+
+    Double_t dpt = Pt1 - Pt2;
+    fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt);
+    fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt);
+    fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1);
+    fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2);
+
+    fHistDeltaPtvsDistance->Fill(d, dpt);
+    fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
+    fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
+
+    fHistDeltaPtvsArea1->Fill(A1, dpt);
+    fHistDeltaPtvsArea2->Fill(A2, dpt);
+
+    Double_t darea = A1 - A2;
+    fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
+
+    fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
+
+    if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
+      Double_t dcorrpt = CorrPt1 - CorrPt2;
+      fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
+      fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
+      fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1);
+      fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2);
+      fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
+      fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
+      fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
+      fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt);
+      fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt);
+      fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
+      fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2);
+    }
+
+    if (fIsEmbedded) {
+      Double_t dmcpt = MCPt1 - Pt2;
+      fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt);
+      fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt);
+      fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1);
+      fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2);
+      fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
+      fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
+      fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
+      fHistDeltaMCPtvsArea1->Fill(A1, dmcpt);
+      fHistDeltaMCPtvsArea2->Fill(A2, dmcpt);
+      fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
+      fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2);
+    }
+  }
+}
+
+//________________________________________________________________________
 Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
 {   
   // Return true if jet is accepted.
@@ -760,6 +1122,8 @@ Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
     return kFALSE;
   if (jet->Area() < fJetAreaCut)
     return kFALSE;
+  if (jet->MCPt() < fMinJetMCPt)
+    return kFALSE;
 
   return kTRUE;
 }
@@ -870,6 +1234,27 @@ void AliJetResponseMaker::ExecOnce()
   if (fJet2MaxPhi == -999)
     fJet2MaxPhi = fJetMaxPhi + fJetRadius;
 
+  if (fMatching == kMCLabel && (!fAreCollections2MC || fAreCollections1MC)) {
+    if (fAreCollections2MC == fAreCollections1MC) {
+      AliWarning("Changing matching type from MC label to same collection...");
+      fMatching = kSameCollections;
+    }
+    else {
+      AliWarning("Changing matching type from MC label to geometrical...");
+      fMatching = kGeometrical;
+    }
+  }
+  else if (fMatching == kSameCollections && fAreCollections2MC != fAreCollections1MC) {
+    if (fAreCollections2MC && !fAreCollections1MC) {
+      AliWarning("Changing matching type from same collection to MC label...");
+      fMatching = kMCLabel;
+    }
+    else {
+      AliWarning("Changing matching type from same collection to geometrical...");
+      fMatching = kGeometrical;
+    }
+  }
+
   AliAnalysisTaskEmcalJet::ExecOnce();
 }
 
@@ -986,9 +1371,6 @@ void AliJetResponseMaker::DoJetLoop(Bool_t order)
     if (!AcceptJet(jet1))
       continue;
 
-    if (jet1->MCPt() < fMinJetMCPt)
-      continue;
-
     if (order) {
      if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
        continue;
@@ -1389,16 +1771,13 @@ Bool_t AliJetResponseMaker::FillHistograms()
 
   static Int_t indexes[9999] = {-1};
 
-  GetSortedArray(indexes, fJets2, fRho2Val);
-
-  const Int_t nJets2 = fJets2->GetEntriesFast();
+  const Int_t nJets2 = GetSortedArray(indexes, fJets2, fRho2Val);
 
   Int_t naccJets2 = 0;
   Int_t naccJets2Acceptance = 0;
 
-  Double_t totalPt2 = 0;
-
   for (Int_t i = 0; i < nJets2; i++) {
+    AliDebug(2,Form("Processing jet (2) %d", indexes[i]));
 
     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(indexes[i]));
 
@@ -1407,93 +1786,54 @@ Bool_t AliJetResponseMaker::FillHistograms()
       continue;
     }
 
+    // Jet 2 cuts
     if (!AcceptJet(jet2))
       continue;
-
-    if (AcceptBiasJet(jet2) &&
-       (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
-      
-      fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
-      fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
-
-      totalPt2 += jet2->Pt();
-      
-      if (naccJets2Acceptance < fNLeadingJets)
-       fHistLeadingJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
-      
-      if (!fRho2Name.IsNull()) {
-       fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
-       if (naccJets2Acceptance < fNLeadingJets)
-         fHistLeadingJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
-      }
-
-      if (fTracks2) {
-       for (Int_t it = 0; it < jet2->GetNumberOfTracks(); it++) {
-         AliVParticle *track2 = jet2->TrackAt(it, fTracks2);
-         if (track2) 
-           fHistJets2ZvsPt->Fill(track2->Pt() / jet2->Pt(), jet2->Pt());
-       }
-      }
-
-      if (fCaloClusters2) {
-       for (Int_t ic = 0; ic < jet2->GetNumberOfClusters(); ic++) {
-         AliVCluster *cluster2 = jet2->ClusterAt(ic, fCaloClusters2);
-         
-         if (cluster2) {
-           TLorentzVector nPart2;
-           cluster2->GetMomentum(nPart2, fVertex);
-           fHistJets2ZvsPt->Fill(nPart2.Et() / jet2->Pt(), jet2->Pt());
-         }
-       }
-      }
-
-      fHistJets2NEFvsPt->Fill(jet2->NEF(), jet2->Pt());
-      fHistJets2CEFvsCEFPt->Fill(1-jet2->NEF(), (1-jet2->NEF())*jet2->Pt());
-      
-      naccJets2Acceptance++;
-    }
-
     if (!AcceptBiasJet2(jet2))
       continue;
-
     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
       continue;
-
     if (jet2->MaxTrackPt() > fMaxTrackPt2 || jet2->MaxClusterPt() > fMaxClusterPt2)
       continue;
-    
-    fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
-    fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
-    
-    if (naccJets2 < fNLeadingJets)
-      fHistLeadingJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 
-    if (!fRho2Name.IsNull()) {
-      fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
-      if (naccJets2 < fNLeadingJets)
+    if (naccJets2 < fNLeadingJets) {
+      fHistLeadingJets2PtArea->Fill(jet2->Area(), jet2->Pt());
+      if (!fRho2Name.IsNull())
        fHistLeadingJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
     }
 
+    if (fDoJet2Histogram)
+      FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), jet2->Pt() - fRho2Val * jet2->Area(), jet2->MCPt(), 2);
+
     naccJets2++;
 
+    // Verify also jet cuts 1 on jet 2
+    if (AcceptBiasJet(jet2) &&
+       (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
+      
+      if (naccJets2Acceptance < fNLeadingJets) {
+       fHistLeadingJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
+       if (!fRho2Name.IsNull()) 
+         fHistLeadingJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
+      }
+      
+      FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), jet2->Pt() - fRho2Val * jet2->Area(), jet2->MCPt(), 3);
+      naccJets2Acceptance++;
+    }
+
     if (jet2->MatchedJet()) {
 
-      if (!AcceptJet(jet2->MatchedJet()) ||
-         jet2->MatchedJet()->Eta() < fJetMinEta || jet2->MatchedJet()->Eta() > fJetMaxEta || 
-         jet2->MatchedJet()->Phi() < fJetMinPhi || jet2->MatchedJet()->Phi() > fJetMaxPhi ||
-         !AcceptBiasJet(jet2->MatchedJet()) || 
-         jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt) {
-       fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
-      }
-      else {
-       if (jet2->MatchedJet()->Pt() > fMaxBinPt)
-         fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
+      if (AcceptJet(jet2->MatchedJet()) &&
+         jet2->MatchedJet()->Eta() > fJetMinEta && jet2->MatchedJet()->Eta() < fJetMaxEta &&
+         jet2->MatchedJet()->Phi() > fJetMinPhi && jet2->MatchedJet()->Phi() < fJetMaxPhi &&
+         AcceptBiasJet(jet2->MatchedJet()) &&
+         jet2->MatchedJet()->MaxTrackPt() < fMaxTrackPt && jet2->MatchedJet()->MaxClusterPt() < fMaxClusterPt) {
 
        Double_t d=-1, ce1=-1, ce2=-1;
        if (jet2->GetMatchingType() == kGeometrical) {
          if (fAreCollections2MC && !fAreCollections1MC) // the other way around is not supported
            GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, ce1, ce2);
-         else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
+         else if (fAreCollections1MC == fAreCollections2MC)
            GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, ce1, ce2);
 
          d = jet2->ClosestJetDistance();
@@ -1505,163 +1845,45 @@ Bool_t AliJetResponseMaker::FillHistograms()
          ce2 = jet2->ClosestJetDistance();
        }
 
-       fHistCommonEnergy1vsJet1Pt->Fill(ce1, jet2->MatchedJet()->Pt());
-       fHistCommonEnergy2vsJet2Pt->Fill(ce2, jet2->Pt());
-
-       fHistDistancevsJet1Pt->Fill(d, jet2->MatchedJet()->Pt());
-       fHistDistancevsJet2Pt->Fill(d, jet2->Pt());
-
-       fHistDistancevsCommonEnergy1->Fill(d, ce1);
-       fHistDistancevsCommonEnergy2->Fill(d, ce2);
-       fHistCommonEnergy1vsCommonEnergy2->Fill(ce1, ce2);
-
-       fHistJet2PtOverJet1PtvsJet2Pt->Fill(jet2->Pt(), jet2->Pt() / jet2->MatchedJet()->Pt());
-       fHistJet1PtOverJet2PtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), jet2->MatchedJet()->Pt() / jet2->Pt());
-
-       Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
-       Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
-       fHistDeltaEtaPhi->Fill(deta, dphi, jet2->Pt());
-
-       Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
-       fHistDeltaPtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dpt);
-       fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
-       fHistDeltaPtOverJet1PtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dpt/jet2->MatchedJet()->Pt());
-       fHistDeltaPtOverJet2PtvsJet2Pt->Fill(jet2->Pt(), dpt/jet2->Pt());
-
-       fHistDeltaPtvsDistance->Fill(d, dpt);
-       fHistDeltaPtvsCommonEnergy1->Fill(ce1, dpt);
-       fHistDeltaPtvsCommonEnergy2->Fill(ce2, dpt);
-
-       fHistDeltaPtvsArea1->Fill(jet2->MatchedJet()->Area(), dpt);
-       fHistDeltaPtvsArea2->Fill(jet2->Area(), dpt);
-
-       Double_t darea = jet2->MatchedJet()->Area() - jet2->Area();
-       fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
-
-       fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
-
-       if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
-         Double_t corrpt1 = jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area();
-         Double_t corrpt2 = jet2->Pt() - fRho2Val * jet2->Area();
-         Double_t dcorrpt = corrpt1 - corrpt2;
-         fHistDeltaCorrPtvsJet1CorrPt->Fill(corrpt1, dcorrpt);
-         fHistDeltaCorrPtvsJet2CorrPt->Fill(corrpt2, dcorrpt);
-         fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(corrpt1, dcorrpt/corrpt1);
-         fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(corrpt2, dcorrpt/corrpt2);
-         fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
-         fHistDeltaCorrPtvsCommonEnergy1->Fill(ce1, dcorrpt);
-         fHistDeltaCorrPtvsCommonEnergy2->Fill(ce2, dcorrpt);
-         fHistDeltaCorrPtvsArea1->Fill(jet2->MatchedJet()->Area(), dcorrpt);
-         fHistDeltaCorrPtvsArea2->Fill(jet2->Area(), dcorrpt);
-         fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
-         fHistJet1CorrPtvsJet2CorrPt->Fill(corrpt1, corrpt2);
-       }
+       Double_t corrpt1 = jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area();
+       Double_t corrpt2 = jet2->Pt() - fRho2Val * jet2->Area();
 
-       if (fIsEmbedded) {
-         Double_t dmcpt = jet2->MatchedJet()->MCPt() - jet2->Pt();
-         fHistDeltaMCPtvsJet1MCPt->Fill(jet2->MatchedJet()->MCPt(), dmcpt);
-         fHistDeltaMCPtvsJet2Pt->Fill(jet2->Pt(), dmcpt);
-         fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(jet2->MatchedJet()->MCPt(), dmcpt/jet2->MatchedJet()->MCPt());
-         fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(jet2->Pt(), dmcpt/jet2->Pt());
-         fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
-         fHistDeltaMCPtvsCommonEnergy1->Fill(ce1, dmcpt);
-         fHistDeltaMCPtvsCommonEnergy2->Fill(ce2, dmcpt);
-         fHistDeltaMCPtvsArea1->Fill(jet2->MatchedJet()->Area(), dmcpt);
-         fHistDeltaMCPtvsArea2->Fill(jet2->Area(), dmcpt);
-         fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
-         fHistJet1MCPtvsJet2Pt->Fill(jet2->MatchedJet()->MCPt(), jet2->Pt());
-       }
+       FillMatchingHistos(jet2->MatchedJet()->Pt(), jet2->Pt(), jet2->MatchedJet()->Eta(), jet2->Eta(), jet2->MatchedJet()->Phi(), jet2->Phi(), jet2->MatchedJet()->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet2->MatchedJet()->MCPt());
       }
     }
-    else {
-      fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
-      fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
-
-      if (!fRho2Name.IsNull())
-       fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
-    }
   }
 
-  GetSortedArray(indexes, fJets, fRhoVal);
+  const Int_t nJets1 = GetSortedArray(indexes, fJets, fRhoVal);
 
-  const Int_t nJets1 = fJets->GetEntriesFast();
   Int_t naccJets1 = 0;
-  Double_t totalMCPt1 = 0;
 
   for (Int_t i = 0; i < nJets1; i++) {
-
-    AliDebug(2,Form("Processing jet %d", i));
+    AliDebug(2,Form("Processing jet (1) %d", indexes[i]));
     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(indexes[i]));
 
     if (!jet1) {
       AliError(Form("Could not receive jet %d", i));
       continue;
     }  
-    
+
     if (!AcceptJet(jet1))
       continue;
-
     if (!AcceptBiasJet(jet1))
       continue;
-
     if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
       continue;
-
     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
       continue;
 
-    if (jet1->MCPt() < fMinJetMCPt)
-      continue;
-
-    if (!jet1->MatchedJet()) {
-      fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
-      if (!fRhoName.IsNull())
-       fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
-    }
-
-    fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
-    fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
-
-    totalMCPt1 += jet1->MCPt();
-
-    if (naccJets1 < fNLeadingJets)
+    if (naccJets1 < fNLeadingJets){
       fHistLeadingJets1PtArea->Fill(jet1->Area(), jet1->Pt());
-
-    if (!fRhoName.IsNull()) {
-      fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
-
-      if (naccJets1 < fNLeadingJets)
+      if (!fRhoName.IsNull())
        fHistLeadingJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
     }
 
-    if (fTracks) {
-      for (Int_t it = 0; it < jet1->GetNumberOfTracks(); it++) {
-       AliVParticle *track1 = jet1->TrackAt(it, fTracks);
-       if (track1) 
-         fHistJets1ZvsPt->Fill(track1->Pt() / jet1->Pt(), jet1->Pt());
-      }
-    }
-    
-    if (fCaloClusters) {
-      for (Int_t ic = 0; ic < jet1->GetNumberOfClusters(); ic++) {
-       AliVCluster *cluster1 = jet1->ClusterAt(ic, fCaloClusters);
-       
-       if (cluster1) {
-         TLorentzVector nPart1;
-         cluster1->GetMomentum(nPart1, fVertex);
-         fHistJets2ZvsPt->Fill(nPart1.Et() / jet1->Pt(), jet1->Pt());
-       }
-      }
-    }
-    
-    fHistJets1NEFvsPt->Fill(jet1->NEF(), jet1->Pt());
-    fHistJets1CEFvsCEFPt->Fill(1-jet1->NEF(), (1-jet1->NEF())*jet1->Pt());
-
+    FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), jet1->MaxPartPt()/jet1->Pt(), jet1->Pt() - fRhoVal * jet1->Area(), jet1->MCPt(), 1);
     naccJets1++;
   }
 
-  if (fIsEmbedded) 
-    fMCEnergy1vsEnergy2->Fill(totalMCPt1, totalPt2);
-
   return kTRUE;
 }
index 9f5d6f7..e95d64f 100644 (file)
@@ -5,6 +5,7 @@
 
 class TClonesArray;
 class TH2;
+class THnSparse;
 class AliNamedArrayI;
 
 #include "AliEmcalJet.h"
@@ -42,6 +43,10 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   void                        SetMinJetMCPt(Float_t pt)                                       { fMinJetMCPt        = pt        ; }
   void                        SetMaxClusterPt2(Float_t b)                                     { fMaxClusterPt2     = b         ; }
   void                        SetMaxTrackPt2(Float_t b)                                       { fMaxTrackPt2       = b         ; }
+  void                        SetHistoType(Int_t b)                                           { fHistoType         = b         ; }
+  void                        SetDeltaPtAxis(Int_t b)                                         { fDeltaPtAxis       = b         ; }
+  void                        SetDeltaEtaDeltaPhiAxis(Int_t b)                                { fDeltaEtaDeltaPhiAxis= b       ; }
+  void                        SetDoJet2Histogram(Int_t b)                                     { fDoJet2Histogram   = b         ; }
 
  protected:
   Bool_t                      AcceptJet(AliEmcalJet* jet) const;
@@ -56,6 +61,10 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   void                        GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const;
   void                        GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const;
   void                        GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const;
+  void                        FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, Double_t MCPt1);
+  void                        FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t Z, Double_t CorrPt, Double_t MCPt, Int_t Set);
+  void                        AllocateTH2();
+  void                        AllocateTHnSparse();
 
   TString                     fTracks2Name;                   // name of second track collection
   TString                     fCalo2Name;                     // name of second cluster collection
@@ -78,6 +87,10 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   Double_t                    fMatchingPar2;                  // matching parameter for jet2-jet1 matching
   Bool_t                      fUseCellsToMatch;               // use cells instead of clusters to match jets (slower but sometimes needed)
   Double_t                    fMinJetMCPt;                    // minimum jet MC pt
+  Int_t                       fHistoType;                     // histogram type (0=TH2, 1=THnSparse)
+  Int_t                       fDeltaPtAxis;                   // add delta pt axis in THnSparse (default=0)
+  Int_t                       fDeltaEtaDeltaPhiAxis;          // add delta eta and delta phi axes in THnSparse (default=0)
+  Int_t                       fDoJet2Histogram;               // add unbiased jet2 histogram (potentially memory consuming if on particle level)
   TClonesArray               *fTracks2;                       //!Tracks 2
   TClonesArray               *fCaloClusters2;                 //!Clusters 2
   TClonesArray               *fJets2;                         //!Jets 2
@@ -85,37 +98,40 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   Double_t                    fRho2Val;                       //!Event rho 2 value 
   AliNamedArrayI             *fTracks2Map;                    //!MC particle map
 
-  // General histograms
-  TH2                        *fMCEnergy1vsEnergy2;                     //!total MC energy jet 1 vs total energy jet 2
+  // General histos
+  TH2                        *fHistLeadingJets1PtArea;                 //!leading jet pt vs. area histogram 1
+  TH2                        *fHistLeadingJets1CorrPtArea;             //!leading jet corr pt vs. area histogram 1
+  TH2                        *fHistLeadingJets2PtArea;                 //!leading jet pt vs. area histogram 2
+  TH2                        *fHistLeadingJets2CorrPtArea;             //!leading jet corr pt vs. area histogram 2
+  TH2                        *fHistLeadingJets2PtAreaAcceptance;       //!leading jet pt vs. area histogram 2 using jet 1 cuts (acceptance, leading hadron bias, ...)
+  TH2                        *fHistLeadingJets2CorrPtAreaAcceptance;   //!leading jet corr pt vs. area histogram 2 using jet 1 cuts (acceptance, leading hadron bias, ...)
+
+  // THnSparse
+  THnSparse                  *fHistJets1;                              //!jet1 THnSparse
+  THnSparse                  *fHistJets2;                              //!jet2 THnSparse
+  THnSparse                  *fHistJets2Acceptance;                    //!jet2 acceptance THnSparse
+  THnSparse                  *fHistMatching;                           //!matching THnSparse
+
   // Jets 1
   TH2                        *fHistJets1PhiEta;                        //!phi-eta distribution of jets 1
-  TH2                        *fHistJets1PtArea;                        //!inclusive jet pt vs area histogram 1
+  TH2                        *fHistJets1PtArea;                        //!inclusive jet pt vs. area histogram 1
   TH2                        *fHistJets1CorrPtArea;                    //!inclusive jet pt vs. area histogram 1
-  TH2                        *fHistLeadingJets1PtArea;                 //!leading jet pt vs area histogram 1
-  TH2                        *fHistLeadingJets1CorrPtArea;             //!leading jet pt vs. area histogram 1
   TH2                        *fHistJets1NEFvsPt;                       //!Jet neutral energy fraction vs. jet pt 1
   TH2                        *fHistJets1CEFvsCEFPt;                    //!Jet charged energy fraction vs. charged jet pt 1
   TH2                        *fHistJets1ZvsPt;                         //!Constituent Pt over Jet Pt ratio vs. jet pt 1
+
   // Jets 2
   TH2                        *fHistJets2PhiEta;                        //!phi-eta distribution of jets 2
   TH2                        *fHistJets2PtArea;                        //!inclusive jet pt vs. area histogram 2
   TH2                        *fHistJets2CorrPtArea;                    //!inclusive jet pt vs. area histogram 2
-  TH2                        *fHistLeadingJets2PtArea;                 //!leading jet pt vs. area histogram 2
-  TH2                        *fHistLeadingJets2CorrPtArea;             //!leading jet pt vs. area histogram 2
   TH2                        *fHistJets2PhiEtaAcceptance;              //!phi-eta distribution of jets 2 using jet 1 cuts (acceptance, leading hadron bias, ...)
   TH2                        *fHistJets2PtAreaAcceptance;              //!inclusive jet pt vs. area histogram 2 using jet 1 cuts (acceptance, leading hadron bias, ...)
   TH2                        *fHistJets2CorrPtAreaAcceptance;          //!inclusive jet pt vs. area histogram 2 using jet 1 cuts (acceptance, leading hadron bias, ...)
-  TH2                        *fHistLeadingJets2PtAreaAcceptance;       //!leading jet pt vs. area histogram 2 using jet 1 cuts (acceptance, leading hadron bias, ...)
-  TH2                        *fHistLeadingJets2CorrPtAreaAcceptance;   //!leading jet pt vs. area histogram 2 using jet 1 cuts (acceptance, leading hadron bias, ...)
   TH2                        *fHistJets2NEFvsPt;                       //!Jet neutral energy fraction vs. jet pt 2
   TH2                        *fHistJets2CEFvsCEFPt;                    //!Jet charged energy fraction vs. charged jet pt 2
   TH2                        *fHistJets2ZvsPt;                         //!Constituent Pt over Jet Pt ratio vs. jet pt 2
+
   // Jet1-Jet2 matching
-  TH2                        *fHistNonMatchedJets1PtArea;              //!non-matched jet 1 pt distribution
-  TH2                        *fHistNonMatchedJets2PtArea;              //!non-matched jet 2 pt distribution
-  TH2                        *fHistNonMatchedJets1CorrPtArea;          //!non-matched jet pt distribution
-  TH2                        *fHistNonMatchedJets2CorrPtArea;          //!non-matched jet pt distribution
-  TH2                        *fHistMissedJets2PtArea;                  //!jets 2 not found in jet 1 collection
   TH2                        *fHistCommonEnergy1vsJet1Pt;              //!common energy 1 (%) vs jet 1 pt
   TH2                        *fHistCommonEnergy2vsJet2Pt;              //!common energy 2 (%) vs jet 2 pt
   TH2                        *fHistDistancevsJet1Pt;                   //!distance vs jet 1 pt
@@ -123,9 +139,10 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   TH2                        *fHistDistancevsCommonEnergy1;            //!distance vs common energy 1 (%)
   TH2                        *fHistDistancevsCommonEnergy2;            //!distance vs common energy 2 (%)
   TH2                        *fHistCommonEnergy1vsCommonEnergy2;       //!common energy 1 (%) vs common energy 2 (%)
+  TH2                        *fHistDeltaEtaDeltaPhi;                   //!delta eta vs delta phi of matched jets
+
   TH2                        *fHistJet2PtOverJet1PtvsJet2Pt;           //!jet 2 pt over jet 1 pt vs jet 2 pt
   TH2                        *fHistJet1PtOverJet2PtvsJet1Pt;           //!jet 1 pt over jet 2 pt vs jet 1 pt
-  TH2                        *fHistDeltaEtaPhi;                        //!delta eta-phi between matched jets
   TH2                        *fHistDeltaPtvsJet1Pt;                    //!delta pt between matched jets vs jet 1 pt
   TH2                        *fHistDeltaPtvsJet2Pt;                    //!delta pt between matched jets vs jet 2 pt
   TH2                        *fHistDeltaPtOverJet1PtvsJet1Pt;          //!delta pt / jet 1 pt between matched jets vs jet 1 pt
@@ -137,10 +154,11 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   TH2                        *fHistDeltaPtvsArea2;                     //!delta pt between matched jets vs jet 2 area
   TH2                        *fHistDeltaPtvsDeltaArea;                 //!delta pt between matched jets vs delta area
   TH2                        *fHistJet1PtvsJet2Pt;                     //!correlation jet 1 pt vs jet 2 pt
-  TH2                        *fHistDeltaCorrPtvsJet1CorrPt;            //!delta pt corr / jet 1 corr pt between matched jets vs jet 1 corr pt
-  TH2                        *fHistDeltaCorrPtvsJet2CorrPt;            //!delta pt corr / jet 2 corr pt between matched jets vs jet 2 corr pt
-  TH2                        *fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt;//!delta pt corr between matched jets vs jet 1 corr pt
-  TH2                        *fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt;//!delta pt corr between matched jets vs jet 2 corr pt
+
+  TH2                        *fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt;//!delta pt corr / jet 1 corr pt between matched jets vs jet 1 corr pt
+  TH2                        *fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt;//!delta pt corr / jet 2 corr pt between matched jets vs jet 2 corr pt
+  TH2                        *fHistDeltaCorrPtvsJet1CorrPt;            //!delta pt corr between matched jets vs jet 1 corr pt 
+  TH2                        *fHistDeltaCorrPtvsJet2CorrPt;            //!delta pt corr between matched jets vs jet 2 corr pt
   TH2                        *fHistDeltaCorrPtvsDistance;              //!delta pt corr between matched jets vs distance
   TH2                        *fHistDeltaCorrPtvsCommonEnergy1;         //!delta pt corr between matched jets vs common energy 1 (%)
   TH2                        *fHistDeltaCorrPtvsCommonEnergy2;         //!delta pt corr between matched jets vs common energy 2 (%)
@@ -148,10 +166,11 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   TH2                        *fHistDeltaCorrPtvsArea2;                 //!delta pt corr between matched jets vs jet 2 area
   TH2                        *fHistDeltaCorrPtvsDeltaArea;             //!delta pt corr between matched jets vs delta area
   TH2                        *fHistJet1CorrPtvsJet2CorrPt;             //!correlation jet 1 corr pt vs jet 2 corr pt
-  TH2                        *fHistDeltaMCPtvsJet1MCPt;                //!jet 1 MC pt - jet2 pt vs jet 1 MC pt
-  TH2                        *fHistDeltaMCPtvsJet2Pt;                  //!jet 1 MC pt - jet2 pt vs jet 2 pt
+
   TH2                        *fHistDeltaMCPtOverJet1MCPtvsJet1MCPt;    //!jet 1 MC pt - jet2 pt / jet 1 MC pt vs jet 1 pt
   TH2                        *fHistDeltaMCPtOverJet2PtvsJet2Pt;        //!jet 1 MC pt - jet2 pt / jet 2 pt vs jet 2 pt
+  TH2                        *fHistDeltaMCPtvsJet1MCPt;                //!jet 1 MC pt - jet2 pt vs jet 1 MC pt
+  TH2                        *fHistDeltaMCPtvsJet2Pt;                  //!jet 1 MC pt - jet2 pt vs jet 2 pt
   TH2                        *fHistDeltaMCPtvsDistance;                //!jet 1 MC pt - jet2 pt vs distance
   TH2                        *fHistDeltaMCPtvsCommonEnergy1;           //!jet 1 MC pt - jet2 pt vs common energy 1 (%)
   TH2                        *fHistDeltaMCPtvsCommonEnergy2;           //!jet 1 MC pt - jet2 pt vs common energy 2 (%)
@@ -164,6 +183,6 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   AliJetResponseMaker(const AliJetResponseMaker&);            // not implemented
   AliJetResponseMaker &operator=(const AliJetResponseMaker&); // not implemented
 
-  ClassDef(AliJetResponseMaker, 19) // Jet response matrix producing task
+  ClassDef(AliJetResponseMaker, 21) // Jet response matrix producing task
 };
 #endif
index 3fc149d..5843f8d 100644 (file)
@@ -5,15 +5,13 @@
 // Author: S.Aiola
 
 #include <TClonesArray.h>
-#include <TH1F.h>
 #include <TH2F.h>
-#include <TH3F.h>
+#include <THnSparse.h>
 #include <TList.h>
 #include <TLorentzVector.h>
 
 #include "AliVCluster.h"
 #include "AliVParticle.h"
-#include "AliVTrack.h"
 #include "AliEmcalJet.h"
 #include "AliRhoParameter.h"
 #include "AliLog.h"
@@ -25,32 +23,17 @@ ClassImp(AliAnalysisTaskSAJF)
 //________________________________________________________________________
 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() : 
   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
-  fNjetsVsCent(0)
+  fHistEventObservables(0),
+  fHistJetObservables(0)
 
 {
   // Default constructor.
 
   for (Int_t i = 0; i < 4; i++) {
-    fHistLeadingJetPhiEta[i] = 0;
-    fHistLeadingJetPtArea[i] = 0;
-    fHistLeadingJetCorrPtArea[i] = 0;
-    fHistLeadingJetMCPtArea[i] = 0;
-    fHistRhoVSleadJetPt[i] = 0;
-    fHistJetPhiEta[i] = 0;
-    fHistJetsPtArea[i] = 0;
-    fHistJetsCorrPtArea[i] = 0;
-    fHistJetPtvsJetCorrPt[i] = 0;
-    fHistJetsMCPtArea[i] = 0;
-    fHistJetPtvsJetMCPt[i] = 0;
-    fHistJetsNEFvsPt[i] = 0;
-    fHistJetsCEFvsCEFPt[i] = 0;
-    fHistJetsZvsPt[i] = 0;
-    fHistConstituents[i] = 0;
     fHistTracksJetPt[i] = 0;
     fHistClustersJetPt[i] = 0;
     fHistTracksPtDist[i] = 0;
     fHistClustersPtDist[i] = 0;
-    fHistJetNconstVsPt[i] = 0;
   }
 
   SetMakeGeneralHistograms(kTRUE);
@@ -59,248 +42,189 @@ AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
 //________________________________________________________________________
 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) : 
   AliAnalysisTaskEmcalJet(name, kTRUE),
-  fNjetsVsCent(0)
+  fHistEventObservables(0),
+  fHistJetObservables(0)
 {
   // Standard constructor.
 
   for (Int_t i = 0; i < 4; i++) {
-    fHistLeadingJetPhiEta[i] = 0;
-    fHistLeadingJetPtArea[i] = 0;
-    fHistLeadingJetCorrPtArea[i] = 0;
-    fHistLeadingJetMCPtArea[i] = 0;
-    fHistRhoVSleadJetPt[i] = 0;
-    fHistJetPhiEta[i] = 0;
-    fHistJetsPtArea[i] = 0;
-    fHistJetsCorrPtArea[i] = 0;
-    fHistJetPtvsJetCorrPt[i] = 0;
-    fHistJetsMCPtArea[i] = 0;
-    fHistJetPtvsJetMCPt[i] = 0;
-    fHistJetsNEFvsPt[i] = 0;
-    fHistJetsCEFvsCEFPt[i] = 0;
-    fHistJetsZvsPt[i] = 0;
-    fHistConstituents[i] = 0;
     fHistTracksJetPt[i] = 0;
     fHistClustersJetPt[i] = 0;
     fHistTracksPtDist[i] = 0;
     fHistClustersPtDist[i] = 0;
-    fHistJetNconstVsPt[i] = 0;
   }
 
   SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
-AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
-{
-  // Destructor.
-}
-
-//________________________________________________________________________
 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
 {
   // Create user output.
 
   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 
-  fNjetsVsCent = new TH2F("fNjetsVsCent","fNjetsVsCent", 100, 0, 100, 150, -0.5, 149.5);
-  fNjetsVsCent->GetXaxis()->SetTitle("Centrality (%)");
-  fNjetsVsCent->GetYaxis()->SetTitle("No. of jets");
-  fOutput->Add(fNjetsVsCent);
+  TString title[20]= {""};
+  Int_t nbins[20]  = {0};
+  Double_t min[20] = {0.};
+  Double_t max[20] = {0.};
+  Int_t dim = 0;
+
+  if (fForceBeamType != kpp) {
+    title[dim] = "Centrality (%)";
+    nbins[dim] = 101;
+    min[dim] = 0;
+    max[dim] = 101;
+    dim++;
+    
+    title[dim] = "#psi_{RP} (rad)";
+    nbins[dim] = 200;
+    min[dim] = -TMath::Pi();
+    max[dim] = TMath::Pi();
+    dim++;
+  }
 
-  const Int_t nbinsZ = 12;
-  Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
+  title[dim] = "No. of jets";
+  nbins[dim] = 150;
+  min[dim] = -0.5;
+  max[dim] = 149.5;
+  dim++;
+
+  title[dim] = "p_{T,lead jet} (GeV/c)";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "A_{lead jet}";
+  nbins[dim] = 100;
+  min[dim] = 0;
+  max[dim] = 1;
+  dim++;
+
+  if (!fRhoName.IsNull()) {
+    title[dim] = "#rho (GeV/c)";
+    nbins[dim] = fNbins;
+    min[dim] = 0;
+    max[dim] = 500;
+    dim++;
+
+    title[dim] = "p_{T,lead jet}^{corr} (GeV/c)";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
 
-  Float_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
-  Float_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
-  Float_t *binsArea     = GenerateFixedBinArray(50, 0, 2);
-  Float_t *binsEta      = GenerateFixedBinArray(50,-1, 1);
-  Float_t *binsPhi      = GenerateFixedBinArray(101, 0, TMath::Pi() * 2.02);
-  Float_t *bins120      = GenerateFixedBinArray(120, 0, 1.2);
-  Float_t *bins150      = GenerateFixedBinArray(150, -0.5, 149.5);
+  fHistEventObservables = new THnSparseD("fHistEventObservables","fHistEventObservables",dim,nbins,min,max);
+  for (Int_t i = 0; i < dim; i++)
+    fHistEventObservables->GetAxis(i)->SetTitle(title[i]);
 
-  TString histname;
+  dim = 0;
 
-  for (Int_t i = 0; i < 4; i++) {
-    histname = "fHistLeadingJetPhiEta_";
-    histname += i;
-    fHistLeadingJetPhiEta[i] = new TH3F(histname.Data(), histname.Data(), 
-                                       50, binsEta, 
-                                       101, binsPhi, 
-                                       nbinsZ, binsZ);
-    fHistLeadingJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
-    fHistLeadingJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
-    fHistLeadingJetPhiEta[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistLeadingJetPhiEta[i]);
-
-    histname = "fHistLeadingJetPtArea_";
-    histname += i;
-    fHistLeadingJetPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
-                                       fNbins, binsPt, 
-                                       50, binsArea,
-                                       nbinsZ, binsZ);
-    fHistLeadingJetPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
-    fHistLeadingJetPtArea[i]->GetYaxis()->SetTitle("area");
-    fHistLeadingJetPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistLeadingJetPtArea[i]);
-
-    if (fIsEmbedded) {
-      histname = "fHistLeadingJetMCPtArea_";
-      histname += i;
-      fHistLeadingJetMCPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
-                                           fNbins, binsPt, 
-                                           50, binsArea,
-                                           nbinsZ, binsZ);
-      fHistLeadingJetMCPtArea[i]->GetXaxis()->SetTitle("p_{T,MC} (GeV/c)");
-      fHistLeadingJetMCPtArea[i]->GetYaxis()->SetTitle("area");
-      fHistLeadingJetMCPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-      fOutput->Add(fHistLeadingJetMCPtArea[i]);
-    }
+  if (fForceBeamType != kpp) {
+    title[dim] = "Centrality (%)";
+    nbins[dim] = 101;
+    min[dim] = 0;
+    max[dim] = 101;
+    dim++;
+    
+    title[dim] = "#psi_{RP} (rad)";
+    nbins[dim] = 200;
+    min[dim] = -TMath::Pi();
+    max[dim] = TMath::Pi();
+    dim++;
+  }
 
-    if (!fRhoName.IsNull()) {
-      histname = "fHistLeadingJetCorrPtArea_";
-      histname += i;
-      fHistLeadingJetCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
-                                             fNbins * 2, binsCorrPt, 
-                                             50, binsArea,
-                                             nbinsZ, binsZ);
-      fHistLeadingJetCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
-      fHistLeadingJetCorrPtArea[i]->GetYaxis()->SetTitle("area");
-      fHistLeadingJetCorrPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-      fOutput->Add(fHistLeadingJetCorrPtArea[i]);
-      
-      histname = "fHistRhoVSleadJetPt_";
-      histname += i;
-      fHistRhoVSleadJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
-      fHistRhoVSleadJetPt[i]->GetXaxis()->SetTitle("#rho * area (GeV/c)");
-      fHistRhoVSleadJetPt[i]->GetYaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
-      fOutput->Add(fHistRhoVSleadJetPt[i]);
-    }
+  title[dim] = "#phi_{jet} (rad)";
+  nbins[dim] = 201;
+  min[dim] = 0;
+  max[dim] = TMath::Pi()*201/100;
+  dim++;
+
+  title[dim] = "#eta";
+  nbins[dim] = 200;
+  min[dim] = -1;
+  max[dim] = 1;
+  dim++;
+
+  title[dim] = "p_{T} (GeV/c)";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "A_{jet}";
+  nbins[dim] = 100;
+  min[dim] = 0;
+  max[dim] = 1;
+  dim++;
+
+  if (fIsEmbedded) {
+    title[dim] = "p_{T}^{MC} (GeV/c)";
+    nbins[dim] = fNbins;
+    min[dim] = 0;
+    max[dim] = 250;
+    dim++;
+  }
 
-    histname = "fHistJetPhiEta_";
-    histname += i;
-    fHistJetPhiEta[i] = new TH3F(histname.Data(), histname.Data(), 
-                                50, binsEta, 
-                                101, binsPhi, 
-                                nbinsZ, binsZ);
-    fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
-    fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
-    fHistJetPhiEta[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistJetPhiEta[i]);
-    
-    histname = "fHistJetsPtArea_";
-    histname += i;
-    fHistJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
-                                 fNbins, binsPt, 
-                                 50, binsArea,
-                                 nbinsZ, binsZ);
-    fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
-    fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
-    fHistJetsPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistJetsPtArea[i]);
-
-    if (!fRhoName.IsNull()) {
-      histname = "fHistJetsCorrPtArea_";
-      histname += i;
-      fHistJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
-                                       fNbins * 2, binsCorrPt, 
-                                       50, binsArea,
-                                       nbinsZ, binsZ);
-      fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
-      fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
-      fHistJetsCorrPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-      fOutput->Add(fHistJetsCorrPtArea[i]);
-
-      histname = "fHistJetPtvsJetCorrPt_";
-      histname += i;
-      fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
-                                         fNbins, fMinBinPt, fMaxBinPt, 
-                                         fNbins * 2, -fMaxBinPt, fMaxBinPt);
-      fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
-      fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
-      fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
-      fOutput->Add(fHistJetPtvsJetCorrPt[i]);
-    }
+  if (!fRhoName.IsNull()) {
+    title[dim] = "p_{T}^{corr} (GeV/c)";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
 
-    if (fIsEmbedded) {
-      histname = "fHistJetsMCPtArea_";
+  title[dim] = "No. of constituents";
+  nbins[dim] = 250;
+  min[dim] = -0.5;
+  max[dim] = 249.5;
+  dim++;
+
+  title[dim] = "NEF";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "Z";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "p_{T,particle}^{leading} (GeV/c)";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 120;
+  dim++;
+
+  fHistJetObservables = new THnSparseD("fHistJetObservables","fHistJetObservables",dim,nbins,min,max);
+  for (Int_t i = 0; i < dim; i++)
+    fHistJetObservables->GetAxis(i)->SetTitle(title[i]);
+
+  for (Int_t i = 0; i < 4; i++) {
+    TString histname;
+
+    if (!fTracksName.IsNull()) {
+      histname = "fHistTracksJetPt_";
       histname += i;
-      fHistJetsMCPtArea[i] = new TH3F(histname.Data(), histname.Data(), 
-                                     fNbins, binsPt, 
-                                     50, binsArea,
-                                     nbinsZ, binsZ);
-      fHistJetsMCPtArea[i]->GetXaxis()->SetTitle("p_{T,MC} (GeV/c)");
-      fHistJetsMCPtArea[i]->GetYaxis()->SetTitle("area");
-      fHistJetsMCPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-      fOutput->Add(fHistJetsMCPtArea[i]);
-
-      histname = "fHistJetPtvsJetMCPt_";
+      fHistTracksJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
+      fHistTracksJetPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
+      fHistTracksJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
+      fHistTracksJetPt[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistTracksJetPt[i]);
+      
+      histname = "fHistTracksPtDist_";
       histname += i;
-      fHistJetPtvsJetMCPt[i] = new TH2F(histname.Data(), histname.Data(),
-                                       fNbins, fMinBinPt, fMaxBinPt, 
-                                       fNbins, fMinBinPt, fMaxBinPt);
-      fHistJetPtvsJetMCPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
-      fHistJetPtvsJetMCPt[i]->GetYaxis()->SetTitle("p_{T,MC} (GeV/c)");
-      fHistJetPtvsJetMCPt[i]->GetZaxis()->SetTitle("counts");
-      fOutput->Add(fHistJetPtvsJetMCPt[i]);
+      fHistTracksPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
+      fHistTracksPtDist[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
+      fHistTracksPtDist[i]->GetYaxis()->SetTitle("d");
+      fHistTracksPtDist[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistTracksPtDist[i]);
     }
 
-    histname = "fHistJetsZvsPt_";
-    histname += i;
-    fHistJetsZvsPt[i] = new TH3F(histname.Data(), histname.Data(), 
-                                120, bins120, 
-                                fNbins, binsPt,
-                                nbinsZ, binsZ);
-    fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
-    fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
-    fHistJetsZvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistJetsZvsPt[i]);
-
-    histname = "fHistJetsNEFvsPt_";
-    histname += i;
-    fHistJetsNEFvsPt[i] = new TH3F(histname.Data(), histname.Data(), 
-                                  120, bins120, 
-                                  fNbins, binsPt,
-                                  nbinsZ, binsZ);
-    fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
-    fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
-    fHistJetsNEFvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistJetsNEFvsPt[i]);
-    
-    histname = "fHistJetsCEFvsCEFPt_";
-    histname += i;
-    fHistJetsCEFvsCEFPt[i] = new TH3F(histname.Data(), histname.Data(), 
-                                     120, bins120, 
-                                     fNbins, binsPt,
-                                     nbinsZ, binsZ);
-    fHistJetsCEFvsCEFPt[i]->GetXaxis()->SetTitle("1-NEF");
-    fHistJetsCEFvsCEFPt[i]->GetYaxis()->SetTitle("(1-NEF)*p_{T}^{raw} (GeV/c)");
-    fHistJetsCEFvsCEFPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistJetsCEFvsCEFPt[i]);
-
-    histname = "fHistConstituents_";
-    histname += i;
-    fHistConstituents[i] = new TH2F(histname.Data(), histname.Data(), 100, 1, 101, 100, -0.5, 99.5);
-    fHistConstituents[i]->GetXaxis()->SetTitle("p_{T,part} (GeV/c)");
-    fHistConstituents[i]->GetYaxis()->SetTitle("no. of particles");
-    fHistConstituents[i]->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistConstituents[i]);
-
-    histname = "fHistTracksJetPt_";
-    histname += i;
-    fHistTracksJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
-    fHistTracksJetPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
-    fHistTracksJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
-    fHistTracksJetPt[i]->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistTracksJetPt[i]);
-
-    histname = "fHistTracksPtDist_";
-    histname += i;
-    fHistTracksPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
-    fHistTracksPtDist[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
-    fHistTracksPtDist[i]->GetYaxis()->SetTitle("d");
-    fHistTracksPtDist[i]->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistTracksPtDist[i]);
-
     if (!fCaloName.IsNull()) {
       histname = "fHistClustersJetPt_";
       histname += i;
@@ -318,24 +242,10 @@ void AliAnalysisTaskSAJF::UserCreateOutputObjects()
       fHistClustersPtDist[i]->GetZaxis()->SetTitle("counts");
       fOutput->Add(fHistClustersPtDist[i]);
     }
-
-    histname = "fHistJetNconstVsPt_";
-    histname += i;
-    fHistJetNconstVsPt[i] = new TH3F(histname.Data(), histname.Data(), 150, bins150, fNbins, binsPt, nbinsZ, binsZ);
-    fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("# of constituents");
-    fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
-    fHistJetNconstVsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
-    fOutput->Add(fHistJetNconstVsPt[i]);
   }
 
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 
-  delete[] binsPt;
-  delete[] binsCorrPt;
-  delete[] binsArea;
-  delete[] binsEta;
-  delete[] binsPhi;
-  delete[] bins120;
 }
 
 //________________________________________________________________________
@@ -348,104 +258,55 @@ Bool_t AliAnalysisTaskSAJF::FillHistograms()
     return kFALSE;
   }
 
-  if (fJets->GetEntriesFast() < 1) // no jets in array, skipping
-    return kTRUE;
-
   static Int_t sortedJets[9999] = {-1};
-  Bool_t r = GetSortedArray(sortedJets, fJets, fRhoVal);
-  
-  if (!r) // no accepted jets, skipping
-    return kTRUE;
-
-  // OK, event accepted
+  Int_t nacc = GetSortedArray(sortedJets, fJets, fRhoVal);
 
-  for (Int_t i = 0; i < fNLeadingJets && i < fJets->GetEntriesFast(); i++) {
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[i]));
+  Float_t corrPt = 0;
 
-    if (!jet) {
-      AliError(Form("Could not receive jet %d", sortedJets[i]));
-      continue;
-    }  
-
-    if (!AcceptJet(jet))
-      continue;
-
-    Float_t ptLeading = GetLeadingHadronPt(jet);
-
-    fHistLeadingJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi(), ptLeading);
-    fHistLeadingJetPtArea[fCentBin]->Fill(jet->Pt(), jet->Area(), ptLeading);
-
-    if (fIsEmbedded) 
-      fHistLeadingJetMCPtArea[fCentBin]->Fill(jet->MCPt(), jet->Area(), ptLeading);      
-
-    Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
-
-    if (fHistLeadingJetCorrPtArea[fCentBin])
-      fHistLeadingJetCorrPtArea[fCentBin]->Fill(corrPt, jet->Area(), ptLeading);
-
-    if (i==0 && fHistRhoVSleadJetPt[fCentBin]) 
-      fHistRhoVSleadJetPt[fCentBin]->Fill(fRhoVal, jet->Pt());
+  if (nacc > 0) {
+    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
+    
+    if (jet) {
+      // Don't need to do AcceptJet since it was done already in GetSortedArray
+      corrPt = jet->Pt() - fRhoVal * jet->Area();
+      
+      DoJetLoop(nacc, sortedJets);
+    } 
+    else {
+      AliError(Form("Could not receive jet %d", sortedJets[0]));
+    }
   }
 
-  Int_t njets = DoJetLoop();
-
-  fNjetsVsCent->Fill(fCent, njets);
+  // Fill THnSparse
 
   return kTRUE;
 }
 
 //________________________________________________________________________
-Int_t AliAnalysisTaskSAJF::DoJetLoop()
+void AliAnalysisTaskSAJF::DoJetLoop(const Int_t nacc, const Int_t* sortedJets)
 {
   // Do the jet loop.
 
-  if (!fJets)
-    return 0;
+  for (Int_t ij = 0; ij < nacc; ij++) {
 
-  const Int_t njets = fJets->GetEntriesFast();
-  Int_t nAccJets = 0;
-
-  TH1F constituents("constituents", "constituents", 
-                   fHistConstituents[0]->GetNbinsX(), fHistConstituents[0]->GetXaxis()->GetXmin(), fHistConstituents[0]->GetXaxis()->GetXmax()); 
-
-  for (Int_t ij = 0; ij < njets; ij++) {
-
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
+    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[ij]));
 
     if (!jet) {
       AliError(Form("Could not receive jet %d", ij));
       continue;
-    }  
-
-    if (!AcceptJet(jet))
-      continue;
+    }
 
-    Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
+    // Don't need to do AcceptJet since it was done already in GetSortedArray
 
     Float_t ptLeading = GetLeadingHadronPt(jet);
+    Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
 
-    fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi(), ptLeading);
-    fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area(), ptLeading);
-    if (fHistJetsCorrPtArea[fCentBin])
-      fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area(), ptLeading);
-    if (fHistJetPtvsJetCorrPt[fCentBin])
-      fHistJetPtvsJetCorrPt[fCentBin]->Fill(jet->Pt(), corrPt);
-    fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt(), ptLeading);
-
-    if (fIsEmbedded) {
-      fHistJetsMCPtArea[fCentBin]->Fill(jet->MCPt(), jet->Area(), ptLeading);
-      fHistJetPtvsJetMCPt[fCentBin]->Fill(jet->Pt(), jet->MCPt());
-    }
-
-    fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt(), ptLeading);
-    fHistJetsCEFvsCEFPt[fCentBin]->Fill(1-jet->NEF(), (1-jet->NEF())*jet->Pt(), ptLeading);
+    // Fill THnSparse
 
     if (fTracks) {
       for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
        AliVParticle *track = jet->TrackAt(it, fTracks);
        if (track) {
-         fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt(), ptLeading);
-         constituents.Fill(track->Pt());
          fHistTracksJetPt[fCentBin]->Fill(track->Pt(), jet->Pt());
          Double_t dist = TMath::Sqrt((track->Eta() - jet->Eta()) * (track->Eta() - jet->Eta()) + (track->Phi() - jet->Phi()) * (track->Phi() - jet->Phi()));
          fHistTracksPtDist[fCentBin]->Fill(track->Pt(), dist);
@@ -460,22 +321,12 @@ Int_t AliAnalysisTaskSAJF::DoJetLoop()
        if (cluster) {
          TLorentzVector nPart;
          cluster->GetMomentum(nPart, fVertex);
-         fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt(), ptLeading);
-         constituents.Fill(nPart.Pt());
+
          fHistClustersJetPt[fCentBin]->Fill(nPart.Pt(), jet->Pt());
          Double_t dist = TMath::Sqrt((nPart.Eta() - jet->Eta()) * (nPart.Eta() - jet->Eta()) + (nPart.Phi() - jet->Phi()) * (nPart.Phi() - jet->Phi()));
          fHistClustersPtDist[fCentBin]->Fill(nPart.Pt(), dist);
        }
       }
     }
-
-    for (Int_t i = 1; i <= constituents.GetNbinsX(); i++) {
-      fHistConstituents[fCentBin]->Fill(constituents.GetBinCenter(i), constituents.GetBinContent(i));
-    }
-
-    constituents.Reset();
-    nAccJets++;
   } //jet loop 
-
-  return nAccJets;
 }
index 5a01308..e68609a 100644 (file)
@@ -3,13 +3,11 @@
 
 // $Id$
 
-class TClonesArray;
-class TString;
-class TH1F;
-class TH2F;
-class TH3F;
-class AliRhoParameter;
 #include <TH3F.h>
+
+class TH2;
+class THnSparse;
+
 #include "AliAnalysisTaskEmcalJet.h"
 
 class AliAnalysisTaskSAJF : public AliAnalysisTaskEmcalJet {
@@ -17,38 +15,25 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskEmcalJet {
 
   AliAnalysisTaskSAJF();
   AliAnalysisTaskSAJF(const char *name);
-  virtual ~AliAnalysisTaskSAJF();
+  virtual ~AliAnalysisTaskSAJF() {;}
 
   void                        UserCreateOutputObjects();
 
  protected:
   Bool_t                      FillHistograms()                                              ;
-  Int_t                       DoJetLoop()                                                   ;
+  void                        DoJetLoop(const Int_t nacc, const Int_t* sortedJets)          ;
 
   // General histograms
-  TH3F                       *fHistLeadingJetPhiEta[4];    //!Leading jet phi-eta
-  TH3F                       *fHistLeadingJetPtArea[4];    //!Leading jet pt spectrum vs. area
-  TH3F                       *fHistLeadingJetCorrPtArea[4];//!Corrected leading jet pt spectrum vs. area
-  TH3F                       *fHistLeadingJetMCPtArea[4];  //!Leading jet MC pt spectrum vs. area
-  TH2F                       *fHistRhoVSleadJetPt[4];      //!Area(leadjet) * rho vs. leading jet pt
-  TH2F                       *fNjetsVsCent;                //!No. of jets vs. centrality
+  THnSparse                  *fHistEventObservables;       //!Event-wise observables (centrality, event plane, rho, number of jets, leading jet pt, leading jet corr pt, 
+                                                           // leading jet area)
 
   // Inclusive jets histograms
-  TH3F                       *fHistJetPhiEta[4];           //!Phi-Eta distribution of jets
-  TH3F                       *fHistJetsPtArea[4];          //!Jet pt vs. area
-  TH3F                       *fHistJetsCorrPtArea[4];      //!Jet corr pt vs. area
-  TH2F                       *fHistJetPtvsJetCorrPt[4];    //!Jet pt vs jet corr pt
-  TH3F                       *fHistJetsMCPtArea[4];        //!Jet pt vs. area
-  TH2F                       *fHistJetPtvsJetMCPt[4];      //!Jet pt vs jet corr pt
-  TH3F                       *fHistJetsNEFvsPt[4];         //!Jet neutral energy fraction vs. jet pt
-  TH3F                       *fHistJetsCEFvsCEFPt[4];      //!Jet charged energy fraction vs. charged jet pt
-  TH3F                       *fHistJetsZvsPt[4];           //!Constituent Pt over Jet Pt ratio vs. jet pt
-  TH2F                       *fHistConstituents[4];        //!x axis = constituents pt; y axis = no. of constituents
-  TH2F                       *fHistTracksJetPt[4];         //!Track pt vs. jet pt
-  TH2F                       *fHistClustersJetPt[4];       //!Cluster pt vs. jet pt
-  TH2F                       *fHistTracksPtDist[4];        //!Track pt vs. distance form jet axis
-  TH2F                       *fHistClustersPtDist[4];      //!Cluster pt vs. distance form jet axis
-  TH3F                       *fHistJetNconstVsPt[4];       //!Jet no. of constituents vs. pt
+  THnSparse                  *fHistJetObservables;         //!Jet-wise observables (centrality, event plane, phi, eta, pt, MCPt, area, corr pt, NEF, Z, no. of constituents
+                                                           // leading hadron pt)
+  TH2                        *fHistTracksJetPt[4];         //!Track pt vs. jet pt
+  TH2                        *fHistClustersJetPt[4];       //!Cluster pt vs. jet pt
+  TH2                        *fHistTracksPtDist[4];        //!Track pt vs. distance form jet axis
+  TH2                        *fHistClustersPtDist[4];      //!Cluster pt vs. distance form jet axis
 
  private:
   AliAnalysisTaskSAJF(const AliAnalysisTaskSAJF&);            // not implemented
index d25d2b6..f599a9e 100644 (file)
@@ -8,6 +8,7 @@
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TH3F.h>
+#include <THnSparse.h>
 #include <TList.h>
 #include <TLorentzVector.h>
 
@@ -25,6 +26,8 @@
 #include "AliEMCALGeometry.h"
 #include "AliEMCALGeoParams.h"
 #include "AliPicoTrack.h"
+#include "AliVVZERO.h"
+#include "AliESDUtils.h"
 
 #include "AliAnalysisTaskSAQA.h"
 
@@ -36,16 +39,15 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
   fCellEnergyCut(0.1),
   fParticleLevel(kFALSE),
   fIsMC(kFALSE),
-  fNclusters(0),
-  fNtracks(0),
-  fNjets(0),
-  fHistTracksCent(0),
-  fHistClusCent(0),
-  fHistJetsCent(0),
-  fHistClusTracks(0),
-  fHistJetsParts(0),
-  fHistCellsCent(0),
-  fHistCellsTracks(0),
+  fCentMethod2(""),
+  fCentMethod3(""),
+  fDoV0QA(0),
+  fCent2(0),
+  fCent3(0),
+  fVZERO(0),
+  fV0ATotMult(0),
+  fV0CTotMult(0),
+  fHistEventQA(0),
   fHistTrNegativeLabels(0),
   fHistTrZeroLabels(0),
   fHistNCellsEnergy(0),
@@ -60,16 +62,18 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
 
   for (Int_t i = 0; i < 4; i++) {
     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+    fHistTrPhiEtaZeroLab[i] = 0;
+    fHistTrPtZeroLab[i] = 0;
     fHistTrEmcPhiEta[i] = 0;
     fHistTrEmcPt[i] = 0;
     fHistTrPhiEtaNonProp[i] = 0;
+    fHistTrPtNonProp[i] = 0;
     fHistDeltaEtaPt[i] = 0;
     fHistDeltaPhiPt[i] = 0;
     fHistDeltaPtvsPt[i] = 0;
-    fHistTrPhiEtaPtZeroLab[i] = 0;
     fHistClusPhiEtaEnergy[i] = 0;
     fHistClusMCEnergyFraction[i] = 0;
-    fHistJetsPhiEtaPt[i] = 0;
+    fHistJetsPhiEta[i] = 0;
     fHistJetsPtArea[i] = 0;
   }
 
@@ -82,16 +86,15 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
   fCellEnergyCut(0.1),
   fParticleLevel(kFALSE),
   fIsMC(kFALSE),
-  fNclusters(0),
-  fNtracks(0),
-  fNjets(0),
-  fHistTracksCent(0),
-  fHistClusCent(0),
-  fHistJetsCent(0),
-  fHistClusTracks(0),
-  fHistJetsParts(0),
-  fHistCellsCent(0),
-  fHistCellsTracks(0),
+  fCentMethod2(""),
+  fCentMethod3(""),
+  fDoV0QA(0),
+  fCent2(0),
+  fCent3(0),
+  fVZERO(0),
+  fV0ATotMult(0),
+  fV0CTotMult(0),
+  fHistEventQA(0),
   fHistTrNegativeLabels(0),
   fHistTrZeroLabels(0),
   fHistNCellsEnergy(0),
@@ -106,16 +109,18 @@ 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;
+    fHistTrPhiEtaZeroLab[i] = 0;
+    fHistTrPtZeroLab[i] = 0;
     fHistTrEmcPhiEta[i] = 0;
     fHistTrEmcPt[i] = 0;
     fHistTrPhiEtaNonProp[i] = 0;
+    fHistTrPtNonProp[i] = 0;
     fHistDeltaEtaPt[i] = 0;
     fHistDeltaPhiPt[i] = 0;
     fHistDeltaPtvsPt[i] = 0;
-    fHistTrPhiEtaPtZeroLab[i] = 0;
     fHistClusPhiEtaEnergy[i] = 0;
     fHistClusMCEnergyFraction[i] = 0;
-    fHistJetsPhiEtaPt[i] = 0;
+    fHistJetsPhiEta[i] = 0;
     fHistJetsPtArea[i] = 0;
   }
 
@@ -148,11 +153,6 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
       fOutput->Add(fHistTrZeroLabels);
     }
 
-    fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
-    fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
-    fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
-    fOutput->Add(fHistTracksCent);
-
     TString histname;
 
     Int_t nlabels = 4;
@@ -171,12 +171,18 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 
       if (!fParticleLevel) {
        if (fIsMC) {
-         histname = Form("fHistTrPhiEtaPtZeroLab_%d",i);
-         fHistTrPhiEtaPtZeroLab[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
-         fHistTrPhiEtaPtZeroLab[i]->GetXaxis()->SetTitle("#eta");
-         fHistTrPhiEtaPtZeroLab[i]->GetYaxis()->SetTitle("#phi");
-         fHistTrPhiEtaPtZeroLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
-         fOutput->Add(fHistTrPhiEtaPtZeroLab[i]);
+         histname = Form("fHistTrPhiEtaZeroLab_%d",i);
+         fHistTrPhiEtaZeroLab[i] = new TH2F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+         fHistTrPhiEtaZeroLab[i]->GetXaxis()->SetTitle("#eta");
+         fHistTrPhiEtaZeroLab[i]->GetYaxis()->SetTitle("#phi");
+         fHistTrPhiEtaZeroLab[i]->GetZaxis()->SetTitle("counts");
+         fOutput->Add(fHistTrPhiEtaZeroLab[i]);
+
+         histname = Form("fHistTrPtZeroLab_%d",i);
+         fHistTrPtZeroLab[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
+         fHistTrPtZeroLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+         fHistTrPtZeroLab[i]->GetYaxis()->SetTitle("counts");
+         fOutput->Add(fHistTrPtZeroLab[i]);
        }
        
        histname = Form("fHistTrEmcPhiEta_%d",i);
@@ -196,6 +202,12 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
        fHistTrPhiEtaNonProp[i]->GetXaxis()->SetTitle("#eta");
        fHistTrPhiEtaNonProp[i]->GetYaxis()->SetTitle("#phi");
        fOutput->Add(fHistTrPhiEtaNonProp[i]);
+
+       histname = Form("fHistTrPtNonProp_%d",i);
+       fHistTrPtNonProp[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
+       fHistTrPtNonProp[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+       fHistTrPtNonProp[i]->GetYaxis()->SetTitle("counts");
+       fOutput->Add(fHistTrPtNonProp[i]);
        
        histname = Form("fHistDeltaEtaPt_%d",i);
        fHistDeltaEtaPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
@@ -219,26 +231,6 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
   }
 
   if (!fCaloName.IsNull()) {
-    fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
-    fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
-    fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
-    fOutput->Add(fHistClusCent);
-
-    fHistClusTracks = new TH2F("fHistClusTracks","Clusters vs. tracks", fNbins, 0, 4000, fNbins, 0, 2000);
-    fHistClusTracks->GetXaxis()->SetTitle("No. of tracks");
-    fHistClusTracks->GetYaxis()->SetTitle("No. of clusters");
-    fOutput->Add(fHistClusTracks);
-
-    fHistCellsCent = new TH2F("fHistCellsCent","Cells vs. centrality", 100, 0, 100, fNbins, 0, 6000);
-    fHistCellsCent->GetXaxis()->SetTitle("Centrality (%)");
-    fHistCellsCent->GetYaxis()->SetTitle("No. of EMCal cells");
-    fOutput->Add(fHistCellsCent);
-
-    fHistCellsTracks = new TH2F("fHistCellsTracks","Cells vs. tracks", fNbins, 0, 4000, fNbins, 0, 6000);
-    fHistCellsTracks->GetXaxis()->SetTitle("No. of tracks");
-    fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
-    fOutput->Add(fHistCellsTracks);
-
     TString histname;
 
     for (Int_t i = 0; i < fNcentBins; i++) {
@@ -301,26 +293,17 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
   }
        
   if (!fJetsName.IsNull()) {
-    fHistJetsCent = new TH2F("fHistJetsCent","Jets vs. centrality", 100, 0, 100, 150, -0.5, 149.5);
-    fHistJetsCent->GetXaxis()->SetTitle("Centrality (%)");
-    fHistJetsCent->GetYaxis()->SetTitle("No. of jets");
-    fOutput->Add(fHistJetsCent);
-    
-    fHistJetsParts = new TH2F("fHistJetsParts","Jets vs. centrality", fNbins, 0, 6000, 150, -0.5, 149.5);
-    fHistJetsParts->GetXaxis()->SetTitle("No. of particles");
-    fHistJetsParts->GetYaxis()->SetTitle("No. of jets");
-    fOutput->Add(fHistJetsParts);
 
     TString histname;
 
     for (Int_t i = 0; i < fNcentBins; i++) {
-      histname = "fHistJetsPhiEtaPt_";
+      histname = "fHistJetsPhiEta_";
       histname += i;
-      fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01,  (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
-      fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
-      fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
-      fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
-      fOutput->Add(fHistJetsPhiEtaPt[i]);
+      fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01);
+      fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
+      fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
+      fHistJetsPhiEta[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistJetsPhiEta[i]);
 
       histname = "fHistJetsPtArea_";
       histname += i;
@@ -330,11 +313,113 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
       fOutput->Add(fHistJetsPtArea[i]);
     }
   }
+  
+  Int_t dim = 0;
+  TString title[10];
+  Int_t nbins[10] = {0};
+  Double_t min[10] = {0};
+  Double_t max[10] = {0};
+  
+  if (fForceBeamType != AliAnalysisTaskEmcal::kpp) {
+    title[dim] = "Centrality %";
+    nbins[dim] = 101;
+    min[dim] = 0;
+    max[dim] = 101;
+    dim++;
+
+    if (!fCentMethod2.IsNull()) {
+      title[dim] = Form("Centrality %s %%", fCentMethod2.Data());
+      nbins[dim] = 101;
+      min[dim] = 0;
+      max[dim] = 101;
+      dim++;
+    }
+
+    if (!fCentMethod3.IsNull()) {
+      title[dim] = Form("Centrality %s %%", fCentMethod3.Data());
+      nbins[dim] = 101;
+      min[dim] = 0;
+      max[dim] = 101;
+      dim++;
+    }
+    
+    if (fDoV0QA) {
+      title[dim] = "V0A total multiplicity";
+      nbins[dim] = 200;
+      min[dim] = 0;
+      max[dim] = 20000;
+      dim++;
+
+      title[dim] = "V0C total multiplicity";
+      nbins[dim] = 200;
+      min[dim] = 0;
+      max[dim] = 20000;
+      dim++;
+    }
+
+    if (!fRhoName.IsNull()) {
+      title[dim] = "#rho (GeV/c)";
+      nbins[dim] = fNbins*4;
+      min[dim] = 0;
+      max[dim] = 400;
+      dim++;
+    }
+  }
+
+  if (!fTracksName.IsNull()) {
+    title[dim] = "No. of tracks";
+    nbins[dim] = 6000;
+    min[dim] = -0.5;
+    max[dim] = 6000-0.5;
+    dim++;
+  }
+
+  if (!fCaloName.IsNull()) {
+    title[dim] = "No. of clusters";
+    nbins[dim] = 4000;
+    min[dim] = 0;
+    max[dim] = 4000-0.5;
+    dim++;
+  }
+
+  if (!fCaloCellsName.IsNull()) {
+    title[dim] = "No. of cells";
+    nbins[dim] = 6000;
+    min[dim] = 0;
+    max[dim] = 6000-0.5;
+    dim++;
+  }
+
+  if (!fJetsName.IsNull()) {
+    title[dim] = "No. of jets";
+    nbins[dim] = 200;
+    min[dim] = 0;
+    max[dim] = 200-0.5;
+    dim++;
+  }
+
+  fHistEventQA = new THnSparseF("fHistEventQA","fHistEventQA",dim,nbins,min,max);
+  for (Int_t i = 0; i < dim; i++)
+    fHistEventQA->GetAxis(i)->SetTitle(title[i]);
+  fOutput->Add(fHistEventQA);
 
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 }
 
 //________________________________________________________________________
+void AliAnalysisTaskSAQA::ExecOnce()
+{
+  AliAnalysisTaskEmcalJet::ExecOnce();
+  
+  if (fDoV0QA) {
+    fVZERO = InputEvent()->GetVZEROData();
+    if (!fVZERO) {
+      AliError("AliVVZERO not available");
+    }
+  }
+}
+
+//________________________________________________________________________
 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
 {
   // Retrieve event objects.
@@ -342,9 +427,22 @@ Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
     return kFALSE;
 
-  fNclusters = 0;
-  fNtracks = 0;
-  fNjets = 0;
+  if (!fCentMethod2.IsNull() || !fCentMethod3.IsNull()) {
+    if (fBeamType == kAA || fBeamType == kpA ) {
+      AliCentrality *aliCent = InputEvent()->GetCentrality();
+      if (aliCent) {
+       if (!fCentMethod2.IsNull()) 
+         fCent2 = aliCent->GetCentralityPercentile(fCentMethod2); 
+       if (!fCentMethod3.IsNull()) 
+         fCent3 = aliCent->GetCentralityPercentile(fCentMethod3);
+      }
+    }
+  }
+
+  if (fVZERO) {
+    fV0ATotMult = AliESDUtils::GetCorrV0A(fVZERO->GetMTotV0A(),fVertex[2]);
+    fV0CTotMult = AliESDUtils::GetCorrV0C(fVZERO->GetMTotV0C(),fVertex[2]);
+  }
 
   return kTRUE;
 }
@@ -355,46 +453,81 @@ Bool_t AliAnalysisTaskSAQA::FillHistograms()
 {
   // Fill histograms.
 
-  Float_t clusSum = 0;
   Float_t trackSum = 0;
+  Float_t clusSum = 0;
+  Float_t cellSum = 0;
+  Float_t cellCutSum = 0;
 
+  Int_t ntracks = 0;
+  Int_t nclusters = 0;
+  Int_t ncells = 0;
+  Int_t njets = 0;
+    
   if (fTracks) {
-    trackSum = DoTrackLoop();
-    AliDebug(2,Form("%d tracks found in the event", fTracks->GetEntriesFast()));
-    fHistTracksCent->Fill(fCent, fNtracks);
+    ntracks = DoTrackLoop(trackSum);
+    AliDebug(2,Form("%d tracks found in the event", ntracks));
   } 
 
   if (fCaloClusters) {
-    clusSum = DoClusterLoop();
+    nclusters = DoClusterLoop(clusSum);
+    AliDebug(2,Form("%d clusters found in the event", nclusters));
 
-    fHistClusCent->Fill(fCent, fNclusters);
-    fHistClusTracks->Fill(fNtracks, fNclusters);
-
-    Float_t cellSum = 0, cellCutSum = 0;
-    
-    Int_t ncells = DoCellLoop(cellSum, cellCutSum);
-
-    if (fTracks)
-      fHistCellsTracks->Fill(fNtracks, ncells);
-
-    fHistCellsCent->Fill(fCent, ncells);
+    fHistChVSneClus->Fill(clusSum, trackSum);
+  }
+  
+  if (fCaloCells) {
+    ncells = DoCellLoop(cellSum, cellCutSum);
+    AliDebug(2,Form("%d cells found in the event", ncells));
     
     fHistChVSneCells->Fill(cellSum, trackSum);
-    fHistChVSneClus->Fill(clusSum, trackSum);
     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
   }
 
   if (fJets) {
-    DoJetLoop();
-    
-    fHistJetsCent->Fill(fCent, fNjets);
-    fHistJetsParts->Fill(fNtracks + fNclusters, fNjets);
+    njets = DoJetLoop();
+    AliDebug(2,Form("%d jets found in the event", njets));
   }
 
+  FillEventQAHisto(fCent, fCent2, fCent3, fV0ATotMult, fV0CTotMult, fRhoVal, ntracks, nclusters, ncells, njets);
+
   return kTRUE;
 }
 
 //________________________________________________________________________
+void AliAnalysisTaskSAQA::FillEventQAHisto(Float_t cent, Float_t cent2, Float_t cent3, Float_t v0a, Float_t v0c, Float_t rho, Int_t ntracks, Int_t nclusters, Int_t ncells, Int_t njets)
+{
+  Double_t contents[10]={0};
+
+  for (Int_t i = 0; i < fHistEventQA->GetNdimensions(); i++) {
+    TString title(fHistEventQA->GetAxis(i)->GetTitle());
+    if (title=="Centrality %")
+      contents[i] = cent;
+    else if (title==Form("Centrality %s %%", fCentMethod2.Data()))
+      contents[i] = cent2;
+    else if (title==Form("Centrality %s %%", fCentMethod3.Data()))
+      contents[i] = cent3;
+    else if (title=="V0A total multiplicity")
+      contents[i] = v0a;
+    else if (title=="V0C total multiplicity")
+      contents[i] = v0c;
+    else if (title=="#rho (GeV/c)")
+      contents[i] = rho;
+    else if (title=="No. of tracks")
+       contents[i] = ntracks;
+    else if (title=="No. of clusters")
+      contents[i] = nclusters;
+    else if (title=="No. of cells")
+      contents[i] = ncells;
+    else if (title=="No. of jets")
+      contents[i] = njets;
+    else 
+      AliWarning(Form("Unable to fill dimension %s!",title.Data()));
+  }
+
+  fHistEventQA->Fill(contents);
+}
+
+//________________________________________________________________________
 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 {
   // Do cell loop.
@@ -479,19 +612,21 @@ Double_t AliAnalysisTaskSAQA::GetFcross(AliVCluster *cluster, AliVCaloCells *cel
 }
 
 //________________________________________________________________________
-Float_t AliAnalysisTaskSAQA::DoClusterLoop()
+Int_t AliAnalysisTaskSAQA::DoClusterLoop(Float_t &sum)
 {
   // Do cluster loop.
 
   if (!fCaloClusters)
     return 0;
 
+  Int_t nAccClusters = 0;
+
   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 
-  Float_t sum = 0;
+  sum = 0;
 
   // Cluster loop
-  Int_t nclusters =  fCaloClusters->GetEntriesFast();
+  const Int_t nclusters = fCaloClusters->GetEntriesFast();
 
   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
@@ -519,23 +654,25 @@ Float_t AliAnalysisTaskSAQA::DoClusterLoop()
     if (fHistClusMCEnergyFraction[fCentBin])
       fHistClusMCEnergyFraction[fCentBin]->Fill(cluster->GetMCEnergyFraction());
 
-    fNclusters++;
+    nAccClusters++;
   }
 
-  return sum;
+  return nAccClusters;
 }
 
 //________________________________________________________________________
-Float_t AliAnalysisTaskSAQA::DoTrackLoop()
+Int_t AliAnalysisTaskSAQA::DoTrackLoop(Float_t &sum)
 {
   // Do track loop.
 
   if (!fTracks)
     return 0;
 
-  Float_t sum = 0;
+  Int_t nAccTracks = 0;
+
+  sum = 0;
 
-  Int_t ntracks = fTracks->GetEntriesFast();
+  const Int_t ntracks = fTracks->GetEntriesFast();
   Int_t neg = 0;
   Int_t zero = 0;
 
@@ -551,7 +688,7 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
     if (!AcceptTrack(track)) 
       continue;
 
-    fNtracks++;
+    nAccTracks++;
 
     sum += track->P();
 
@@ -562,8 +699,10 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
       fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
       if (track->GetLabel() == 0) {
        zero++;
-       if (fHistTrPhiEtaPtZeroLab[fCentBin])
-         fHistTrPhiEtaPtZeroLab[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
+       if (fHistTrPhiEtaZeroLab[fCentBin]) {
+         fHistTrPhiEtaZeroLab[fCentBin]->Fill(track->Eta(), track->Phi());
+         fHistTrPtZeroLab[fCentBin]->Fill(track->Pt());
+       }
       }
 
       if (track->GetLabel() < 0)
@@ -586,8 +725,10 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
     if (!vtrack)
       continue;
 
-    if ((vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999) && fHistTrPhiEtaNonProp[fCentBin])
+    if ((vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999) && fHistTrPhiEtaNonProp[fCentBin]) {
       fHistTrPhiEtaNonProp[fCentBin]->Fill(vtrack->Eta(), vtrack->Phi());
+      fHistTrPtNonProp[fCentBin]->Fill(vtrack->Pt());
+    }
 
     if (fHistTrEmcPhiEta[fCentBin])
       fHistTrEmcPhiEta[fCentBin]->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());   
@@ -607,16 +748,18 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
   if (fHistTrZeroLabels)
     fHistTrZeroLabels->Fill(1. * zero / ntracks);
 
-  return sum;
+  return nAccTracks;
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskSAQA::DoJetLoop()
+Int_t AliAnalysisTaskSAQA::DoJetLoop()
 {
   // Do jet loop.
 
   if (!fJets)
-    return;
+    return 0;
+
+  Int_t nAccJets = 0;
 
   Int_t njets = fJets->GetEntriesFast();
 
@@ -632,16 +775,11 @@ void AliAnalysisTaskSAQA::DoJetLoop()
     if (!AcceptJet(jet))
       continue;
 
-    fNjets++;
+    nAccJets++;
 
-    fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
+    fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
   }
-}
-
 
-//________________________________________________________________________
-void AliAnalysisTaskSAQA::Terminate(Option_t *) 
-{
-  // Called once at the end of the analysis.
+  return nAccJets;
 }
index 0ca0ef1..7cf92ee 100644 (file)
@@ -6,6 +6,8 @@
 class TH1;
 class TH2;
 class TH3;
+class THnSparse;
+class AliVVZERO;
 
 #include "AliAnalysisTaskEmcalJet.h"
 
@@ -16,46 +18,50 @@ class AliAnalysisTaskSAQA : public AliAnalysisTaskEmcalJet {
   virtual ~AliAnalysisTaskSAQA();
 
   void                        UserCreateOutputObjects();
-  void                        Terminate(Option_t *option);
 
   void                        SetCellEnergyCut(Float_t cut)                        { fCellEnergyCut      = cut        ; }
   void                        SetParticleLevel(Bool_t s)                           { fParticleLevel      = s          ; }
   void                        SetMC(Bool_t m)                                      { fIsMC               = m          ; }
+  void                        SetAdditionalCentEst(const char* meth2, const char* meth3="") { fCentMethod2 = meth2; fCentMethod3 = meth3; }
+  void                        SetDoV0QA(Int_t b)                                   { fDoV0QA             = b          ; }
 
  protected:
 
+  void                        ExecOnce()                                                    ;
   Bool_t                      FillHistograms()                                              ;
+  void                        FillEventQAHisto(Float_t cent, Float_t cent2, Float_t cent3, Float_t v0a, Float_t v0c, Float_t rho, Int_t ntracks, Int_t nclusters, Int_t ncells, Int_t njets);
   Bool_t                      RetrieveEventObjects()                                        ;
   Int_t                       DoCellLoop(Float_t &sum, Float_t &sum_cut)                    ;
-  Float_t                     DoTrackLoop()                                                 ;
-  Float_t                     DoClusterLoop()                                               ;
-  void                        DoJetLoop()                                                   ;
+  Int_t                       DoTrackLoop(Float_t &sum)                                     ;
+  Int_t                       DoClusterLoop(Float_t &sum)                                   ;
+  Int_t                       DoJetLoop()                                                   ;
   Double_t                    GetFcross(AliVCluster *cluster, AliVCaloCells *cells)         ;
 
   Float_t                     fCellEnergyCut;            // Energy cell cut
   Bool_t                      fParticleLevel;            // Set particle level analysis
   Bool_t                      fIsMC;                     // Trigger, MC analysis
-  Int_t                       fNclusters;                //!Number of accepted clusters in the event
-  Int_t                       fNtracks;                  //!Number of accepted tracks in the event
-  Int_t                       fNjets;                    //!Number of accepted jets in the event
+  TString                     fCentMethod2;              // Centrality method 2
+  TString                     fCentMethod3;              // Centrality method 3
+  Int_t                       fDoV0QA;                   // Add V0 QA histograms
+  Double_t                    fCent2;                    //!Event centrality with method 2
+  Double_t                    fCent3;                    //!Event centrality with method 3
+  AliVVZERO                  *fVZERO;                    //!Event V0 object
+  Double_t                    fV0ATotMult;               //!Event V0A total multiplicity
+  Double_t                    fV0CTotMult;               //!Event V0C total multiplicity
  
   // General histograms
-  TH2                        *fHistTracksCent;           //!Number of tracks vs. centrality
-  TH2                        *fHistClusCent;             //!Number of clusters vs. centrality
-  TH2                        *fHistJetsCent;             //!Number of jets vs. centrality
-  TH2                        *fHistClusTracks;           //!Number of clusters vs. number of tracks
-  TH2                        *fHistJetsParts;            //!Number of jets vs. number of particles (tracks+clusters)
-  TH2                        *fHistCellsCent;            //!Number of cells vs. centrality
-  TH2                        *fHistCellsTracks;          //!Number of cells vs. number of tracks
+  THnSparse                  *fHistEventQA;              //!Event-wise QA observables (cent[M], cent[A], cent[C], rho, # tracks, # clusters, # cells, # jets)
 
   // Tracks
   TH1                        *fHistTrNegativeLabels;     //!Percentage of negative label tracks
   TH1                        *fHistTrZeroLabels;         //!Percentage of tracks with label=0
   TH3                        *fHistTrPhiEtaPt[4][4];     //!Phi-Eta-Pt distribution of tracks
-  TH3                        *fHistTrPhiEtaPtZeroLab[4]; //!Phi-Eta-Pt distribution of tracks with label=0
+  TH2                        *fHistTrPhiEtaZeroLab[4];   //!Phi-Eta distribution of tracks with label=0
+  TH1                        *fHistTrPtZeroLab[4];       //!Pt distribution of tracks with label=0
   TH2                        *fHistTrEmcPhiEta[4];       //!Phi-Eta emcal propagated distribution of tracks
-  TH1                        *fHistTrEmcPt[4];           //!Phi-Eta emcal propagated distribution of tracks
+  TH1                        *fHistTrEmcPt[4];           //!Pt emcal propagated distribution of tracks
   TH2                        *fHistTrPhiEtaNonProp[4];   //!Phi-Eta distribution of non emcal propagated tracks
+  TH1                        *fHistTrPtNonProp[4];       //!Pt distribution of non emcal propagated tracks
   TH2                        *fHistDeltaEtaPt[4];        //!Eta-EtaProp vs. Pt
   TH2                        *fHistDeltaPhiPt[4];        //!Phi-PhiProp vs. Pt
   TH2                        *fHistDeltaPtvsPt[4];       //!Pt-PtProp vs. Pt
@@ -68,7 +74,7 @@ class AliAnalysisTaskSAQA : public AliAnalysisTaskEmcalJet {
   TH1                        *fHistClusMCEnergyFraction[4];//!MC energy fraction (embedding)
 
   // Jets
-  TH3                        *fHistJetsPhiEtaPt[4];      //!Phi-Eta distribution of jets
+  TH2                        *fHistJetsPhiEta[4];        //!Phi-Eta distribution of jets
   TH2                        *fHistJetsPtArea[4];        //!Pt vs. area of jets
 
   // EMCAL Cells
@@ -83,6 +89,6 @@ class AliAnalysisTaskSAQA : public AliAnalysisTaskEmcalJet {
   AliAnalysisTaskSAQA(const AliAnalysisTaskSAQA&);            // not implemented
   AliAnalysisTaskSAQA &operator=(const AliAnalysisTaskSAQA&); // not implemented
 
-  ClassDef(AliAnalysisTaskSAQA, 18) // Quality task for Emcal analysis
+  ClassDef(AliAnalysisTaskSAQA, 19) // Quality task for Emcal analysis
 };
 #endif
index 82238bf..f32971c 100644 (file)
@@ -3,7 +3,9 @@
 AliAnalysisTaskSAQA* AddTaskSAQA(
   const char *ntracks            = "Tracks",
   const char *nclusters          = "CaloClusters",
+  const char *ncells             = "EMCALCells",
   const char *njets              = "Jets",
+  const char *nrho               = "Rho",
   Double_t    jetradius          = 0.2,
   Double_t    jetptcut           = 1,
   Double_t    jetareacut         = 0.557,
@@ -57,7 +59,9 @@ AliAnalysisTaskSAQA* AddTaskSAQA(
   AliAnalysisTaskSAQA* qaTask = new AliAnalysisTaskSAQA(name);
   qaTask->SetTracksName(ntracks);
   qaTask->SetClusName(nclusters);
+  qaTask->SetCaloCellsName(ncells);
   qaTask->SetJetsName(njets);
+  qaTask->SetRhoName(nrho);
   qaTask->SetJetRadius(jetradius);
   qaTask->SetJetPtCut(jetptcut);
   qaTask->SetPercAreaCut(jetareacut);