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