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