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