Add sigma2 jet shape and fill constituent info. for subtracted jets

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJet.cxx b/PWGJE/EMCALJetTasks/AliEmcalJet.cxx
index e98a457eb8630d11762dfc600dfd634af55f7d15..77e4a94002fa65be8c0e51a67142687962051dce 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJet.cxx
+++ b/PWGJE/EMCALJetTasks/AliEmcalJet.cxx
@@ -61,6 +61,10 @@ AliEmcalJet::AliEmcalJet() :
   fJetShapeCircularitySecondDer(0),
   fJetShapeCircularityFirstSub(0),
   fJetShapeCircularitySecondSub(0),
+  fJetShapeSigma2FirstDer(0),
+  fJetShapeSigma2SecondDer(0),
+  fJetShapeSigma2FirstSub(0),
+  fJetShapeSigma2SecondSub(0),
   fJetShapeConstituentFirstDer(0),
   fJetShapeConstituentSecondDer(0),
   fJetShapeConstituentFirstSub(0),
@@ -128,11 +132,15 @@ AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
   fJetShapeCircularitySecondDer(0),
   fJetShapeCircularityFirstSub(0),
   fJetShapeCircularitySecondSub(0),
+  fJetShapeSigma2FirstDer(0),
+  fJetShapeSigma2SecondDer(0),
+  fJetShapeSigma2FirstSub(0),
+  fJetShapeSigma2SecondSub(0),
   fJetShapeConstituentFirstDer(0),
   fJetShapeConstituentSecondDer(0),
   fJetShapeConstituentFirstSub(0),
   fJetShapeConstituentSecondSub(0),
-   fJetShapeLeSubFirstDer(0),
+  fJetShapeLeSubFirstDer(0),
   fJetShapeLeSubSecondDer(0),
   fJetShapeLeSubFirstSub(0),
   fJetShapeLeSubSecondSub(0)
@@ -202,6 +210,10 @@ AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
   fJetShapeCircularitySecondDer(0),
   fJetShapeCircularityFirstSub(0),
   fJetShapeCircularitySecondSub(0),
+  fJetShapeSigma2FirstDer(0),
+  fJetShapeSigma2SecondDer(0),
+  fJetShapeSigma2FirstSub(0),
+  fJetShapeSigma2SecondSub(0),
   fJetShapeConstituentFirstDer(0),
   fJetShapeConstituentSecondDer(0),
   fJetShapeConstituentFirstSub(0),
@@ -274,6 +286,10 @@ AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
   fJetShapeCircularitySecondDer(jet.fJetShapeCircularitySecondDer),
   fJetShapeCircularityFirstSub(jet.fJetShapeCircularityFirstSub),
   fJetShapeCircularitySecondSub(jet.fJetShapeCircularitySecondSub),
+  fJetShapeSigma2FirstDer(jet.fJetShapeSigma2FirstDer),
+  fJetShapeSigma2SecondDer(jet.fJetShapeSigma2SecondDer),
+  fJetShapeSigma2FirstSub(jet.fJetShapeSigma2FirstSub),
+  fJetShapeSigma2SecondSub(jet.fJetShapeSigma2SecondSub),
   fJetShapeConstituentFirstDer(jet.fJetShapeConstituentFirstDer),
   fJetShapeConstituentSecondDer(jet.fJetShapeConstituentSecondDer),
   fJetShapeConstituentFirstSub(jet.fJetShapeConstituentFirstSub),
@@ -348,6 +364,10 @@ AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
     fJetShapeCircularitySecondDer = jet.fJetShapeCircularitySecondDer;
     fJetShapeCircularityFirstSub  = jet.fJetShapeCircularityFirstSub;
     fJetShapeCircularitySecondSub = jet.fJetShapeCircularitySecondSub;
