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