updates to AliEmcalJet, AliEmcalJetTask and AliFJWrapper to include subtraction of...
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliEmcalJet.cxx
CommitLineData
04f6e9de 1//
2// Emcal jet class.
3//
cd231d42 4// Author: C.Loizides
7efbea04 5
914d486c 6#include "AliEmcalJet.h"
7df864a3 7
59ebe296 8#include "AliLog.h"
9#include "Riostream.h"
10
914d486c 11ClassImp(AliEmcalJet)
7efbea04 12
13//__________________________________________________________________________________________________
96919f91 14AliEmcalJet::AliEmcalJet() :
15 AliVParticle(),
16 fPt(0),
17 fEta(0),
18 fPhi(0),
19 fM(0),
20 fNEF(0),
21 fArea(0),
db8a0184 22 fAreaEta(0),
23 fAreaPhi(0),
96919f91 24 fAreaEmc(-1),
25 fAxisInEmcal(0),
3c9775d9 26 fFlavourTagging(0),
96919f91 27 fMaxCPt(0),
28 fMaxNPt(0),
29 fMCPt(0),
30 fNn(0),
31 fNch(0),
8895e61b 32 fPtEmc(0),
33 fNEmc(0),
96919f91 34 fClusterIDs(),
35 fTrackIDs(),
8895e61b 36 fMatched(2),
cfc2ac24 37 fMatchingType(0),
254370e1 38 fTaggedJet(0x0),
39 fTagStatus(-1),
db8a0184 40 fPtSub(0),
85d48773 41 fPtVectSub(0),
8082a80b 42 fTriggers(0),
43 fJetShapeMassFirstDer(0),
44 fJetShapeMassSecondDer(0),
45 fJetShapeMassFirstSub(0),
46 fJetShapeMassSecondSub(0),
726c9488 47 fLabel(-1),
48 fGRNumerator(0),
49 fGRDenominator(0),
50 fGRNumeratorSub(0),
51 fGRDenominatorSub(0)
96919f91 52{
53 // Constructor.
54
55 fClosestJets[0] = 0;
56 fClosestJets[1] = 0;
57 fClosestJetsDist[0] = 999;
58 fClosestJetsDist[1] = 999;
96919f91 59}
60
61//__________________________________________________________________________________________________
101cefde 62AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
63 AliVParticle(),
64 fPt(TMath::Sqrt(px*px+py*py)),
65 fEta(TMath::ASinH(pz/fPt)),
66 fPhi(0),
67 fM(0),
68 fNEF(0),
69 fArea(0),
db8a0184 70 fAreaEta(0),
71 fAreaPhi(0),
101cefde 72 fAreaEmc(-1),
73 fAxisInEmcal(0),
3c9775d9 74 fFlavourTagging(0),
101cefde 75 fMaxCPt(0),
76 fMaxNPt(0),
96919f91 77 fMCPt(0),
78 fNn(0),
79 fNch(0),
8895e61b 80 fPtEmc(0),
81 fNEmc(0),
101cefde 82 fClusterIDs(),
96919f91 83 fTrackIDs(),
8895e61b 84 fMatched(2),
cfc2ac24 85 fMatchingType(0),
254370e1 86 fTaggedJet(0x0),
87 fTagStatus(-1),
db8a0184 88 fPtSub(0),
85d48773 89 fPtVectSub(0),
8082a80b 90 fTriggers(0),
91 fJetShapeMassFirstDer(0),
92 fJetShapeMassSecondDer(0),
93 fJetShapeMassFirstSub(0),
94 fJetShapeMassSecondSub(0),
726c9488 95 fLabel(-1),
96 fGRNumerator(0),
97 fGRDenominator(0),
98 fGRNumeratorSub(0),
99 fGRDenominatorSub(0)
7efbea04 100{
101 // Constructor.
102
103 if (fPt != 0) {
104 fPhi = TMath::ATan2(py, px);
105 if (fPhi<0.)
106 fPhi += 2. * TMath::Pi();
107 }
96919f91 108
109 fClosestJets[0] = 0;
110 fClosestJets[1] = 0;
111 fClosestJetsDist[0] = 999;
112 fClosestJetsDist[1] = 999;
7efbea04 113}
114
115//_________________________________________________________________________________________________
914d486c 116AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
9cffd2d0 117 AliVParticle(),
f8087a81 118 fPt(pt),
119 fEta(eta),
120 fPhi(phi),
121 fM(m),
122 fNEF(0),
123 fArea(0),
db8a0184 124 fAreaEta(0),
125 fAreaPhi(0),
101cefde 126 fAreaEmc(-1),
127 fAxisInEmcal(0),
3c9775d9 128 fFlavourTagging(0),
f8087a81 129 fMaxCPt(0),
96919f91 130 fMaxNPt(0),
131 fMCPt(0),
132 fNn(0),
133 fNch(0),
8895e61b 134 fPtEmc(0),
135 fNEmc(0),
f8087a81 136 fClusterIDs(),
96919f91 137 fTrackIDs(),
8895e61b 138 fMatched(2),
cfc2ac24 139 fMatchingType(0),
254370e1 140 fTaggedJet(0x0),
141 fTagStatus(-1),
db8a0184 142 fPtSub(0),
85d48773 143 fPtVectSub(0),
8082a80b 144 fTriggers(0),
145 fJetShapeMassFirstDer(0),
146 fJetShapeMassSecondDer(0),
147 fJetShapeMassFirstSub(0),
148 fJetShapeMassSecondSub(0),
726c9488 149 fLabel(-1),
150 fGRNumerator(0),
151 fGRDenominator(0),
152 fGRNumeratorSub(0),
153 fGRDenominatorSub(0)
9cffd2d0 154{
155 // Constructor.
156
4f149872 157 if (fPhi<0.)
158 fPhi += TMath::TwoPi();
96919f91 159
160 fClosestJets[0] = 0;
161 fClosestJets[1] = 0;
162 fClosestJetsDist[0] = 999;
163 fClosestJetsDist[1] = 999;
9cffd2d0 164}
165
166//_________________________________________________________________________________________________
914d486c 167AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
7efbea04 168 AliVParticle(jet),
f8087a81 169 fPt(jet.fPt),
170 fEta(jet.fEta),
171 fPhi(jet.fPhi),
172 fM(jet.fM),
173 fNEF(jet.fNEF),
174 fArea(jet.fArea),
db8a0184 175 fAreaEta(jet.fAreaEta),
176 fAreaPhi(jet.fAreaPhi),
101cefde 177 fAreaEmc(jet.fAreaEmc),
178 fAxisInEmcal(jet.fAxisInEmcal),
3c9775d9 179 fFlavourTagging(jet.fFlavourTagging),
f8087a81 180 fMaxCPt(jet.fMaxCPt),
181 fMaxNPt(jet.fMaxNPt),
96919f91 182 fMCPt(jet.fMCPt),
183 fNn(jet.fNn),
184 fNch(jet.fNch),
8895e61b 185 fPtEmc(jet.fPtEmc),
186 fNEmc(jet.fNEmc),
f8087a81 187 fClusterIDs(jet.fClusterIDs),
4f149872 188 fTrackIDs(jet.fTrackIDs),
8895e61b 189 fMatched(jet.fMatched),
cfc2ac24 190 fMatchingType(jet.fMatchingType),
254370e1 191 fTaggedJet(jet.fTaggedJet),
192 fTagStatus(jet.fTagStatus),
db8a0184 193 fPtSub(jet.fPtSub),
85d48773 194 fPtVectSub(jet.fPtVectSub),
8082a80b 195 fTriggers(jet.fTriggers),
196 fJetShapeMassFirstDer(jet.fJetShapeMassFirstDer),
197 fJetShapeMassSecondDer(jet.fJetShapeMassSecondDer),
198 fJetShapeMassFirstSub(jet.fJetShapeMassFirstSub),
199 fJetShapeMassSecondSub(jet.fJetShapeMassSecondSub),
726c9488 200 fLabel(jet.fLabel),
201 fGRNumerator(jet.fGRNumerator),
202 fGRDenominator(jet.fGRDenominator),
203 fGRNumeratorSub(jet.fGRNumeratorSub),
204 fGRDenominatorSub(jet.fGRDenominatorSub)
7df864a3 205{
101cefde 206 // Copy constructor.
4f149872 207
208 fClosestJets[0] = jet.fClosestJets[0];
209 fClosestJets[1] = jet.fClosestJets[1];
210 fClosestJetsDist[0] = jet.fClosestJetsDist[0];
211 fClosestJetsDist[1] = jet.fClosestJetsDist[1];
7efbea04 212}
213
214//_________________________________________________________________________________________________
914d486c 215AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
7efbea04 216{
217 // Assignment operator.
218
219 if (this!=&jet) {
220 AliVParticle::operator=(jet);
4f149872 221 fPt = jet.fPt;
222 fEta = jet.fEta;
223 fPhi = jet.fPhi;
224 fM = jet.fM;
225 fNEF = jet.fNEF;
226 fArea = jet.fArea;
db8a0184 227 fAreaEta = jet.fAreaEta;
228 fAreaPhi = jet.fAreaPhi;
4f149872 229 fAreaEmc = jet.fAreaEmc;
230 fAxisInEmcal = jet.fAxisInEmcal;
3c9775d9 231 fFlavourTagging = jet.fFlavourTagging;
4f149872 232 fMaxCPt = jet.fMaxCPt;
233 fMaxNPt = jet.fMaxNPt;
234 fMCPt = jet.fMCPt;
235 fNn = jet.fNn;
236 fNch = jet.fNch;
8895e61b 237 fPtEmc = jet.fPtEmc;
238 fNEmc = jet.fNEmc;
4f149872 239 fClusterIDs = jet.fClusterIDs;
240 fTrackIDs = jet.fTrackIDs;
241 fClosestJets[0] = jet.fClosestJets[0];
242 fClosestJets[1] = jet.fClosestJets[1];
243 fClosestJetsDist[0] = jet.fClosestJetsDist[0];
244 fClosestJetsDist[1] = jet.fClosestJetsDist[1];
245 fMatched = jet.fMatched;
254370e1 246 fTaggedJet = jet.fTaggedJet;
247 fTagStatus = jet.fTagStatus;
8895e61b 248 fPtSub = jet.fPtSub;
f660c2d6 249 fPtVectSub = jet.fPtVectSub;
85d48773 250 fTriggers = jet.fTriggers;
8082a80b 251 fJetShapeMassFirstDer = jet.fJetShapeMassFirstDer;
252 fJetShapeMassSecondDer = jet.fJetShapeMassSecondDer;
253 fJetShapeMassFirstSub = jet.fJetShapeMassFirstSub;
254 fJetShapeMassSecondSub = jet.fJetShapeMassSecondSub;
255 fLabel = jet.fLabel;
726c9488 256 fGRNumerator = jet.fGRNumerator;
257 fGRDenominator = jet.fGRDenominator;
258 fGRNumeratorSub = jet.fGRNumeratorSub;
259 fGRDenominatorSub = jet.fGRDenominatorSub;
7efbea04 260 }
261
262 return *this;
7df864a3 263}
7efbea04 264
2103dc6a 265//_________________________________________________________________________________________________
db8a0184 266Int_t AliEmcalJet::Compare(const TObject* obj) const
267{
268 //Return -1 if this is smaller than obj, 0 if objects are equal and 1 if this is larger than obj.
269
270 const AliEmcalJet *jet = static_cast<const AliEmcalJet *>(obj);
271 if (!obj)
272 return 0;
273 if (Pt()>jet->Pt())
274 return -1;
275 return 1;
276}
277
7efbea04 278//__________________________________________________________________________________________________
914d486c 279void AliEmcalJet::GetMom(TLorentzVector &vec) const
7efbea04 280{
281 // Return momentum as four vector.
282
283 Double_t p = fPt *TMath::CosH(fEta);
284 vec.SetPtEtaPhiE(fPt,fEta,fPhi,TMath::Sqrt(p*p+fM*fM));
285}
286
287//__________________________________________________________________________________________________
914d486c 288void AliEmcalJet::Print(Option_t* /*option*/) const
7efbea04 289{
290 // Print jet information.
291
59ebe296 292 Printf("Jet pt=%.2f, eta=%.2f, phi=%.2f, area=%.2f, NEF=%.2f", fPt, fEta, fPhi, fArea, fNEF);
7efbea04 293}
629e03fc 294
295//__________________________________________________________________________________________________
db8a0184 296Double_t AliEmcalJet::PtSubVect(Double_t rho) const
297{
298 // Return vectorial subtracted transverse momentum.
299
300 Double_t dx = Px() - rho * fArea * TMath::Cos(fAreaPhi);
301 Double_t dy = Py() - rho * fArea * TMath::Sin(fAreaPhi);
302 //Double_t dz = Pz() - rho * fArea * TMath::SinH(fAreaEta);
303 return TMath::Sqrt(dx*dx+dy*dy);
304}
305
306//__________________________________________________________________________________________________
629e03fc 307void AliEmcalJet::SortConstituents()
308{
f8087a81 309 // Sort constituent by index (increasing).
310
311 std::sort(fClusterIDs.GetArray(), fClusterIDs.GetArray() + fClusterIDs.GetSize());
312 std::sort(fTrackIDs.GetArray(), fTrackIDs.GetArray() + fTrackIDs.GetSize());
629e03fc 313}
f660c2d6 314
315//__________________________________________________________________________________________________
d1f0108e 316Double_t AliEmcalJet::DeltaR (const AliVParticle* part) const
317{ // Helper function to calculate the distance between two jets or a jet and a particle
318 Double_t dPhi = this->Phi() - part->Phi();
319 Double_t dEta = this->Eta() - part->Eta();
320 dPhi = TVector2::Phi_mpi_pi ( dPhi );
321
322 return TMath::Sqrt ( dPhi * dPhi + dEta * dEta );
323}
324
325
326//__________________________________________________________________________________________________
327std::vector<int> AliEmcalJet::SortConstituentsPt( TClonesArray *tracks ) const
328{ //___________________________________________
329 // Sorting by p_T (decreasing) jet constituents
330 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
331 typedef std::pair<Double_t, Int_t> ptidx_pair;
332
333 // Create vector for Pt sorting
334 std::vector<ptidx_pair> pair_list ;
335
336 for ( Int_t i_entry = 0; i_entry < GetNumberOfTracks(); i_entry++ )
337 {
338 AliVParticle *track = TrackAt(i_entry, tracks);
339 if (!track)
340 {
341 AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)", i_entry, tracks->GetName(), TrackAt(i_entry), tracks->GetEntriesFast()));
342 continue;
343 }
344
345 pair_list.push_back( std::make_pair ( track->Pt(), i_entry ) );
346 }
347
348 std::stable_sort( pair_list.begin() , pair_list.end() , sort_descend() );
349
350 // return an vector of indexes of constituents (sorted descending by pt)
351 std::vector <int> index_sorted_list;
352
353 for ( std::vector< std::pair<Double_t,Int_t> >::iterator it = pair_list.begin(); it != pair_list.end(); ++it)
354 { index_sorted_list.push_back( (*it).second ); } // populating the return object with indexes of sorted tracks
355
356 return index_sorted_list;
357}
358
359//__________________________________________________________________________________________________
f660c2d6 360AliVParticle* AliEmcalJet::GetLeadingTrack(TClonesArray *tracks) const
361{
362 AliVParticle* maxTrack = 0;
363 for (Int_t i = 0; i < GetNumberOfTracks(); i++) {
364 AliVParticle *track = TrackAt(i, tracks);
59ebe296 365 if (!track) {
3c9775d9 366 AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)",
367 i,tracks->GetName(),TrackAt(i),tracks->GetEntriesFast()));
59ebe296 368 continue;
369 }
f660c2d6 370 if (!maxTrack || track->Pt() > maxTrack->Pt())
371 maxTrack = track;
372 }
373
374 return maxTrack;
375}
376
377//__________________________________________________________________________________________________
378AliVCluster* AliEmcalJet::GetLeadingCluster(TClonesArray *clusters) const
379{
380 AliVCluster* maxCluster = 0;
381 for (Int_t i = 0; i < GetNumberOfClusters(); i++) {
382 AliVCluster *cluster = ClusterAt(i, clusters);
59ebe296 383 if (!cluster) {
3c9775d9 384 AliError(Form("Unable to find jet cluster %d in collection %s (pos in collection %d, max %d)",
385 i,clusters->GetName(),ClusterAt(i),clusters->GetEntriesFast()));
59ebe296 386 continue;
387 }
f660c2d6 388 if (!maxCluster || cluster->E() > maxCluster->E())
389 maxCluster = cluster;
390 }
391
392 return maxCluster;
393}
2103dc6a 394
395//__________________________________________________________________________________________________
396void AliEmcalJet::ResetMatching()
397{
398 fClosestJets[0] = 0;
399 fClosestJets[1] = 0;
400 fClosestJetsDist[0] = 999;
401 fClosestJetsDist[1] = 999;
402 fMatched = 2;
403}
726c9488 404
405//__________________________________________________________________________________________________
406void AliEmcalJet::PrintGR() {
407 for(Int_t i = 0; i<fGRNumerator.GetSize(); i++) {
408 Printf("num[%d] = %f",i,fGRNumerator.At(i));
409 }
410}