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