]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Adding new set of jet shapes
authorlcunquei <lcunquei@cern.ch>
Tue, 19 Aug 2014 15:47:12 +0000 (17:47 +0200)
committermvl <marco.van.leeuwen@cern.ch>
Wed, 20 Aug 2014 09:59:10 +0000 (11:59 +0200)
PWGJE/EMCALJetTasks/AliEmcalJet.cxx
PWGJE/EMCALJetTasks/AliEmcalJet.h
PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
PWGJE/EMCALJetTasks/AliEmcalJetTask.h
PWGJE/EMCALJetTasks/AliFJWrapper.h
PWGJE/EMCALJetTasks/AliJetShape.h

index e299a77968454e99295796b1c89ead2000c1a886..9f99aeaa271cc33da423229667d1463f57224b44 100644 (file)
@@ -48,7 +48,27 @@ AliEmcalJet::AliEmcalJet() :
   fGRNumerator(0),
   fGRDenominator(0),
   fGRNumeratorSub(0),
-  fGRDenominatorSub(0)
+  fGRDenominatorSub(0),
+  fJetShapeAngularityFirstDer(0),
+  fJetShapeAngularitySecondDer(0),
+  fJetShapeAngularityFirstSub(0),
+  fJetShapeAngularitySecondSub(0),
+  fJetShapepTDFirstDer(0),
+  fJetShapepTDSecondDer(0),
+  fJetShapepTDFirstSub(0),
+  fJetShapepTDSecondSub(0),
+  fJetShapeCircularityFirstDer(0),
+  fJetShapeCircularitySecondDer(0),
+  fJetShapeCircularityFirstSub(0),
+  fJetShapeCircularitySecondSub(0),
+  fJetShapeConstituentFirstDer(0),
+  fJetShapeConstituentSecondDer(0),
+  fJetShapeConstituentFirstSub(0),
+  fJetShapeConstituentSecondSub(0),
+  fJetShapeLeSubFirstDer(0),
+  fJetShapeLeSubSecondDer(0),
+  fJetShapeLeSubFirstSub(0),
+  fJetShapeLeSubSecondSub(0)
 {
   // Constructor.
 
@@ -96,7 +116,27 @@ AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
   fGRNumerator(0),
   fGRDenominator(0),
   fGRNumeratorSub(0),
-  fGRDenominatorSub(0)
+  fGRDenominatorSub(0),
+    fJetShapeAngularityFirstDer(0),
+  fJetShapeAngularitySecondDer(0),
+  fJetShapeAngularityFirstSub(0),
+  fJetShapeAngularitySecondSub(0),
+   fJetShapepTDFirstDer(0),
+  fJetShapepTDSecondDer(0),
+  fJetShapepTDFirstSub(0),
+  fJetShapepTDSecondSub(0),
+  fJetShapeCircularityFirstDer(0),
+  fJetShapeCircularitySecondDer(0),
+  fJetShapeCircularityFirstSub(0),
+  fJetShapeCircularitySecondSub(0),
+  fJetShapeConstituentFirstDer(0),
+  fJetShapeConstituentSecondDer(0),
+  fJetShapeConstituentFirstSub(0),
+  fJetShapeConstituentSecondSub(0),
+   fJetShapeLeSubFirstDer(0),
+  fJetShapeLeSubSecondDer(0),
+  fJetShapeLeSubFirstSub(0),
+  fJetShapeLeSubSecondSub(0)
 {    
   // Constructor.
 
@@ -150,7 +190,28 @@ AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
   fGRNumerator(0),
   fGRDenominator(0),
   fGRNumeratorSub(0),
-  fGRDenominatorSub(0)
+  fGRDenominatorSub(0),
+  fJetShapeAngularityFirstDer(0),
+  fJetShapeAngularitySecondDer(0),
+  fJetShapeAngularityFirstSub(0),
+  fJetShapeAngularitySecondSub(0),
+  fJetShapepTDFirstDer(0),
+  fJetShapepTDSecondDer(0),
+  fJetShapepTDFirstSub(0),
+  fJetShapepTDSecondSub(0),
+  fJetShapeCircularityFirstDer(0),
+  fJetShapeCircularitySecondDer(0),
+  fJetShapeCircularityFirstSub(0),
+  fJetShapeCircularitySecondSub(0),
+  fJetShapeConstituentFirstDer(0),
+  fJetShapeConstituentSecondDer(0),
+  fJetShapeConstituentFirstSub(0),
+  fJetShapeConstituentSecondSub(0),
+  fJetShapeLeSubFirstDer(0),
+  fJetShapeLeSubSecondDer(0),
+  fJetShapeLeSubFirstSub(0),
+  fJetShapeLeSubSecondSub(0)
+
 {
   // Constructor.
 
@@ -201,7 +262,31 @@ AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
   fGRNumerator(jet.fGRNumerator),
   fGRDenominator(jet.fGRDenominator),
   fGRNumeratorSub(jet.fGRNumeratorSub),
-  fGRDenominatorSub(jet.fGRDenominatorSub)
+  fGRDenominatorSub(jet.fGRDenominatorSub),
+   fJetShapeAngularityFirstDer(jet.fJetShapeAngularityFirstDer),
+  fJetShapeAngularitySecondDer(jet.fJetShapeAngularitySecondDer),
+  fJetShapeAngularityFirstSub(jet.fJetShapeAngularityFirstSub),
+  fJetShapeAngularitySecondSub(jet.fJetShapeAngularitySecondSub),
+
+  fJetShapepTDFirstDer(jet.fJetShapepTDFirstDer),
+  fJetShapepTDSecondDer(jet.fJetShapepTDSecondDer),
+  fJetShapepTDFirstSub(jet.fJetShapepTDFirstSub),
+  fJetShapepTDSecondSub(jet.fJetShapepTDSecondSub),
+
+  fJetShapeCircularityFirstDer(jet.fJetShapeCircularityFirstDer),
+  fJetShapeCircularitySecondDer(jet.fJetShapeCircularitySecondDer),
+  fJetShapeCircularityFirstSub(jet.fJetShapeCircularityFirstSub),
+  fJetShapeCircularitySecondSub(jet.fJetShapeCircularitySecondSub),
+
+  fJetShapeConstituentFirstDer(jet.fJetShapeConstituentFirstDer),
+  fJetShapeConstituentSecondDer(jet.fJetShapeConstituentSecondDer),
+  fJetShapeConstituentFirstSub(jet.fJetShapeConstituentFirstSub),
+  fJetShapeConstituentSecondSub(jet.fJetShapeConstituentSecondSub),
+
+  fJetShapeLeSubFirstDer(jet.fJetShapeLeSubFirstDer),
+  fJetShapeLeSubSecondDer(jet.fJetShapeLeSubSecondDer),
+  fJetShapeLeSubFirstSub(jet.fJetShapeLeSubFirstSub),
+  fJetShapeLeSubSecondSub(jet.fJetShapeLeSubSecondSub)
 {
   // Copy constructor.
 
@@ -257,6 +342,30 @@ AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
     fGRDenominator      = jet.fGRDenominator;
     fGRNumeratorSub     = jet.fGRNumeratorSub;
     fGRDenominatorSub   = jet.fGRDenominatorSub;
+           fJetShapeAngularityFirstDer  = jet.fJetShapeAngularityFirstDer;
+    fJetShapeAngularitySecondDer = jet.fJetShapeAngularitySecondDer;
+    fJetShapeAngularityFirstSub  = jet.fJetShapeAngularityFirstSub;
+    fJetShapeAngularitySecondSub = jet.fJetShapeAngularitySecondSub;
+
+    fJetShapepTDFirstDer  = jet.fJetShapepTDFirstDer;
+    fJetShapepTDSecondDer = jet.fJetShapepTDSecondDer;
+    fJetShapepTDFirstSub  = jet.fJetShapepTDFirstSub;
+    fJetShapepTDSecondSub = jet.fJetShapepTDSecondSub;
+
+    fJetShapeCircularityFirstDer  = jet.fJetShapeCircularityFirstDer;
+    fJetShapeCircularitySecondDer = jet.fJetShapeCircularitySecondDer;
+    fJetShapeCircularityFirstSub  = jet.fJetShapeCircularityFirstSub;
+    fJetShapeCircularitySecondSub = jet.fJetShapeCircularitySecondSub;
+
+    fJetShapeConstituentFirstDer  = jet.fJetShapeConstituentFirstDer;
+    fJetShapeConstituentSecondDer = jet.fJetShapeConstituentSecondDer;
+    fJetShapeConstituentFirstSub  = jet.fJetShapeConstituentFirstSub;
+    fJetShapeConstituentSecondSub = jet.fJetShapeConstituentSecondSub;
+    fJetShapeLeSubFirstDer  = jet.fJetShapeLeSubFirstDer;
+    fJetShapeLeSubSecondDer = jet.fJetShapeLeSubSecondDer;
+    fJetShapeLeSubFirstSub  = jet.fJetShapeLeSubFirstSub;
+    fJetShapeLeSubSecondSub = jet.fJetShapeLeSubSecondSub;
   }
 
   return *this;
index ff2d80bd14695d6b32487ab0b355134b3ba96c4c..05e939a7eceaa79df5d7799ff6bef388e9ef5c16 100644 (file)
@@ -173,6 +173,56 @@ class AliEmcalJet : public AliVParticle
   void              SetGRNumSubSize(UInt_t s) {fGRNumeratorSub.Set(s); }
   void              SetGRDenSubSize(UInt_t s) {fGRDenominatorSub.Set(s); }
   void              PrintGR();
+  void              SetFirstDerivativeAngularity(Double_t d)                  { fJetShapeAngularityFirstDer = d                       ; }
+  void              SetSecondDerivativeAngularity(Double_t d)                 { fJetShapeAngularitySecondDer = d                      ; }
+  void              SetFirstOrderSubtractedAngularity(Double_t d)             { fJetShapeAngularityFirstSub = d                       ; }
+  void              SetSecondOrderSubtractedAngularity(Double_t d)            { fJetShapeAngularitySecondSub = d                      ; }
+  Double_t          GetFirstDerivativeAngularity()                      const { return fJetShapeAngularityFirstDer                    ; }
+  Double_t          GetSecondDerivativeAngularity()                     const { return fJetShapeAngularitySecondDer                   ; }
+  Double_t          GetFirstOrderSubtractedAngularity()                 const { return fJetShapeAngularityFirstSub                    ; }
+  Double_t          GetSecondOrderSubtractedAngularity()                const { return fJetShapeAngularitySecondSub                   ; }
+
+ void              SetFirstDerivativepTD(Double_t d)                  { fJetShapepTDFirstDer = d                       ; }
+  void              SetSecondDerivativepTD(Double_t d)                 { fJetShapepTDSecondDer = d                      ; }
+  void              SetFirstOrderSubtractedpTD(Double_t d)             { fJetShapepTDFirstSub = d                       ; }
+  void              SetSecondOrderSubtractedpTD(Double_t d)            { fJetShapepTDSecondSub = d                      ; }
+  Double_t          GetFirstDerivativepTD()                      const { return fJetShapepTDFirstDer                    ; }
+  Double_t          GetSecondDerivativepTD()                     const { return fJetShapepTDSecondDer                   ; }
+  Double_t          GetFirstOrderSubtractedpTD()                 const { return fJetShapepTDFirstSub                    ; }
+  Double_t          GetSecondOrderSubtractedpTD()                const { return fJetShapepTDSecondSub                   ; }
+
+ void              SetFirstDerivativeCircularity(Double_t d)                  { fJetShapeCircularityFirstDer = d                       ; }
+  void              SetSecondDerivativeCircularity(Double_t d)                 { fJetShapeCircularitySecondDer = d                      ; }
+  void              SetFirstOrderSubtractedCircularity(Double_t d)             { fJetShapeCircularityFirstSub = d                       ; }
+  void              SetSecondOrderSubtractedCircularity(Double_t d)            { fJetShapeCircularitySecondSub = d                      ; }
+  Double_t          GetFirstDerivativeCircularity()                      const { return fJetShapeCircularityFirstDer                    ; }
+  Double_t          GetSecondDerivativeCircularity()                     const { return fJetShapeCircularitySecondDer                   ; }
+  Double_t          GetFirstOrderSubtractedCircularity()                 const { return fJetShapeCircularityFirstSub                    ; }
+  Double_t          GetSecondOrderSubtractedCircularity()                const { return fJetShapeCircularitySecondSub                   ; }
+
+
+
+
+ void              SetFirstDerivativeConstituent(Double_t d)                  { fJetShapeConstituentFirstDer = d                       ; }
+  void              SetSecondDerivativeConstituent(Double_t d)                 { fJetShapeConstituentSecondDer = d                      ; }
+  void              SetFirstOrderSubtractedConstituent(Double_t d)             { fJetShapeConstituentFirstSub = d                       ; }
+  void              SetSecondOrderSubtractedConstituent(Double_t d)            { fJetShapeConstituentSecondSub = d                      ; }
+  Double_t          GetFirstDerivativeConstituent()                      const { return fJetShapeConstituentFirstDer                    ; }
+  Double_t          GetSecondDerivativeConstituent()                     const { return fJetShapeConstituentSecondDer                   ; }
+  Double_t          GetFirstOrderSubtractedConstituent()                 const { return fJetShapeConstituentFirstSub                    ; }
+  Double_t          GetSecondOrderSubtractedConstituent()                const { return fJetShapeConstituentSecondSub                   ; }
+
+ void              SetFirstDerivativeLeSub(Double_t d)                  { fJetShapeLeSubFirstDer = d                       ; }
+  void              SetSecondDerivativeLeSub(Double_t d)                 { fJetShapeLeSubSecondDer = d                      ; }
+  void              SetFirstOrderSubtractedLeSub(Double_t d)             { fJetShapeLeSubFirstSub = d                       ; }
+  void              SetSecondOrderSubtractedLeSub(Double_t d)            { fJetShapeLeSubSecondSub = d                      ; }
+  Double_t          GetFirstDerivativeLeSub()                      const { return fJetShapeLeSubFirstDer                    ; }
+  Double_t          GetSecondDerivativeLeSub()                     const { return fJetShapeLeSubSecondDer                   ; }
+  Double_t          GetFirstOrderSubtractedLeSub()                 const { return fJetShapeLeSubFirstSub                    ; }
+  Double_t          GetSecondOrderSubtractedLeSub()                const { return fJetShapeLeSubSecondSub                   ; }
+
 
  protected:
   Double32_t        fPt;                  //[0,0,12]   pt 
@@ -214,6 +264,31 @@ class AliEmcalJet : public AliVParticle
   TArrayF           fGRDenominator;       //!           array with angular structure function denominator
   TArrayF           fGRNumeratorSub;      //!           array with angular structure function numerator
   TArrayF           fGRDenominatorSub;    //!           array with angular structure function denominator
+ Double_t          fJetShapeAngularityFirstDer;  //!        result from shape derivatives for jet Angularity: 1st derivative
+  Double_t          fJetShapeAngularitySecondDer; //!        result from shape derivatives for jet Angularity: 2nd derivative
+  Double_t          fJetShapeAngularityFirstSub;  //!        result from shape derivatives for jet Angularity: 1st order subtracted
+  Double_t          fJetShapeAngularitySecondSub; //!        result from shape derivatives for jet Angularity: 2nd order subtracted
+
+  Double_t          fJetShapepTDFirstDer;  //!        result from shape derivatives for jet pTD: 1st derivative
+  Double_t          fJetShapepTDSecondDer; //!        result from shape derivatives for jet pTD: 2nd derivative
+  Double_t          fJetShapepTDFirstSub;  //!        result from shape derivatives for jet pTD: 1st order subtracted
+  Double_t          fJetShapepTDSecondSub; //!        result from shape derivatives for jet pTD: 2nd order subtracted
+  Double_t          fJetShapeCircularityFirstDer;  //!        result from shape derivatives for jet circularity: 1st derivative
+  Double_t          fJetShapeCircularitySecondDer; //!        result from shape derivatives for jet circularity: 2nd derivative
+  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          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
+  Double_t          fJetShapeConstituentSecondSub; //!        result from shape derivatives for jet const: 2nd order subtracted
+
+    Double_t          fJetShapeLeSubFirstDer;  //!        result from shape derivatives for jet LeSub: 1st derivative
+  Double_t          fJetShapeLeSubSecondDer; //!        result from shape derivatives for jet LeSub: 2nd derivative
+  Double_t          fJetShapeLeSubFirstSub;  //!        result from shape derivatives for jet LeSub: 1st order subtracted
+  Double_t          fJetShapeLeSubSecondSub; //!        result from shape derivatives for jet LeSub: 2nd order subtracted 
 
   private:
     struct sort_descend
index dddc0c76ee228e1b56510608362e2ede998b0d6d..ff7515d0021bbdc0ad75976c367360a807f060c2 100644 (file)
@@ -70,6 +70,7 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fCodeDebug(kFALSE),
   fDoGenericSubtractionJetMass(kFALSE),
   fDoGenericSubtractionGR(kFALSE),
+  fDoGenericSubtractionExtraJetShapes(kFALSE),
   fDoConstituentSubtraction(kFALSE),
   fUseExternalBkg(kFALSE),
   fRhoName(""),
@@ -129,6 +130,7 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fCodeDebug(kFALSE),
   fDoGenericSubtractionJetMass(kFALSE),
   fDoGenericSubtractionGR(kFALSE),
+  fDoGenericSubtractionExtraJetShapes(kFALSE),
   fDoConstituentSubtraction(kFALSE),
   fUseExternalBkg(kFALSE),
   fRhoName(""),
@@ -407,6 +409,17 @@ void AliEmcalJetTask::FindJets()
     fjw.DoGenericSubtractionJetMass();
   }
 
+     if(fDoGenericSubtractionExtraJetShapes){
+
+    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+    fjw.DoGenericSubtractionJetAngularity();
+    fjw.DoGenericSubtractionJetpTD();
+    fjw.DoGenericSubtractionJetCircularity();
+    fjw.DoGenericSubtractionJetConstituent();
+    fjw.DoGenericSubtractionJetLeSub();
+}
+
+
   //run constituent subtractor
   if(fDoConstituentSubtraction) {
     fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
@@ -482,6 +495,67 @@ void AliEmcalJetTask::FindJets()
         jet->AddGRDenSubAt(dens[g],g);
       }
     }
+
+   if(fDoGenericSubtractionExtraJetShapes){
+     
+      std::vector<fastjet::contrib::GenericSubtractorInfo> jetAngularityInfo = fjw.GetGenSubtractorInfoJetAngularity();
+     
+      
+      Int_t na = (Int_t)jetAngularityInfo.size();
+      if(na>ij && na>0) {
+       jet->SetFirstDerivativeAngularity(jetAngularityInfo[ij].first_derivative());
+       jet->SetSecondDerivativeAngularity(jetAngularityInfo[ij].second_derivative());
+       jet->SetFirstOrderSubtractedAngularity(jetAngularityInfo[ij].first_order_subtracted());
+       jet->SetSecondOrderSubtractedAngularity(jetAngularityInfo[ij].second_order_subtracted());
+
+    }
+
+      std::vector<fastjet::contrib::GenericSubtractorInfo> jetpTDInfo = fjw.GetGenSubtractorInfoJetpTD();
+    
+      Int_t np = (Int_t)jetpTDInfo.size();
+    
+      if(np>ij && np>0) {
+       jet->SetFirstDerivativepTD(jetpTDInfo[ij].first_derivative());
+       jet->SetSecondDerivativepTD(jetpTDInfo[ij].second_derivative());
+       jet->SetFirstOrderSubtractedpTD(jetpTDInfo[ij].first_order_subtracted());
+       jet->SetSecondOrderSubtractedpTD(jetpTDInfo[ij].second_order_subtracted());
+        
+      }
+      
+      std::vector<fastjet::contrib::GenericSubtractorInfo> jetCircularityInfo = fjw.GetGenSubtractorInfoJetCircularity();
+   
+       Int_t nc = (Int_t)jetCircularityInfo.size();
+      if(nc>ij && nc>0) {
+       jet->SetFirstDerivativeCircularity(jetCircularityInfo[ij].first_derivative());
+       jet->SetSecondDerivativeCircularity(jetCircularityInfo[ij].second_derivative());
+       jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted());
+       jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted());
+
+       
+     }
+
+  
+      std::vector<fastjet::contrib::GenericSubtractorInfo> jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent();
+          Int_t nco= (Int_t)jetConstituentInfo.size();
+      if(nco>ij && nco>0) {
+       jet->SetFirstDerivativeConstituent(jetConstituentInfo[ij].first_derivative());
+       jet->SetSecondDerivativeConstituent(jetConstituentInfo[ij].second_derivative());
+       jet->SetFirstOrderSubtractedConstituent(jetConstituentInfo[ij].first_order_subtracted());
+       jet->SetSecondOrderSubtractedConstituent(jetConstituentInfo[ij].second_order_subtracted());
+        
+      }
+   std::vector<fastjet::contrib::GenericSubtractorInfo> jetLeSubInfo = fjw.GetGenSubtractorInfoJetLeSub();
+          Int_t nlsub= (Int_t)jetLeSubInfo.size();
+      if(nlsub>ij && nlsub>0) {
+       jet->SetFirstDerivativeLeSub(jetLeSubInfo[ij].first_derivative());
+       jet->SetSecondDerivativeLeSub(jetLeSubInfo[ij].second_derivative());
+       jet->SetFirstOrderSubtractedLeSub(jetLeSubInfo[ij].first_order_subtracted());
+       jet->SetSecondOrderSubtractedLeSub(jetLeSubInfo[ij].second_order_subtracted());
+        
+      }
+
+   }
+
 #endif
 
     // loop over constituents
