]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliJetShape.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add jet mass to tagged jet correlation
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetShape.h
CommitLineData
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"
14
15using namespace std;
16
17/*
18const Int_t knbins = 130;
19const Double_t kxmin = -0.8;
20const Double_t kxmax = 0.5;
21static Int_t GetBin(Double_t x) {Double_t mdx = (kxmax-kxmin)/(double)knbins; Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
22static Double_t gfcn[knbins] = {0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000001, 0.000005, 0.000020, 0.000069, 0.000226, 0.000696, 0.002014, 0.005473, 0.013977, 0.033528, 0.075556, 0.159953, 0.318105, 0.594298, 1.043025, 1.719657, 2.663457, 3.875307, 5.296916, 6.801375, 8.204024, 9.296377, 9.895942, 9.895942, 9.296377, 8.204024, 6.801375, 5.296916, 3.875307, 2.663457, 1.719657, 1.043025, 0.594298, 0.318105, 0.159953, 0.075556, 0.033528, 0.013977, 0.005473, 0.002014, 0.000696, 0.000226, 0.000069, 0.000020, 0.000005, 0.000001, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000};
23
24static Double_t efcn[knbins] = {0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000001, 0.000002, 0.000006, 0.000019, 0.000053, 0.000144, 0.000369, 0.000889, 0.002020, 0.004332, 0.008774, 0.016793, 0.030396, 0.052081, 0.084566, 0.130295, 0.190787, 0.265986, 0.353830, 0.450262, 0.549738, 0.646170, 0.734014, 0.809213, 0.869705, 0.915434, 0.947919, 0.969604, 0.983207, 0.991226, 0.995668, 0.997980, 0.999111, 0.999631, 0.999856, 0.999947, 0.999981, 0.999994, 0.999998, 0.999999, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000};
25*/
26
8082a80b 27#ifdef FASTJET_VERSION
28class AliJetShapeMass : public fastjet::FunctionOfPseudoJet<Double32_t>
29{
30 public:
31 virtual std::string description() const{return "jet mass";}
32 Double32_t result(const fastjet::PseudoJet &jet) const{ return jet.m();}
33};
874eb409 34
35class AliJetShapeGRNum : public fastjet::FunctionOfPseudoJet<Double32_t>
36{
37 public:
38 // default ctor
39 AliJetShapeGRNum(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
40 virtual std::string description() const{return "Numerator angular structure function";}
41 // static Int_t GetBin(Double_t x) {Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
42 Double32_t result(const fastjet::PseudoJet &jet) const {
43 if (!jet.has_constituents()) {
44 Printf("Angular structure can only be applied on jets for which the constituents are known.");
45 return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");
46 }
47 Double_t A = 0.;
48 std::vector<fastjet::PseudoJet> constits = jet.constituents();
49 for(UInt_t ic = 0; ic < constits.size(); ++ic) {
50 /* Int_t uid = constits[ic].user_index(); */
51 /* if (uid == -1) //skip ghost particle */
52 /* continue; */
53 for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
54 /* Int_t uid = constits[jc].user_index(); */
55 /* if (uid == -1) //skip ghost particle */
56 /* continue; */
57 Double_t dphi = constits[ic].phi()-constits[jc].phi();
58 if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
59 if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
60 Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
61 if(dr2>0.) {
62 Double_t dr = TMath::Sqrt(dr2);
63 Double_t x = fR-dr;
64 // Printf("Calculating numerator between %d-%d for r: %f dr: %f",ic,jc,fR,dr);
65 //noisy function
66 Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
67 /* Int_t bin = GetBin(x); */
68 /* Double_t noise = 0; */
69 /* if(bin<0) noise = gfcn[0]; */
70 /* else if(bin>=knbins) noise = gfcn[knbins-1]; */
71 /* else noise = gfcn[bin]; */
72 /* Printf("fR: %f dr: %f x=%f bin=%d noise: %f %f",fR,dr,x,bin,noise,TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep)); */
73 A += constits[ic].perp()*constits[jc].perp()*dr2*noise;
74 }
75 }
76 }
77 return A;
78 }
79
80 protected:
81 Double_t fR;
82 Double_t fDRStep;
83};
84
85class AliJetShapeGRDen : public fastjet::FunctionOfPseudoJet<Double32_t>
86{
87 public:
88 // default ctor
89 AliJetShapeGRDen(Double_t r = 0.2, Double_t wr = 0.04) : fR(r),fDRStep(wr){}
90 virtual std::string description() const{return "Denominator angular structure function";}
91 Double32_t result(const fastjet::PseudoJet &jet) const {
92 if (!jet.has_constituents())
93 return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");
94
95 Double_t A = 0.;
96 std::vector<fastjet::PseudoJet> constits = jet.constituents();
97 for(UInt_t ic = 0; ic < constits.size(); ++ic) {
98 /* Int_t uid = constits[ic].user_index(); */
99 /* if (uid == -1) //skip ghost particle */
100 /* continue; */
101 for(UInt_t jc = ic+1; jc < constits.size(); ++jc) {
102 /* Int_t uid = constits[jc].user_index(); */
103 /* if (uid == -1) //skip ghost particle */
104 /* continue; */
105 Double_t dphi = constits[ic].phi()-constits[jc].phi();
106 if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
107 if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
108 Double_t dr2 = (constits[ic].eta()-constits[jc].eta())*(constits[ic].eta()-constits[jc].eta()) + dphi*dphi;
109 if(dr2>0.) {
110 Double_t dr = TMath::Sqrt(dr2);
111 Double_t x = fR-dr;
112 //error function
113 Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
114 /* Int_t bin = GetBin(x); */
115 /* Double_t erf = 0; */
116 /* if(bin<0) erf = efcn[0]; */
117 /* else if(bin>=knbins) erf = efcn[knbins-1]; */
118 /* else erf = efcn[bin]; */
119 A += constits[ic].perp()*constits[jc].perp()*dr2*erf;
120 }
121 }
122 }
123 return A;
124 }
125
126 protected:
127 Double_t fR;
128 Double_t fDRStep;
129};
130
8082a80b 131#endif
132#endif