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