index 3772cd18d65a01f3301272ba2c6ff7a09e4c1124..eed2b5502bd9d9803cf3d112b3dc62c0ecf7bd02 100644 (file)
@@ -96,6 +96,7 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   void                   SetGenericSubtractionJetMass(Bool_t b) { fDoGenericSubtractionJetMass = b; }
   void                   SetGenericSubtractionGR(Bool_t b, Double_t rmax = 2., Double_t dr = 0.04, Double_t ptmin = 0.) { fDoGenericSubtractionGR = b; fRMax=rmax; fDRStep=dr; fPtMinGR=ptmin;}
   void                   SetConstituentSubtraction(Bool_t b)    { fDoConstituentSubtraction = b; }
+  void                   SetGenericSubtractionExtraJetShapes(Bool_t b)        {fDoGenericSubtractionExtraJetShapes =b;}
   void                   SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
 
   UInt_t                 GetJetType()                     { return fJetType; }
@@ -172,6 +173,7 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
 
   Bool_t                 fDoGenericSubtractionJetMass; // calculate generic subtraction
   Bool_t                 fDoGenericSubtractionGR;      // calculate generic subtraction for angular structure function GR
+  Bool_t                 fDoGenericSubtractionExtraJetShapes;    //calculate generic subtraction for other jet shapes like radialmoment,pTD etc
   Bool_t                 fDoConstituentSubtraction;    // calculate constituent subtraction
   Bool_t                 fUseExternalBkg;         // use external background for generic subtractor
   TString                fRhoName;                // name of rho
