updates to AliEmcalJet, AliEmcalJetTask and AliFJWrapper to include subtraction of...
authormverweij <marta.verweij@cern.ch>
Wed, 30 Jul 2014 15:10:50 +0000 (11:10 -0400)
committermverweij <marta.verweij@cern.ch>
Wed, 30 Jul 2014 15:48:51 +0000 (11:48 -0400)
fix warnings

PWGJE/EMCALJetTasks/AliEmcalJet.cxx
PWGJE/EMCALJetTasks/AliEmcalJet.h
PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
PWGJE/EMCALJetTasks/AliEmcalJetTask.h
PWGJE/EMCALJetTasks/AliFJWrapper.h

index 75d97c2..e299a77 100644 (file)
@@ -44,7 +44,11 @@ AliEmcalJet::AliEmcalJet() :
   fJetShapeMassSecondDer(0),
   fJetShapeMassFirstSub(0),
   fJetShapeMassSecondSub(0),
-  fLabel(-1)
+  fLabel(-1),
+  fGRNumerator(0),
+  fGRDenominator(0),
+  fGRNumeratorSub(0),
+  fGRDenominatorSub(0)
 {
   // Constructor.
 
@@ -88,7 +92,11 @@ AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
   fJetShapeMassSecondDer(0),
   fJetShapeMassFirstSub(0),
   fJetShapeMassSecondSub(0),
-  fLabel(-1)
+  fLabel(-1),
+  fGRNumerator(0),
+  fGRDenominator(0),
+  fGRNumeratorSub(0),
+  fGRDenominatorSub(0)
 {    
   // Constructor.
 
@@ -138,7 +146,11 @@ AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
   fJetShapeMassSecondDer(0),
   fJetShapeMassFirstSub(0),
   fJetShapeMassSecondSub(0),
-  fLabel(-1)
+  fLabel(-1),
+  fGRNumerator(0),
+  fGRDenominator(0),
+  fGRNumeratorSub(0),
+  fGRDenominatorSub(0)
 {
   // Constructor.
 
@@ -185,7 +197,11 @@ AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
   fJetShapeMassSecondDer(jet.fJetShapeMassSecondDer),
   fJetShapeMassFirstSub(jet.fJetShapeMassFirstSub),
   fJetShapeMassSecondSub(jet.fJetShapeMassSecondSub),
-  fLabel(jet.fLabel)
+  fLabel(jet.fLabel),
+  fGRNumerator(jet.fGRNumerator),
+  fGRDenominator(jet.fGRDenominator),
+  fGRNumeratorSub(jet.fGRNumeratorSub),
+  fGRDenominatorSub(jet.fGRDenominatorSub)
 {
   // Copy constructor.
 
@@ -237,6 +253,10 @@ AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
     fJetShapeMassFirstSub  = jet.fJetShapeMassFirstSub;
     fJetShapeMassSecondSub = jet.fJetShapeMassSecondSub;
     fLabel              = jet.fLabel;
+    fGRNumerator        = jet.fGRNumerator;
+    fGRDenominator      = jet.fGRDenominator;
+    fGRNumeratorSub     = jet.fGRNumeratorSub;
+    fGRDenominatorSub   = jet.fGRDenominatorSub;
   }
 
   return *this;
@@ -381,3 +401,10 @@ void AliEmcalJet::ResetMatching()
   fClosestJetsDist[1] = 999; 
   fMatched = 2;
 }
+
+//__________________________________________________________________________________________________
+void AliEmcalJet::PrintGR() {
+  for(Int_t i = 0; i<fGRNumerator.GetSize(); i++) {
+    Printf("num[%d] = %f",i,fGRNumerator.At(i));
+  }
+}
index f42b9d7..ff2d80b 100644 (file)
@@ -158,6 +158,21 @@ class AliEmcalJet : public AliVParticle
   Double_t          GetSecondDerivative()                     const { return fJetShapeMassSecondDer                   ; }
   Double_t          GetFirstOrderSubtracted()                 const { return fJetShapeMassFirstSub                    ; }
   Double_t          GetSecondOrderSubtracted()                const { return fJetShapeMassSecondSub                   ; }
+  
+  TArrayF GetGRNumerator()                                       const { return fGRNumerator                             ; }
+  TArrayF GetGRDenominator()                                     const { return fGRDenominator                           ; }
+  TArrayF GetGRNumeratorSub()                                    const { return fGRNumeratorSub                          ; }
+  TArrayF GetGRDenominatorSub()                                  const { return fGRDenominatorSub                        ; }
+  void              AddGRNumAt(Float_t num, Int_t idx)                 { fGRNumerator.AddAt(num, idx);                     }
+  void              AddGRDenAt(Float_t den, Int_t idx)                 { fGRDenominator.AddAt(den, idx);                   }
+  void              SetGRNumSize(UInt_t s) {fGRNumerator.Set(s); }
+  void              SetGRDenSize(UInt_t s) {fGRDenominator.Set(s); }
+
+  void              AddGRNumSubAt(Float_t num, Int_t idx)                 { fGRNumeratorSub.AddAt(num, idx);                     }
+  void              AddGRDenSubAt(Float_t den, Int_t idx)                 { fGRDenominatorSub.AddAt(den, idx);                   }
+  void              SetGRNumSubSize(UInt_t s) {fGRNumeratorSub.Set(s); }
+  void              SetGRDenSubSize(UInt_t s) {fGRDenominatorSub.Set(s); }
+  void              PrintGR();
 
  protected:
   Double32_t        fPt;                  //[0,0,12]   pt 
@@ -195,6 +210,11 @@ class AliEmcalJet : public AliVParticle
   Double_t          fJetShapeMassSecondSub; //!        result from shape derivatives for jet mass: 2nd order subtracted
   Int_t             fLabel;               //           label to inclusive jet for constituent subtracted jet    
 
+  TArrayF           fGRNumerator;         //!           array with angular structure function numerator
+  TArrayF           fGRDenominator;       //!           array with angular structure function denominator
+  TArrayF           fGRNumeratorSub;      //!           array with angular structure function numerator
+  TArrayF           fGRDenominatorSub;    //!           array with angular structure function denominator
+
   private:
     struct sort_descend
         { // sort in decreasing order
@@ -202,6 +222,6 @@ class AliEmcalJet : public AliVParticle
         bool operator () (const std::pair<Double_t, Int_t>& p1, const std::pair<Double_t, Int_t>& p2)  { return p1.first > p2.first ; }
         };
 
-  ClassDef(AliEmcalJet,14) // Emcal jet class in cylindrical coordinates
+  ClassDef(AliEmcalJet,15) // Emcal jet class in cylindrical coordinates
 };
 #endif
index 545ca91..dddc0c7 100644 (file)
@@ -68,13 +68,17 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fIsEmcPart(0),
   fLegacyMode(kFALSE),
   fCodeDebug(kFALSE),
-  fDoGenericSubtraction(kFALSE),
+  fDoGenericSubtractionJetMass(kFALSE),
+  fDoGenericSubtractionGR(kFALSE),
   fDoConstituentSubtraction(kFALSE),
   fUseExternalBkg(kFALSE),
   fRhoName(""),
   fRhomName(""),
   fRho(0),
   fRhom(0),
+  fRMax(0.4),
+  fDRStep(0.04),
+  fPtMinGR(40.),
   fJets(0),
   fJetsSub(0),
   fEvent(0),
@@ -123,13 +127,17 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fIsEmcPart(0),
   fLegacyMode(kFALSE),
   fCodeDebug(kFALSE),
-  fDoGenericSubtraction(kFALSE),
+  fDoGenericSubtractionJetMass(kFALSE),
+  fDoGenericSubtractionGR(kFALSE),
   fDoConstituentSubtraction(kFALSE),
   fUseExternalBkg(kFALSE),
   fRhoName(""),
   fRhomName(""),
   fRho(0),
   fRhom(0),
+  fRMax(0.4),
+  fDRStep(0.04),
+  fPtMinGR(40.),
   fJets(0),
   fJetsSub(0),
   fEvent(0),
@@ -394,7 +402,7 @@ void AliEmcalJetTask::FindJets()
   fjw.Run();
 
   //run generic subtractor
-  if(fDoGenericSubtraction) {
+  if(fDoGenericSubtractionJetMass) {
     fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
     fjw.DoGenericSubtractionJetMass();
   }
@@ -436,18 +444,45 @@ void AliEmcalJetTask::FindJets()
     jet->SetLabel(ij);
 
     //do generic subtraction if requested
-    if(fDoGenericSubtraction) {
 #ifdef FASTJET_VERSION
+    if(fDoGenericSubtractionJetMass) {
       std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
       Int_t n = (Int_t)jetMassInfo.size();
       if(n>ij && n>0) {
-       jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
-       jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
-       jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
-       jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
+        jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
+        jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
+        jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
+        jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
       }
-#endif
     }
+    //here do generic subtraction for angular structure function
+    Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
+    fRMax = fRadius+0.2;
+    if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
+      fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+      fjw.SetRMaxAndStep(fRMax,fDRStep);
+      fjw.DoGenericSubtractionGR(ij);
+      std::vector<double> num = fjw.GetGRNumerator();
+      std::vector<double> den = fjw.GetGRDenominator();
+      std::vector<double> nums = fjw.GetGRNumeratorSub();
+      std::vector<double> dens = fjw.GetGRDenominatorSub();
+      //pass this to AliEmcalJet
+      jet->SetGRNumSize(num.size());
+      jet->SetGRDenSize(den.size());
+      jet->SetGRNumSubSize(nums.size());
+      jet->SetGRDenSubSize(dens.size());
+      Int_t nsize = (Int_t)num.size();
+      for(Int_t g = 0; g<nsize; ++g) {
+        jet->AddGRNumAt(num[g],g);
+        jet->AddGRNumSubAt(nums[g],g);
+      }
+      Int_t dsize = (Int_t)den.size();
+      for(Int_t g = 0; g<dsize; ++g) {
+        jet->AddGRDenAt(den[g],g);
+        jet->AddGRDenSubAt(dens[g],g);
+      }
+    }
+#endif
 
     // loop over constituents
     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
index 18f8cf6..3772cd1 100644 (file)
@@ -93,7 +93,8 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
 
   void                   SetRhoName(const char *n)              { fRhoName      = n            ; }
   void                   SetRhomName(const char *n)             { fRhomName     = n            ; }
-  void                   SetGenericSubtraction(Bool_t b)        { fDoGenericSubtraction     = b; }
+  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                   SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
 
@@ -169,13 +170,18 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   Bool_t                 fLegacyMode;             //! if true, enable FJ 2.x behavior
   Bool_t                 fCodeDebug;              // use nontested code changes 
 
-  Bool_t                 fDoGenericSubtraction;   // calculate generic subtraction
-  Bool_t                 fDoConstituentSubtraction; // calculate constituent subtraction
+  Bool_t                 fDoGenericSubtractionJetMass; // calculate generic subtraction
+  Bool_t                 fDoGenericSubtractionGR;      // calculate generic subtraction for angular structure function GR
+  Bool_t                 fDoConstituentSubtraction;    // calculate constituent subtraction
   Bool_t                 fUseExternalBkg;         // use external background for generic subtractor
   TString                fRhoName;                // name of rho
   TString                fRhomName;               // name of rhom
   Double_t               fRho;                    // pT background density
   Double_t               fRhom;                   // mT background density
+  Double_t               fRMax;                   // R max for GR calculation
+  Double_t               fDRStep;                 // step width for GR calculation
+  Double_t               fPtMinGR;                // min pT for GR calculation
 
   TClonesArray          *fJets;                   //!jet collection
   TClonesArray          *fJetsSub;                //!subtracted jet collection
@@ -189,6 +195,6 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   AliEmcalJetTask(const AliEmcalJetTask&);            // not implemented
   AliEmcalJetTask &operator=(const AliEmcalJetTask&); // not implemented
 
-  ClassDef(AliEmcalJetTask, 14) // Jet producing task
+  ClassDef(AliEmcalJetTask, 15) // Jet producing task
 };
 #endif
index 533ad0b..4e102bf 100644 (file)
@@ -40,9 +40,14 @@ class AliFJWrapper
   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetMass() const {return fGenSubtractorInfoJetMass;}
   const std::vector<fastjet::PseudoJet>   GetConstituentSubtrJets() const { return fConstituentSubtrJets;  }
 #endif
+  virtual std::vector<double>             GetGRNumerator()     const { return fGRNumerator;                }
+  virtual std::vector<double>             GetGRDenominator()   const { return fGRDenominator;              }
+  virtual std::vector<double>             GetGRNumeratorSub()  const { return fGRNumeratorSub;             }
+  virtual std::vector<double>             GetGRDenominatorSub()const { return fGRDenominatorSub;           }
 
   virtual Int_t Run();
   virtual Int_t DoGenericSubtractionJetMass();
+  virtual Int_t DoGenericSubtractionGR(Int_t ijet);
   virtual Int_t DoConstituentSubtraction();
 
   void SetStrategy(const fastjet::Strategy &strat)                 { fStrategy = strat;  }
@@ -65,6 +70,7 @@ class AliFJWrapper
   void SetLegacyMode (Bool_t mode)      { fLegacyMode ^= mode; }
   void SetLegacyFJ();
   void SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
+  void SetRMaxAndStep(Double_t rmax, Double_t dr) {fRMax = rmax; fDRStep = dr; }
 
  protected:
   TString                                fName;               //!
@@ -97,15 +103,24 @@ class AliFJWrapper
   Double_t                               fMedUsedForBgSub;    //!
   Bool_t                                 fUseArea4Vector;     //!
 #ifdef FASTJET_VERSION
-  fastjet::JetMedianBackgroundEstimator   *fBkrdEstimator;        //!
+  fastjet::JetMedianBackgroundEstimator   *fBkrdEstimator;    //!
   //from contrib package
-  fastjet::contrib::GenericSubtractor     *fGenSubtractor;        //!
+  fastjet::contrib::GenericSubtractor     *fGenSubtractor;    //!
   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetMass;    //!
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRNum;      //!
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRDen;      //!
+  
 #endif
   Bool_t                                   fLegacyMode;           //!
   Bool_t                                   fUseExternalBkg;       //!
   Double_t                                 fRho;                  //  pT background density
   Double_t                                 fRhom;                 //  mT background density
+  Double_t                                 fRMax;             //!
+  Double_t                                 fDRStep;           //!
+  std::vector<double>                      fGRNumerator;      //!
+  std::vector<double>                      fGRDenominator;    //!
+  std::vector<double>                      fGRNumeratorSub;   //!
+  std::vector<double>                      fGRDenominatorSub; //!
 
   virtual void   SubtractBackground(const Double_t median_pt = -1);
 
@@ -165,6 +180,14 @@ AliFJWrapper::AliFJWrapper(const char *name, const char *title)
   , fUseExternalBkg    (false)
   , fRho               (0)
   , fRhom              (0)
+  , fRMax(2.)
+  , fDRStep(0.04)
+  , fGenSubtractorInfoGRNum ( )
+  , fGenSubtractorInfoGRDen ( )
+  , fGRNumerator()
+  , fGRDenominator()
+  , fGRNumeratorSub()
+  , fGRDenominatorSub()
 {
   // Constructor.
 }
@@ -549,6 +572,45 @@ Int_t AliFJWrapper::DoGenericSubtractionJetMass() {
 }
 
 //_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionGR(Int_t ijet) {
+  //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);
+
+  if(ijet>fInclusiveJets.size()) return 0;
+
+  fGRNumerator.clear();
+  fGRDenominator.clear();
+  fGRNumeratorSub.clear();
+  fGRDenominatorSub.clear();
+
+  // Define jet shape
+  for(Double_t r = 0.; r<fRMax; r+=fDRStep) {
+    AliJetShapeGRNum shapeGRNum(r,fDRStep);
+    AliJetShapeGRDen shapeGRDen(r,fDRStep);
+
+    // clear the generic subtractor info vector
+    fGenSubtractorInfoGRNum.clear();
+    fGenSubtractorInfoGRDen.clear();
+    fj::contrib::GenericSubtractorInfo infoNum;
+    fj::contrib::GenericSubtractorInfo infoDen;
+    if(fInclusiveJets[ijet].perp()>0.) {
+      double sub_num = (*fGenSubtractor)(shapeGRNum, fInclusiveJets[ijet], infoNum);
+      double sub_den = (*fGenSubtractor)(shapeGRDen, fInclusiveJets[ijet], infoDen);
+    }
+    fGenSubtractorInfoGRNum.push_back(infoNum);
+    fGenSubtractorInfoGRDen.push_back(infoDen);
+    fGRNumerator.push_back(infoNum.unsubtracted());
+    fGRDenominator.push_back(infoDen.unsubtracted());
+    fGRNumeratorSub.push_back(infoNum.second_order_subtracted());
+    fGRDenominatorSub.push_back(infoDen.second_order_subtracted());
+  }
+#endif
+  return 0;
+}
+
+//_________________________________________________________________________________________________
 Int_t AliFJWrapper::DoConstituentSubtraction() {
   //Do constituent subtraction
 #ifdef FASTJET_VERSION