+    fJetShapeSigma2FirstDer  = jet.fJetShapeSigma2FirstDer;
+    fJetShapeSigma2SecondDer = jet.fJetShapeSigma2SecondDer;
+    fJetShapeSigma2FirstSub  = jet.fJetShapeSigma2FirstSub;
+    fJetShapeSigma2SecondSub = jet.fJetShapeSigma2SecondSub;
     fJetShapeConstituentFirstDer  = jet.fJetShapeConstituentFirstDer;
     fJetShapeConstituentSecondDer = jet.fJetShapeConstituentSecondDer;
     fJetShapeConstituentFirstSub  = jet.fJetShapeConstituentFirstSub;

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJet.h b/PWGJE/EMCALJetTasks/AliEmcalJet.h
index 589579e46a752f99325e5433f7e8a90663c8e3e0..c6b2fb678404e47f42139e04f27d68244ec3500b 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJet.h
+++ b/PWGJE/EMCALJetTasks/AliEmcalJet.h
@@ -210,6 +210,17 @@ class AliEmcalJet : public AliVParticle
   Double_t          GetFirstOrderSubtractedCircularity()      const { return fJetShapeCircularityFirstSub             ; }
   Double_t          GetSecondOrderSubtractedCircularity()     const { return fJetShapeCircularitySecondSub            ; }
 
+  //Sigma2
+  void              SetFirstDerivativeSigma2(Double_t d)       { fJetShapeSigma2FirstDer = d                ; }
+  void              SetSecondDerivativeSigma2(Double_t d)      { fJetShapeSigma2SecondDer = d               ; }
+  void              SetFirstOrderSubtractedSigma2(Double_t d)  { fJetShapeSigma2FirstSub = d                ; }
+  void              SetSecondOrderSubtractedSigma2(Double_t d) { fJetShapeSigma2SecondSub = d               ; }
+  Double_t          GetFirstDerivativeSigma2()           const { return fJetShapeSigma2FirstDer             ; }
+  Double_t          GetSecondDerivativeSigma2()          const { return fJetShapeSigma2SecondDer            ; }
+  Double_t          GetFirstOrderSubtractedSigma2()      const { return fJetShapeSigma2FirstSub             ; }
+  Double_t          GetSecondOrderSubtractedSigma2()     const { return fJetShapeSigma2SecondSub            ; }
+
+
   //number of contituents
   void              SetFirstDerivativeConstituent(Double_t d)       { fJetShapeConstituentFirstDer = d                ; }
   void              SetSecondDerivativeConstituent(Double_t d)      { fJetShapeConstituentSecondDer = d               ; }
@@ -287,6 +298,11 @@ class AliEmcalJet : public AliVParticle
   Double_t          fJetShapeCircularityFirstSub;  //!   result from shape derivatives for jet circularity: 1st order subtracted
   Double_t          fJetShapeCircularitySecondSub; //!   result from shape derivatives for jetcircularity: 2nd order subtracted
 
+  Double_t          fJetShapeSigma2FirstDer;  //!   result from shape derivatives for jet sigma2: 1st derivative
+  Double_t          fJetShapeSigma2SecondDer; //!   result from shape derivatives for jet sigma2: 2nd derivative
+  Double_t          fJetShapeSigma2FirstSub;  //!   result from shape derivatives for jet sigma2: 1st order subtracted
+  Double_t          fJetShapeSigma2SecondSub; //!   result from shape derivatives for jetsigma2: 2nd order subtracted
+
   Double_t          fJetShapeConstituentFirstDer;  //!   result from shape derivatives for jet const: 1st derivative
   Double_t          fJetShapeConstituentSecondDer; //!   result from shape derivatives for jet const: 2nd derivative
   Double_t          fJetShapeConstituentFirstSub;  //!   result from shape derivatives for jet const: 1st order subtracted

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
index afc5f8afaad70a30db66145b112b9082426e66e7..ddf2de02c4c2a67e3ce7f77f9b41757cf5904595 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
+++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
@@ -419,6 +419,7 @@ void AliEmcalJetTask::FindJets()
     fjw.DoGenericSubtractionJetAngularity();
     fjw.DoGenericSubtractionJetpTD();
     fjw.DoGenericSubtractionJetCircularity();
+    fjw.DoGenericSubtractionJetSigma2(); 
     fjw.DoGenericSubtractionJetConstituent();
     fjw.DoGenericSubtractionJetLeSub();
   }
@@ -527,6 +528,16 @@ void AliEmcalJetTask::FindJets()
        jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted());
        jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted());
       }
