Named ArrayI
[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
b16bb001 9#include <TFile.h>
10#include <TTree.h>
11#include <TClonesArray.h>
12#include <TObjArray.h>
13#include <TObjString.h>
14#include <TGrid.h>
15#include <TH2C.h>
cd6431de 16#include <TList.h>
17#include <TStreamerInfo.h>
6a20534a 18#include <TRandom.h>
b16bb001 19
20#include "AliVEvent.h"
21#include "AliAODTrack.h"
22#include "AliESDtrack.h"
23#include "AliPicoTrack.h"
24#include "AliVTrack.h"
25#include "AliVCluster.h"
26#include "AliVCaloCells.h"
787a3c4f 27#include "AliAODMCParticle.h"
b16bb001 28#include "AliCentrality.h"
29#include "AliVHeader.h"
30#include "AliVVertex.h"
31#include "AliAODHeader.h"
32#include "AliLog.h"
33#include "AliInputEventHandler.h"
34
35ClassImp(AliJetEmbeddingFromAODTask)
36
37//________________________________________________________________________
38AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
39 AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
40 fFileList(0),
6a20534a 41 fRandomAccess(kFALSE),
b16bb001 42 fAODTreeName(),
43 fAODHeaderName(),
44 fAODVertexName(),
45 fAODTrackName(),
46 fAODClusName(),
47 fAODCellsName(),
787a3c4f 48 fAODMCParticlesName(),
b16bb001 49 fMinCentrality(0),
50 fMaxCentrality(10),
51 fTriggerMask(AliVEvent::kAny),
52 fZVertexCut(10),
53 fIncludeNoITS(kTRUE),
5ce8ae64 54 fUseNegativeLabels(kTRUE),
55 fTrackEfficiency(1),
56 fIsMC(kFALSE),
6a20534a 57 fTotalFiles(2050),
58 fAttempts(5),
b16bb001 59 fEsdTreeMode(kFALSE),
60 fCurrentFileID(0),
6a20534a 61 fCurrentAODFileID(0),
b16bb001 62 fCurrentAODFile(0),
cd6431de 63 fPicoTrackVersion(0),
6a20534a 64 fCurrentAODTree(0),
b16bb001 65 fAODHeader(0),
66 fAODVertex(0),
67 fAODTracks(0),
68 fAODClusters(0),
69 fAODCaloCells(0),
787a3c4f 70 fAODMCParticles(0),
6a20534a 71 fCurrentAODEntry(-1),
72 fHistFileMatching(0),
73 fHistAODFileError(0),
74 fHistNotEmbedded(0)
b16bb001 75{
76 // Default constructor.
77 SetSuffix("AODEmbedding");
7f76e479 78 SetMarkMC(0);
b16bb001 79 fAODfilterBits[0] = -1;
80 fAODfilterBits[1] = -1;
81}
82
83//________________________________________________________________________
84AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) :
85 AliJetModelBaseTask(name, drawqa),
86 fFileList(0),
6a20534a 87 fRandomAccess(kFALSE),
b16bb001 88 fAODTreeName("aodTree"),
89 fAODHeaderName("header"),
90 fAODVertexName("vertices"),
91 fAODTrackName("tracks"),
787a3c4f 92 fAODClusName(""),
b16bb001 93 fAODCellsName("emcalCells"),
787a3c4f 94 fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
b16bb001 95 fMinCentrality(0),
96 fMaxCentrality(10),
97 fTriggerMask(AliVEvent::kAny),
98 fZVertexCut(10),
99 fIncludeNoITS(kTRUE),
5ce8ae64 100 fUseNegativeLabels(kTRUE),
101 fTrackEfficiency(1),
102 fIsMC(kFALSE),
6a20534a 103 fTotalFiles(2050),
104 fAttempts(5),
b16bb001 105 fEsdTreeMode(kFALSE),
106 fCurrentFileID(0),
6a20534a 107 fCurrentAODFileID(0),
b16bb001 108 fCurrentAODFile(0),
cd6431de 109 fPicoTrackVersion(0),
6a20534a 110 fCurrentAODTree(0),
b16bb001 111 fAODHeader(0),
112 fAODVertex(0),
113 fAODTracks(0),
114 fAODClusters(0),
787a3c4f 115 fAODCaloCells(0),
116 fAODMCParticles(0),
6a20534a 117 fCurrentAODEntry(-1),
118 fHistFileMatching(0),
119 fHistAODFileError(0),
120 fHistNotEmbedded(0)
b16bb001 121{
122 // Standard constructor.
123 SetSuffix("AODEmbedding");
7f76e479 124 SetMarkMC(0);
b16bb001 125 fAODfilterBits[0] = -1;
126 fAODfilterBits[1] = -1;
127}
128
129//________________________________________________________________________
130AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
131{
132 // Destructor
133
134 if (fCurrentAODFile) {
135 fCurrentAODFile->Close();
136 delete fCurrentAODFile;
137 }
138}
139
140//________________________________________________________________________
141void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
142{
143 if (!fQAhistos)
144 return;
145
146 AliJetModelBaseTask::UserCreateOutputObjects();
147
6a20534a 148 fHistFileMatching = new TH2C("fHistFileMatching", "fHistFileMatching", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
149 fHistFileMatching->GetXaxis()->SetTitle("File no. (PYTHIA)");
150 fHistFileMatching->GetYaxis()->SetTitle("File no. (Embedded AOD)");
151 fHistFileMatching->GetZaxis()->SetTitle("counts");
152 fOutput->Add(fHistFileMatching);
153
154 fHistAODFileError = new TH1C("fHistAODFileError", "fHistAODFileError", fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
155 fHistAODFileError->GetXaxis()->SetTitle("File no. (Embedded AOD)");
156 fHistAODFileError->GetYaxis()->SetTitle("counts");
157 fOutput->Add(fHistAODFileError);
158
159 fHistNotEmbedded = new TH1C("fHistNotEmbedded", "fHistNotEmbedded", fTotalFiles, -0.5, fTotalFiles - 0.5);
160 fHistNotEmbedded->GetXaxis()->SetTitle("File no. (PYTHIA)");
161 fHistNotEmbedded->GetYaxis()->SetTitle("counts");
162 fOutput->Add(fHistNotEmbedded);
b16bb001 163
164 PostData(1, fOutput);
165}
166
167//________________________________________________________________________
168Bool_t AliJetEmbeddingFromAODTask::UserNotify()
169{
170 TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
171 if (path.EndsWith("/"))
172 path.Remove(path.Length()-1);
173 path.Remove(path.Last('/'));
174 path.Remove(0, path.Last('/')+1);
175 if (!path.IsDec()) {
176 AliWarning(Form("Could not extract file number from path %s", path.Data()));
177 return kTRUE;
178 }
179 fCurrentFileID = path.Atoi();
6a20534a 180 if (!fRandomAccess) {
181 fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles-1;
182 AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
183 }
b16bb001 184 return kTRUE;
185}
186
187//________________________________________________________________________
188Bool_t AliJetEmbeddingFromAODTask::ExecOnce()
189{
190 if (fAODTreeName.Contains("aod"))
191 fEsdTreeMode = kFALSE;
192 else
193 fEsdTreeMode = kTRUE;
194
195 return AliJetModelBaseTask::ExecOnce();
196}
197
198//________________________________________________________________________
199Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
200{
201 if (fCurrentAODFile) {
202 fCurrentAODFile->Close();
203 delete fCurrentAODFile;
204 fCurrentAODFile = 0;
205 }
206
6a20534a 207 Int_t i = 0;
208 TString fileName;
b16bb001 209
6a20534a 210 do {
211 if (i>0) {
7030f36f 212 AliDebug(3,Form("Failed to open file %s...", fileName.Data()));
213 if (fHistAODFileError)
6a20534a 214 fHistAODFileError->Fill(fCurrentAODFileID);
215 }
216
217 fileName = GetNextFileName();
218
219 if (fileName.IsNull())
220 break;
221
222 if (fileName.BeginsWith("alien://") && !gGrid) {
b16bb001 223 AliInfo("Trying to connect to AliEn ...");
224 TGrid::Connect("alien://");
6a20534a 225 }
226
7030f36f 227 AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
6a20534a 228 fCurrentAODFile = TFile::Open(fileName);
b16bb001 229
7030f36f 230 i++;
231
6a20534a 232 } while ((!fCurrentAODFile || fCurrentAODFile->IsZombie()) && i < fAttempts);
b16bb001 233
6a20534a 234 if (!fCurrentAODFile || fCurrentAODFile->IsZombie())
b16bb001 235 return kFALSE;
b16bb001 236
cd6431de 237 const TList *clist = fCurrentAODFile->GetStreamerInfoCache();
238 if(clist) {
239 TStreamerInfo *cinfo = static_cast<TStreamerInfo*>(clist->FindObject("AliPicoTrack"));
240 if(cinfo)
241 fPicoTrackVersion = cinfo->GetClassVersion();
242 else
243 fPicoTrackVersion = 0;
244 }
245
6a20534a 246 fCurrentAODTree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
247 if (!fCurrentAODTree)
248 return kFALSE;
249
250 if (!fAODHeaderName.IsNull())
251 fCurrentAODTree->SetBranchAddress(fAODHeaderName, &fAODHeader);
252
253 if (!fAODVertexName.IsNull())
254 fCurrentAODTree->SetBranchAddress(fAODVertexName, &fAODVertex);
255
256 if (!fAODTrackName.IsNull())
257 fCurrentAODTree->SetBranchAddress(fAODTrackName, &fAODTracks);
258
259 if (!fAODClusName.IsNull())
260 fCurrentAODTree->SetBranchAddress(fAODClusName, &fAODClusters);
261
262 if (!fAODCellsName.IsNull())
263 fCurrentAODTree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
264
265 if (!fAODMCParticlesName.IsNull())
266 fCurrentAODTree->SetBranchAddress(fAODMCParticlesName, &fAODMCParticles);
267
268 if (fRandomAccess)
269 fCurrentAODEntry = TMath::Nint(gRandom->Rndm()*fCurrentAODTree->GetEntries());
270 else
271 fCurrentAODEntry = 0;
7030f36f 272
273 AliDebug(2,Form("Will start embedding from entry %d", fCurrentAODEntry));
6a20534a 274
7030f36f 275 if (fHistFileMatching)
6a20534a 276 fHistFileMatching->Fill(fCurrentFileID, fCurrentAODFileID-1);
b16bb001 277
b16bb001 278 return kTRUE;
279}
280
281//________________________________________________________________________
6a20534a 282TString AliJetEmbeddingFromAODTask::GetNextFileName()
b16bb001 283{
6a20534a 284 if (fRandomAccess)
285 fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*fFileList->GetEntriesFast());
286 else
287 fCurrentAODFileID++;
b16bb001 288
6a20534a 289 if (fCurrentAODFileID >= fFileList->GetEntriesFast())
290 return "";
291
292 TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
293 return objFileName->GetString();
294}
295
296//________________________________________________________________________
297Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
298{
b16bb001 299 Int_t attempts = 0;
300
301 do {
6a20534a 302 if (!fCurrentAODFile || !fCurrentAODTree || fCurrentAODEntry >= fCurrentAODTree->GetEntries()) {
b16bb001 303 if (!OpenNextFile())
304 return kFALSE;
b16bb001 305 }
306
6a20534a 307 fCurrentAODTree->GetEntry(fCurrentAODEntry);
308 fCurrentAODEntry++;
b16bb001 309 attempts++;
310
311 if (attempts == 1000)
312 AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
313
6a20534a 314 } while (!IsAODEventSelected() && fCurrentAODTree);
b16bb001 315
6a20534a 316 return (fCurrentAODTree!=0);
b16bb001 317}
318
319//________________________________________________________________________
320Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
321{
322 if (!fEsdTreeMode && fAODHeader) {
323 AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
324
6a20534a 325 if (fTriggerMask != AliVEvent::kAny) {
326 UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
327
328 if ((offlineTrigger & fTriggerMask) == 0) {
329 AliDebug(2,Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask));
330 return kFALSE;
331 }
b16bb001 332 }
333
6a20534a 334 if (fMinCentrality >= 0) {
335 AliCentrality *cent = aodHeader->GetCentralityP();
336 Float_t centVal = cent->GetCentralityPercentile("V0M");
337 if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
338 AliDebug(2,Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality));
339 return kFALSE;
340 }
b16bb001 341 }
342 }
343
344 if (fAODVertex) {
345 Double_t vert[3]={0};
346 ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
347 if (TMath::Abs(vert[2]) > fZVertexCut) {
348 AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f", vert[2], fZVertexCut));
349 return kFALSE;
350 }
351 }
352
353 return kTRUE;
354}
355
356//________________________________________________________________________
357void AliJetEmbeddingFromAODTask::Run()
358{
359 if (!GetNextEntry()) {
7030f36f 360 if (fHistNotEmbedded)
361 fHistNotEmbedded->Fill(fCurrentFileID);
6a20534a 362 AliError("Unable to get the AOD event to embed. Nothing will be embedded.");
b16bb001 363 return;
364 }
365
787a3c4f 366 if (fOutMCParticles) {
367
368 if (fCopyArray && fMCParticles)
369 CopyMCParticles();
370
371 if (fAODMCParticles) {
372 for (Int_t i = 0; i < fAODMCParticles->GetEntriesFast(); i++) {
373 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fAODMCParticles->At(i));
374 if (!part) {
375 AliError(Form("Could not find MC particle %d in branch %s of tree %s!", i, fAODMCParticlesName.Data(), fAODTreeName.Data()));
376 continue;
377 }
378
379 if (!part->IsPhysicalPrimary())
380 continue;
381
ca5c29fa 382 AliDebug(3, Form("Embedding MC particle with pT = %f, eta = %f, phi = %f", part->Pt(), part->Eta(), part->Phi()));
787a3c4f 383 AddMCParticle(part, i);
384 }
385 }
386 }
387
388 if (fOutTracks) {
b16bb001 389
787a3c4f 390 if (fCopyArray && fTracks)
b16bb001 391 CopyTracks();
392
393 if (fAODTracks) {
394 AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
395 for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
396 AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
397 if (!track) {
398 AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
399 continue;
400 }
401
402 Int_t type = 0;
403 Bool_t isEmc = kFALSE;
404 if (!fEsdTreeMode) {
405 AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
406 if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
407 type = 0;
408 }
409 else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
410 if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
411 if (fIncludeNoITS)
412 type = 2;
413 else
414 continue;
415 }
416 else {
417 type = 1;
418 }
419 }
420 else { /*not a good track*/
421 continue;
422 }
423
424 if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
425 aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
426 aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
427 isEmc = kTRUE;
428 }
cd6431de 429 else if (fPicoTrackVersion > 0) { /*not AOD mode, let's see if it is PicoTrack*/
b16bb001 430 AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
431 if (ptrack) {
cd6431de 432 if (fPicoTrackVersion >= 3)
433 type = ptrack->GetTrackType();
434 else
435 type = ptrack->GetLabel();
b16bb001 436 isEmc = ptrack->IsEMCAL();
437 }
438 }
439
5ce8ae64 440 if (fTrackEfficiency < 1) {
441 Double_t r = gRandom->Rndm();
442 if (fTrackEfficiency < r)
443 continue;
444 }
445
7030f36f 446 Int_t label = 0;
5ce8ae64 447 if (fIsMC) {
448 if (fUseNegativeLabels)
449 label = track->GetLabel();
450 else
451 label = TMath::Abs(track->GetLabel());
452
453 if (label == 0) {
454 AliWarning(Form("%s: Track %d with label==0", GetName(), i));
7030f36f 455 label = 99999;
5ce8ae64 456 }
7030f36f 457 }
458
ca5c29fa 459 AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f, label = %d", track->Pt(), track->Eta(), track->Phi(), label));
7030f36f 460 AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, label);
b16bb001 461 }
462 }
463 }
464
787a3c4f 465 if (fOutClusters) {
b16bb001 466
787a3c4f 467 if (fCopyArray && fClusters)
b16bb001 468 CopyClusters();
469
470 if (fAODClusters) {
471 for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
472 AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
473 if (!clus) {
474 AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
475 continue;
476 }
477 TLorentzVector vect;
478 Double_t vert[3] = {0,0,0};
479 clus->GetMomentum(vect,vert);
787a3c4f 480 AddCluster(clus->E(), vect.Eta(), vect.Phi(), clus->GetLabel());
b16bb001 481 }
482 }
483 }
484
787a3c4f 485 if (fOutCaloCells) {
486
487 Int_t totalCells = 0;
b16bb001 488
787a3c4f 489 if (fCaloCells)
490 totalCells += fCaloCells->GetNumberOfCells();
b16bb001 491 if (fAODCaloCells)
492 totalCells += fAODCaloCells->GetNumberOfCells();
493
494 SetNumberOfOutCells(totalCells);
495
787a3c4f 496 if (fCaloCells)
497 CopyCells();
b16bb001 498
499 if (fAODCaloCells) {
500 for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
60d77596 501 Int_t mclabel = 0;
b16bb001 502 Double_t efrac = 0.;
503 Double_t time = -1;
504 Short_t cellNum = -1;
505 Double_t amp = -1;
506
507 fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
508 AliDebug(3,Form("Adding cell with amplitude %f, absolute ID %d, time %f", amp, cellNum, time));
787a3c4f 509 AddCell(amp, cellNum, time, mclabel);
b16bb001 510 }
511 }
cd6431de 512
513 AliDebug(2,Form("Added cells = %d, total cells = %d", fAddedCells, totalCells));
b16bb001 514 }
515}