Introduction of PartonInfo object to trace the original hard partons +
authordcaffarr <davide.caffarri@cern.ch>
Wed, 8 Oct 2014 16:58:31 +0000 (18:58 +0200)
committermverweij <marta.verweij@cern.ch>
Thu, 9 Oct 2014 08:40:37 +0000 (10:40 +0200)
analysis classed for quark-gluon tagging from Davide and Leticia
testing the patch

15 files changed:
PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJet.cxx
PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
PWGJE/EMCALJetTasks/AliJetContainer.cxx
PWGJE/EMCALJetTasks/AliJetContainer.h
PWGJE/EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx
PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
PWGJE/EMCALJetTasks/AliJetModelBaseTask.h
PWGJE/EMCALJetTasks/AliStackPartonInfo.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliStackPartonInfo.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/macros/AddTaskEmcalQGTagging.C [new file with mode: 0644]
PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromGen.C
PWGJE/PWGJEEMCALJetTasksLinkDef.h

index dd6591d..8ac6d64 100644 (file)
@@ -51,6 +51,7 @@ set ( SRCS
  EMCALJetTasks/AliRhoParameter.cxx
  EMCALJetTasks/AliLocalRhoParameter.cxx
  EMCALJetTasks/AliJetTriggerSelectionTask.cxx
+ EMCALJetTasks/AliStackPartonInfo.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskCLQA.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskDcalDijetPerf.cxx
@@ -82,6 +83,7 @@ set ( SRCS
  EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMassBase.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMassScale.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskSOH.cxx
index ca2b25f..14a56f5 100644 (file)
@@ -175,7 +175,8 @@ void AliAnalysisTaskEmcalJet::ExecOnce()
     cont->SetRunNumber(InputEvent()->GetRunNumber());
     cont->SetEMCALGeometry();
     cont->SetArray(InputEvent());
-    cont->LoadRho(InputEvent());
+    cont->LoadRho(InputEvent()); 
+    cont->LoadPartonsInfo(InputEvent());
   }
 
   //Get Jets, cuts and rho for first jet container
index 9efd53c..2d251ca 100644 (file)
@@ -622,7 +622,7 @@ void AliEmcalJetTask::FindJets()
        AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount]) 
          AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
        jet_sub->SetLabel(ijet);
-        // Fill constituent info
+        // Fill constituent info
        std::vector<fastjet::PseudoJet> constituents_sub(fjw.GetJetConstituents(ijet));
        jet_sub->SetNumberOfTracks(constituents_sub.size());
        jet_sub->SetNumberOfClusters(constituents_sub.size());
@@ -875,7 +875,7 @@ void  AliEmcalJetTask::FillJetConstituents(std::vector<fastjet::PseudoJet>& cons
       } else {
         AliError(Form("%s: No logical way to end up here.", GetName()));
         continue;
-      } 
+      }
     }
 }
 //______________________________________________________________________________________
index c753fd0..cb56c9f 100644 (file)
@@ -13,6 +13,7 @@
 #include "AliParticleContainer.h"
 #include "AliClusterContainer.h"
 #include "AliLocalRhoParameter.h"
+#include "AliStackPartonInfo.h"
 
 #include "AliJetContainer.h"
 
@@ -26,6 +27,7 @@ AliJetContainer::AliJetContainer():
   fRhoName(),
   fLocalRhoName(),
   fRhoMassName(),
+  fPartonInfoName(),
   fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
@@ -52,6 +54,7 @@ AliJetContainer::AliJetContainer():
   fRho(0),
   fLocalRho(0),
   fRhoMass(0),
+  fPartonsInfo(0),
   fGeom(0),
   fRunNumber(0)
 {
@@ -68,6 +71,7 @@ AliJetContainer::AliJetContainer(const char *name):
   fRhoName(),
   fLocalRhoName(),
   fRhoMassName(),
+  fPartonInfoName(),
   fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
@@ -94,6 +98,7 @@ AliJetContainer::AliJetContainer(const char *name):
   fRho(0),
   fLocalRho(0),
   fRhoMass(0),
+  fPartonsInfo(0),
   fGeom(0),
   fRunNumber(0)
 {
@@ -170,6 +175,19 @@ void AliJetContainer::LoadRhoMass(AliVEvent *event)
     }
   }
 }
+//________________________________________________________________________
+void AliJetContainer::LoadPartonsInfo(AliVEvent *event)
+{
+    // Load parton info
+    
+    if (!fPartonInfoName.IsNull() && !fPartonsInfo) {
+        fPartonsInfo = dynamic_cast<AliStackPartonInfo*>(event->FindListObject(fPartonInfoName));
+        if (!fPartonsInfo) {
+           AliError(Form("%s: Could not retrieve parton infos! %s!", GetName(), fPartonInfoName.Data()));            return;
+        }
+    }
+}
+
 
 
 //________________________________________________________________________
index ec7b790..62eb137 100644 (file)
@@ -12,6 +12,7 @@ class AliVEvent;
 class AliParticleContainer;
 class AliClusterContainer;
 class AliLocalRhoParameter;
+class AliStackPartonInfo;
 
 #include "AliRhoParameter.h"
 