+
+      std::vector<fastjet::contrib::GenericSubtractorInfo> jetSigma2Info = fjw.GetGenSubtractorInfoJetSigma2();
+      Int_t ns = (Int_t)jetSigma2Info.size();
+      if(ns>ij && ns>0) {
+       jet->SetFirstDerivativeSigma2(jetSigma2Info[ij].first_derivative());
+       jet->SetSecondDerivativeSigma2(jetSigma2Info[ij].second_derivative());
+       jet->SetFirstOrderSubtractedSigma2(jetSigma2Info[ij].first_order_subtracted());
+       jet->SetSecondOrderSubtractedSigma2(jetSigma2Info[ij].second_order_subtracted());
+      } 
+
       
       std::vector<fastjet::contrib::GenericSubtractorInfo> jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent();
       Int_t nco= (Int_t)jetConstituentInfo.size();
@@ -548,7 +559,7 @@ void AliEmcalJetTask::FindJets()
    }
 #endif
 
-    // loop over constituents
+    // Fill constituent info
     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
     jet->SetNumberOfTracks(constituents.size());
     jet->SetNumberOfClusters(constituents.size());
@@ -566,114 +577,7 @@ void AliEmcalJetTask::FindJets()
     Double_t mcpt       = 0;
     Double_t emcpt      = 0;
 
-    for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
-      Int_t uid = constituents[ic].user_index();
-      AliDebug(3,Form("Processing constituent %d", uid));
-      if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
-        ++gall;
-        Double_t gphi = constituents[ic].phi();
-        if (gphi<0) 
-          gphi += TMath::TwoPi();
-        gphi *= TMath::RadToDeg();
-        Double_t geta = constituents[ic].eta();
-        if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
-            (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
-          ++gemc; 
-      }        else if ((uid > 0) && fTracks) { // track constituent
-       Int_t tid = uid - 100;
-        AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
-        if (!t) {
-         AliError(Form("Could not find track %d",tid));
-          continue;
-       }
-       if (jetCount < fMarkConst) {
-         AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
-         t->SetBit(fJetType);
-       }
-        Double_t cEta = t->Eta();
-        Double_t cPhi = t->Phi();
-        Double_t cPt  = t->Pt();
-        Double_t cP   = t->P();
-        if (t->Charge() == 0) {
-          neutralE += cP;
-          ++nneutral;
-          if (cPt > maxNe)
-            maxNe = cPt;
-        } else {
-          ++ncharged;
-          if (cPt > maxCh)
-            maxCh = cPt;
-        }
-        if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
-          mcpt += cPt;
-
-        if (cPhi<0) 
-          cPhi += TMath::TwoPi();
-        cPhi *= TMath::RadToDeg();
-        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
-            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
-          emcpt += cPt;
-          ++cemc;
-        }
-
-        jet->AddTrackAt(tid, nt);
-        ++nt;
-      } else if (fClus) { // cluster constituent
-       Int_t cid = -uid - 100;
-       AliVCluster *c = 0;
-        Double_t cEta=0,cPhi=0,cPt=0,cP=0;
-        if (fIsEmcPart) {
-          AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
-          if (!ep)
-            continue;
-          c = ep->GetCluster();
-          if (!c)
-            continue;
-         if (jetCount < fMarkConst)
-           ep->SetBit(fJetType);
-          cEta = ep->Eta();
-          cPhi = ep->Phi();
-          cPt  = ep->Pt();
-          cP   = ep->P();
-        } else {
-          c = static_cast<AliVCluster*>(fClus->At(cid));
-          if (!c)
-            continue;
-         if (jetCount < fMarkConst)
-           c->SetBit(fJetType);
-          TLorentzVector nP;
-          c->GetMomentum(nP, vertex);
-          cEta = nP.Eta();
-          cPhi = nP.Phi();
-          cPt  = nP.Pt();
-          cP   = nP.P();
-        }
-
-        neutralE += cP;
-        if (cPt > maxNe)
-          maxNe = cPt;
-
-        if (c->GetLabel() > fMinMCLabel) // MC particle
-          mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
-
-        if (cPhi<0) 
-          cPhi += TMath::TwoPi();
-        cPhi *= TMath::RadToDeg();
-        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
-            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
-          emcpt += cPt;
-          ++cemc;
-        }
-
-        jet->AddClusterAt(cid, nc);
-        ++nc;
-        ++nneutral;
-      } else {
-        AliError(Form("%s: No logical way to end up here.", GetName()));
-        continue;
-      }
-    } /* end of constituent loop */
-
+    FillJetConstituents(constituents,jet,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc);
     jet->SetNumberOfTracks(nt);
     jet->SetNumberOfClusters(nc);
     jet->SortConstituents();
