PWGJE/EMCALJetTasks/AliJetShape.cxx

   1 // $Id$
   2 //
   3 // Author: M. Verweij
   4
   5 #define AliJetShape_CXX
   6 #include "AliJetShape.h"
   7
   8 /*
   9 Double32_t AliJetShapeGRNum::result(const fastjet::PseudoJet &jet) const {
  10
  11 #ifdef FASTJET_VERSION
  12   if (!jet.has_constituents())
  13     throw Error("Angular structure can only be applied on jets for which the constituents are known.");
  14
  15   Double_t A = 0.;
  16   vector<PseudoJet> constits = jet.constituents();
  17   for (vector<PseudoJet>::iterator ci = constits.begin(); ci!=constits.end(); ci++){
  18     for (vector<PseudoJet>::iterator cj = ci+1; cj!=constits.end(); cj++){
  19       Double_t dphi = ci.phi()-cj.phi();
  20       if(dphi<0.) dphi+=TMath::TwoPi();
  21       if(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
  22       Double_t dr2 = (ci.eta()-cj.eta())*(ci.eta()-cj.eta()) + dphi*dphi;
  23       if(dr2>0.) {
  24         Double_t dr = TMath::Sqrt(dr2);
  25         Double_t x = fR-dr;
  26         //noisy function
  27         Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
  28         //error function
  29         // Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
  30         A += ci.perp()*cj.perp()*dr2*noise;
  31       }
  32     }
  33   }
  34   return A;
  35 #endif
  36   return 0;
  37 }
  38 */
  39
  40 Double32_t AliJetShapeGRDen::result(const fastjet::PseudoJet &jet) const {
  41
  42 #ifdef FASTJET_VERSION
  43   if (!jet.has_constituents())
  44     throw Error("Angular structure can only be applied on jets for which the constituents are known.");
  45
  46   Double_t A = 0.;
  47   vector<PseudoJet> constits = jet.constituents();
  48   for (vector<PseudoJet>::iterator ci = constits.begin(); ci!=constits.end(); ci++){
  49     for (vector<PseudoJet>::iterator cj = ci+1; cj!=constits.end(); cj++){
  50       Double_t dphi = ci.phi()-cj.phi();
  51       if(dphi<0.) dphi+=TMath::TwoPi();
  52       if(dphi>TMath::TwoPi()) dphi-=TMath::TwoPi();
  53       Double_t dr2 = (ci.eta()-cj.eta())*(ci.eta()-cj.eta()) + dphi*dphi;
  54       if(dr2>0.) {
  55         Double_t dr = TMath::Sqrt(dr2);
  56         Double_t x = fR-dr;
  57         //error function
  58          Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
  59         A += ci.perp()*cj.perp()*dr2*erf;
  60       }
  61     }
  62   }
  63   return A;
  64 #endif
  65   return 0;
  66 }