@@ -36,6 +37,7 @@ class AliJetContainer : public AliEmcalContainer {
   void LoadRho(AliVEvent *event);
   void LoadLocalRho(AliVEvent *event);
   void LoadRhoMass(AliVEvent *event);
+  void LoadPartonsInfo(AliVEvent *event);
 
   void                        SetJetAcceptanceType(JetAcceptanceType type)         { fJetAcceptanceType          = type ; }
   void                        PrintCuts();
@@ -67,6 +69,8 @@ class AliJetContainer : public AliEmcalContainer {
   virtual void                SetRhoName(const char *n)                            { fRhoName        = n                ; }
   virtual void                SetLocalRhoName(const char *n)                       { fLocalRhoName   = n                ; }
   virtual void                SetRhoMassName(const char *n)                        { fRhoMassName    = n                ; }
+  virtual void                SetPartonInfoName(const char *n)                      { fPartonInfoName    = n; }
+    
   void                        ConnectParticleContainer(AliParticleContainer *c)    { fParticleContainer = c             ; }
   void                        ConnectClusterContainer(AliClusterContainer *c)      { fClusterContainer  = c             ; }
 
@@ -95,6 +99,8 @@ class AliJetContainer : public AliEmcalContainer {
   AliRhoParameter            *GetRhoMassParameter()                          {return fRhoMass;}
   Double_t                    GetRhoMassVal()                       const    {if (fRhoMass) return fRhoMass->GetVal(); else return 0;}
   const TString&              GetRhoMassName()                      const    {return fRhoMassName;}
+  const TString&              GetPartonInfoName()                   const    {return fPartonInfoName;}
+  AliStackPartonInfo         *GetStackPartonInfo()                  const    {return fPartonsInfo;}
   Double_t                    GetJetPtCorr(Int_t i)                 const;
   Double_t                    GetJetPtCorrLocal(Int_t i)            const;
   Float_t                     GetJetRadius()                        const    {return fJetRadius;}
@@ -107,6 +113,7 @@ class AliJetContainer : public AliEmcalContainer {
   AliParticleContainer       *GetParticleContainer()                         {return fParticleContainer;}
   AliClusterContainer        *GetClusterContainer()                          {return fClusterContainer;}
   Double_t                    GetFractionSharedPt(AliEmcalJet *jet) const;
 
  protected:
   JetAcceptanceType           fJetAcceptanceType;    //  acceptance type
@@ -114,6 +121,7 @@ class AliJetContainer : public AliEmcalContainer {
   TString                     fRhoName;              //  Name of rho object
   TString                     fLocalRhoName;         //  Name of local rho object
   TString                     fRhoMassName;          //  Name of rho mass object
+  TString                     fPartonInfoName;       //  Name of parton info object
   Int_t                       fFlavourSelection;     //  selection on jet flavour
   Float_t                     fPtBiasJetTrack;       //  select jets with a minimum pt track
   Float_t                     fPtBiasJetClus;        //  select jets with a minimum pt cluster
@@ -140,6 +148,7 @@ class AliJetContainer : public AliEmcalContainer {
   AliRhoParameter            *fRho;                  //! event rho for these jets
   AliLocalRhoParameter       *fLocalRho;             //! event local rho for these jets
   AliRhoParameter            *fRhoMass;              //! event rho mass for these jets
+  AliStackPartonInfo         *fPartonsInfo;          //! event parton info
   AliEMCALGeometry           *fGeom;                 //! emcal geometry
   Int_t                       fRunNumber;            //! run number
 
index beb9614..70e0624 100644 (file)
@@ -28,6 +28,7 @@
 #include "AliVCluster.h"
 #include "AliVEvent.h"
 #include "AliGenPythiaEventHeader.h"
+#include "AliStackPartonInfo.h"
 
 ClassImp(AliJetEmbeddingFromGenTask)
 
@@ -82,7 +83,7 @@ void AliJetEmbeddingFromGenTask::UserCreateOutputObjects()
   fHistPt = new TH1F("fHistpt","fHistPt;#it{p}_{T};N",100,0.,100.);
   fOutput->Add(fHistPt);
 
-  fHistEtaPhi = new TH2F("fHistEtaPhi","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
+  fHistEtaPhi = new TH2F("fHistEtapHI","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
   fOutput->Add(fHistEtaPhi);
 
   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
@@ -126,7 +127,12 @@ Bool_t AliJetEmbeddingFromGenTask::ExecOnce()
     InputEvent()->AddObject(fOutTracks);
     fNTracks = 0;
   }
-  
+
+  if (!(InputEvent()->FindListObject(fPartonInfoName))) {
+  fStackPartonInfo = new AliStackPartonInfo("PartonsInfo");
+  fStackPartonInfo->SetName(fPartonInfoName);
+  InputEvent()->AddObject(fStackPartonInfo);
+  }
   return kTRUE;
 }
 
@@ -142,6 +148,19 @@ void AliJetEmbeddingFromGenTask::Run()
   stack->Reset();
   fGen->Generate();
   const Int_t nprim = stack->GetNprimary();
+  TParticle *part6 = stack->Particle(6);
+  TParticle *part7 = stack->Particle(7);
+
+  fStackPartonInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
+  fStackPartonInfo->SetPartonPt6(part6->Pt());
+  fStackPartonInfo->SetPartonEta6(part6->Eta());
+  fStackPartonInfo->SetPartonPhi6(part6->Phi());
+  
+  fStackPartonInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
+  fStackPartonInfo->SetPartonPt7(part7->Pt());
+  fStackPartonInfo->SetPartonEta7(part7->Eta());
+  fStackPartonInfo->SetPartonPhi7(part7->Phi());
+  
   for (Int_t i=0;i<nprim;++i) {
     if (!stack->IsPhysicalPrimary(i))
       continue;
index abdb744..d73df93 100644 (file)
@@ -78,7 +78,8 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fOutMCParticlesMap(0),
   fMCLabelShift(0),
   fEsdMode(kFALSE),
-  fOutput(0)
+  fOutput(0),
+  fStackPartonInfo(0x0)
 {
   // Default constructor.
 
@@ -135,7 +136,8 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fOutMCParticlesMap(0),
   fMCLabelShift(0),
   fEsdMode(kFALSE),
-  fOutput(0)
+  fOutput(0),
+  fStackPartonInfo(0x0)
 {
   // Standard constructor.
 
@@ -218,9 +220,9 @@ void AliJetModelBaseTask::UserExec(Option_t *)
 
   if (fPtPhiEvPlDistribution || fAddV2)
     fPsi = gRandom->Rndm() * TMath::Pi();
-
+  
   Run();
-
+  
   if (fCaloCells && !fCopyArray) {
     delete fCaloCells;
     fCaloCells = tempCaloCells;
index 85fbf5a..31dea7c 100644 (file)
@@ -11,6 +11,7 @@ class AliVCaloCells;
 class AliAODMCParticle;
 class AliNamedArrayI;
 class TF2;
+class AliStackPartonInfo;
 
 #include <TH1F.h>
 #include <TF1.h>
@@ -26,11 +27,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                   SetPtSpectrum(TH1F *f)                { fPtSpectrum   = f;    }
+  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(TH1F *f)           { fDensitySpectrum = f;    }
+  void                   SetDensitySpectrum(TH1 *f)            { fDensitySpectrum = f;    }
   void                   SetDensitySpectrum(TF1 *f)            { fDensitySpectrum = new TH1F("densitypectrum","densitypectrum",1000,f->GetXmin(),f->GetXmax()); 
                                                                  fDensitySpectrum->Add(f); }
   void                   SetDifferentialV2(TF1* f)             { fDifferentialV2 = f;  }
@@ -42,6 +43,7 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   void                   SetClusName(const char *n)            { fCaloName     = n;    }
   void                   SetCellsName(const char *n)           { fCellsName    = n;    }
   void                   SetMCParticlesName(const char *n)     { fMCParticlesName = n; }
+    void                 SetPartonInfoName(const char *n)      { fPartonInfoName = n; }
   void                   SetSuffix(const char *s)              { fSuffix       = s;    }
   void                   SetGeometryName(const char *n)        { fGeomName     = n;    }
   void                   SetMarkMC(Int_t m)                    { fMarkMC       = m;    }
@@ -84,6 +86,7 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   TString                fOutCellsName;           // name of output cells collection
   TString                fMCParticlesName;        // name of MC particle collection
   TString                fOutMCParticlesName;     // name of output MC particle collection
+  TString                fPartonInfoName;         // name of partons info
   Bool_t                 fIsMC;                   // whether the current event is MC or not
   TString                fSuffix;                 // suffix to add in the name of new collections
   Float_t                fEtaMin;                 // eta minimum value
@@ -97,9 +100,9 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   Int_t                  fNCells;                 // how many cells are being processed
   Int_t                  fNTracks;                // how many tracks are being processed
   Int_t                  fMarkMC;                 // which MC label is to be used (default=100)
-  TH1F                  *fPtSpectrum;             // pt spectrum to extract random pt values
+  TH1                   *fPtSpectrum;             // pt spectrum to extract random pt values
   TF2                   *fPtPhiEvPlDistribution;  // pt vs. (phi-psi) distribution to extract random pt/phi values
-  TH1F                  *fDensitySpectrum;        // particle density spectrum to extract random density values
+  TH1                   *fDensitySpectrum;        // particle density spectrum to extract random density values
   TF1                   *fDifferentialV2;         // v2 as function of pt
   Bool_t                 fAddV2;                  // add v2 sampled from a tf1
   Bool_t                 fFlowFluctuations;       // introduce gaussian flow fluctuation 
@@ -122,11 +125,12 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   Int_t                  fMCLabelShift;           //!MC label shift
   Bool_t                 fEsdMode;                //!ESD/AOD mode
   TList                 *fOutput;                 //!output list for QA histograms
+  AliStackPartonInfo    *fStackPartonInfo;     //!Info on original partons:PDG,pt, eta, phi 
 
  private:
   AliJetModelBaseTask(const AliJetModelBaseTask&);            // not implemented
   AliJetModelBaseTask &operator=(const AliJetModelBaseTask&); // not implemented
 
-  ClassDef(AliJetModelBaseTask, 12) // Jet modelling task
+  ClassDef(AliJetModelBaseTask, 11) // Jet modelling task
 };
 #endif
diff --git a/PWGJE/EMCALJetTasks/AliStackPartonInfo.cxx b/PWGJE/EMCALJetTasks/AliStackPartonInfo.cxx
new file mode 100644 (file)
index 0000000..f12cccb
--- /dev/null
@@ -0,0 +1,34 @@
+#include "AliStackPartonInfo.h" 
+
+ClassImp(AliStackPartonInfo)
+
+AliStackPartonInfo::AliStackPartonInfo() :
+  TNamed(),
+  fPartonFlag6(0),
+  fPartonPt6(0),
+  fPartonEta6(0),
+  fPartonPhi6(0),
+  fPartonFlag7(0),
+  fPartonPt7(0),
+  fPartonEta7(0),
+  fPartonPhi7(0)
+{
+  
+}
+
+//_______________________________________________
+
+AliStackPartonInfo::AliStackPartonInfo(const char* name) :
+  TNamed(),
+  fPartonFlag6(0),
+  fPartonPt6(0),
+  fPartonEta6(0),
+  fPartonPhi6(0),
+  fPartonFlag7(0),
+  fPartonPt7(0),
+  fPartonEta7(0),
+  fPartonPhi7(0)
+{
+  
+}
+
diff --git a/PWGJE/EMCALJetTasks/AliStackPartonInfo.h b/PWGJE/EMCALJetTasks/AliStackPartonInfo.h
new file mode 100644 (file)
index 0000000..72b96ed
--- /dev/null
@@ -0,0 +1,51 @@
+#ifndef ALISTACKPARTONINFO_H
+#define ALISTACKPARTONINFO_H
+
+#include <TMath.h>
+#include <TNamed.h>
+
+class AliStackPartonInfo : public TNamed{
+
+ public:
+  AliStackPartonInfo();
+  AliStackPartonInfo(const char* name); 
+
+  void SetPartonFlag6(Int_t flag6) {fPartonFlag6 = flag6;}
+  void SetPartonPt6(Float_t pt6) {fPartonPt6 = pt6;}
+  void SetPartonEta6(Float_t eta6) {fPartonEta6 = eta6;}
+  void SetPartonPhi6(Float_t phi6) {fPartonPhi6 = phi6;}
+
+  void SetPartonFlag7(Int_t flag7) {fPartonFlag7 = flag7;}
+  void SetPartonPt7(Float_t pt7) {fPartonPt7 = pt7;}
+  void SetPartonEta7(Float_t eta7) {fPartonEta7 = eta7;}
+  void SetPartonPhi7(Float_t phi7) {fPartonPhi7 = phi7;}
+
+  
+  Int_t GetPartonFlag6() {return fPartonFlag6;}
+  Float_t GetPartonPt6() {return fPartonPt6;}
+  Float_t GetPartonEta6() {return fPartonEta6;}
+  Float_t GetPartonPhi6() {return fPartonPhi6;}
+
+  Int_t GetPartonFlag7() {return fPartonFlag7;}
+  Float_t GetPartonPt7() {return fPartonPt7;}
+  Float_t GetPartonEta7() {return fPartonEta7;}
+  Float_t GetPartonPhi7() {return fPartonPhi7;}
+
+ private: 
+  Int_t fPartonFlag6; //! parton flag 
+  Float_t fPartonPt6; //! pT parton 
+  Float_t fPartonEta6; //!eta parton 
+  Float_t fPartonPhi6; //! phi parton
+
+  Int_t fPartonFlag7; //! parton flag 
+  Float_t fPartonPt7; //! pT parton 
+  Float_t fPartonEta7; //!eta parton 
+  Float_t fPartonPhi7; //! phi parton
+    
+  AliStackPartonInfo(const AliStackPartonInfo&);
+  AliStackPartonInfo& operator=(const AliStackPartonInfo&);
+  
+  ClassDef(AliStackPartonInfo, 1);
+
+};
+#endif
diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.cxx
new file mode 100644 (file)
index 0000000..456f3dc
--- /dev/null
@@ -0,0 +1,543 @@
+//
+// Jet QG tagging analysis task.
+//
+// Author: D. Caffarri, L. Cunqueiro 
+
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <THnSparse.h>
+#include <TTree.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TProfile.h>
+#include <TChain.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TKey.h>
+#include "TMatrixD.h"
+#include "TMatrixDSym.h"
+#include "TMatrixDSymEigen.h"
+#include "TVector3.h"
+#include "TVector2.h"
+
+#include "AliVCluster.h"
+#include "AliVTrack.h"
+#include "AliEmcalJet.h"
+#include "AliRhoParameter.h"
+#include "AliLog.h"
+#include "AliEmcalParticle.h"
+#include "AliMCEvent.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliAODMCHeader.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliJetContainer.h"
+#include "AliParticleContainer.h"
+#include "AliStackPartonInfo.h"
+
+
+#include "AliAODEvent.h"
+
+#include "AliAnalysisTaskEmcalQGTagging.h"
+
+ClassImp(AliAnalysisTaskEmcalQGTagging)
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging() : 
+  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalQGTagging", kTRUE),
+  fContainer(0),
+  fMinFractionShared(0),
+  fJetShapeType(kRaw),
+  fShapesVar(0),
+  fIsMC(kFALSE),
+  fIsEmbedding(kFALSE),
+  fIsConstSub(kFALSE),
+  fRMatching(0.3),
+  fPhiJetCorr6(0x0), 
+  fPhiJetCorr7(0x0),
+  fEtaJetCorr6(0x0),
+  fEtaJetCorr7(0x0),
+  fPtJetCorr(0x0),
+  fPtJet(0x0),
+  fTreeObservableTagging(0)
+{
+  SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging(const char *name) : 
+  AliAnalysisTaskEmcalJet(name, kTRUE),
+  fContainer(0),
+  fMinFractionShared(0),
+  fJetShapeType(kRaw),
+  fShapesVar(0),
+  fIsMC(kFALSE),
+  fIsEmbedding(kFALSE),
+  fIsConstSub(kFALSE),
+  fRMatching(0.3),
+  fPhiJetCorr6(0x0), 
+  fPhiJetCorr7(0x0),
+  fEtaJetCorr6(0x0),
+  fEtaJetCorr7(0x0),
+  fPtJetCorr(0x0),
+  fPtJet(0x0),
+  fTreeObservableTagging(0)
+{
+  // Standard constructor.
+  
+  SetMakeGeneralHistograms(kTRUE);
+
+  DefineOutput(1, TTree::Class());
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalQGTagging::~AliAnalysisTaskEmcalQGTagging()
+{
+  // Destructor.
+}
+
+//________________________________________________________________________
+ void AliAnalysisTaskEmcalQGTagging::UserCreateOutputObjects()
+{
+  // Create user output.
+
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
+  Int_t nVar = 9; 
+  fShapesVar = new Float_t [nVar]; 
+  TString *fShapesVarNames = new TString [nVar];
+
+  fShapesVarNames[0] = "partonCode"; 
+  fShapesVarNames[1] = "ptJet"; 
+  fShapesVarNames[2] = "ptDJet"; 
+  fShapesVarNames[3] = "mJet";
+  fShapesVarNames[4] = "nbOfConst";
+  fShapesVarNames[5] = "angularity";
+  fShapesVarNames[6] = "circularity";
+  fShapesVarNames[7] = "lesub";
+  fShapesVarNames[8] = "sigma2";
+  
+  for(Int_t ivar=0; ivar < nVar; ivar++){
+    cout<<"looping over variables"<<endl;
+    fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));
+
+    //if( ivar == 4 )  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/I", fShapesVarNames[ivar].Data()));
+
+  }
+  
+  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
+  fOutput->Add(fPhiJetCorr6);
+  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
+  fOutput->Add(fEtaJetCorr6);
+  
+  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
+  fOutput->Add(fPhiJetCorr7);
+  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
+  fOutput->Add(fEtaJetCorr7);
+  
+  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200,  100, 0, 200);
+  fOutput->Add(fPtJetCorr);
+  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
+  fOutput->Add(fPtJet); 
+
+
+  fOutput->Add(fTreeObservableTagging);
+  TH1::AddDirectory(oldStatus);
+  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
+
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalQGTagging::Run()
+{
+  // Run analysis code here, if needed. It will be executed before FillHistograms().
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalQGTagging::FillHistograms()
+{
+  // Fill histograms.
+  //cout<<"base container"<<endl;
+  AliEmcalJet* jet1 = NULL;
+  AliJetContainer *jetCont = GetJetContainer(0);
+  
+  if(jetCont) {
+    jetCont->ResetCurrentID();
+    while((jet1 = jetCont->GetNextAcceptJet())) {
+      if (!jet1) continue;
+      fPtJet->Fill(jet1->Pt());
+      if (fIsMC) {
+       AliStackPartonInfo *partonsInfo = 0x0;
+       AliEmcalJet* jet2 = 0x0;
+       if (fIsEmbedding){ 
+         AliJetContainer *jetContTrue = GetJetContainer(1);
+         jet2 = jet1->ClosestJet();
+         if (!jet2) {
+           Printf("jet2 not exists, returning");
+           continue;
+         }
+         
+         Double_t fraction = jetCont->GetFractionSharedPt(jet1);
+         if(fraction<fMinFractionShared) continue;
+         partonsInfo = (AliStackPartonInfo*) jetContTrue->GetStackPartonInfo();     
+         
+       }
+       else {
+         partonsInfo = (AliStackPartonInfo*) jetCont->GetStackPartonInfo(); 
+         jet2=jet1;
+       }
+       
+       Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6()); 
+       Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
+       
+       if (jp1< -1*TMath::Pi()) jp1 = (-2*TMath::Pi())-jp1;
+       else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
+       Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
+       fEtaJetCorr6->Fill(jet2->Eta(), partonsInfo->GetPartonEta6());
+       fPhiJetCorr6->Fill(jet2->Phi(), partonsInfo->GetPartonPhi6());
+       if(dRp1 < fRMatching) {
+         fShapesVar[0] = partonsInfo->GetPartonFlag6();
+         fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet2->Pt());
+       }
+       else {
+         jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi7());
+         detap1=(jet2->Eta())-(partonsInfo->GetPartonEta7());
+         if (jp1< -1*TMath::Pi()) jp1= (-2*TMath::Pi())-jp1;
+         else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
+         dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
+         fEtaJetCorr7->Fill(jet2->Eta(), partonsInfo->GetPartonEta7());
+         fPhiJetCorr7->Fill(jet2->Phi(), partonsInfo->GetPartonPhi7());
+         if(dRp1 < fRMatching) {
+           fShapesVar[0] = partonsInfo->GetPartonFlag7();
+           fPtJetCorr ->Fill(partonsInfo->GetPartonPt7(), jet2->Pt());
+         }
+         else continue;
+       }
+      }
+      else
+       fShapesVar[0] = 0.;
+     
+
+      if ((fJetShapeType==AliAnalysisTaskEmcalQGTagging::kRaw && fIsConstSub==kFALSE) || (fJetShapeType==AliAnalysisTaskEmcalQGTagging::kDeriv)) fShapesVar[1] = jet1->Pt() - GetRhoVal(0)*jet1->Area();
+      else fShapesVar[1] = jet1->Pt(); 
+     
+      fShapesVar[2] = GetJetpTD(jet1);
+      fShapesVar[3] = GetJetMass(jet1);
+      fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1);
+      fShapesVar[5] = GetJetAngularity(jet1);
+      fShapesVar[6] = GetJetCircularity(jet1);
+      fShapesVar[7] = GetJetLeSub(jet1);
+      fShapesVar[8] = GetSigma2(jet1);
+      fTreeObservableTagging->Fill();
+       
+    }
+    
+  } 
+  
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet) {
+  //calc subtracted jet mass
+  if(fJetShapeType==kDeriv)
+    return jet->GetSecondOrderSubtracted();
+  else 
+    return jet->M();
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet){
+
+  AliJetContainer *jetCont = GetJetContainer(0);
+  if (!jet->GetNumberOfTracks())
+      return 0; 
+    Double_t den=0.;
+    Double_t num = 0.;
+    AliVParticle *vp1 = 0x0;
+    for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+      vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
+      Double_t dphi = vp1->Phi()-jet->Phi();
+      if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
+      if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
+      Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
+      Double_t dr = TMath::Sqrt(dr2);
+      num=num+vp1->Pt()*dr;
+      den=den+vp1->Pt();
+    }
+    return num/den;
+} 
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet) {
+
+  if(fJetShapeType==kDeriv)
+    return jet->GetSecondOrderSubtractedAngularity();
+  else
+    return Angularity(jet);
+}
+
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet){
+
+  AliJetContainer *jetCont = GetJetContainer(0);
+  if (!jet->GetNumberOfTracks())
+      return 0; 
+    Double_t den=0.;
+    Double_t num = 0.;
+    AliVParticle *vp1 = 0x0;
+    for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+      vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
+      num=num+vp1->Pt()*vp1->Pt();
+      den=den+vp1->Pt();
+    }
+    return TMath::Sqrt(num)/den;
+} 
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet) {
+  //calc subtracted jet mass
+  if(fJetShapeType==kDeriv)
+    return jet->GetSecondOrderSubtractedpTD();
+  else
+    return PTD(jet);
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet){
+
+  AliJetContainer *jetCont = GetJetContainer(0);
+  if (!jet->GetNumberOfTracks())
+    return 0; 
+  Double_t mxx    = 0.;
+  Double_t myy    = 0.;
+  Double_t mxy    = 0.;
+  int  nc     = 0;
+  Double_t sump2  = 0.;
+  Double_t pxjet=jet->Px();
+  Double_t pyjet=jet->Py();
+  Double_t pzjet=jet->Pz();
+  
+  
+  //2 general normalized vectors perpendicular to the jet
+  TVector3  ppJ1(pxjet, pyjet, pzjet);
+  TVector3  ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
+  ppJ3.SetMag(1.);
+  TVector3  ppJ2(-pyjet, pxjet, 0);
+  ppJ2.SetMag(1.);
+  AliVParticle *vp1 = 0x0;
+  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+    vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
+    
+    
+    TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
+   
+    //local frame
+    TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
+    TVector3 pPerp = pp - pLong;
+    //projection onto the two perpendicular vectors defined above
+    
+    Float_t ppjX = pPerp.Dot(ppJ2);
+    Float_t ppjY = pPerp.Dot(ppJ3);
+    Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
+    if(ppjT<=0) return 0;
+    
+    mxx += (ppjX * ppjX / ppjT);
+    myy += (ppjY * ppjY / ppjT);
+    mxy += (ppjX * ppjY / ppjT);
+    nc++;
+    sump2 += ppjT;}
+  
+  if(nc<2) return 0;
+  if(sump2==0) return 0;
+  // Sphericity Matrix
+  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
+  TMatrixDSym m0(2,ele);
+  
+  // Find eigenvectors
+  TMatrixDSymEigen m(m0);
+  TVectorD eval(2);
+  TMatrixD evecm = m.GetEigenVectors();
+  eval  = m.GetEigenValues();
+  // Largest eigenvector
+  int jev = 0;
+  //  cout<<eval[0]<<" "<<eval[1]<<endl;
+  if (eval[0] < eval[1]) jev = 1;
+  TVectorD evec0(2);
+  // Principle axis
+  evec0 = TMatrixDColumn(evecm, jev);
+  Double_t compx=evec0[0];
+  Double_t compy=evec0[1];
+  TVector2 evec(compx, compy);
+  Double_t circ=0;
+  if(jev==1) circ=2*eval[0];
+  if(jev==0) circ=2*eval[1];
+  
+  return circ;
+  
+  
+  
+}
+
+
+
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet) {
+  //calc subtracted jet mass
+  if(fJetShapeType==kDeriv)
+    return jet->GetSecondOrderSubtractedCircularity();
+  else
+    return Circularity(jet);
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet){
+
+  AliJetContainer *jetCont = GetJetContainer(0);
+  if (!jet->GetNumberOfTracks())
+      return 0; 
+    Double_t den=0.;
+    Double_t num = 0.;
+    AliVParticle *vp1 = 0x0;
+    AliVParticle *vp2 = 0x0;
+    std::vector<int> ordindex;
+    ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
+    
+   vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));  
+   vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));  
+     
+  num=vp1->Pt();
+  den=vp2->Pt();
+  
+return num-den;
+} 
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet) {
+  //calc subtracted jet mass
+  if(fJetShapeType==kDeriv)
+    return jet->GetSecondOrderSubtractedLeSub();
+  else
+    return LeSub(jet);
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet) {
+  //calc subtracted jet mass
+  
+  if(fJetShapeType==kDeriv)
+    return jet->GetSecondOrderSubtractedConstituent();
+  else
+    return jet->GetNumberOfTracks();
+}
+   
+
+//______________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet){
+
+  AliJetContainer *jetCont = GetJetContainer(0);
+  if (!jet->GetNumberOfTracks())
+      return 0; 
+      Double_t mxx    = 0.;
+      Double_t myy    = 0.;
+      Double_t mxy    = 0.;
+      int  nc     = 0;
+      Double_t sump2  = 0.;
+       
+     AliVParticle *vp1 = 0x0;
+     for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+       vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));  
+       Double_t ppt=vp1->Pt();
+       Double_t dphi = vp1->Phi()-jet->Phi();
+       if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
+       if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
+       Double_t deta = vp1->Eta()-jet->Eta();
+       mxx += ppt*ppt*deta*deta;
+       myy += ppt*ppt*dphi*dphi;
+       mxy -= ppt*ppt*deta*dphi;
+       nc++;
+       sump2 += ppt*ppt;
+       
+     }  
+     if(nc<2) return 0;
+     if(sump2==0) return 0;
+     // Sphericity Matrix
+     Double_t ele[4] = {mxx , mxy , mxy , myy };
+     TMatrixDSym m0(2,ele);
+     
+     // Find eigenvectors
+     TMatrixDSymEigen m(m0);
+     TVectorD eval(2);
+     TMatrixD evecm = m.GetEigenVectors();
+     eval  = m.GetEigenValues();
+     // Largest eigenvector
+     int jev = 0;
+     //  cout<<eval[0]<<" "<<eval[1]<<endl;
+     if (eval[0] < eval[1]) jev = 1;
+     TVectorD evec0(2);
+     // Principle axis
+     evec0 = TMatrixDColumn(evecm, jev);
+     Double_t compx=evec0[0];
+     Double_t compy=evec0[1];
+     TVector2 evec(compx, compy);
+     Double_t sig=0;
+     if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
+     if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
+     
+     return sig;
+     
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet) {
+  //calc subtracted jet mass
+  if(fJetShapeType==kDeriv)
+    return jet->GetSecondOrderSubtractedSigma2();
+  else
+    return Sigma2(jet);
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() {
+  //
+  // retrieve event objects
+  //
+  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
+    return kFALSE;
+
+  return kTRUE;
+}
+
+//_______________________________________________________________________
+void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *) 
+{
+  // Called once at the end of the analysis.
+
+  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
+  // if (!fTreeObservableTagging){
+  //   Printf("ERROR: fTreeObservableTagging not available"); 
+  //   return;
+  // }
+
+}
+
diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.h
new file mode 100644 (file)
index 0000000..9951de5
--- /dev/null
@@ -0,0 +1,89 @@
+#ifndef ALIANALYSISTASKEMCALQGTAGGING_H
+#define ALIANALYSISTASKEMCALQGTAGGING_H
+
+class TH1;
+class TH2;
+class TH3;
+class TH3F;
+class TTree;
+class THnSparse;
+class TClonesArray;
+class TArrayI;
+class AliAnalysisManager;
+class AliJetContainer;
+
+#include "AliAnalysisTaskEmcalJet.h"
+
+
+
+class AliAnalysisTaskEmcalQGTagging : public AliAnalysisTaskEmcalJet {
+ public:
+  enum JetShapeType {
+    kRaw   = 0,  //mass form anti-kt 4-vector
+    kConstSub = 1,   //constituent subtracted jetshape
+    kTrue  = 2, 
+    kDeriv = 3 //area based subtracted jet mass
+  };
+
+  AliAnalysisTaskEmcalQGTagging();
+  AliAnalysisTaskEmcalQGTagging(const char *name);
+  virtual ~AliAnalysisTaskEmcalQGTagging();
+
+  void                                UserCreateOutputObjects();
+  void                                Terminate(Option_t *option);
+
+  //Setters
+  void SetJetContainer(Int_t c)                             { fContainer     = c   ; }
+  void SetMinFractionShared(Double_t f)                     { fMinFractionShared = f   ; }
+  void SetJetShapeType(JetShapeType t)                      { fJetShapeType       = t   ; }
+  void SetMCTask(Int_t f)                                   { fIsMC       = f   ; }
+  void SetEmbeddingTask(Int_t f)                            { fIsEmbedding       = f   ; }
+  void SetIsConstSub(Bool_t f)                              { fIsConstSub     = f   ; }
+  void SetRMatching(Float_t f)                              { fRMatching = f ;}
+
+ protected:
+  Bool_t                              RetrieveEventObjects();
+  Bool_t                              Run();
+  Bool_t                              FillHistograms();
+
+  Float_t                            GetJetMass(AliEmcalJet *jet);
+  Float_t                            Angularity(AliEmcalJet *jet);
+  Float_t                            GetJetAngularity(AliEmcalJet *jet);
+  Float_t                            PTD(AliEmcalJet *jet);
+  Float_t                            GetJetpTD(AliEmcalJet *jet);
+  Float_t                            Circularity(AliEmcalJet *jet); 
+  Float_t                            GetJetCircularity(AliEmcalJet *jet);
+  Float_t                            LeSub(AliEmcalJet *jet);
+  Float_t                            GetJetLeSub(AliEmcalJet *jet);
+  Float_t                            GetJetNumberOfConstituents(AliEmcalJet *jet);
+  Float_t                            GetSigma2(AliEmcalJet *jet);
+  Float_t                            Sigma2(AliEmcalJet *jet);
+
+
+  Int_t                               fContainer;              // jets to be analyzed 0 for Base, 1 for subtracted. 
+  Float_t                            fMinFractionShared;          // only fill histos for jets if shared fraction larger than X
+  JetShapeType                        fJetShapeType;                // jet mass type to be used
+  Float_t                            *fShapesVar;                  // jet shapes used for the tagging
+  Int_t                               fIsMC;
+  Int_t                               fIsEmbedding;
+  Bool_t                              fIsConstSub;
+  Float_t                             fRMatching;
+  
+  TH2F                                *fPhiJetCorr6;
+  TH2F                                *fPhiJetCorr7;
+  TH2F                                *fEtaJetCorr6;
+  TH2F                                *fEtaJetCorr7;
+  TH2F                                *fPtJetCorr;
+  TH1F                                *fPtJet;
+  
+
+  TTree           *fTreeObservableTagging;  //Tree with tagging variables subtracted MC or true MC or raw 
+
+ private:
+  AliAnalysisTaskEmcalQGTagging(const AliAnalysisTaskEmcalQGTagging&);            // not implemented
+  AliAnalysisTaskEmcalQGTagging &operator=(const AliAnalysisTaskEmcalQGTagging&); // not implemented
+
+  ClassDef(AliAnalysisTaskEmcalQGTagging, 1)
+};
+#endif
+
diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskEmcalQGTagging.C b/PWGJE/EMCALJetTasks/macros/AddTaskEmcalQGTagging.C
new file mode 100644 (file)
index 0000000..13eefe9
--- /dev/null
@@ -0,0 +1,95 @@
+AliAnalysisTaskEmcalQGTagging* AddTaskEmcalQGTagging(const char * njetsBase,
+                                                    const char * njetsTrue,
+                                                    const Double_t R,
+                                                    const char * nrhoBase, 
+                                                    const char * ntracks, 
+                                                    const char * nclusters,
+                                                    const char * ntracksTrue,
+                                                    const char *type,                                
+                                                    const char *CentEst,
+                                                    Int_t       pSel,
+                                                    TString     trigClass      = "",
+                                                    TString     kEmcalTriggers = "",
+                                                    TString     tag            = "",
+                                                    AliAnalysisTaskEmcalQGTagging::JetShapeType jetShapeType, Int_t isEmbedding = 0) {
+  
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+    {
+      Error("AddTaskEmcalQGTagging","No analysis manager found.");
+      return 0;
+    }
+  Bool_t ismc=kFALSE;
+  ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
+
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+    {
+      ::Error("AddTaskEmcalQGTagging", "This task requires an input event handler");
+      return NULL;
+    }
+
+  TString wagonName = Form("JetQGTaggings_%s_TC%s%s",njetsBase,trigClass.Data(),tag.Data());
+
+  //Configure jet tagger task
+  AliAnalysisTaskEmcalQGTagging *task = new AliAnalysisTaskEmcalQGTagging(wagonName.Data());
+
+  //task->SetNCentBins(4);
+  task->SetJetShapeType(jetShapeType);
+  TString thename(njetsBase);
+  if(thename.Contains("Sub")) task->SetIsConstSub(kTRUE);
+  //task->SetVzRange(-10.,10.);
+
+  AliParticleContainer *trackCont  = task->AddParticleContainer(ntracks);
+  AliParticleContainer *trackContTrue  = task->AddParticleContainer(ntracksTrue);
+  AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
+
+  task->SetJetContainer(0);
+  
+  TString strType(type);
+  AliJetContainer *jetContBase = task->AddJetContainer(njetsBase,strType,R);
+  if(jetContBase) {
+    jetContBase->SetRhoName(nrhoBase);
+    jetContBase->ConnectParticleContainer(trackCont);
+    jetContBase->ConnectClusterContainer(clusterCont);
+    jetContBase->SetPercAreaCut(0.6);
+    if(jetShapeType == AliAnalysisTaskEmcalQGTagging::kTrue){ jetContBase->SetPartonInfoName("PartonsInfo");}
+  }
+
+
+  if(isEmbedding){
+    task->SetJetContainer(1);
+    
+    AliJetContainer *jetContTrue = task->AddJetContainer(njetsTrue,strType,R);
+    if(jetContTrue) {
+      jetContTrue->SetRhoName(nrhoBase);
+      jetContTrue->ConnectParticleContainer(trackContTrue);
+      jetContTrue->SetPercAreaCut(0.6); 
+      jetContTrue->SetPartonInfoName("PartonsInfo");
+    }
+  } 
+    //  task->SetJetContainer(1);
+
+  task->SetCaloTriggerPatchInfoName(kEmcalTriggers.Data());
+  task->SetCentralityEstimator(CentEst);
+  task->SelectCollisionCandidates(pSel);
+  task->SetUseAliAnaUtils(kFALSE);
+
+  mgr->AddTask(task);
+
+  //Connnect input
+  mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
+
+  //Connect output
+  TString contName1(wagonName);
+  
+  TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName1.Data(), TTree::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+    
+  mgr->ConnectOutput(task,1,coutput1);
+
+  return task;  
+
+}
+
index 59a9f39..0c7de45 100644 (file)
@@ -3,6 +3,7 @@
 AliJetEmbeddingFromGenTask* AddTaskJetEmbeddingFromGen(
   AliGenerator   *genGen,
   const char     *tracksName   = "GenParticles",
+  const char     *partonInfoName = "PartonsInfo",
   const char     *taskName     = "JetEmbeddingFromGenTask",
   const Double_t  minPt        = 10,
   const Double_t  maxPt        = 10,
@@ -40,6 +41,7 @@ AliJetEmbeddingFromGenTask* AddTaskJetEmbeddingFromGen(
   if(TString(genGen->IsA()->GetName()).EqualTo("AliGenPythia")) genGen->AliGenPythia::SetEventListRange(0, 0);
   jetEmb->SetGen(genGen);
   jetEmb->SetTracksName(tracksName);
+  jetEmb->SetPartonInfoName(partonInfoName);
   jetEmb->SetEtaRange(minEta, maxEta);
   jetEmb->SetPhiRange(minPhi, maxPhi);
   jetEmb->SetPtRange(minPt, maxPt);
index e43cc4f..0b889e8 100644 (file)
@@ -30,6 +30,7 @@
 #pragma link C++ class AliRhoParameter+;
 #pragma link C++ class AliLocalRhoParameter+;
 #pragma link C++ class AliJetTriggerSelectionTask+;
+#pragma link C++ class AliStackPartonInfo+;
 
 // user tasks
 #pragma link C++ class AliAnalysisTaskCLQA+;
@@ -64,6 +65,7 @@
 #pragma link C++ class AliAnalysisTaskRhoMass+;
 #pragma link C++ class AliAnalysisTaskRhoMassBase+;
 #pragma link C++ class AliAnalysisTaskRhoMassScale+;
+#pragma link C++ class AliAnalysisTaskEmcalQGTagging+;
 #pragma link C++ class AliAnalysisTaskSAJF+;
 #pragma link C++ class AliAnalysisTaskSAQA+;
 #pragma link C++ class AliAnalysisTaskSOH+;