@@ -718,6 +622,30 @@ 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
+       std::vector<fastjet::PseudoJet> constituents_sub(fjw.GetJetConstituents(ijet));
+       jet_sub->SetNumberOfTracks(constituents_sub.size());
+       jet_sub->SetNumberOfClusters(constituents_sub.size());
+       Int_t nt            = 0;
+       Int_t nc            = 0;
+       Double_t neutralE   = 0;
+       Double_t maxCh      = 0;
+       Double_t maxNe      = 0;
+       Int_t gall          = 0;
+       Int_t gemc       =0;
+       Int_t cemc          = 0;
+       Int_t ncharged      = 0;
+       Int_t nneutral      = 0;
+       Double_t mcpt       = 0;
+       Double_t emcpt      = 0;
+
+       FillJetConstituents(constituents_sub,jet_sub,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc);
+      jet_sub->SetNumberOfTracks(nt);
+      jet_sub->SetNumberOfClusters(nc);
+      jet_sub->SortConstituents();
+       
+
+
        fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
        jet_sub->SetArea(area.perp());
        jet_sub->SetAreaEta(area.eta());
@@ -727,6 +655,11 @@ void AliEmcalJetTask::FindJets()
       }
     }
 #endif
+
+
+
+
+
   } //constituent subtraction
 }
 
@@ -828,3 +761,134 @@ Bool_t AliEmcalJetTask::DoInit()
   
   return 1;
 }
