]> 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
Adapt add macro and particle/cluster containers for the analysis on jets
[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();
b9ccc456 422 fjw.DoGenericSubtractionJetSigma2();
0d13a63c 423 fjw.DoGenericSubtractionJetConstituent();
424 fjw.DoGenericSubtractionJetLeSub();
3ec5da8e 425 }
426
8082a80b 427 //run constituent subtractor
428 if(fDoConstituentSubtraction) {
429 fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
430 fjw.DoConstituentSubtraction();
431 }
432
f5f3c8e9 433 // get geometry
434 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
435 if (!geom) {
02c7e943 436 AliFatal(Form("%s: Can not create geometry", GetName()));
f5f3c8e9 437 return;
438 }
04905231 439
440 // loop over fastjet jets
48d874ff 441 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
ca5c29fa 442 // sort jets according to jet pt
443 static Int_t indexes[9999] = {-1};
444 GetSortedArray(indexes, jets_incl);
445
446 AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
447 for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
448 Int_t ij = indexes[ijet];
449 AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
450
f5f3c8e9 451 if (jets_incl[ij].perp()<fMinJetPt)
452 continue;
453 if (fjw.GetJetArea(ij)<fMinJetArea)
48d874ff 454 continue;
c9ad772c 455 if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
456 (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
457 continue;
04905231 458
914d486c 459 AliEmcalJet *jet = new ((*fJets)[jetCount])
460 AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
8082a80b 461 jet->SetLabel(ij);
462
463 //do generic subtraction if requested
93daf3f7 464#ifdef FASTJET_VERSION
726c9488 465 if(fDoGenericSubtractionJetMass) {
8082a80b 466 std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
1fd2a48e 467 Int_t n = (Int_t)jetMassInfo.size();
8082a80b 468 if(n>ij && n>0) {
726c9488 469 jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
470 jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
471 jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
472 jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
8082a80b 473 }
8082a80b 474 }
3ec5da8e 475
726c9488 476 //here do generic subtraction for angular structure function
477 Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
478 fRMax = fRadius+0.2;
479 if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
480 fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
481 fjw.SetRMaxAndStep(fRMax,fDRStep);
482 fjw.DoGenericSubtractionGR(ij);
483 std::vector<double> num = fjw.GetGRNumerator();
484 std::vector<double> den = fjw.GetGRDenominator();
485 std::vector<double> nums = fjw.GetGRNumeratorSub();
486 std::vector<double> dens = fjw.GetGRDenominatorSub();
487 //pass this to AliEmcalJet
488 jet->SetGRNumSize(num.size());
489 jet->SetGRDenSize(den.size());
490 jet->SetGRNumSubSize(nums.size());
491 jet->SetGRDenSubSize(dens.size());
492 Int_t nsize = (Int_t)num.size();
493 for(Int_t g = 0; g<nsize; ++g) {
494 jet->AddGRNumAt(num[g],g);
495 jet->AddGRNumSubAt(nums[g],g);
496 }
497 Int_t dsize = (Int_t)den.size();
498 for(Int_t g = 0; g<dsize; ++g) {
499 jet->AddGRDenAt(den[g],g);
500 jet->AddGRDenSubAt(dens[g],g);
501 }
502 }
0d13a63c 503
504 if(fDoGenericSubtractionExtraJetShapes){
0d13a63c 505 std::vector<fastjet::contrib::GenericSubtractorInfo> jetAngularityInfo = fjw.GetGenSubtractorInfoJetAngularity();
0d13a63c 506 Int_t na = (Int_t)jetAngularityInfo.size();
507 if(na>ij && na>0) {
508 jet->SetFirstDerivativeAngularity(jetAngularityInfo[ij].first_derivative());
509 jet->SetSecondDerivativeAngularity(jetAngularityInfo[ij].second_derivative());
510 jet->SetFirstOrderSubtractedAngularity(jetAngularityInfo[ij].first_order_subtracted());
511 jet->SetSecondOrderSubtractedAngularity(jetAngularityInfo[ij].second_order_subtracted());
3ec5da8e 512 }
0d13a63c 513
514 std::vector<fastjet::contrib::GenericSubtractorInfo> jetpTDInfo = fjw.GetGenSubtractorInfoJetpTD();
0d13a63c 515 Int_t np = (Int_t)jetpTDInfo.size();
0d13a63c 516 if(np>ij && np>0) {
517 jet->SetFirstDerivativepTD(jetpTDInfo[ij].first_derivative());
518 jet->SetSecondDerivativepTD(jetpTDInfo[ij].second_derivative());
519 jet->SetFirstOrderSubtractedpTD(jetpTDInfo[ij].first_order_subtracted());
520 jet->SetSecondOrderSubtractedpTD(jetpTDInfo[ij].second_order_subtracted());
0d13a63c 521 }
522
523 std::vector<fastjet::contrib::GenericSubtractorInfo> jetCircularityInfo = fjw.GetGenSubtractorInfoJetCircularity();
3ec5da8e 524 Int_t nc = (Int_t)jetCircularityInfo.size();
0d13a63c 525 if(nc>ij && nc>0) {
526 jet->SetFirstDerivativeCircularity(jetCircularityInfo[ij].first_derivative());
527 jet->SetSecondDerivativeCircularity(jetCircularityInfo[ij].second_derivative());
528 jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted());
529 jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted());
3ec5da8e 530 }
b9ccc456 531
532 std::vector<fastjet::contrib::GenericSubtractorInfo> jetSigma2Info = fjw.GetGenSubtractorInfoJetSigma2();
533 Int_t ns = (Int_t)jetSigma2Info.size();
534 if(ns>ij && ns>0) {
535 jet->SetFirstDerivativeSigma2(jetSigma2Info[ij].first_derivative());
536 jet->SetSecondDerivativeSigma2(jetSigma2Info[ij].second_derivative());
537 jet->SetFirstOrderSubtractedSigma2(jetSigma2Info[ij].first_order_subtracted());
538 jet->SetSecondOrderSubtractedSigma2(jetSigma2Info[ij].second_order_subtracted());
539 }
540
3ec5da8e 541
0d13a63c 542 std::vector<fastjet::contrib::GenericSubtractorInfo> jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent();
3ec5da8e 543 Int_t nco= (Int_t)jetConstituentInfo.size();
0d13a63c 544 if(nco>ij && nco>0) {
545 jet->SetFirstDerivativeConstituent(jetConstituentInfo[ij].first_derivative());
546 jet->SetSecondDerivativeConstituent(jetConstituentInfo[ij].second_derivative());
547 jet->SetFirstOrderSubtractedConstituent(jetConstituentInfo[ij].first_order_subtracted());
548 jet->SetSecondOrderSubtractedConstituent(jetConstituentInfo[ij].second_order_subtracted());
0d13a63c 549 }
3ec5da8e 550
551 std::vector<fastjet::contrib::GenericSubtractorInfo> jetLeSubInfo = fjw.GetGenSubtractorInfoJetLeSub();
552 Int_t nlsub= (Int_t)jetLeSubInfo.size();
0d13a63c 553 if(nlsub>ij && nlsub>0) {
554 jet->SetFirstDerivativeLeSub(jetLeSubInfo[ij].first_derivative());
555 jet->SetSecondDerivativeLeSub(jetLeSubInfo[ij].second_derivative());
556 jet->SetFirstOrderSubtractedLeSub(jetLeSubInfo[ij].first_order_subtracted());
557 jet->SetSecondOrderSubtractedLeSub(jetLeSubInfo[ij].second_order_subtracted());
0d13a63c 558 }
0d13a63c 559 }
726c9488 560#endif
04905231 561
b9ccc456 562 // Fill constituent info
c64cb1f6 563 std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
04905231 564 jet->SetNumberOfTracks(constituents.size());
565 jet->SetNumberOfClusters(constituents.size());
566
b65fac7a 567 Int_t nt = 0;
568 Int_t nc = 0;
f5f3c8e9 569 Double_t neutralE = 0;
111318e8 570 Double_t maxCh = 0;
571 Double_t maxNe = 0;
f5f3c8e9 572 Int_t gall = 0;
573 Int_t gemc = 0;
04905231 574 Int_t cemc = 0;
111318e8 575 Int_t ncharged = 0;
576 Int_t nneutral = 0;
02c7e943 577 Double_t mcpt = 0;
04905231 578 Double_t emcpt = 0;
b65fac7a 579
b00cebf9 580 FillJetConstituents(constituents,jet,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents,0);
35789a2d 581 jet->SetNumberOfTracks(nt);
582 jet->SetNumberOfClusters(nc);
f5f3c8e9 583 jet->SortConstituents();
111318e8 584 jet->SetMaxNeutralPt(maxNe);
585 jet->SetMaxChargedPt(maxCh);
b65fac7a 586 jet->SetNEF(neutralE / jet->E());
a88b5c99 587 fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
588 jet->SetArea(area.perp());
589 jet->SetAreaEta(area.eta());
590 jet->SetAreaPhi(area.phi());
111318e8 591 jet->SetNumberOfCharged(ncharged);
592 jet->SetNumberOfNeutrals(nneutral);
02c7e943 593 jet->SetMCPt(mcpt);
04905231 594 jet->SetNEmc(cemc);
595 jet->SetPtEmc(emcpt);
596
b65fac7a 597 if (gall > 0)
598 jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
f5f3c8e9 599 else
600 jet->SetAreaEmc(-1);
b65fac7a 601 if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
1772fe65 602 (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
603 (jet->Eta() > geom->GetArm1EtaMin()) &&
604 (jet->Eta() < geom->GetArm1EtaMax()))
f5f3c8e9 605 jet->SetAxisInEmcal(kTRUE);
ca5c29fa 606
607 AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
48d874ff 608 jetCount++;
609 }
7030f36f 610 //fJets->Sort();
8082a80b 611
612 //run constituent subtractor if requested
613 if(fDoConstituentSubtraction) {
93daf3f7 614#ifdef FASTJET_VERSION
8082a80b 615 if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
616 else {
617 std::vector<fastjet::PseudoJet> jets_sub;
618 jets_sub = fjw.GetConstituentSubtrJets();
619 AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
620 for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
621 //Only storing 4-vector and jet area of unsubtracted jet
622 AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount])
623 AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
624 jet_sub->SetLabel(ijet);
cfa7b4df 625 // Fill constituent info
b00cebf9 626 std::vector<fastjet::PseudoJet> constituents_unsub(fjw.GetJetConstituents(ijet));
627 std::vector<fastjet::PseudoJet> constituents_sub = jets_sub[ijet].constituents();
628 jet_sub->SetNumberOfTracks(constituents_sub.size());
df20822b 629 jet_sub->SetNumberOfClusters(constituents_sub.size());
630 Int_t nt = 0;
631 Int_t nc = 0;
632 Double_t neutralE = 0;
633 Double_t maxCh = 0;
634 Double_t maxNe = 0;
635 Int_t gall = 0;
636 Int_t gemc = 0;
637 Int_t cemc = 0;
638 Int_t ncharged = 0;
639 Int_t nneutral = 0;
640 Double_t mcpt = 0;
641 Double_t emcpt = 0;
b00cebf9 642
643 FillJetConstituents(constituents_sub,jet_sub,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents_unsub,1);
df20822b 644 jet_sub->SetNumberOfTracks(nt);
645 jet_sub->SetNumberOfClusters(nc);
646 jet_sub->SortConstituents();
647
8082a80b 648 fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
649 jet_sub->SetArea(area.perp());
650 jet_sub->SetAreaEta(area.eta());
651 jet_sub->SetAreaPhi(area.phi());
652 jet_sub->SetAreaEmc(area.perp());
653 jetCount++;
654 }
655 }
656#endif
657 } //constituent subtraction
48d874ff 658}
04905231 659
ca5c29fa 660//________________________________________________________________________
661Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
662{
663 // Get the leading jets.
664
665 static Float_t pt[9999] = {0};
666
667 const Int_t n = (Int_t)array.size();
668
669 if (n < 1)
670 return kFALSE;
671
672 for (Int_t i = 0; i < n; i++)
673 pt[i] = array[i].perp();
674
675 TMath::Sort(n, pt, indexes);
676
677 return kTRUE;
678}
679
04905231 680//________________________________________________________________________
681Bool_t AliEmcalJetTask::DoInit()
682{
683 // Init. Return true if successful.
684
6421eeb0 685 if (fTrackEfficiency < 1.) {
686 if (gRandom) delete gRandom;
687 gRandom = new TRandom3(0);
688 }
689
04905231 690 // get input collections
691 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
692
693 // get the event
694 fEvent = InputEvent();
695 if (!fEvent) {
696 AliError(Form("%s: Could not retrieve event! Returning", GetName()));
697 return 0;
698 }
699
8082a80b 700 if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
701 fEvent->AddObject(fJetsSub);
702
da67b88a 703 if (fTracksName == "Tracks")
704 am->LoadBranch("Tracks");
705 if (!fTracks && !fTracksName.IsNull()) {
706 fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
707 if (!fTracks) {
708 AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
709 return 0;
04905231 710 }
da67b88a 711 }
712 if (fTracks) {
713 TClass cls(fTracks->GetClass()->GetName());
5be3857d 714 if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
04905231 715 fIsMcPart = 1;
716 }
da67b88a 717
718 if (fCaloName == "CaloClusters")
719 am->LoadBranch("CaloClusters");
720 if (!fClus && !fCaloName.IsNull()) {
721 fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
722 if (!fClus) {
723 AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
724 return 0;
04905231 725 }
da67b88a 726 }
727 if (fClus) {
728 TClass cls(fClus->GetClass()->GetName());
04905231 729 if (cls.InheritsFrom("AliEmcalParticle"))
730 fIsEmcPart = 1;
731 }
732
8082a80b 733 if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
734 fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
735 if (!fRhoParam) {
736 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
737 return 0;
738 }
739 }
740 if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
741 fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
742 if (!fRhomParam) {
743 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
744 return 0;
745 }
746 }
cbbb66b0 747
748 // add jets to event if not yet there
749 if (!(fEvent->FindListObject(fJetsName)))
750 fEvent->AddObject(fJets);
751 else {
752 AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
753 return 0;
754 }
8082a80b 755
04905231 756 return 1;
757}
b9ccc456 758
759//___________________________________________________________________________________________________________________
b00cebf9 760void AliEmcalJetTask::FillJetConstituents(std::vector<fastjet::PseudoJet>& constituents,AliEmcalJet *jet,Double_t vertex[3],UInt_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc,std::vector<fastjet::PseudoJet>& constituentsunsub,Int_t flagsub){
b9ccc456 761 nt = 0;
762 nc = 0;
763 neutralE = 0;
764 maxCh = 0;
765 maxNe = 0;
766 gall = 0;
767 gemc =0;
768 cemc = 0;
769 ncharged = 0;
770 nneutral = 0;
771 mcpt = 0;
772 emcpt = 0;
b00cebf9 773 Int_t uid = -1;
b9ccc456 774 AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
775 for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
b00cebf9 776 if(flagsub==0) uid = constituents[ic].user_index();
777 if(flagsub!=0) uid=GetIndexSub(constituents[ic].perp(),constituentsunsub);
b9ccc456 778 AliDebug(3,Form("Processing constituent %d", uid));
779 if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
780 ++gall;
781 Double_t gphi = constituents[ic].phi();
782 if (gphi<0)
783 gphi += TMath::TwoPi();
784 gphi *= TMath::RadToDeg();
785 Double_t geta = constituents[ic].eta();
786 if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
787 (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
788 ++gemc;
789 } else if ((uid > 0) && fTracks) { // track constituent
790 Int_t tid = uid - 100;
791 AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
792 if (!t) {
793 AliError(Form("Could not find track %d",tid));
794 continue;
795 }
796 if (jetCount < fMarkConst) {
797 AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
798 t->SetBit(fJetType);
799 }
800 Double_t cEta = t->Eta();
801 Double_t cPhi = t->Phi();
802 Double_t cPt = t->Pt();
803 Double_t cP = t->P();
804 if (t->Charge() == 0) {
805 neutralE += cP;
806 ++nneutral;
807 if (cPt > maxNe)
808 maxNe = cPt;
809 } else {
810 ++ncharged;
811 if (cPt > maxCh)
812 maxCh = cPt;
813 }
814 if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
815 mcpt += cPt;
816
817 if (cPhi<0)
818 cPhi += TMath::TwoPi();
819 cPhi *= TMath::RadToDeg();
820 if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
821 (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
822 emcpt += cPt;
823 ++cemc;
824 }
b00cebf9 825 // cout<<"Adding the track"<<endl;
b9ccc456 826 jet->AddTrackAt(tid, nt);
827 ++nt;
828 } else if (fClus) { // cluster constituent
829 Int_t cid = -uid - 100;
830 AliVCluster *c = 0;
831 Double_t cEta=0,cPhi=0,cPt=0,cP=0;
832 if (fIsEmcPart) {
833 AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
834 if (!ep)
835 continue;
836 c = ep->GetCluster();
837 if (!c)
838 continue;
839 if (jetCount < fMarkConst)
840 ep->SetBit(fJetType);
841 cEta = ep->Eta();
842 cPhi = ep->Phi();
843 cPt = ep->Pt();
844 cP = ep->P();
845 } else {
846 c = static_cast<AliVCluster*>(fClus->At(cid));
847 if (!c)
848 continue;
849 if (jetCount < fMarkConst)
850 c->SetBit(fJetType);
851 TLorentzVector nP;
852 c->GetMomentum(nP, vertex);
853 cEta = nP.Eta();
854 cPhi = nP.Phi();
855 cPt = nP.Pt();
856 cP = nP.P();
857 }
858
859 neutralE += cP;
860 if (cPt > maxNe)
861 maxNe = cPt;
862
863 if (c->GetLabel() > fMinMCLabel) // MC particle
864 mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
865
866 if (cPhi<0)
867 cPhi += TMath::TwoPi();
868 cPhi *= TMath::RadToDeg();
869 if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
870 (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
871 emcpt += cPt;
872 ++cemc;
873 }
874
875 jet->AddClusterAt(cid, nc);
876 ++nc;
877 ++nneutral;
878 } else {
879 AliError(Form("%s: No logical way to end up here.", GetName()));
880 continue;
cfa7b4df 881 }
b9ccc456 882 }
b9ccc456 883}
884//______________________________________________________________________________________
b00cebf9 885Int_t AliEmcalJetTask::GetIndexSub(Double_t ptsub,std::vector<fastjet::PseudoJet>& constituentsunsub){
886
887 for(UInt_t ii=0;ii<constituentsunsub.size();ii++){
888 Double_t ptunsub=constituentsunsub[ii].perp();
889 //cout<<ptunsub<<" "<<ptsub<<endl;
890 if(ptsub==ptunsub) return constituentsunsub[ii].user_index();
891
892 } return -1;}
893//______________________________________________________________________________________