AliPhysics  58f3d52 (58f3d52)
 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  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(1,0,0.2,1.0,0.4,jet,0);
224  fFastjetWrapper->Clear();
225  delete fFastjetWrapper;
226  return Result;
227 }
228 
229 
230 Double32_t AliJetShape2subjettiness_kt::result(const fastjet::PseudoJet &jet) const {
231  if (!jet.has_constituents())
232  return 0;
233  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
234  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,0,0.2,1.0,0.4,jet,0);
235  fFastjetWrapper->Clear();
236  delete fFastjetWrapper;
237  return Result;
238 }
239 
240 Double32_t AliJetShape3subjettiness_kt::result(const fastjet::PseudoJet &jet) const {
241  if (!jet.has_constituents())
242  return 0;
243  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
244  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(3,0,0.2,1.0,0.4,jet,0);
245  fFastjetWrapper->Clear();
246  delete fFastjetWrapper;
247  return Result;
248 }
249 
250 Double32_t AliJetShapeOpeningAngle_kt::result(const fastjet::PseudoJet &jet) const {
251  if (!jet.has_constituents())
252  return 0;
253  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
254  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,0,0.2,1.0,0.4,jet,1);
255  fFastjetWrapper->Clear();
256  delete fFastjetWrapper;
257  return Result;
258 }
259 
260 
261 Double32_t AliJetShape1subjettiness_ca::result(const fastjet::PseudoJet &jet) const {
262  if (!jet.has_constituents())
263  return 0;
264  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
265  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(1,1,0.2,1.0,0.4,jet,0);
266  fFastjetWrapper->Clear();
267  delete fFastjetWrapper;
268  return Result;
269 }
270 
271 
272 Double32_t AliJetShape2subjettiness_ca::result(const fastjet::PseudoJet &jet) const {
273  if (!jet.has_constituents())
274  return 0;
275  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
276  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,1,0.2,1.0,0.4,jet,0);
277  fFastjetWrapper->Clear();
278  delete fFastjetWrapper;
279  return Result;
280 }
281 
282 Double32_t AliJetShapeOpeningAngle_ca::result(const fastjet::PseudoJet &jet) const {
283  if (!jet.has_constituents())
284  return 0;
285  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
286  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,1,0.2,1.0,0.4,jet,1);
287  fFastjetWrapper->Clear();
288  delete fFastjetWrapper;
289  return Result;
290 }
291 
292 
293 Double32_t AliJetShape1subjettiness_akt02::result(const fastjet::PseudoJet &jet) const {
294  if (!jet.has_constituents())
295  return 0;
296  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
297  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(1,2,0.2,1.0,0.4,jet,0);
298  fFastjetWrapper->Clear();
299  delete fFastjetWrapper;
300  return Result;
301 }
302 
303 
304 Double32_t AliJetShape2subjettiness_akt02::result(const fastjet::PseudoJet &jet) const {
305  if (!jet.has_constituents())
306  return 0;
307  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
308  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,2,0.2,1.0,0.4,jet,0);
309  fFastjetWrapper->Clear();
310  delete fFastjetWrapper;
311  return Result;
312 }
313 
314 
315 Double32_t AliJetShapeOpeningAngle_akt02::result(const fastjet::PseudoJet &jet) const {
316  if (!jet.has_constituents())
317  return 0;
318  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
319  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,2,0.2,1.0,0.4,jet,1);
320  fFastjetWrapper->Clear();
321  delete fFastjetWrapper;
322  return Result;
323 }
324 
325 
326 Double32_t AliJetShape1subjettiness_casd::result(const fastjet::PseudoJet &jet) const {
327  if (!jet.has_constituents())
328  return 0;
329  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
330  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(1,1,0.2,1.0,0.4,jet,0,0,0.0,0.1,1);
331  fFastjetWrapper->Clear();
332  delete fFastjetWrapper;
333  return Result;
334 }
335 
336 
337 Double32_t AliJetShape2subjettiness_casd::result(const fastjet::PseudoJet &jet) const {
338  if (!jet.has_constituents())
339  return 0;
340  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
341  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,1,0.2,1.0,0.4,jet,0,0,0.0,0.1,1);
342  fFastjetWrapper->Clear();
343  delete fFastjetWrapper;
344  return Result;
345 }
346 
347 
348 Double32_t AliJetShapeOpeningAngle_casd::result(const fastjet::PseudoJet &jet) const {
349  if (!jet.has_constituents())
350  return 0;
351  AliFJWrapper *fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
352  Double32_t Result = fFastjetWrapper->AliFJWrapper::NSubjettinessDerivativeSub(2,1,0.2,1.0,0.4,jet,1,0,0.0,0.1,1);
353  fFastjetWrapper->Clear();
354  delete fFastjetWrapper;
355  return Result;
356 }
357 
358 
359 #endif
360 
double Double_t
Definition: External.C:58
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
virtual void Clear(const Option_t *="")