index dfbb9f2fcb527101c4ca724d4f966c68d275c1ab..25fee67c4258bac7b9aefbce84ddd0135726a700 100644 (file)
@@ -38,6 +38,11 @@ class AliFJWrapper
   Bool_t                                  GetLegacyMode()            { return fLegacyMode; }
 #ifdef FASTJET_VERSION
   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetMass() const {return fGenSubtractorInfoJetMass;}
+  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> GetGenSubtractorInfoJetConstituent() const {return fGenSubtractorInfoJetConstituent;} 
+  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetLeSub() const {return fGenSubtractorInfoJetLeSub;}
   const std::vector<fastjet::PseudoJet>   GetConstituentSubtrJets() const { return fConstituentSubtrJets;  }
 #endif
   virtual std::vector<double>             GetGRNumerator()     const { return fGRNumerator;                }
@@ -48,6 +53,11 @@ class AliFJWrapper
   virtual Int_t Run();
   virtual Int_t DoGenericSubtractionJetMass();
   virtual Int_t DoGenericSubtractionGR(Int_t ijet);
+  virtual Int_t DoGenericSubtractionJetAngularity();
+  virtual Int_t DoGenericSubtractionJetpTD();
+  virtual Int_t DoGenericSubtractionJetCircularity();
+  virtual Int_t DoGenericSubtractionJetConstituent();
+  virtual Int_t DoGenericSubtractionJetLeSub();
   virtual Int_t DoConstituentSubtraction();
 
   void SetStrategy(const fastjet::Strategy &strat)                 { fStrategy = strat;  }
