]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
pp PbPb transparent, sparse bin changes, set up cuts on tracks for aods and esds...
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliEmcalJetTask.cxx
CommitLineData
542ac102 1//
2// Emcal jet finder task.
3//
8082a80b 4// Authors: C.Loizides, S.Aiola, M. Verweij
48d874ff 5
c64cb1f6 6#include <vector>
04905231 7#include "AliEmcalJetTask.h"
8
6ea168d0 9#include <TChain.h>
48d874ff 10#include <TClonesArray.h>
6ea168d0 11#include <TList.h>
12#include <TLorentzVector.h>
ca5c29fa 13#include <TMath.h>
6421eeb0 14#include <TRandom3.h>
b65fac7a 15
48d874ff 16#include "AliAnalysisManager.h"
6ea168d0 17#include "AliCentrality.h"
f5f3c8e9 18#include "AliEMCALGeometry.h"
04905231 19#include "AliESDEvent.h"
f5f3c8e9 20#include "AliEmcalJet.h"
04905231 21#include "AliEmcalParticle.h"
63206072 22#include "AliFJWrapper.h"
04905231 23#include "AliMCEvent.h"
35789a2d 24#include "AliVCluster.h"
b65fac7a 25#include "AliVEvent.h"
04905231 26#include "AliVParticle.h"
8082a80b 27#include "AliRhoParameter.h"
c2aad3ae 28using std::cout;
29using std::endl;
30using std::cerr;
48d874ff 31
914d486c 32ClassImp(AliEmcalJetTask)
48d874ff 33
6ea168d0 34//________________________________________________________________________
914d486c 35AliEmcalJetTask::AliEmcalJetTask() :
36 AliAnalysisTaskSE("AliEmcalJetTask"),
6ea168d0 37 fTracksName("Tracks"),
38 fCaloName("CaloClusters"),
39 fJetsName("Jets"),
8082a80b 40 fJetsSubName(""),
6a20534a 41 fJetType(kNone),
d1f2e0d9 42 fConstSel(0),
43 fMCConstSel(0),
6a20534a 44 fMarkConst(kFALSE),
6ea168d0 45 fRadius(0.4),
914d486c 46 fMinJetTrackPt(0.15),
47 fMinJetClusPt(0.15),
04905231 48 fPhiMin(0),
49 fPhiMax(TMath::TwoPi()),
83086e6c 50 fEtaMin(-0.9),
51 fEtaMax(+0.9),
cba7b46e 52 fMinJetArea(0.005),
f5f3c8e9 53 fMinJetPt(1.0),
c9ad772c 54 fJetPhiMin(-10),
83086e6c 55 fJetPhiMax(+10),
c9ad772c 56 fJetEtaMin(-1),
83086e6c 57 fJetEtaMax(+1),
cba7b46e 58 fGhostArea(0.005),
a7477843 59 fMinMCLabel(0),
751cb2dc 60 fRecombScheme(fastjet::pt_scheme),
6421eeb0 61 fTrackEfficiency(1.),
bf9072ca 62 fMCFlag(0),
63 fGeneratorIndex(-1),
04905231 64 fIsInit(0),
83c097da 65 fLocked(0),
d14d6e2c 66 fIsPSelSet(0),
04905231 67 fIsMcPart(0),
68 fIsEmcPart(0),
5d0a5007 69 fLegacyMode(kFALSE),
ffdca6aa 70 fCodeDebug(kFALSE),
8082a80b 71 fDoGenericSubtraction(kFALSE),
72 fDoConstituentSubtraction(kFALSE),
73 fUseExternalBkg(kFALSE),
74 fRhoName(""),
75 fRhomName(""),
76 fRho(0),
77 fRhom(0),
b65fac7a 78 fJets(0),
8082a80b 79 fJetsSub(0),
04905231 80 fEvent(0),
81 fTracks(0),
8082a80b 82 fClus(0),
83 fRhoParam(0),
84 fRhomParam(0)
6ea168d0 85{
86 // Default constructor.
6ea168d0 87}
88
48d874ff 89//________________________________________________________________________
914d486c 90AliEmcalJetTask::AliEmcalJetTask(const char *name) :
63206072 91 AliAnalysisTaskSE(name),
7efbea04 92 fTracksName("Tracks"),
48d874ff 93 fCaloName("CaloClusters"),
7efbea04 94 fJetsName("Jets"),
8082a80b 95 fJetsSubName(""),
6a20534a 96 fJetType(kAKT|kFullJet|kRX1Jet),
d1f2e0d9 97 fConstSel(0),
98 fMCConstSel(0),
6a20534a 99 fMarkConst(kFALSE),
48d874ff 100 fRadius(0.4),
6ea168d0 101 fMinJetTrackPt(0.15),
914d486c 102 fMinJetClusPt(0.15),
04905231 103 fPhiMin(0),
104 fPhiMax(TMath::TwoPi()),
83086e6c 105 fEtaMin(-0.9),
106 fEtaMax(+0.9),
cba7b46e 107 fMinJetArea(0.001),
f5f3c8e9 108 fMinJetPt(1.0),
c9ad772c 109 fJetPhiMin(-10),
83086e6c 110 fJetPhiMax(+10),
c9ad772c 111 fJetEtaMin(-1),
83086e6c 112 fJetEtaMax(+1),
cba7b46e 113 fGhostArea(0.005),
a7477843 114 fMinMCLabel(0),
751cb2dc 115 fRecombScheme(fastjet::pt_scheme),
6421eeb0 116 fTrackEfficiency(1.),
bf9072ca 117 fMCFlag(0),
118 fGeneratorIndex(-1),
04905231 119 fIsInit(0),
83c097da 120 fLocked(0),
d14d6e2c 121 fIsPSelSet(0),
04905231 122 fIsMcPart(0),
123 fIsEmcPart(0),
31bab07e 124 fLegacyMode(kFALSE),
ffdca6aa 125 fCodeDebug(kFALSE),
8082a80b 126 fDoGenericSubtraction(kFALSE),
127 fDoConstituentSubtraction(kFALSE),
128 fUseExternalBkg(kFALSE),
129 fRhoName(""),
130 fRhomName(""),
131 fRho(0),
132 fRhom(0),
b65fac7a 133 fJets(0),
8082a80b 134 fJetsSub(0),
04905231 135 fEvent(0),
136 fTracks(0),
8082a80b 137 fClus(0),
138 fRhoParam(0),
139 fRhomParam(0)
48d874ff 140{
141 // Standard constructor.
7efbea04 142
7efbea04 143 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
48d874ff 144}
145
146//________________________________________________________________________
914d486c 147AliEmcalJetTask::~AliEmcalJetTask()
48d874ff 148{
149 // Destructor
150}
151
152//________________________________________________________________________
914d486c 153void AliEmcalJetTask::UserCreateOutputObjects()
48d874ff 154{
155 // Create user objects.
156
914d486c 157 fJets = new TClonesArray("AliEmcalJet");
48d874ff 158 fJets->SetName(fJetsName);
8082a80b 159
160 if(!fJetsSubName.IsNull()) {
161 fJetsSub = new TClonesArray("AliEmcalJet");
162 fJetsSub->SetName(fJetsSubName);
163 }
48d874ff 164}
165
166//________________________________________________________________________
914d486c 167void AliEmcalJetTask::UserExec(Option_t *)
48d874ff 168{
169 // Main loop, called for each event.
04905231 170 if (!fIsInit) {
171 if (!DoInit())
172 return;
173 fIsInit = kTRUE;
1772fe65 174 }
b65fac7a 175
918dd8c8 176 // clear the jet array (normally a null operation)
177 fJets->Delete();
8082a80b 178 if(fJetsSub) fJetsSub->Delete();
918dd8c8 179
04905231 180 FindJets();
48d874ff 181}
182
183//________________________________________________________________________
914d486c 184void AliEmcalJetTask::Terminate(Option_t *)
48d874ff 185{
186 // Called once at the end of the analysis.
48d874ff 187}
188
189//________________________________________________________________________
04905231 190void AliEmcalJetTask::FindJets()
48d874ff 191{
192 // Find jets.
3c5aff5d 193 if (!fTracks && !fClus){
194 cout << "WARNING NO TRACKS OR CLUSTERS:" <<endl;
d03084cd 195 return;
3c5aff5d 196 }
197
8082a80b 198 if (fRhoParam)
199 fRho = fRhoParam->GetVal();
200 if (fRhomParam)
201 fRhom = fRhomParam->GetVal();
202
48d874ff 203 TString name("kt");
204 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
6a20534a 205 if ((fJetType & kAKT) != 0) {
48d874ff 206 name = "antikt";
207 jalgo = fastjet::antikt_algorithm;
ca5c29fa 208 AliDebug(1,"Using AKT algorithm");
209 }
210 else {
211 AliDebug(1,"Using KT algorithm");
48d874ff 212 }
213
6a20534a 214 if ((fJetType & kR020Jet) != 0)
215 fRadius = 0.2;
216 else if ((fJetType & kR030Jet) != 0)
217 fRadius = 0.3;
218 else if ((fJetType & kR040Jet) != 0)
219 fRadius = 0.4;
220
04905231 221 // setup fj wrapper
48d874ff 222 AliFJWrapper fjw(name, name);
f5f3c8e9 223 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
1f9c287f 224 fjw.SetGhostArea(fGhostArea);
04905231 225 fjw.SetR(fRadius);
7071e730 226 fjw.SetAlgorithm(jalgo);
31bab07e 227 fjw.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
83086e6c 228 fjw.SetMaxRap(fEtaMax);
48d874ff 229 fjw.Clear();
230
04905231 231 // get primary vertex
232 Double_t vertex[3] = {0, 0, 0};
816bb2c8 233 if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex);
04905231 234
ca5c29fa 235 AliDebug(2,Form("Jet type = %d", fJetType));
6a20534a 236
237 if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) {
04905231 238 const Int_t Ntracks = fTracks->GetEntries();
48d874ff 239 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
04905231 240 AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
48d874ff 241 if (!t)
242 continue;
6a20534a 243 if (fIsMcPart) {
5be3857d 244 if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
245 AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
6a20534a 246 continue;
5be3857d 247 }
248 if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
249 AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
6a20534a 250 continue;
5be3857d 251 }
6a20534a 252 }
a7477843 253 if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) {
d1f2e0d9 254 if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
255 AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
6a20534a 256 continue;
257 }
d1f2e0d9 258 else {
259 AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
ca5c29fa 260 }
6a20534a 261 }
262 else {
d1f2e0d9 263 if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
264 AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
6a20534a 265 continue;
266 }
d1f2e0d9 267 else {
268 AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
ca5c29fa 269 }
bf9072ca 270 }
271 if ((t->GetFlag() & fMCFlag) != fMCFlag) {
272 AliDebug(2,Form("Skipping track %d because it does not match the MC flags", iTracks));
273 continue;
274 }
275 if (fGeneratorIndex >= 0 && t->GetGeneratorIndex() != fGeneratorIndex) {
276 AliDebug(2,Form("Skipping track %d because it does not match the MC generator index", iTracks));
277 continue;
278 }
1772fe65 279 if (t->Pt() < fMinJetTrackPt)
914d486c 280 continue;
04905231 281 Double_t eta = t->Eta();
282 Double_t phi = t->Phi();
ddf74a88 283 if ((eta<fEtaMin) || (eta>fEtaMax) ||
284 (phi<fPhiMin) || (phi>fPhiMax))
04905231 285 continue;
da67b88a 286
6421eeb0 287 // artificial inefficiency
288 if (fTrackEfficiency < 1.) {
289 Double_t rnd = gRandom->Rndm();
290 if (fTrackEfficiency < rnd) {
291 AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", iTracks));
292 continue;
293 }
294 }
295
a3c8c8c8 296 // offset of 100 for consistency with cluster ids
ca5c29fa 297 AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
751cb2dc 298 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
48d874ff 299 }
300 }
6ea168d0 301
6a20534a 302 if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) {
da67b88a 303 const Int_t Nclus = fClus->GetEntries();
304 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
305 AliVCluster *c = 0;
306 Double_t cEta=0,cPhi=0,cPt=0;
307 Double_t cPx=0,cPy=0,cPz=0;
308 if (fIsEmcPart) {
309 AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
310 if (!ep)
311 continue;
6a20534a 312
da67b88a 313 c = ep->GetCluster();
314 if (!c)
6a20534a 315 continue;
316
a7477843 317 if (c->GetLabel() > fMinMCLabel) {
d1f2e0d9 318 if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
319 AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
6a20534a 320 continue;
ca5c29fa 321 }
d1f2e0d9 322 else {
323 AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
ca5c29fa 324 }
6a20534a 325 }
326 else {
d1f2e0d9 327 if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
328 AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
6a20534a 329 continue;
ca5c29fa 330 }
d1f2e0d9 331 else {
332 AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
ca5c29fa 333 }
6a20534a 334 }
335
da67b88a 336 cEta = ep->Eta();
337 cPhi = ep->Phi();
338 cPt = ep->Pt();
339 cPx = ep->Px();
340 cPy = ep->Py();
341 cPz = ep->Pz();
342 } else {
343 c = static_cast<AliVCluster*>(fClus->At(iClus));
344 if (!c)
345 continue;
6a20534a 346
a7477843 347 if (c->GetLabel() > fMinMCLabel) {
d1f2e0d9 348 if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
349 AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
6a20534a 350 continue;
ca5c29fa 351 }
d1f2e0d9 352 else {
353 AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
ca5c29fa 354 }
6a20534a 355 }
356 else {
d1f2e0d9 357 if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
358 AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
6a20534a 359 continue;
ca5c29fa 360 }
d1f2e0d9 361 else {
362 AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
ca5c29fa 363 }
6a20534a 364 }
365
da67b88a 366 TLorentzVector nP;
367 c->GetMomentum(nP, vertex);
368 cEta = nP.Eta();
369 cPhi = nP.Phi();
370 cPt = nP.Pt();
371 cPx = nP.Px();
372 cPy = nP.Py();
373 cPz = nP.Pz();
04905231 374 }
da67b88a 375 if (!c->IsEMCAL())
376 continue;
377 if (cPt < fMinJetClusPt)
378 continue;
ddf74a88 379 if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
380 (cPhi<fPhiMin) || (cPhi>fPhiMax))
da67b88a 381 continue;
382 // offset of 100 to skip ghost particles uid = -1
ca5c29fa 383 AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
f660c2d6 384 fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
48d874ff 385 }
386 }
5d0a5007 387
388 // setting legacy mode
389 if (fLegacyMode) {
390 fjw.SetLegacyMode(kTRUE);
391 }
392
7efbea04 393 // run jet finder
48d874ff 394 fjw.Run();
f5f3c8e9 395
8082a80b 396 //run generic subtractor
397 if(fDoGenericSubtraction) {
398 fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
399 fjw.DoGenericSubtractionJetMass();
400 }
401
402 //run constituent subtractor
403 if(fDoConstituentSubtraction) {
404 fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
405 fjw.DoConstituentSubtraction();
406 }
407
f5f3c8e9 408 // get geometry
409 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
410 if (!geom) {
02c7e943 411 AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9 412 return;
413 }
04905231 414
415 // loop over fastjet jets
48d874ff 416 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
ca5c29fa 417 // sort jets according to jet pt
418 static Int_t indexes[9999] = {-1};
419 GetSortedArray(indexes, jets_incl);
420
421 AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
422 for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
423 Int_t ij = indexes[ijet];
424 AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
425
f5f3c8e9 426 if (jets_incl[ij].perp()<fMinJetPt)
427 continue;
428 if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff 429 continue;
c9ad772c 430 if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
431 (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
432 continue;
04905231 433
914d486c 434 AliEmcalJet *jet = new ((*fJets)[jetCount])
435 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
8082a80b 436 jet->SetLabel(ij);
437
438 //do generic subtraction if requested
439 if(fDoGenericSubtraction) {
289a9376 440#ifdef FASTJET_CONTRIB
8082a80b 441 std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
442 UInt_t n = (UInt_t)jetMassInfo.size();
443 if(n>ij && n>0) {
444 jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
445 jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
446 jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
447 jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
448 }
449#endif
450 }
04905231 451
452 // loop over constituents
c64cb1f6 453 std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
04905231 454 jet->SetNumberOfTracks(constituents.size());
455 jet->SetNumberOfClusters(constituents.size());
456
b65fac7a 457 Int_t nt = 0;
458 Int_t nc = 0;
f5f3c8e9 459 Double_t neutralE = 0;
111318e8 460 Double_t maxCh = 0;
461 Double_t maxNe = 0;
f5f3c8e9 462 Int_t gall = 0;
463 Int_t gemc = 0;
04905231 464 Int_t cemc = 0;
111318e8 465 Int_t ncharged = 0;
466 Int_t nneutral = 0;
02c7e943 467 Double_t mcpt = 0;
04905231 468 Double_t emcpt = 0;
b65fac7a 469
b65fac7a 470 for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
48d874ff 471 Int_t uid = constituents[ic].user_index();
5be3857d 472 AliDebug(3,Form("Processing constituent %d", uid));
04905231 473 if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
f5f3c8e9 474 ++gall;
04905231 475 Double_t gphi = constituents[ic].phi();
476 if (gphi<0)
477 gphi += TMath::TwoPi();
478 gphi *= TMath::RadToDeg();
f5f3c8e9 479 Double_t geta = constituents[ic].eta();
b65fac7a 480 if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
481 (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
04905231 482 ++gemc;
483 } else if ((uid > 0) && fTracks) { // track constituent
111318e8 484 Int_t tid = uid - 100;
04905231 485 AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
ca5c29fa 486 if (!t) {
487 AliError(Form("Could not find track %d",tid));
04905231 488 continue;
ca5c29fa 489 }
490 if (jetCount < fMarkConst) {
491 AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
6a20534a 492 t->SetBit(fJetType);
ca5c29fa 493 }
04905231 494 Double_t cEta = t->Eta();
495 Double_t cPhi = t->Phi();
496 Double_t cPt = t->Pt();
497 Double_t cP = t->P();
498 if (t->Charge() == 0) {
499 neutralE += cP;
500 ++nneutral;
501 if (cPt > maxNe)
502 maxNe = cPt;
503 } else {
504 ++ncharged;
505 if (cPt > maxCh)
506 maxCh = cPt;
507 }
a7477843 508 if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
04905231 509 mcpt += cPt;
510
511 if (cPhi<0)
512 cPhi += TMath::TwoPi();
513 cPhi *= TMath::RadToDeg();
514 if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
515 (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
516 emcpt += cPt;
517 ++cemc;
d03084cd 518 }
04905231 519
520 jet->AddTrackAt(tid, nt);
521 ++nt;
522 } else if (fClus) { // cluster constituent
111318e8 523 Int_t cid = -uid - 100;
04905231 524 AliVCluster *c = 0;
525 Double_t cEta=0,cPhi=0,cPt=0,cP=0;
526 if (fIsEmcPart) {
527 AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
528 if (!ep)
529 continue;
530 c = ep->GetCluster();
531 if (!c)
532 continue;
7030f36f 533 if (jetCount < fMarkConst)
6a20534a 534 ep->SetBit(fJetType);
04905231 535 cEta = ep->Eta();
536 cPhi = ep->Phi();
537 cPt = ep->Pt();
538 cP = ep->P();
539 } else {
540 c = static_cast<AliVCluster*>(fClus->At(cid));
541 if (!c)
542 continue;
7030f36f 543 if (jetCount < fMarkConst)
6a20534a 544 c->SetBit(fJetType);
04905231 545 TLorentzVector nP;
546 c->GetMomentum(nP, vertex);
547 cEta = nP.Eta();
548 cPhi = nP.Phi();
549 cPt = nP.Pt();
550 cP = nP.P();
551 }
552
553 neutralE += cP;
554 if (cPt > maxNe)
555 maxNe = cPt;
556
a7477843 557 if (c->GetLabel() > fMinMCLabel) // MC particle
43032ce2 558 mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
04905231 559
560 if (cPhi<0)
561 cPhi += TMath::TwoPi();
562 cPhi *= TMath::RadToDeg();
563 if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
564 (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
565 emcpt += cPt;
566 ++cemc;
567 }
568
569 jet->AddClusterAt(cid, nc);
570 ++nc;
571 ++nneutral;
02c7e943 572 } else {
573 AliError(Form("%s: No logical way to end up here.", GetName()));
574 continue;
48d874ff 575 }
04905231 576 } /* end of constituent loop */
577
35789a2d 578 jet->SetNumberOfTracks(nt);
579 jet->SetNumberOfClusters(nc);
f5f3c8e9 580 jet->SortConstituents();
111318e8 581 jet->SetMaxNeutralPt(maxNe);
582 jet->SetMaxChargedPt(maxCh);
b65fac7a 583 jet->SetNEF(neutralE / jet->E());
a88b5c99 584 fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
585 jet->SetArea(area.perp());
586 jet->SetAreaEta(area.eta());
587 jet->SetAreaPhi(area.phi());
111318e8 588 jet->SetNumberOfCharged(ncharged);
589 jet->SetNumberOfNeutrals(nneutral);
02c7e943 590 jet->SetMCPt(mcpt);
04905231 591 jet->SetNEmc(cemc);
592 jet->SetPtEmc(emcpt);
593
b65fac7a 594 if (gall > 0)
595 jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9 596 else
597 jet->SetAreaEmc(-1);
b65fac7a 598 if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65 599 (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
600 (jet->Eta() > geom->GetArm1EtaMin()) &&
601 (jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9 602 jet->SetAxisInEmcal(kTRUE);
ca5c29fa 603
604 AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
48d874ff 605 jetCount++;
606 }
7030f36f 607 //fJets->Sort();
8082a80b 608
609 //run constituent subtractor if requested
610 if(fDoConstituentSubtraction) {
289a9376 611#ifdef FASTJET_CONTRIB
8082a80b 612 if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
613 else {
614 std::vector<fastjet::PseudoJet> jets_sub;
615 jets_sub = fjw.GetConstituentSubtrJets();
616 AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
617 for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
618 //Only storing 4-vector and jet area of unsubtracted jet
619 AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount])
620 AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
621 jet_sub->SetLabel(ijet);
622 fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
623 jet_sub->SetArea(area.perp());
624 jet_sub->SetAreaEta(area.eta());
625 jet_sub->SetAreaPhi(area.phi());
626 jet_sub->SetAreaEmc(area.perp());
627 jetCount++;
628 }
629 }
630#endif
631 } //constituent subtraction
48d874ff 632}
04905231 633
ca5c29fa 634//________________________________________________________________________
635Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
636{
637 // Get the leading jets.
638
639 static Float_t pt[9999] = {0};
640
641 const Int_t n = (Int_t)array.size();
642
643 if (n < 1)
644 return kFALSE;
645
646 for (Int_t i = 0; i < n; i++)
647 pt[i] = array[i].perp();
648
649 TMath::Sort(n, pt, indexes);
650
651 return kTRUE;
652}
653
04905231 654//________________________________________________________________________
655Bool_t AliEmcalJetTask::DoInit()
656{
657 // Init. Return true if successful.
658
6421eeb0 659 if (fTrackEfficiency < 1.) {
660 if (gRandom) delete gRandom;
661 gRandom = new TRandom3(0);
662 }
663
04905231 664 // get input collections
665 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
666
667 // get the event
668 fEvent = InputEvent();
669 if (!fEvent) {
670 AliError(Form("%s: Could not retrieve event! Returning", GetName()));
671 return 0;
672 }
673
674 // add jets to event if not yet there
675 if (!(fEvent->FindListObject(fJetsName)))
676 fEvent->AddObject(fJets);
677 else {
f4b49fa6 678 AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
04905231 679 return 0;
680 }
681
8082a80b 682 if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
683 fEvent->AddObject(fJetsSub);
684
da67b88a 685 if (fTracksName == "Tracks")
686 am->LoadBranch("Tracks");
687 if (!fTracks && !fTracksName.IsNull()) {
688 fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
689 if (!fTracks) {
690 AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
691 return 0;
04905231 692 }
da67b88a 693 }
694 if (fTracks) {
695 TClass cls(fTracks->GetClass()->GetName());
5be3857d 696 if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
04905231 697 fIsMcPart = 1;
698 }
da67b88a 699
700 if (fCaloName == "CaloClusters")
701 am->LoadBranch("CaloClusters");
702 if (!fClus && !fCaloName.IsNull()) {
703 fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
704 if (!fClus) {
705 AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
706 return 0;
04905231 707 }
da67b88a 708 }
709 if (fClus) {
710 TClass cls(fClus->GetClass()->GetName());
04905231 711 if (cls.InheritsFrom("AliEmcalParticle"))
712 fIsEmcPart = 1;
713 }
714
8082a80b 715 if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
716 fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
717 if (!fRhoParam) {
718 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
719 return 0;
720 }
721 }
722 if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
723 fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
724 if (!fRhomParam) {
725 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
726 return 0;
727 }
728 }
729
04905231 730 return 1;
731}