PWGJE/EMCALJetTasks/AliJetShape.h

   1 #ifndef AliJetShape_H
   2 #define AliJetShape_H
   3
   4 // $Id$
   5
   6 #include <vector>
   7 #include <TString.h>
   8 #include "fastjet/PseudoJet.hh"
   9 #ifdef FASTJET_VERSION
  10 #include "fastjet/FunctionOfPseudoJet.hh"
  11 #endif
  12
  13 #include "TMath.h"
  14
  15 using namespace std;
  16
  17 /*
  18 const Int_t knbins = 130;
  19 const Double_t kxmin = -0.8;
  20 const Double_t kxmax = 0.5;
  21 static Int_t GetBin(Double_t x) {Double_t mdx = (kxmax-kxmin)/(double)knbins; Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
  22 static Double_t gfcn[knbins] = {0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000001, 0.000005, 0.000020, 0.000069, 0.000226, 0.000696, 0.002014, 0.005473, 0.013977, 0.033528, 0.075556, 0.159953, 0.318105, 0.594298, 1.043025, 1.719657, 2.663457, 3.875307, 5.296916, 6.801375, 8.204024, 9.296377, 9.895942, 9.895942, 9.296377, 8.204024, 6.801375, 5.296916, 3.875307, 2.663457, 1.719657, 1.043025, 0.594298, 0.318105, 0.159953, 0.075556, 0.033528, 0.013977, 0.005473, 0.002014, 0.000696, 0.000226, 0.000069, 0.000020, 0.000005, 0.000001, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000};
  23
  24 static Double_t efcn[knbins] = {0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000001, 0.000002, 0.000006, 0.000019, 0.000053, 0.000144, 0.000369, 0.000889, 0.002020, 0.004332, 0.008774, 0.016793, 0.030396, 0.052081, 0.084566, 0.130295, 0.190787, 0.265986, 0.353830, 0.450262, 0.549738, 0.646170, 0.734014, 0.809213, 0.869705, 0.915434, 0.947919, 0.969604, 0.983207, 0.991226, 0.995668, 0.997980, 0.999111, 0.999631, 0.999856, 0.999947, 0.999981, 0.999994, 0.999998, 0.999999, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000};
  25 */
  26
  27 #ifdef FASTJET_VERSION
  28 class AliJetShapeMass : public fastjet::FunctionOfPseudoJet<Double32_t>
  29 {
  30  public:
  31   virtual std::string description() const{return "jet mass";}
  32   Double32_t result(const fastjet::PseudoJet &jet) const{ return jet.m();}
  33 };
  34
  35 class AliJetShapeGRNum : public fastjet::FunctionOfPseudoJet<Double32_t>
  36 {
  37  public:
  38   // default ctor
  39   AliJetShapeGRNum(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
  40   virtual std::string description() const{return "Numerator angular structure function";}
  41   //  static Int_t GetBin(Double_t x) {Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
  42   Double32_t result(const fastjet::PseudoJet &jet) const {
  43     if (!jet.has_constituents()) {
  44       Printf("Angular structure can only be applied on jets for which the constituents are known.");
  45       return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");
  46     }
  47   Double_t A = 0.;
  48   std::vector<fastjet::PseudoJet> constits = jet.constituents();
  49   for(UInt_t ic = 0; ic < constits.size(); ++ic) {
  50     /* Int_t uid = constits[ic].user_index(); */
  51     /* if (uid == -1) //skip ghost particle */
  52     /*   continue; */
  53     for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
  54       /* Int_t uid = constits[jc].user_index(); */
  55       /* if (uid == -1) //skip ghost particle */
  56       /*        continue; */
  57       Double_t dphi = constits[ic].phi()-constits[jc].phi();
  58       if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
  59       if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
  60       Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
  61       if(dr2>0.) {
  62         Double_t dr = TMath::Sqrt(dr2);
  63         Double_t x = fR-dr;
  64         //      Printf("Calculating numerator between %d-%d for r: %f dr: %f",ic,jc,fR,dr);
  65         //noisy function
  66         Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
  67         /* Int_t bin = GetBin(x); */
  68         /* Double_t noise = 0; */
  69         /* if(bin<0) noise = gfcn[0]; */
  70         /* else if(bin>=knbins) noise = gfcn[knbins-1]; */
  71         /* else  noise = gfcn[bin]; */
  72         /* Printf("fR: %f dr: %f x=%f bin=%d noise: %f  %f",fR,dr,x,bin,noise,TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep)); */
  73         A += constits[ic].perp()*constits[jc].perp()*dr2*noise;
  74       }
  75     }
  76   }
  77   return A;
  78   }
  79
  80  protected:
  81   Double_t fR;
  82   Double_t fDRStep;
  83 };
  84
  85 class AliJetShapeGRDen : public fastjet::FunctionOfPseudoJet<Double32_t>
  86 {
  87  public:
  88   // default ctor
  89   AliJetShapeGRDen(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
  90   virtual std::string description() const{return "Denominator angular structure function";}
  91   Double32_t result(const fastjet::PseudoJet &jet) const {
  92   if (!jet.has_constituents())
  93     return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");
  94
  95   Double_t A = 0.;
  96   std::vector<fastjet::PseudoJet> constits = jet.constituents();
  97   for(UInt_t ic = 0; ic < constits.size(); ++ic) {
  98     /* Int_t uid = constits[ic].user_index(); */
  99     /* if (uid == -1) //skip ghost particle */
 100     /*   continue; */
 101     for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
 102       /* Int_t uid = constits[jc].user_index(); */
 103       /* if (uid == -1) //skip ghost particle */
 104       /*        continue; */
 105       Double_t dphi = constits[ic].phi()-constits[jc].phi();
 106       if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
 107       if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 108       Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
 109       if(dr2>0.) {
 110         Double_t dr = TMath::Sqrt(dr2);
 111         Double_t x = fR-dr;
 112         //error function
 113         Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
 114         /* Int_t bin = GetBin(x); */
 115         /* Double_t erf = 0; */
 116         /* if(bin<0) erf = efcn[0]; */
 117         /* else if(bin>=knbins) erf = efcn[knbins-1]; */
 118         /* else  erf = efcn[bin]; */
 119         A += constits[ic].perp()*constits[jc].perp()*dr2*erf;
 120       }
 121     }
 122   }
 123   return A;
 124   }
 125
 126  protected:
 127   Double_t fR;
 128   Double_t fDRStep;
 129 };
 130
 131 #endif
 132 #endif