]>
Commit | Line | Data |
---|---|---|
4d3b366f | 1 | // $Id$ |
2 | // | |
3 | // Calculation of rho mass from a collection of jets. | |
4 | // If scale function is given the scaled rho will be exported | |
5 | // with the name as "fOutRhoMassName".Apppend("_Scaled"). | |
6 | // | |
7 | // Authors: M. Verweij | |
8 | ||
9 | #include "AliAnalysisTaskRhoMass.h" | |
10 | ||
11 | #include <TClonesArray.h> | |
12 | #include <TMath.h> | |
13 | ||
14 | #include "AliAnalysisManager.h" | |
15 | #include "AliEmcalJet.h" | |
16 | #include "AliLog.h" | |
17 | #include "AliRhoParameter.h" | |
18 | ||
19 | ClassImp(AliAnalysisTaskRhoMass) | |
20 | ||
21 | //________________________________________________________________________ | |
22 | AliAnalysisTaskRhoMass::AliAnalysisTaskRhoMass() : | |
23 | AliAnalysisTaskRhoMassBase("AliAnalysisTaskRhoMass"), | |
24 | fNExclLeadJets(0), | |
25 | fJetRhoMassType(kMd), | |
26 | fHistMdAreavsCent(0) | |
27 | { | |
28 | // Constructor. | |
29 | } | |
30 | ||
31 | //________________________________________________________________________ | |
32 | AliAnalysisTaskRhoMass::AliAnalysisTaskRhoMass(const char *name, Bool_t histo) : | |
33 | AliAnalysisTaskRhoMassBase(name, histo), | |
34 | fNExclLeadJets(0), | |
35 | fJetRhoMassType(kMd), | |
36 | fHistMdAreavsCent(0) | |
37 | { | |
38 | // Constructor. | |
39 | } | |
40 | ||
41 | //________________________________________________________________________ | |
42 | void AliAnalysisTaskRhoMass::UserCreateOutputObjects() | |
43 | { | |
44 | // User create output objects, called at the beginning of the analysis. | |
45 | ||
46 | if (!fCreateHisto) | |
47 | return; | |
48 | ||
49 | AliAnalysisTaskRhoMassBase::UserCreateOutputObjects(); | |
50 | ||
51 | fHistMdAreavsCent = new TH2F("fHistMdAreavsCent", "fHistMdAreavsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt/2.); | |
52 | fHistMdAreavsCent->GetXaxis()->SetTitle("Centrality (%)"); | |
53 | fHistMdAreavsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})"); | |
54 | fOutput->Add(fHistMdAreavsCent); | |
55 | } | |
56 | ||
57 | ||
58 | //________________________________________________________________________ | |
59 | Bool_t AliAnalysisTaskRhoMass::Run() | |
60 | { | |
61 | // Run the analysis. | |
62 | ||
63 | fOutRhoMass->SetVal(0); | |
64 | if (fOutRhoMassScaled) | |
65 | fOutRhoMassScaled->SetVal(0); | |
66 | ||
67 | if (!fJets) | |
68 | return kFALSE; | |
69 | ||
70 | const Int_t Njets = fJets->GetEntries(); | |
71 | ||
72 | Int_t maxJetIds[] = {-1, -1}; | |
73 | Float_t maxJetPts[] = { 0, 0}; | |
74 | ||
75 | if (fNExclLeadJets > 0) { | |
76 | for (Int_t ij = 0; ij < Njets; ++ij) { | |
77 | AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij)); | |
78 | if (!jet) { | |
79 | AliError(Form("%s: Could not receive jet %d", GetName(), ij)); | |
80 | continue; | |
81 | } | |
82 | ||
83 | if (!AcceptJet(jet)) | |
84 | continue; | |
85 | ||
86 | if (jet->Pt() > maxJetPts[0]) { | |
87 | maxJetPts[1] = maxJetPts[0]; | |
88 | maxJetIds[1] = maxJetIds[0]; | |
89 | maxJetPts[0] = jet->Pt(); | |
90 | maxJetIds[0] = ij; | |
91 | } else if (jet->Pt() > maxJetPts[1]) { | |
92 | maxJetPts[1] = jet->Pt(); | |
93 | maxJetIds[1] = ij; | |
94 | } | |
95 | } | |
96 | if (fNExclLeadJets < 2) { | |
97 | maxJetIds[1] = -1; | |
98 | maxJetPts[1] = 0; | |
99 | } | |
100 | } | |
101 | ||
102 | static Double_t rhomvec[999]; | |
103 | static Double_t Evec[999]; | |
104 | static Double_t Mvec[999]; | |
105 | Int_t NjetAcc = 0; | |
106 | ||
107 | // push all jets within selected acceptance into stack | |
108 | for (Int_t iJets = 0; iJets < Njets; ++iJets) { | |
109 | ||
110 | // exlcuding lead jets | |
111 | if (iJets == maxJetIds[0] || iJets == maxJetIds[1]) | |
112 | continue; | |
113 | ||
114 | AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets)); | |
115 | if (!jet) { | |
116 | AliError(Form("%s: Could not receive jet %d", GetName(), iJets)); | |
117 | continue; | |
118 | } | |
119 | ||
120 | if (!AcceptJet(jet)) | |
121 | continue; | |
122 | ||
123 | // Double_t sumM = GetSumMConstituents(jet); | |
124 | // Double_t sumPt = GetSumPtConstituents(jet); | |
125 | if(jet->Area()>0.) {// && (jet->M()*jet->M() + jet->Pt()*jet->Pt())>0.) { | |
126 | //rhomvec[NjetAcc] = (TMath::Sqrt(sumM*sumM + sumPt*sumPt) - sumPt ) / jet->Area(); | |
127 | // rhomvec[NjetAcc] = (TMath::Sqrt(jet->M()*jet->M() + jet->Pt()*jet->Pt()) - jet->Pt() ) / jet->Area(); | |
128 | rhomvec[NjetAcc] = GetMd(jet) / jet->Area(); | |
129 | fHistMdAreavsCent->Fill(fCent,rhomvec[NjetAcc]); | |
130 | Evec[NjetAcc] = jet->E(); | |
131 | Mvec[NjetAcc] = jet->M(); | |
132 | ++NjetAcc; | |
133 | } | |
134 | } | |
135 | ||
136 | if (NjetAcc > 0) { | |
137 | //find median value | |
138 | Double_t rhom = TMath::Median(NjetAcc, rhomvec); | |
139 | fOutRhoMass->SetVal(rhom); | |
140 | ||
141 | Int_t Ntracks = fTracks->GetEntries(); | |
142 | Double_t meanM = TMath::Mean(NjetAcc, Mvec); | |
143 | Double_t meanE = TMath::Mean(NjetAcc, Evec); | |
144 | Double_t gamma = 0.; | |
145 | if(meanM>0.) gamma = meanE/meanM; | |
146 | fHistGammaVsNtrack->Fill(Ntracks,gamma); | |
147 | ||
148 | if (fOutRhoMassScaled) { | |
149 | Double_t rhomScaled = rhom * GetScaleFactor(fCent); | |
150 | fOutRhoMassScaled->SetVal(rhomScaled); | |
151 | } | |
152 | } | |
153 | ||
154 | return kTRUE; | |
155 | } | |
156 | ||
157 | //________________________________________________________________________ | |
158 | Double_t AliAnalysisTaskRhoMass::GetSumMConstituents(AliEmcalJet *jet) { | |
159 | ||
160 | Double_t sum = 0.; | |
161 | ||
162 | AliVParticle *vp; | |
163 | for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) { | |
164 | vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks)); | |
165 | if(!vp) continue; | |
166 | sum+=vp->M(); | |
167 | } | |
168 | return sum; | |
169 | } | |
170 | ||
171 | //________________________________________________________________________ | |
172 | Double_t AliAnalysisTaskRhoMass::GetSumPtConstituents(AliEmcalJet *jet) { | |
173 | ||
174 | Double_t sum = 0.; | |
175 | ||
176 | AliVParticle *vp; | |
177 | for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) { | |
178 | vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks)); | |
179 | if(!vp) continue; | |
180 | sum+=vp->Pt(); | |
181 | } | |
182 | return sum; | |
183 | } | |
184 | ||
185 | //________________________________________________________________________ | |
186 | Double_t AliAnalysisTaskRhoMass::GetMd(AliEmcalJet *jet) { | |
187 | //get md as defined in http://arxiv.org/pdf/1211.2811.pdf | |
188 | Double_t sum = 0.; | |
189 | Double_t px = 0.; | |
190 | Double_t py = 0.; | |
191 | Double_t pz = 0.; | |
192 | Double_t E = 0.; | |
193 | AliVParticle *vp; | |
194 | for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) { | |
195 | vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks)); | |
196 | if(!vp) continue; | |
197 | if(fJetRhoMassType==kMd) sum += TMath::Sqrt(vp->M()*vp->M() + vp->Pt()*vp->Pt()) - vp->Pt(); //sqrt(E^2-P^2+pt^2)=sqrt(E^2-pz^2) | |
198 | else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()*vp->M() + vp->P()*vp->P()) - vp->P(); | |
199 | else if(fJetRhoMassType==kMd4) { | |
200 | px+=vp->Px(); | |
201 | py+=vp->Py(); | |
202 | pz+=vp->Pz(); | |
203 | E+=vp->E(); | |
204 | } | |
205 | } | |
206 | if(fJetRhoMassType==kMd4) { | |
207 | Double_t pt = TMath::Sqrt(px*px + py*py); | |
208 | Double_t m2 = E*E - pt*pt - pz*pz; | |
209 | sum = TMath::Sqrt(m2 + pt*pt) - pt; | |
210 | } | |
211 | return sum; | |
212 | } |