@@ -109,7 +119,11 @@ class AliFJWrapper
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetMass;    //!
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRNum;      //!
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRDen;      //!
-  
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetAngularity;  
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetpTD;  
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetCircularity;
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetConstituent;  
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetLeSub;
 #endif
   Bool_t                                   fLegacyMode;           //!
   Bool_t                                   fUseExternalBkg;       //!
@@ -177,6 +191,11 @@ AliFJWrapper::AliFJWrapper(const char *name, const char *title)
   , fGenSubtractorInfoJetMass ( )
   , fGenSubtractorInfoGRNum ( )
   , fGenSubtractorInfoGRDen ( )
+  , fGenSubtractorInfoJetAngularity ( )
+  , fGenSubtractorInfoJetpTD ( )
+  , fGenSubtractorInfoJetCircularity( )
+  , fGenSubtractorInfoJetConstituent ( )
+  , fGenSubtractorInfoJetLeSub ( )
 #endif
   , fLegacyMode        (false)
   , fUseExternalBkg    (false)
@@ -609,6 +628,119 @@ Int_t AliFJWrapper::DoGenericSubtractionGR(Int_t ijet) {
 #endif
   return 0;
 }
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionJetAngularity() {
+  //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
+  AliJetShapeAngularity shapeAngularity;
+
+  // clear the generic subtractor info vector
+  fGenSubtractorInfoJetAngularity.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::contrib::GenericSubtractorInfo infoAng;
+    if(fInclusiveJets[i].perp()>0.)
+      double subtracted_shape = (*fGenSubtractor)(shapeAngularity, fInclusiveJets[i], infoAng);
+    fGenSubtractorInfoJetAngularity.push_back(infoAng);
+  }
+  
+#endif
+  return 0;
+}
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionJetpTD() {
+  //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
+  AliJetShapepTD shapepTD;
+
+  // clear the generic subtractor info vector
+  fGenSubtractorInfoJetpTD.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::contrib::GenericSubtractorInfo infopTD;
+    if(fInclusiveJets[i].perp()>0.)
+      double subtracted_shape = (*fGenSubtractor)(shapepTD, fInclusiveJets[i], infopTD);
+    fGenSubtractorInfoJetpTD.push_back(infopTD);
+  }
+  
+#endif
+  return 0;
+}
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionJetCircularity() {
+  //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
+  AliJetShapeCircularity shapecircularity;
+
+  // clear the generic subtractor info vector
+  fGenSubtractorInfoJetCircularity.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::contrib::GenericSubtractorInfo infoCirc;
+    if(fInclusiveJets[i].perp()>0.)
+      double subtracted_shape = (*fGenSubtractor)(shapecircularity, fInclusiveJets[i], infoCirc);
+    fGenSubtractorInfoJetCircularity.push_back(infoCirc);
+  }
+  
+#endif
+ return 0;
+}
+
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionJetConstituent() {
+  //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
+  AliJetShapeConstituent shapeconst;
+
+  // clear the generic subtractor info vector
+   fGenSubtractorInfoJetConstituent.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::contrib::GenericSubtractorInfo infoConst;
+    if(fInclusiveJets[i].perp()>0.)
+      double subtracted_shape = (*fGenSubtractor)(shapeconst, fInclusiveJets[i], infoConst);
+    fGenSubtractorInfoJetConstituent.push_back(infoConst);
+  }
+  #endif
+  return 0;
+}
+
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionJetLeSub() {
+  //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
+  AliJetShapeLeSub shapeLeSub;
+
+  // clear the generic subtractor info vector
+  fGenSubtractorInfoJetLeSub.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::contrib::GenericSubtractorInfo infoLeSub;
+    if(fInclusiveJets[i].perp()>0.)
+      double subtracted_shape = (*fGenSubtractor)(shapeLeSub, fInclusiveJets[i], infoLeSub);
+    fGenSubtractorInfoJetLeSub.push_back(infoLeSub);
+  }
+  
+#endif
+  return 0;
+}
+//_________________________
+
 
 //_________________________________________________________________________________________________
 Int_t AliFJWrapper::DoConstituentSubtraction() {
index c263db9eb49ae5623b77ec14556351ec897e8913..3702ab405baa3d7ff66952d7819e8a401c326dc1 100644 (file)
 #endif
 
 #include "TMath.h"
-
+#include "TMatrixD.h"
+#include "TMatrixDSym.h"
+#include "TMatrixDSymEigen.h"
+#include "TVector3.h"
+#include "TVector2.h"
 using namespace std;
 
 #ifdef FASTJET_VERSION
@@ -104,6 +108,154 @@ class AliJetShapeGRDen : public fastjet::FunctionOfPseudoJet<Double32_t>
   Double_t fR;
   Double_t fDRStep;
 };
