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