+
+//___________________________________________________________________________________________________________________
+void  AliEmcalJetTask::FillJetConstituents(std::vector<fastjet::PseudoJet>& constituents,AliEmcalJet *jet,Double_t vertex[3],Int_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc){
+    nt            = 0;
+    nc            = 0;
+    neutralE   = 0;
+    maxCh      = 0;
+    maxNe      = 0;
+    gall          = 0;
+    gemc       =0;
+    cemc          = 0;
+    ncharged      = 0;
+    nneutral      = 0;
+    mcpt       = 0;
+    emcpt      = 0;
+   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
+    for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
+      Int_t uid = constituents[ic].user_index();
+      AliDebug(3,Form("Processing constituent %d", uid));
+      if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
+        ++gall;
+        Double_t gphi = constituents[ic].phi();
+        if (gphi<0) 
+          gphi += TMath::TwoPi();
+        gphi *= TMath::RadToDeg();
+        Double_t geta = constituents[ic].eta();
+        if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
+            (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
+          ++gemc; 
+      }        else if ((uid > 0) && fTracks) { // track constituent
+       Int_t tid = uid - 100;
+        AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
+        if (!t) {
+         AliError(Form("Could not find track %d",tid));
+          continue;
+       }
+       if (jetCount < fMarkConst) {
+         AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
+         t->SetBit(fJetType);
+       }
+        Double_t cEta = t->Eta();
+        Double_t cPhi = t->Phi();
+        Double_t cPt  = t->Pt();
+        Double_t cP   = t->P();
+        if (t->Charge() == 0) {
+          neutralE += cP;
+          ++nneutral;
+          if (cPt > maxNe)
+            maxNe = cPt;
+        } else {
+          ++ncharged;
+          if (cPt > maxCh)
+            maxCh = cPt;
+        }
+        if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
+          mcpt += cPt;
+
+        if (cPhi<0) 
+          cPhi += TMath::TwoPi();
+        cPhi *= TMath::RadToDeg();
+        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
+            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+          emcpt += cPt;
+          ++cemc;
+        }
+
+        jet->AddTrackAt(tid, nt);
+        ++nt;
+      } else if (fClus) { // cluster constituent
+       Int_t cid = -uid - 100;
+       AliVCluster *c = 0;
+        Double_t cEta=0,cPhi=0,cPt=0,cP=0;
+        if (fIsEmcPart) {
+          AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
+          if (!ep)
+            continue;
+          c = ep->GetCluster();
+          if (!c)
+            continue;
+         if (jetCount < fMarkConst)
+           ep->SetBit(fJetType);
+          cEta = ep->Eta();
+          cPhi = ep->Phi();
+          cPt  = ep->Pt();
+          cP   = ep->P();
+        } else {
+          c = static_cast<AliVCluster*>(fClus->At(cid));
+          if (!c)
+            continue;
+         if (jetCount < fMarkConst)
+           c->SetBit(fJetType);
+          TLorentzVector nP;
+          c->GetMomentum(nP, vertex);
+          cEta = nP.Eta();
+          cPhi = nP.Phi();
+          cPt  = nP.Pt();
+          cP   = nP.P();
+        }
+
+        neutralE += cP;
+        if (cPt > maxNe)
+          maxNe = cPt;
+
+        if (c->GetLabel() > fMinMCLabel) // MC particle
+          mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
+
+        if (cPhi<0) 
+          cPhi += TMath::TwoPi();
+        cPhi *= TMath::RadToDeg();
+        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
+            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+          emcpt += cPt;
+          ++cemc;
+        }
+
+        jet->AddClusterAt(cid, nc);
+        ++nc;
+        ++nneutral;
+      } else {
+        AliError(Form("%s: No logical way to end up here.", GetName()));
+        continue;
+      } 
+
+
+
+
+    }
+
+     
+}
+//______________________________________________________________________________________

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.h b/PWGJE/EMCALJetTasks/AliEmcalJetTask.h
index 6077913e3687e8a7c011a6e591316ab839c38571..0b9e06217df96557336ee3ed583fe2994e1623bb 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.h
+++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.h
@@ -13,7 +13,7 @@ namespace fastjet {
 #include "AliAnalysisTaskSE.h"
 #include "AliFJWrapper.h"
 #include "AliAODMCParticle.h"
-
+#include "AliEmcalJet.h"
 class AliEmcalJetTask : public AliAnalysisTaskSE {
  public:
 
@@ -40,7 +40,8 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   void                   UserCreateOutputObjects();
   void                   UserExec(Option_t *option);
   void                   Terminate(Option_t *option);
-
+  void                    FillJetConstituents(std::vector<fastjet::PseudoJet>& constituents,AliEmcalJet *jet,Double_t vertex[3],Int_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc); 
+ 
   Bool_t                 IsLocked()                       { if(fLocked) {AliFatal("Jet finder task is locked! Changing properties is not allowed."); return kTRUE;} else return kFALSE;};
   void                   SetLocked()                      { fLocked = kTRUE;}
   void                   SelectConstituents(UInt_t constSel, UInt_t MCconstSel)  { fConstSel = constSel; fMCConstSel = MCconstSel; };

diff --git a/PWGJE/EMCALJetTasks/AliFJWrapper.h b/PWGJE/EMCALJetTasks/AliFJWrapper.h
index a4b904c7d96627b3f08735a3003b573a55f58122..45a2b81e5b84b2b9654281c7124f1718581f0042 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliFJWrapper.h
+++ b/PWGJE/EMCALJetTasks/AliFJWrapper.h
@@ -41,6 +41,7 @@ class AliFJWrapper
   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetAngularity()  const {return fGenSubtractorInfoJetAngularity  ; } 
   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetpTD()         const {return fGenSubtractorInfoJetpTD         ; }
   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetCircularity() const {return fGenSubtractorInfoJetCircularity ; }
+  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetSigma2() const {return fGenSubtractorInfoJetSigma2; }
   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetConstituent() const {return fGenSubtractorInfoJetConstituent ; } 
   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetLeSub()       const {return fGenSubtractorInfoJetLeSub       ; }
   const std::vector<fastjet::PseudoJet>                      GetConstituentSubtrJets()            const {return fConstituentSubtrJets            ; }
@@ -56,6 +57,7 @@ class AliFJWrapper
   virtual Int_t DoGenericSubtractionJetAngularity();
   virtual Int_t DoGenericSubtractionJetpTD();
   virtual Int_t DoGenericSubtractionJetCircularity();
+  virtual Int_t DoGenericSubtractionJetSigma2();
   virtual Int_t DoGenericSubtractionJetConstituent();
   virtual Int_t DoGenericSubtractionJetLeSub();
   virtual Int_t DoConstituentSubtraction();
@@ -122,6 +124,7 @@ class AliFJWrapper
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetAngularity;  //!  
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetpTD;         //!
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetCircularity; //!
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetSigma2; //!
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetConstituent; //!
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetLeSub;       //!
 #endif
@@ -194,6 +197,7 @@ AliFJWrapper::AliFJWrapper(const char *name, const char *title)
   , fGenSubtractorInfoJetAngularity ( )
   , fGenSubtractorInfoJetpTD ( )
   , fGenSubtractorInfoJetCircularity( )
+  , fGenSubtractorInfoJetSigma2()
   , fGenSubtractorInfoJetConstituent ( )
   , fGenSubtractorInfoJetLeSub ( )
 #endif
@@ -690,7 +694,27 @@ Int_t AliFJWrapper::DoGenericSubtractionJetCircularity() {
 #endif
  return 0;
 }
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionJetSigma2() {
+  //Do generic subtraction for jet mass
+#ifdef FASTJET_VERSION
+  if(fUseExternalBkg)   fGenSubtractor     = new fj::contrib::GenericSubtractor(fRho,fRhom);
+  else                  fGenSubtractor     = new fj::contrib::GenericSubtractor(fBkrdEstimator);
+
+  // Define jet shape
+  AliJetShapeSigma2 shapesigma2;
 
+  // clear the generic subtractor info vector
+  fGenSubtractorInfoJetSigma2.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::contrib::GenericSubtractorInfo infoSigma;
+    if(fInclusiveJets[i].perp()>0.)
+      double subtracted_shape = (*fGenSubtractor)(shapesigma2, fInclusiveJets[i], infoSigma);
+    fGenSubtractorInfoJetSigma2.push_back(infoSigma);
+  }
+#endif
+ return 0;
+}
 //_________________________________________________________________________________________________
 Int_t AliFJWrapper::DoGenericSubtractionJetConstituent() {
   //Do generic subtraction for jet mass

diff --git a/PWGJE/EMCALJetTasks/AliJetShape.h b/PWGJE/EMCALJetTasks/AliJetShape.h
index 250a691faffaab028aaa7ea032d8e3dfe2d5ea48..970b608198ca8283158986c4e7d0d925aabee0b2 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliJetShape.h
+++ b/PWGJE/EMCALJetTasks/AliJetShape.h
@@ -229,6 +229,63 @@ class AliJetShapeCircularity : public fastjet::FunctionOfPseudoJet<Double32_t>{
   }
 };
 
+class AliJetShapeSigma2 : public fastjet::FunctionOfPseudoJet<Double32_t>{
+ public:
+  virtual std::string description() const{return "cms sigma2";}
+  Double32_t result(const fastjet::PseudoJet &jet) const {
+    if (!jet.has_constituents())
+      return 0;
+    Double_t mxx    = 0.;
+    Double_t myy    = 0.;
+    Double_t mxy    = 0.;
+    int  nc     = 0;
+    Double_t sump2  = 0.;
+           
+       
+    std::vector<fastjet::PseudoJet> constits = jet.constituents();
+    for(UInt_t ic = 0; ic < constits.size(); ++ic) {
+      Double_t ppt=constits[ic].perp();
+      Double_t dphi = constits[ic].phi()-jet.phi();
+      if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
+      if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
+      Double_t deta = constits[ic].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;
+    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 sigma2=0;
+    if(jev==1) sigma2=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
+    if(jev==0) sigma2=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
+    
+    return sigma2;
+    
+  }
+};
+
+
+
 class AliJetShapeLeSub : public fastjet::FunctionOfPseudoJet<Double32_t>{
  public:
   virtual std::string description() const{return "leading mins subleading";}