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