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