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