+class AliJetShapeAngularity : public fastjet::FunctionOfPseudoJet<Double32_t>{
+public:
+  virtual std::string description() const{return "Angularity:radial moment";}
+  Double32_t result(const fastjet::PseudoJet &jet) const {
+  if (!jet.has_constituents())
+    return 0; 
+  Double_t den=0.;
+  Double_t num = 0.;
+  std::vector<fastjet::PseudoJet> constits = jet.constituents();
+  for(UInt_t ic = 0; ic < constits.size(); ++ic) {
+      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 dr2 = (constits[ic].eta()-jet.eta())*(constits[ic].eta()-jet.eta()) + dphi*dphi;
+      Double_t dr = TMath::Sqrt(dr2);
+      num=num+constits[ic].perp()*dr;
+      den=den+constits[ic].perp();
+      }
+   
+  return num/den;
+  }
+
+};
+
+
+ class AliJetShapepTD : public fastjet::FunctionOfPseudoJet<Double32_t>{
+ public:
+  virtual std::string description() const{return "pTD";}
+  Double32_t result(const fastjet::PseudoJet &jet) const {
+  if (!jet.has_constituents())
+    return 0; 
+  Double_t den=0;
+  Double_t num = 0.;
+  std::vector<fastjet::PseudoJet> constits = jet.constituents();
+  for(UInt_t ic = 0; ic < constits.size(); ++ic) {
+    num=num+constits[ic].perp()*constits[ic].perp();
+    den=den+constits[ic].perp();
+      }
+   
+  return TMath::Sqrt(num)/den;
+  }};
+
+
+ class AliJetShapeConstituent : public fastjet::FunctionOfPseudoJet<Double32_t>{
+ public:
+  virtual std::string description() const{return "constituents";}
+  Double_t result(const fastjet::PseudoJet &jet) const {
+  if (!jet.has_constituents())
+    return 0; 
+
+  Double_t num = 0.;
+  std::vector<fastjet::PseudoJet> constits = jet.constituents();
+  num=1.*constits.size();  
+  return num;
+  }};
+
+
+class AliJetShapeCircularity : public fastjet::FunctionOfPseudoJet<Double32_t>{
+ public:
+  virtual std::string description() const{return "circularity denominator";}
+  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.;
+      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.);
+  
+  std::vector<fastjet::PseudoJet> constits = jet.constituents();
+  for(UInt_t ic = 0; ic < constits.size(); ++ic) {
+   
+    
+    TVector3 pp(constits[ic].px(), constits[ic].py(), constits[ic].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;
+  }};
+
+
+
+
+class AliJetShapeLeSub : public fastjet::FunctionOfPseudoJet<Double32_t>{
+ public:
+  virtual std::string description() const{return "leading mins subleading";}
+  Double32_t result(const fastjet::PseudoJet &jet) const {
+  if (!jet.has_constituents())
+    return 0;
+  std::vector<fastjet::PseudoJet> constits = jet.constituents();
+  std::vector<fastjet::PseudoJet> sortedconstits=sorted_by_pt(constits); 
+  if(sortedconstits.size()<2) return 0;
+  Double_t num=TMath::Abs(sortedconstits[0].perp()-sortedconstits[1].perp());
+  return num;
+  }};
 
 #endif
 #endif