Setter and getter for centrality estimator (Chiara Bianchin)
[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
9//#include <iostream>
10
11#include <TFile.h>
12#include <TTree.h>
13#include <TClonesArray.h>
14#include <TObjArray.h>
15#include <TObjString.h>
16#include <TGrid.h>
17#include <TH2C.h>
18
19#include "AliVEvent.h"
20#include "AliAODTrack.h"
21#include "AliESDtrack.h"
22#include "AliPicoTrack.h"
23#include "AliVTrack.h"
24#include "AliVCluster.h"
25#include "AliVCaloCells.h"
26#include "AliCentrality.h"
27#include "AliVHeader.h"
28#include "AliVVertex.h"
29#include "AliAODHeader.h"
30#include "AliLog.h"
31#include "AliInputEventHandler.h"
32
33ClassImp(AliJetEmbeddingFromAODTask)
34
35//________________________________________________________________________
36AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
37 AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
38 fFileList(0),
39 fAODTreeName(),
40 fAODHeaderName(),
41 fAODVertexName(),
42 fAODTrackName(),
43 fAODClusName(),
44 fAODCellsName(),
45 fMinCentrality(0),
46 fMaxCentrality(10),
47 fTriggerMask(AliVEvent::kAny),
48 fZVertexCut(10),
49 fIncludeNoITS(kTRUE),
50 fTotalFiles(2200),
51 fEsdTreeMode(kFALSE),
52 fCurrentFileID(0),
53 fCurrentAODFileID(-1),
54 fCurrentAODFile(0),
55 fAODHeader(0),
56 fAODVertex(0),
57 fAODTracks(0),
58 fAODClusters(0),
59 fAODCaloCells(0),
60 fHistFileIDs(0)
61{
62 // Default constructor.
63 SetSuffix("AODEmbedding");
64 fMarkMC = kFALSE;
65 fAODfilterBits[0] = -1;
66 fAODfilterBits[1] = -1;
67}
68
69//________________________________________________________________________
70AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) :
71 AliJetModelBaseTask(name, drawqa),
72 fFileList(0),
73 fAODTreeName("aodTree"),
74 fAODHeaderName("header"),
75 fAODVertexName("vertices"),
76 fAODTrackName("tracks"),
77 fAODClusName("caloClusters"),
78 fAODCellsName("emcalCells"),
79 fMinCentrality(0),
80 fMaxCentrality(10),
81 fTriggerMask(AliVEvent::kAny),
82 fZVertexCut(10),
83 fIncludeNoITS(kTRUE),
84 fTotalFiles(2200),
85 fEsdTreeMode(kFALSE),
86 fCurrentFileID(0),
87 fCurrentAODFileID(-1),
88 fCurrentAODFile(0),
89 fAODHeader(0),
90 fAODVertex(0),
91 fAODTracks(0),
92 fAODClusters(0),
93 fAODCaloCells(0),
94 fHistFileIDs(0)
95{
96 // Standard constructor.
97 SetSuffix("AODEmbedding");
98 fMarkMC = kFALSE;
99 fAODfilterBits[0] = -1;
100 fAODfilterBits[1] = -1;
101}
102
103//________________________________________________________________________
104AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
105{
106 // Destructor
107
108 if (fCurrentAODFile) {
109 fCurrentAODFile->Close();
110 delete fCurrentAODFile;
111 }
112}
113
114//________________________________________________________________________
115void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
116{
117 if (!fQAhistos)
118 return;
119
120 AliJetModelBaseTask::UserCreateOutputObjects();
121
122 fHistFileIDs = new TH2C("fHistFileIDs", "fHistFileIDs", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
123 fHistFileIDs->GetXaxis()->SetTitle("File no. (PYTHIA)");
124 fHistFileIDs->GetYaxis()->SetTitle("File no. (Embedded AOD)");
125 fOutput->Add(fHistFileIDs);
126
127 PostData(1, fOutput);
128}
129
130//________________________________________________________________________
131Bool_t AliJetEmbeddingFromAODTask::UserNotify()
132{
133 TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
134 if (path.EndsWith("/"))
135 path.Remove(path.Length()-1);
136 path.Remove(path.Last('/'));
137 path.Remove(0, path.Last('/')+1);
138 if (!path.IsDec()) {
139 AliWarning(Form("Could not extract file number from path %s", path.Data()));
140 return kTRUE;
141 }
142 fCurrentFileID = path.Atoi();
143 fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles;
144 AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
145 return kTRUE;
146}
147
148//________________________________________________________________________
149Bool_t AliJetEmbeddingFromAODTask::ExecOnce()
150{
151 if (fAODTreeName.Contains("aod"))
152 fEsdTreeMode = kFALSE;
153 else
154 fEsdTreeMode = kTRUE;
155
156 return AliJetModelBaseTask::ExecOnce();
157}
158
159//________________________________________________________________________
160Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
161{
162 if (fCurrentAODFile) {
163 fCurrentAODFile->Close();
164 delete fCurrentAODFile;
165 fCurrentAODFile = 0;
166 }
167
168 if (fCurrentAODFileID >= fFileList->GetEntriesFast())
169 return kFALSE;
170
171 TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
172 TString fileName = objFileName->GetString();
173
174 if (fileName.BeginsWith("alien://") && !gGrid) {
175 AliInfo("Trying to connect to AliEn ...");
176 TGrid::Connect("alien://");
177 }
178
179 fCurrentAODFile = TFile::Open(fileName);
180
181 if (!fCurrentAODFile || fCurrentAODFile->IsZombie()) {
182 return kFALSE;
183 }
184
185 if (fQAhistos)
186 fHistFileIDs->Fill(fCurrentFileID, fCurrentAODFileID);
187
188 fCurrentAODFileID++;
189 return kTRUE;
190}
191
192//________________________________________________________________________
193Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
194{
195 static TTree *tree = 0;
196 static Int_t entry = 0;
197
198 Int_t attempts = 0;
199
200 do {
201
202 if (!fCurrentAODFile || !tree || entry >= tree->GetEntries()) {
203 if (!OpenNextFile())
204 return kFALSE;
205
206 tree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
207 if (!tree)
208 return kFALSE;
209
210 if (!fAODHeaderName.IsNull())
211 tree->SetBranchAddress(fAODHeaderName, &fAODHeader);
212
213 if (!fAODVertexName.IsNull())
214 tree->SetBranchAddress(fAODVertexName, &fAODVertex);
215
216 if (!fAODTrackName.IsNull())
217 tree->SetBranchAddress(fAODTrackName, &fAODTracks);
218
219 if (!fAODClusName.IsNull())
220 tree->SetBranchAddress(fAODClusName, &fAODClusters);
221
222 if (!fAODCellsName.IsNull())
223 tree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
224
225 entry = 0;
226 }
227
228 tree->GetEntry(entry);
229 entry++;
230 attempts++;
231
232 if (attempts == 1000)
233 AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
234
235 } while (!IsAODEventSelected() && tree);
236
237 return (tree!=0);
238}
239
240//________________________________________________________________________
241Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
242{
243 if (!fEsdTreeMode && fAODHeader) {
244 AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
245
246 UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
247
248 if ((offlineTrigger & fTriggerMask) == 0) {
249 AliDebug(2, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask));
250 return kFALSE;
251 }
252
253 AliCentrality *cent = aodHeader->GetCentralityP();
254 Float_t centVal = cent->GetCentralityPercentile("V0M");
255 if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
256 AliDebug(2, Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality));
257 return kFALSE;
258 }
259 }
260
261 if (fAODVertex) {
262 Double_t vert[3]={0};
263 ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
264 if (TMath::Abs(vert[2]) > fZVertexCut) {
265 AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f", vert[2], fZVertexCut));
266 return kFALSE;
267 }
268 }
269
270 return kTRUE;
271}
272
273//________________________________________________________________________
274void AliJetEmbeddingFromAODTask::Run()
275{
276 if (!GetNextEntry()) {
277 AliError("Unable to get AOD event to embed. Nothing will be embedded.");
278 return;
279 }
280
281 if (fTracks && fOutTracks) {
282
283 if (fCopyArray)
284 CopyTracks();
285
286 if (fAODTracks) {
287 AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
288 for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
289 AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
290 if (!track) {
291 AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
292 continue;
293 }
294
295 Int_t type = 0;
296 Bool_t isEmc = kFALSE;
297 if (!fEsdTreeMode) {
298 AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
299 if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
300 type = 0;
301 }
302 else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
303 if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
304 if (fIncludeNoITS)
305 type = 2;
306 else
307 continue;
308 }
309 else {
310 type = 1;
311 }
312 }
313 else { /*not a good track*/
314 continue;
315 }
316
317 if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
318 aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
319 aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
320 isEmc = kTRUE;
321 }
322 else { /*not AOD mode, let's see if it is PicoTrack*/
323 AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
324 if (ptrack) {
325 type = ptrack->GetTrackType();
326 isEmc = ptrack->IsEMCAL();
327 }
328 }
329
330 AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f", track->Pt(), track->Eta(), track->Phi()));
331 AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc);
332 }
333 }
334 }
335
336 if (fClusters && fOutClusters) {
337
338 if (fCopyArray)
339 CopyClusters();
340
341 if (fAODClusters) {
342 for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
343 AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
344 if (!clus) {
345 AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
346 continue;
347 }
348 TLorentzVector vect;
349 Double_t vert[3] = {0,0,0};
350 clus->GetMomentum(vect,vert);
351 AddCluster(clus->E(), vect.Eta(), vect.Phi());
352 }
353 }
354 }
355
356 if (fCaloCells && fOutCaloCells) {
357
358 Int_t totalCells = fCaloCells->GetNumberOfCells();
359 if (fAODCaloCells)
360 totalCells += fAODCaloCells->GetNumberOfCells();
361
362 SetNumberOfOutCells(totalCells);
363
364 CopyCells();
365
366 if (fAODCaloCells) {
367 for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
368 Short_t mclabel = -1;
369 Double_t efrac = 0.;
370 Double_t time = -1;
371 Short_t cellNum = -1;
372 Double_t amp = -1;
373
374 fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
375 AliDebug(3,Form("Adding cell with amplitude %f, absolute ID %d, time %f", amp, cellNum, time));
376 AddCell(amp, cellNum, time);
377 }
378 }
379 }
380}