[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetShape.h

#ifndef AliJetShape_H
#define AliJetShape_H

// $Id$

#include <vector>
#include <TString.h>
#include "fastjet/PseudoJet.hh"
#ifdef FASTJET_VERSION
#include "fastjet/FunctionOfPseudoJet.hh"
#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
//________________________________________________________________________
class AliJetShapeMass : public fastjet::FunctionOfPseudoJet<Double32_t>
{
 public:
  virtual std::string description() const{return "jet mass";}
  Double32_t result(const fastjet::PseudoJet &jet) const{ return jet.m();}
};

//________________________________________________________________________
class AliJetShapeGRNum : public fastjet::FunctionOfPseudoJet<Double32_t>
{
 public:
  // default ctor
  AliJetShapeGRNum(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
  virtual std::string description() const{return "Numerator angular structure function";}
  //  static Int_t GetBin(Double_t x) {Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
  Double32_t result(const fastjet::PseudoJet &jet) const {
    if (!jet.has_constituents())
      return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");

  Double_t A = 0.;
  std::vector<fastjet::PseudoJet> constits = jet.constituents();
  for(UInt_t ic = 0; ic < constits.size(); ++ic) {
    /* Int_t uid = constits[ic].user_index(); */
    /* if (uid == -1) //skip ghost particle */
    /*   continue; */
    for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
      /* Int_t uid = constits[jc].user_index(); */
      /* if (uid == -1) //skip ghost particle */
      /* 	continue; */
      Double_t dphi = constits[ic].phi()-constits[jc].phi();
      if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
      if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
      Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
      if(dr2>0.) {
	Double_t dr = TMath::Sqrt(dr2);
	Double_t x = fR-dr;
	//noisy function
	Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
	A += constits[ic].perp()*constits[jc].perp()*dr2*noise;
      }
    }
  }
  return A;
  }

 protected:
  Double_t fR;
  Double_t fDRStep;
};

//________________________________________________________________________
class AliJetShapeGRDen : public fastjet::FunctionOfPseudoJet<Double32_t>
{
 public:
  // default ctor
  AliJetShapeGRDen(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
  virtual std::string description() const{return "Denominator angular structure function";}
  Double32_t result(const fastjet::PseudoJet &jet) const {
  if (!jet.has_constituents())
    return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");

  Double_t A = 0.;
  std::vector<fastjet::PseudoJet> constits = jet.constituents();
  for(UInt_t ic = 0; ic < constits.size(); ++ic) {
    /* Int_t uid = constits[ic].user_index(); */
    /* if (uid == -1) //skip ghost particle */
    /*   continue; */
    for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
      /* Int_t uid = constits[jc].user_index(); */
      /* if (uid == -1) //skip ghost particle */
      /* 	continue; */
      Double_t dphi = constits[ic].phi()-constits[jc].phi();
      if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
      if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
      Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
      if(dr2>0.) {
	Double_t dr = TMath::Sqrt(dr2);
	Double_t x = fR-dr;
	//error function
	Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
	A += constits[ic].perp()*constits[jc].perp()*dr2*erf;
      }
    }
  }
  return A;
  }

 protected:
  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;
    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 AliJetShapeSigma2 : public fastjet::FunctionOfPseudoJet<Double32_t>{
 public:
  virtual std::string description() const{return "cms sigma2";}
  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.;

    std::vector<fastjet::PseudoJet> constits = jet.constituents();
    for(UInt_t ic = 0; ic < constits.size(); ++ic) {
      Double_t ppt=constits[ic].perp();
      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 deta = constits[ic].eta()-jet.eta();
      mxx += ppt*ppt*deta*deta;
      myy += ppt*ppt*dphi*dphi;
      mxy -= ppt*ppt*deta*dphi;
      nc++;
      sump2 += ppt*ppt;
    }
    if(nc<2) return 0;
    if(sump2==0) return 0;
    // Sphericity Matrix
    Double_t ele[4] = {mxx , mxy, mxy, myy };
    TMatrixDSym m0(2,ele);
      
    // Find eigenvectors
    TMatrixDSymEigen m(m0);
    TVectorD eval(2);
    TMatrixD evecm = m.GetEigenVectors();
    eval  = m.GetEigenValues();
    // Largest eigenvector
    int jev = 0;
    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 sigma2=0;
    if(jev==1) sigma2=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
    if(jev==0) sigma2=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
    
    return sigma2;
  }
};

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