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