Added functionality for combining methodes/pid.
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetEmbeddingFromAODTask.cxx
CommitLineData
b16bb001 1// $Id: AliJetEmbeddingFromAODTask.cxx $
2//
3// Jet embedding from AOD task.
4//
5// Author: S.Aiola, C.Loizides
6
7#include "AliJetEmbeddingFromAODTask.h"
8
4358e58a 9// C++ standard library
10#include <vector>
11
12// ROOT
b16bb001 13#include <TFile.h>
14#include <TTree.h>
15#include <TClonesArray.h>
16#include <TObjArray.h>
17#include <TObjString.h>
18#include <TGrid.h>
19#include <TH2C.h>
cd6431de 20#include <TList.h>
21#include <TStreamerInfo.h>
6a20534a 22#include <TRandom.h>
2103dc6a 23#include <TSystem.h>
4358e58a 24#include <TLorentzVector.h>
b16bb001 25
4358e58a 26// AliRoot
b16bb001 27#include "AliVEvent.h"
28#include "AliAODTrack.h"
29#include "AliESDtrack.h"
30#include "AliPicoTrack.h"
31#include "AliVTrack.h"
32#include "AliVCluster.h"
33#include "AliVCaloCells.h"
787a3c4f 34#include "AliAODMCParticle.h"
b16bb001 35#include "AliCentrality.h"
36#include "AliVHeader.h"
37#include "AliVVertex.h"
38#include "AliAODHeader.h"
4358e58a 39#include "AliFJWrapper.h"
b16bb001 40#include "AliLog.h"
41#include "AliInputEventHandler.h"
42
43ClassImp(AliJetEmbeddingFromAODTask)
44
45//________________________________________________________________________
46AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
47 AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
48 fFileList(0),
6a20534a 49 fRandomAccess(kFALSE),
b16bb001 50 fAODTreeName(),
51 fAODHeaderName(),
52 fAODVertexName(),
53 fAODTrackName(),
54 fAODClusName(),
55 fAODCellsName(),
787a3c4f 56 fAODMCParticlesName(),
b16bb001 57 fMinCentrality(0),
58 fMaxCentrality(10),
59 fTriggerMask(AliVEvent::kAny),
60 fZVertexCut(10),
4358e58a 61 fJetMinPt(0),
62 fJetMinEta(-0.5),
63 fJetMaxEta(0.5),
64 fJetMinPhi(-999),
65 fJetMaxPhi(999),
05077f28 66 fJetConstituentMinPt(0),
4358e58a 67 fJetRadius(0.4),
68 fJetType(0),
69 fJetAlgo(1),
70 fJetParticleLevel(kTRUE),
9adcb46d 71 fIncludeNoITS(kFALSE),
72 fCutMaxFractionSharedTPCClusters(0.4),
5ce8ae64 73 fUseNegativeLabels(kTRUE),
74 fTrackEfficiency(1),
4358e58a 75 fIsAODMC(kFALSE),
6a20534a 76 fTotalFiles(2050),
77 fAttempts(5),
b16bb001 78 fEsdTreeMode(kFALSE),
79 fCurrentFileID(0),
6a20534a 80 fCurrentAODFileID(0),
b16bb001 81 fCurrentAODFile(0),
cd6431de 82 fPicoTrackVersion(0),
6a20534a 83 fCurrentAODTree(0),
b16bb001 84 fAODHeader(0),
85 fAODVertex(0),
86 fAODTracks(0),
87 fAODClusters(0),
88 fAODCaloCells(0),
787a3c4f 89 fAODMCParticles(0),
6a20534a 90 fCurrentAODEntry(-1),
91 fHistFileMatching(0),
92 fHistAODFileError(0),
2103dc6a 93 fHistNotEmbedded(0),
94 fHistEmbeddingQA(0)
b16bb001 95{
96 // Default constructor.
97 SetSuffix("AODEmbedding");
b16bb001 98 fAODfilterBits[0] = -1;
99 fAODfilterBits[1] = -1;
5be3857d 100 fEtaMin = -1;
101 fEtaMax = 1;
102 fPhiMin = -10;
103 fPhiMax = 10;
104 fPtMin = 0;
105 fPtMax = 1000;
b16bb001 106}
107
108//________________________________________________________________________
109AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) :
110 AliJetModelBaseTask(name, drawqa),
111 fFileList(0),
6a20534a 112 fRandomAccess(kFALSE),
b16bb001 113 fAODTreeName("aodTree"),
114 fAODHeaderName("header"),
115 fAODVertexName("vertices"),
116 fAODTrackName("tracks"),
787a3c4f 117 fAODClusName(""),
b16bb001 118 fAODCellsName("emcalCells"),
787a3c4f 119 fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
b16bb001 120 fMinCentrality(0),
121 fMaxCentrality(10),
122 fTriggerMask(AliVEvent::kAny),
123 fZVertexCut(10),
4358e58a 124 fJetMinPt(0),
125 fJetMinEta(-0.5),
126 fJetMaxEta(0.5),
127 fJetMinPhi(-999),
128 fJetMaxPhi(999),
05077f28 129 fJetConstituentMinPt(0),
4358e58a 130 fJetRadius(0.4),
131 fJetType(0),
132 fJetAlgo(1),
133 fJetParticleLevel(kTRUE),
9adcb46d 134 fIncludeNoITS(kFALSE),
135 fCutMaxFractionSharedTPCClusters(0.4),
5ce8ae64 136 fUseNegativeLabels(kTRUE),
137 fTrackEfficiency(1),
4358e58a 138 fIsAODMC(kFALSE),
6a20534a 139 fTotalFiles(2050),
140 fAttempts(5),
b16bb001 141 fEsdTreeMode(kFALSE),
142 fCurrentFileID(0),
6a20534a 143 fCurrentAODFileID(0),
b16bb001 144 fCurrentAODFile(0),
cd6431de 145 fPicoTrackVersion(0),
6a20534a 146 fCurrentAODTree(0),
b16bb001 147 fAODHeader(0),
148 fAODVertex(0),
149 fAODTracks(0),
150 fAODClusters(0),
787a3c4f 151 fAODCaloCells(0),
152 fAODMCParticles(0),
6a20534a 153 fCurrentAODEntry(-1),
154 fHistFileMatching(0),
155 fHistAODFileError(0),
2103dc6a 156 fHistNotEmbedded(0),
157 fHistEmbeddingQA(0)
b16bb001 158{
159 // Standard constructor.
160 SetSuffix("AODEmbedding");
b16bb001 161 fAODfilterBits[0] = -1;
162 fAODfilterBits[1] = -1;
5be3857d 163 fEtaMin = -1;
164 fEtaMax = 1;
165 fPhiMin = -10;
166 fPhiMax = 10;
167 fPtMin = 0;
168 fPtMax = 1000;
b16bb001 169}
170
171//________________________________________________________________________
172AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
173{
174 // Destructor
175
176 if (fCurrentAODFile) {
177 fCurrentAODFile->Close();
178 delete fCurrentAODFile;
179 }
180}
181
182//________________________________________________________________________
183void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
184{
185 if (!fQAhistos)
186 return;
187
188 AliJetModelBaseTask::UserCreateOutputObjects();
189
6a20534a 190 fHistFileMatching = new TH2C("fHistFileMatching", "fHistFileMatching", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
191 fHistFileMatching->GetXaxis()->SetTitle("File no. (PYTHIA)");
192 fHistFileMatching->GetYaxis()->SetTitle("File no. (Embedded AOD)");
193 fHistFileMatching->GetZaxis()->SetTitle("counts");
194 fOutput->Add(fHistFileMatching);
195
196 fHistAODFileError = new TH1C("fHistAODFileError", "fHistAODFileError", fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
197 fHistAODFileError->GetXaxis()->SetTitle("File no. (Embedded AOD)");
198 fHistAODFileError->GetYaxis()->SetTitle("counts");
199 fOutput->Add(fHistAODFileError);
200
201 fHistNotEmbedded = new TH1C("fHistNotEmbedded", "fHistNotEmbedded", fTotalFiles, -0.5, fTotalFiles - 0.5);
202 fHistNotEmbedded->GetXaxis()->SetTitle("File no. (PYTHIA)");
203 fHistNotEmbedded->GetYaxis()->SetTitle("counts");
204 fOutput->Add(fHistNotEmbedded);
b16bb001 205
2103dc6a 206 fHistEmbeddingQA = new TH1F("fHistEmbeddingQA", "fHistEmbeddingQA", 2, 0, 2);
207 fHistEmbeddingQA->GetXaxis()->SetTitle("Event state");
208 fHistEmbeddingQA->GetYaxis()->SetTitle("counts");
209 fHistEmbeddingQA->GetXaxis()->SetBinLabel(1, "OK");
210 fHistEmbeddingQA->GetXaxis()->SetBinLabel(2, "Not embedded");
211 fOutput->Add(fHistEmbeddingQA);
212
b16bb001 213 PostData(1, fOutput);
214}
215
216//________________________________________________________________________
217Bool_t AliJetEmbeddingFromAODTask::UserNotify()
218{
219 TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
220 if (path.EndsWith("/"))
221 path.Remove(path.Length()-1);
222 path.Remove(path.Last('/'));
223 path.Remove(0, path.Last('/')+1);
224 if (!path.IsDec()) {
225 AliWarning(Form("Could not extract file number from path %s", path.Data()));
226 return kTRUE;
227 }
228 fCurrentFileID = path.Atoi();
6a20534a 229 if (!fRandomAccess) {
230 fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles-1;
231 AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
232 }
b16bb001 233 return kTRUE;
234}
235
236//________________________________________________________________________
237Bool_t AliJetEmbeddingFromAODTask::ExecOnce()
238{
239 if (fAODTreeName.Contains("aod"))
240 fEsdTreeMode = kFALSE;
241 else
242 fEsdTreeMode = kTRUE;
243
244 return AliJetModelBaseTask::ExecOnce();
245}
246
247//________________________________________________________________________
248Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
249{
250 if (fCurrentAODFile) {
251 fCurrentAODFile->Close();
252 delete fCurrentAODFile;
253 fCurrentAODFile = 0;
254 }
255
6a20534a 256 Int_t i = 0;
b16bb001 257
2103dc6a 258 while ((!fCurrentAODFile || fCurrentAODFile->IsZombie()) && i < fAttempts) {
259 if (i > 0 && fHistAODFileError) {
6a20534a 260 fHistAODFileError->Fill(fCurrentAODFileID);
261 }
262
2103dc6a 263 fCurrentAODFile = GetNextFile();
7030f36f 264 i++;
2103dc6a 265 }
b16bb001 266
6a20534a 267 if (!fCurrentAODFile || fCurrentAODFile->IsZombie())
b16bb001 268 return kFALSE;
b16bb001 269
cd6431de 270 const TList *clist = fCurrentAODFile->GetStreamerInfoCache();
271 if(clist) {
272 TStreamerInfo *cinfo = static_cast<TStreamerInfo*>(clist->FindObject("AliPicoTrack"));
273 if(cinfo)
274 fPicoTrackVersion = cinfo->GetClassVersion();
275 else
276 fPicoTrackVersion = 0;
277 }
278
6a20534a 279 fCurrentAODTree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
280 if (!fCurrentAODTree)
281 return kFALSE;
282
283 if (!fAODHeaderName.IsNull())
284 fCurrentAODTree->SetBranchAddress(fAODHeaderName, &fAODHeader);
285
286 if (!fAODVertexName.IsNull())
287 fCurrentAODTree->SetBranchAddress(fAODVertexName, &fAODVertex);
288
289 if (!fAODTrackName.IsNull())
290 fCurrentAODTree->SetBranchAddress(fAODTrackName, &fAODTracks);
291
292 if (!fAODClusName.IsNull())
293 fCurrentAODTree->SetBranchAddress(fAODClusName, &fAODClusters);
294
295 if (!fAODCellsName.IsNull())
296 fCurrentAODTree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
297
298 if (!fAODMCParticlesName.IsNull())
299 fCurrentAODTree->SetBranchAddress(fAODMCParticlesName, &fAODMCParticles);
300
301 if (fRandomAccess)
302 fCurrentAODEntry = TMath::Nint(gRandom->Rndm()*fCurrentAODTree->GetEntries());
303 else
304 fCurrentAODEntry = 0;
7030f36f 305
4358e58a 306 AliDebug(3,Form("Will start embedding from entry %d", fCurrentAODEntry));
6a20534a 307
7030f36f 308 if (fHistFileMatching)
6a20534a 309 fHistFileMatching->Fill(fCurrentFileID, fCurrentAODFileID-1);
b16bb001 310
b16bb001 311 return kTRUE;
312}
313
314//________________________________________________________________________
2103dc6a 315TFile* AliJetEmbeddingFromAODTask::GetNextFile()
b16bb001 316{
2103dc6a 317 if (fRandomAccess)
318 fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*fFileList->GetEntriesFast());
319 else
320 fCurrentAODFileID++;
321
322 if (fCurrentAODFileID >= fFileList->GetEntriesFast()) {
323 AliError("No more file in the list!");
324 return 0;
325 }
326
327 TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
328 TString fileName(objFileName->GetString());
b16bb001 329
4358e58a 330 if (fileName.BeginsWith("alien://") && !gGrid) {
331 AliInfo("Trying to connect to AliEn ...");
332 TGrid::Connect("alien://");
333 }
334
2103dc6a 335 TString baseFileName(fileName);
336 if (baseFileName.Contains(".zip#")) {
337 Ssiz_t pos = baseFileName.Last('#');
338 baseFileName.Remove(pos);
339 }
340
341 if (gSystem->AccessPathName(baseFileName)) {
342 AliDebug(3,Form("File %s does not exist!", baseFileName.Data()));
343 return 0;
344 }
345
346 AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
347 TFile *file = TFile::Open(fileName);
348
349 if (!file)
350 AliDebug(3,Form("Unable to open file: %s!", fileName.Data()));
351
352 return file;
6a20534a 353}
354
355//________________________________________________________________________
356Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
357{
b16bb001 358 Int_t attempts = 0;
359
360 do {
6a20534a 361 if (!fCurrentAODFile || !fCurrentAODTree || fCurrentAODEntry >= fCurrentAODTree->GetEntries()) {
b16bb001 362 if (!OpenNextFile())
363 return kFALSE;
b16bb001 364 }
365
6a20534a 366 fCurrentAODTree->GetEntry(fCurrentAODEntry);
367 fCurrentAODEntry++;
b16bb001 368 attempts++;
369
370 if (attempts == 1000)
371 AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
372
6a20534a 373 } while (!IsAODEventSelected() && fCurrentAODTree);
b16bb001 374
6a20534a 375 return (fCurrentAODTree!=0);
b16bb001 376}
377
378//________________________________________________________________________
379Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
380{
4358e58a 381 // AOD event selection.
382
b16bb001 383 if (!fEsdTreeMode && fAODHeader) {
384 AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
385
4358e58a 386 // Trigger selection
6a20534a 387 if (fTriggerMask != AliVEvent::kAny) {
388 UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
389
390 if ((offlineTrigger & fTriggerMask) == 0) {
4358e58a 391 AliDebug(2,Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
392 offlineTrigger, fTriggerMask));
6a20534a 393 return kFALSE;
394 }
b16bb001 395 }
396
4358e58a 397 // Centrality selection
6a20534a 398 if (fMinCentrality >= 0) {
399 AliCentrality *cent = aodHeader->GetCentralityP();
400 Float_t centVal = cent->GetCentralityPercentile("V0M");
401 if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
4358e58a 402 AliDebug(2,Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f",
403 centVal, fMinCentrality, fMaxCentrality));
6a20534a 404 return kFALSE;
405 }
b16bb001 406 }
407 }
408
4358e58a 409 // Vertex selection
b16bb001 410 if (fAODVertex) {
411 Double_t vert[3]={0};
412 ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
413 if (TMath::Abs(vert[2]) > fZVertexCut) {
4358e58a 414 AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
415 vert[2], fZVertexCut));
b16bb001 416 return kFALSE;
417 }
418 }
419
4358e58a 420 // Jet selection
421 if (fJetMinPt > 0) {
422 TLorentzVector jet;
423
424 if (fJetParticleLevel) {
425 if (fAODMCParticles)
426 jet = GetLeadingJet(fAODMCParticles);
427 else {
428 AliWarning("Particle level jets selected, but not MC particles found. The jet event selection will be skipped.");
429 return kTRUE;
430 }
431 }
432 else {
433 if (fAODTracks || fAODClusters)
434 jet = GetLeadingJet(fAODTracks, fAODClusters);
435 else {
436 AliWarning("Detector level jets selected, but not tracks or clusters found. The jet event selection will be skipped.");
437 return kTRUE;
438 }
439 }
440
05077f28 441 AliDebug(1, Form("Leading jet pt = %f", jet.Pt()));
4358e58a 442 if (jet.Pt() < fJetMinPt)
443 return kFALSE;
444 }
445
b16bb001 446 return kTRUE;
447}
448
449//________________________________________________________________________
450void AliJetEmbeddingFromAODTask::Run()
451{
452 if (!GetNextEntry()) {
7030f36f 453 if (fHistNotEmbedded)
454 fHistNotEmbedded->Fill(fCurrentFileID);
2103dc6a 455 if (fHistEmbeddingQA)
456 fHistEmbeddingQA->Fill("Not embedded", 1);
6a20534a 457 AliError("Unable to get the AOD event to embed. Nothing will be embedded.");
b16bb001 458 return;
459 }
460
2103dc6a 461 if (fHistEmbeddingQA)
462 fHistEmbeddingQA->Fill("OK", 1);
463
787a3c4f 464 if (fOutMCParticles) {
465
466 if (fCopyArray && fMCParticles)
467 CopyMCParticles();
468
469 if (fAODMCParticles) {
4358e58a 470 AliDebug(3, Form("%d MC particles will be processed for embedding.", fAODMCParticles->GetEntriesFast()));
787a3c4f 471 for (Int_t i = 0; i < fAODMCParticles->GetEntriesFast(); i++) {
472 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fAODMCParticles->At(i));
473 if (!part) {
474 AliError(Form("Could not find MC particle %d in branch %s of tree %s!", i, fAODMCParticlesName.Data(), fAODTreeName.Data()));
475 continue;
476 }
477
5be3857d 478 AliDebug(3, Form("Processing MC particle %d with pT = %f, eta = %f, phi = %f", i, part->Pt(), part->Eta(), part->Phi()));
479
787a3c4f 480 if (!part->IsPhysicalPrimary())
481 continue;
482
5be3857d 483 if (part->Pt() < fPtMin || part->Pt() > fPtMax ||
484 part->Eta() < fEtaMin || part->Eta() > fEtaMax ||
485 part->Phi() < fPhiMin || part->Phi() > fPhiMax)
486 continue;
487
787a3c4f 488 AddMCParticle(part, i);
5be3857d 489 AliDebug(3, "Embedded!");
787a3c4f 490 }
491 }
492 }
493
494 if (fOutTracks) {
b16bb001 495
787a3c4f 496 if (fCopyArray && fTracks)
b16bb001 497 CopyTracks();
498
4358e58a 499 AliDebug(3, Form("Start embedding with %d tracks.", fOutTracks->GetEntriesFast()));
5be3857d 500
b16bb001 501 if (fAODTracks) {
4358e58a 502 AliDebug(3, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
b16bb001 503 for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
504 AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
505 if (!track) {
506 AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
507 continue;
508 }
509
5be3857d 510 AliDebug(3, Form("Processing track %d with pT = %f, eta = %f, phi = %f, label = %d", i, track->Pt(), track->Eta(), track->Phi(), track->GetLabel()));
511
b16bb001 512 Int_t type = 0;
513 Bool_t isEmc = kFALSE;
514 if (!fEsdTreeMode) {
515 AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
516 if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
517 type = 0;
518 }
519 else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
520 if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
521 if (fIncludeNoITS)
522 type = 2;
5be3857d 523 else {
524 AliDebug(3, "Track not embedded because ITS refit failed.");
b16bb001 525 continue;
526 }
5be3857d 527 }
b16bb001 528 else {
529 type = 1;
530 }
531 }
532 else { /*not a good track*/
5be3857d 533 AliDebug(3, "Track not embedded because not an hybrid track.");
b16bb001 534 continue;
535 }
9adcb46d 536 if (fCutMaxFractionSharedTPCClusters > 0) {
537 Double_t frac = Double_t(aodtrack->GetTPCnclsS()) / Double_t(aodtrack->GetTPCncls());
538 if (frac > fCutMaxFractionSharedTPCClusters)
539 continue;
540 }
b16bb001 541 if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
542 aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
543 aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
544 isEmc = kTRUE;
545 }
cd6431de 546 else if (fPicoTrackVersion > 0) { /*not AOD mode, let's see if it is PicoTrack*/
b16bb001 547 AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
548 if (ptrack) {
cd6431de 549 if (fPicoTrackVersion >= 3)
550 type = ptrack->GetTrackType();
551 else
552 type = ptrack->GetLabel();
b16bb001 553 isEmc = ptrack->IsEMCAL();
554 }
555 }
556
5be3857d 557 if (track->Pt() < fPtMin || track->Pt() > fPtMax ||
558 track->Eta() < fEtaMin || track->Eta() > fEtaMax ||
559 track->Phi() < fPhiMin || track->Phi() > fPhiMax) {
560 AliDebug(3, "Track not embedded because out of limits.");
561 continue;
562 }
563
5ce8ae64 564 if (fTrackEfficiency < 1) {
565 Double_t r = gRandom->Rndm();
5be3857d 566 if (fTrackEfficiency < r) {
567 AliDebug(3, "Track not embedded because of artificial inefiiciency.");
5ce8ae64 568 continue;
4358e58a 569 }
5be3857d 570 }
5ce8ae64 571
7030f36f 572 Int_t label = 0;
4358e58a 573 if (fIsAODMC) {
5ce8ae64 574 if (fUseNegativeLabels)
575 label = track->GetLabel();
576 else
577 label = TMath::Abs(track->GetLabel());
578
4358e58a 579 if (label == 0)
cac0e10b 580 AliDebug(1,Form("%s: Track %d with label==0", GetName(), i));
7030f36f 581 }
582
7030f36f 583 AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, label);
5be3857d 584 AliDebug(3, "Track embedded!");
b16bb001 585 }
586 }
587 }
588
787a3c4f 589 if (fOutClusters) {
b16bb001 590
787a3c4f 591 if (fCopyArray && fClusters)
b16bb001 592 CopyClusters();
593
594 if (fAODClusters) {
595 for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
596 AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
597 if (!clus) {
598 AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
599 continue;
600 }
fde82e42 601
602 if (!clus->IsEMCAL())
603 continue;
604
b16bb001 605 TLorentzVector vect;
606 Double_t vert[3] = {0,0,0};
607 clus->GetMomentum(vect,vert);
5be3857d 608
609 if (vect.Pt() < fPtMin || vect.Pt() > fPtMax ||
610 vect.Eta() < fEtaMin || vect.Eta() > fEtaMax ||
611 vect.Phi() < fPhiMin || vect.Phi() > fPhiMax)
612 continue;
fde82e42 613
4358e58a 614 Int_t label = 0;
615 if (fIsAODMC) {
616 label = clus->GetLabel();
617 if (label <= 0)
cac0e10b 618 AliDebug(1,Form("%s: Clus %d with label<=0", GetName(), i));
4358e58a 619 }
507f74bc 620 AddCluster(clus);
b16bb001 621 }
622 }
623 }
624
787a3c4f 625 if (fOutCaloCells) {
626
4358e58a 627 Double_t totalEnergy = 0;
787a3c4f 628 Int_t totalCells = 0;
b16bb001 629
787a3c4f 630 if (fCaloCells)
631 totalCells += fCaloCells->GetNumberOfCells();
b16bb001 632 if (fAODCaloCells)
633 totalCells += fAODCaloCells->GetNumberOfCells();
634
635 SetNumberOfOutCells(totalCells);
636
787a3c4f 637 if (fCaloCells)
638 CopyCells();
b16bb001 639
640 if (fAODCaloCells) {
641 for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
60d77596 642 Int_t mclabel = 0;
b16bb001 643 Double_t efrac = 0.;
644 Double_t time = -1;
645 Short_t cellNum = -1;
646 Double_t amp = -1;
647
648 fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
4358e58a 649
650 if (fIsAODMC) {
651 if (mclabel <= 0)
cac0e10b 652 AliDebug(1,Form("%s: Cell %d with label<=0", GetName(), i));
4358e58a 653 }
654 else {
655 mclabel = 0;
656 }
657
658 AliDebug(2,Form("Adding cell with amplitude %f, absolute ID %d, time %f, mc label %d", amp, cellNum, time, mclabel));
787a3c4f 659 AddCell(amp, cellNum, time, mclabel);
4358e58a 660 totalEnergy += amp;
b16bb001 661 }
662 }
4358e58a 663 AliDebug(2,Form("Added cells = %d (energy = %f), total cells = %d", fAddedCells, totalEnergy, totalCells));
b16bb001 664 }
665}
4358e58a 666
667//________________________________________________________________________
668TLorentzVector AliJetEmbeddingFromAODTask::GetLeadingJet(TClonesArray *tracks, TClonesArray *clusters)
669{
670 TString name("kt");
671 fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
672 if (fJetAlgo == 1) {
673 name = "antikt";
674 jalgo = fastjet::antikt_algorithm;
675 }
676
677 // setup fj wrapper
678 AliFJWrapper fjw(name, name);
679 fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
680 fjw.SetGhostArea(1); // set a very large ghost area to speed up jet finding
681 fjw.SetR(fJetRadius);
682 fjw.SetAlgorithm(jalgo);
683 fjw.SetMaxRap(1);
684 fjw.Clear();
685
686 if (tracks) {
687 const Int_t Ntracks = tracks->GetEntries();
688 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
689 AliVParticle *t = static_cast<AliVParticle*>(tracks->At(iTracks));
690 if (!t)
691 continue;
05077f28 692
693 AliAODMCParticle *aodmcpart = dynamic_cast<AliAODMCParticle*>(t);
694 if (aodmcpart && !aodmcpart->IsPhysicalPrimary())
695 continue;
4358e58a 696
697 if ((fJetType == 1 && t->Charge() == 0) ||
698 (fJetType == 2 && t->Charge() != 0))
699 continue;
700
05077f28 701 if (t->Pt() < fJetConstituentMinPt)
4358e58a 702 continue;
703
704 fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
705 }
706 }
707
708 if (clusters && fJetType != 1) {
709 Double_t vert[3]={0};
710 if (fAODVertex)
711 ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
712
713 const Int_t Nclus = clusters->GetEntries();
714 for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
715 AliVCluster *c = static_cast<AliVCluster*>(clusters->At(iClus));
716 if (!c)
717 continue;
718
719 if (!c->IsEMCAL())
720 continue;
721
722 TLorentzVector nP;
723 c->GetMomentum(nP, vert);
724
05077f28 725 if (nP.Pt() < fJetConstituentMinPt)
4358e58a 726 continue;
727
728 fjw.AddInputVector(nP.Px(), nP.Py(), nP.Pz(), nP.P(), -iClus - 100);
729 }
730 }
731
732 // run jet finder
733 fjw.Run();
734
735 std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
736 AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
737
738 TLorentzVector jet;
739
740 Int_t njets = jets_incl.size();
741
742 if (njets > 0) {
05077f28 743 //std::vector<fastjet::PseudoJet> jets_incl_sorted = fastjet::sorted_by_pt(jets_incl);
4358e58a 744 for (Int_t i = 0; i < njets; i++) {
05077f28 745 if (jet.Pt() >= jets_incl[i].perp())
746 continue;
747 if (jets_incl[i].eta() < fJetMinEta || jets_incl[i].eta() > fJetMaxEta || jets_incl[i].phi() < fJetMinPhi || jets_incl[i].phi() > fJetMaxPhi)
748 continue;
749 jet.SetPxPyPzE(jets_incl[i].px(), jets_incl[i].py(), jets_incl[i].pz(), jets_incl[i].E());
4358e58a 750 }
751 }
752
753 return jet;
754}