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