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 #include "TMatrixD.h"
  15 #include "TMatrixDSym.h"
  16 #include "TMatrixDSymEigen.h"
  17 #include "TVector3.h"
  18 #include "TVector2.h"
  19 using namespace std;
  20
  21 #ifdef FASTJET_VERSION
  22 class AliJetShapeMass : public fastjet::FunctionOfPseudoJet<Double32_t>
  23 {
  24  public:
  25   virtual std::string description() const{return "jet mass";}
  26   Double32_t result(const fastjet::PseudoJet &jet) const{ return jet.m();}
  27 };
  28
  29 class AliJetShapeGRNum : public fastjet::FunctionOfPseudoJet<Double32_t>
  30 {
  31  public:
  32   // default ctor
  33   AliJetShapeGRNum(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
  34   virtual std::string description() const{return "Numerator angular structure function";}
  35   //  static Int_t GetBin(Double_t x) {Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
  36   Double32_t result(const fastjet::PseudoJet &jet) const {
  37     if (!jet.has_constituents())
  38       return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");
  39
  40   Double_t A = 0.;
  41   std::vector<fastjet::PseudoJet> constits = jet.constituents();
  42   for(UInt_t ic = 0; ic < constits.size(); ++ic) {
  43     /* Int_t uid = constits[ic].user_index(); */
  44     /* if (uid == -1) //skip ghost particle */
  45     /*   continue; */
  46     for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
  47       /* Int_t uid = constits[jc].user_index(); */
  48       /* if (uid == -1) //skip ghost particle */
  49       /*        continue; */
  50       Double_t dphi = constits[ic].phi()-constits[jc].phi();
  51       if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
  52       if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
  53       Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
  54       if(dr2>0.) {
  55         Double_t dr = TMath::Sqrt(dr2);
  56         Double_t x = fR-dr;
  57         //noisy function
  58         Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
  59         A += constits[ic].perp()*constits[jc].perp()*dr2*noise;
  60       }
  61     }
  62   }
  63   return A;
  64   }
  65
  66  protected:
  67   Double_t fR;
  68   Double_t fDRStep;
  69 };
  70
  71 class AliJetShapeGRDen : public fastjet::FunctionOfPseudoJet<Double32_t>
  72 {
  73  public:
  74   // default ctor
  75   AliJetShapeGRDen(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
  76   virtual std::string description() const{return "Denominator angular structure function";}
  77   Double32_t result(const fastjet::PseudoJet &jet) const {
  78   if (!jet.has_constituents())
  79     return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");
  80
  81   Double_t A = 0.;
  82   std::vector<fastjet::PseudoJet> constits = jet.constituents();
  83   for(UInt_t ic = 0; ic < constits.size(); ++ic) {
  84     /* Int_t uid = constits[ic].user_index(); */
  85     /* if (uid == -1) //skip ghost particle */
  86     /*   continue; */
  87     for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
  88       /* Int_t uid = constits[jc].user_index(); */
  89       /* if (uid == -1) //skip ghost particle */
  90       /*        continue; */
  91       Double_t dphi = constits[ic].phi()-constits[jc].phi();
  92       if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
  93       if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
  94       Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
  95       if(dr2>0.) {
  96         Double_t dr = TMath::Sqrt(dr2);
  97         Double_t x = fR-dr;
  98         //error function
  99         Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
 100         A += constits[ic].perp()*constits[jc].perp()*dr2*erf;
 101       }
 102     }
 103   }
 104   return A;
 105   }
 106
 107  protected:
 108   Double_t fR;
 109   Double_t fDRStep;
 110 };
 111
 112 class AliJetShapeAngularity : public fastjet::FunctionOfPseudoJet<Double32_t>{
 113 public:
 114   virtual std::string description() const{return "Angularity:radial moment";}
 115   Double32_t result(const fastjet::PseudoJet &jet) const {
 116     if (!jet.has_constituents())
 117       return 0;
 118     Double_t den=0.;
 119     Double_t num = 0.;
 120     std::vector<fastjet::PseudoJet> constits = jet.constituents();
 121     for(UInt_t ic = 0; ic < constits.size(); ++ic) {
 122       Double_t dphi = constits[ic].phi()-jet.phi();
 123       if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
 124       if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
 125       Double_t dr2 = (constits[ic].eta()-jet.eta())*(constits[ic].eta()-jet.eta()) + dphi*dphi;
 126       Double_t dr = TMath::Sqrt(dr2);
 127       num=num+constits[ic].perp()*dr;
 128       den=den+constits[ic].perp();
 129     }
 130     return num/den;
 131   }
 132 };
 133
 134
 135 class AliJetShapepTD : public fastjet::FunctionOfPseudoJet<Double32_t>{
 136  public:
 137   virtual std::string description() const{return "pTD";}
 138   Double32_t result(const fastjet::PseudoJet &jet) const {
 139     if (!jet.has_constituents())
 140       return 0;
 141     Double_t den=0;
 142     Double_t num = 0.;
 143     std::vector<fastjet::PseudoJet> constits = jet.constituents();
 144     for(UInt_t ic = 0; ic < constits.size(); ++ic) {
 145       num=num+constits[ic].perp()*constits[ic].perp();
 146       den=den+constits[ic].perp();
 147     }
 148     return TMath::Sqrt(num)/den;
 149   }
 150 };
 151
 152 class AliJetShapeConstituent : public fastjet::FunctionOfPseudoJet<Double32_t>{
 153  public:
 154   virtual std::string description() const{return "constituents";}
 155   Double_t result(const fastjet::PseudoJet &jet) const {
 156     if (!jet.has_constituents())
 157       return 0;
 158     Double_t num = 0.;
 159     std::vector<fastjet::PseudoJet> constits = jet.constituents();
 160     num=1.*constits.size();
 161     return num;
 162   }
 163 };
 164
 165 class AliJetShapeCircularity : public fastjet::FunctionOfPseudoJet<Double32_t>{
 166  public:
 167   virtual std::string description() const{return "circularity denominator";}
 168   Double32_t result(const fastjet::PseudoJet &jet) const {
 169     if (!jet.has_constituents())
 170       return 0;
 171     Double_t mxx    = 0.;
 172     Double_t myy    = 0.;
 173     Double_t mxy    = 0.;
 174     int  nc     = 0;
 175     Double_t sump2  = 0.;
 176     Double_t pxjet=jet.px();
 177     Double_t pyjet=jet.py();
 178     Double_t pzjet=jet.pz();
 179
 180     //2 general normalized vectors perpendicular to the jet
 181     TVector3  ppJ1(pxjet, pyjet, pzjet);
 182     TVector3  ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
 183     ppJ3.SetMag(1.);
 184     TVector3  ppJ2(-pyjet, pxjet, 0);
 185     ppJ2.SetMag(1.);
 186
 187     std::vector<fastjet::PseudoJet> constits = jet.constituents();
 188     for(UInt_t ic = 0; ic < constits.size(); ++ic) {
 189       TVector3 pp(constits[ic].px(), constits[ic].py(), constits[ic].pz());
 190       //local frame
 191       TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
 192       TVector3 pPerp = pp - pLong;
 193       //projection onto the two perpendicular vectors defined above
 194       Float_t ppjX = pPerp.Dot(ppJ2);
 195       Float_t ppjY = pPerp.Dot(ppJ3);
 196       Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
 197       if(ppjT<=0) return 0;
 198       mxx += (ppjX * ppjX / ppjT);
 199       myy += (ppjY * ppjY / ppjT);
 200       mxy += (ppjX * ppjY / ppjT);
 201       nc++;
 202       sump2 += ppjT;
 203     }
 204     if(nc<2) return 0;
 205     if(sump2==0) return 0;
 206     // Sphericity Matrix
 207     Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
 208     TMatrixDSym m0(2,ele);
 209
 210     // Find eigenvectors
 211     TMatrixDSymEigen m(m0);
 212     TVectorD eval(2);
 213     TMatrixD evecm = m.GetEigenVectors();
 214     eval  = m.GetEigenValues();
 215     // Largest eigenvector
 216     int jev = 0;
 217     if (eval[0] < eval[1]) jev = 1;
 218     TVectorD evec0(2);
 219     // Principle axis
 220     evec0 = TMatrixDColumn(evecm, jev);
 221     Double_t compx=evec0[0];
 222     Double_t compy=evec0[1];
 223     TVector2 evec(compx, compy);
 224     Double_t circ=0;
 225     if(jev==1) circ=2*eval[0];
 226     if(jev==0) circ=2*eval[1];
 227
 228     return circ;
 229   }
 230 };
 231
 232 class AliJetShapeLeSub : public fastjet::FunctionOfPseudoJet<Double32_t>{
 233  public:
 234   virtual std::string description() const{return "leading mins subleading";}
 235   Double32_t result(const fastjet::PseudoJet &jet) const {
 236     if (!jet.has_constituents())
 237       return 0;
 238     std::vector<fastjet::PseudoJet> constits = jet.constituents();
 239     std::vector<fastjet::PseudoJet> sortedconstits=sorted_by_pt(constits);
 240     if(sortedconstits.size()<2) return 0;
 241     Double_t num=TMath::Abs(sortedconstits[0].perp()-sortedconstits[1].perp());
 242     return num;
 243   }
 244 };
 245
 246 #endif
 247 #endif