8 #include "fastjet/PseudoJet.hh"
10 #include "fastjet/FunctionOfPseudoJet.hh"
15 #include "TMatrixDSym.h"
16 #include "TMatrixDSymEigen.h"
21 #ifdef FASTJET_VERSION
22 class AliJetShapeMass : public fastjet::FunctionOfPseudoJet<Double32_t>
25 virtual std::string description() const{return "jet mass";}
26 Double32_t result(const fastjet::PseudoJet &jet) const{ return jet.m();}
29 class AliJetShapeGRNum : public fastjet::FunctionOfPseudoJet<Double32_t>
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.");
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 */
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 */
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;
55 Double_t dr = TMath::Sqrt(dr2);
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;
71 class AliJetShapeGRDen : public fastjet::FunctionOfPseudoJet<Double32_t>
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.");
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 */
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 */
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;
96 Double_t dr = TMath::Sqrt(dr2);
99 Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
100 A += constits[ic].perp()*constits[jc].perp()*dr2*erf;
112 class AliJetShapeAngularity : public fastjet::FunctionOfPseudoJet<Double32_t>{
114 virtual std::string description() const{return "Angularity:radial moment";}
115 Double32_t result(const fastjet::PseudoJet &jet) const {
116 if (!jet.has_constituents())
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();
135 class AliJetShapepTD : public fastjet::FunctionOfPseudoJet<Double32_t>{
137 virtual std::string description() const{return "pTD";}
138 Double32_t result(const fastjet::PseudoJet &jet) const {
139 if (!jet.has_constituents())
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();
148 return TMath::Sqrt(num)/den;
152 class AliJetShapeConstituent : public fastjet::FunctionOfPseudoJet<Double32_t>{
154 virtual std::string description() const{return "constituents";}
155 Double_t result(const fastjet::PseudoJet &jet) const {
156 if (!jet.has_constituents())
159 std::vector<fastjet::PseudoJet> constits = jet.constituents();
160 num=1.*constits.size();
165 class AliJetShapeCircularity : public fastjet::FunctionOfPseudoJet<Double32_t>{
167 virtual std::string description() const{return "circularity denominator";}
168 Double32_t result(const fastjet::PseudoJet &jet) const {
169 if (!jet.has_constituents())
176 Double_t pxjet=jet.px();
177 Double_t pyjet=jet.py();
178 Double_t pzjet=jet.pz();
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);
184 TVector3 ppJ2(-pyjet, pxjet, 0);
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());
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);
205 if(sump2==0) return 0;
207 Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
208 TMatrixDSym m0(2,ele);
211 TMatrixDSymEigen m(m0);
213 TMatrixD evecm = m.GetEigenVectors();
214 eval = m.GetEigenValues();
215 // Largest eigenvector
217 if (eval[0] < eval[1]) jev = 1;
220 evec0 = TMatrixDColumn(evecm, jev);
221 Double_t compx=evec0[0];
222 Double_t compy=evec0[1];
223 TVector2 evec(compx, compy);
225 if(jev==1) circ=2*eval[0];
226 if(jev==0) circ=2*eval[1];
232 class AliJetShapeLeSub : public fastjet::FunctionOfPseudoJet<Double32_t>{
234 virtual std::string description() const{return "leading mins subleading";}
235 Double32_t result(const fastjet::PseudoJet &jet) const {
236 if (!jet.has_constituents())
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());