AliPhysics  c923f52 (c923f52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  if (fHistEmbeddingQA)
583  fHistEmbeddingQA->Fill("OK", 1);
584 
585  fAODFilePath->SetString(fCurrentAODFile->GetName());
586 
587  if (fOutMCParticles) {
588 
589  if (fCopyArray && fMCParticles)
590  CopyMCParticles();
591 
592  if (fAODMCParticles) {
593  AliDebug(3, Form("%d MC particles will be processed for embedding.", fAODMCParticles->GetEntriesFast()));
594  for (Int_t i = 0; i < fAODMCParticles->GetEntriesFast(); i++) {
595  AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fAODMCParticles->At(i));
596  if (!part) {
597  AliError(Form("Could not find MC particle %d in branch %s of tree %s!", i, fAODMCParticlesName.Data(), fAODTreeName.Data()));
598  continue;
599  }
600 
601  AliDebug(3, Form("Processing MC particle %d with pT = %f, eta = %f, phi = %f", i, part->Pt(), part->Eta(), part->Phi()));
602 
603  if (!part->IsPhysicalPrimary())
604  continue;
605 
606  if (part->Pt() < fPtMin || part->Pt() > fPtMax ||
607  part->Eta() < fEtaMin || part->Eta() > fEtaMax ||
608  part->Phi() < fPhiMin || part->Phi() > fPhiMax)
609  continue;
610 
611  AddMCParticle(part, i);
612  AliDebug(3, "Embedded!");
613  }
614  }
615  }
616 
617  if (fOutTracks) {
618 
619  if (fCopyArray && fTracks)
620  CopyTracks();
621 
622  AliDebug(3, Form("Start embedding with %d tracks.", fOutTracks->GetEntriesFast()));
623 
624  if (fAODTracks) {
625  AliDebug(3, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
626  for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
627  AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
628  if (!track) {
629  AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
630  continue;
631  }
632 
633  AliDebug(3, Form("Processing track %d with pT = %f, eta = %f, phi = %f, label = %d", i, track->Pt(), track->Eta(), track->Phi(), track->GetLabel()));
634 
635  Int_t type = 0;
636  Bool_t isEmc = kFALSE;
637  if (!fEsdTreeMode) {
638  AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
639  if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
640  type = 0;
641  }
642  else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
643  if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
644  if (fIncludeNoITS)
645  type = 2;
646  else {
647  AliDebug(3, "Track not embedded because ITS refit failed.");
648  continue;
649  }
650  }
651  else {
652  type = 1;
653  }
654  }
655  else { /*not a good track*/
656  AliDebug(3, "Track not embedded because not an hybrid track.");
657  continue;
658  }
660  Double_t frac = Double_t(aodtrack->GetTPCnclsS()) / Double_t(aodtrack->GetTPCncls());
662  continue;
663  }
664  if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
665  aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
666  aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
667  isEmc = kTRUE;
668  }
669  else if (fPicoTrackVersion > 0) { /*not AOD mode, let's see if it is PicoTrack*/
670  AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
671  if (ptrack) {
672  if (fPicoTrackVersion >= 3)
673  type = ptrack->GetTrackType();
674  else
675  type = ptrack->GetLabel();
676  isEmc = ptrack->IsEMCAL();
677 
678  if (!fIncludeNoITS && type==2) {
679  AliDebug(3, "Track not embedded because ITS refit failed.");
680  continue;
681  }
682  }
683  }
684 
685  if (track->Pt() < fPtMin || track->Pt() > fPtMax ||
686  track->Eta() < fEtaMin || track->Eta() > fEtaMax ||
687  track->Phi() < fPhiMin || track->Phi() > fPhiMax) {
688  AliDebug(3, "Track not embedded because out of limits.");
689  continue;
690  }
691 
692  if (fTrackEfficiency) {
693  Double_t r = gRandom->Rndm();
694  if (fTrackEfficiency->Eval(track->Pt()) < r) {
695  AliDebug(3, "Track not embedded because of artificial inefficiency.");
696  continue;
697  }
698  }
699 
700  Int_t label = 0;
701  if (fIsAODMC) {
702  if (fUseNegativeLabels)
703  label = track->GetLabel();
704  else
705  label = TMath::Abs(track->GetLabel());
706 
707  if (label == 0)
708  AliDebug(1,Form("%s: Track %d with label==0", GetName(), i));
709  }
710 
711  AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), isEmc, label, track->Charge());
712  AliDebug(3, "Track embedded!");
713  }
714  }
715  }
716 
717  if (fOutClusters) {
718 
719  if (fCopyArray && fClusters)
720  CopyClusters();
721 
722  if (fAODClusters) {
723  for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
724  AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
725  if (!clus) {
726  AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
727  continue;
728  }
729 
730  if (!clus->IsEMCAL())
731  continue;
732 
733  TLorentzVector vect;
734  Double_t vert[3] = {0,0,0};
735  clus->GetMomentum(vect,vert);
736 
737  if (vect.Pt() < fPtMin || vect.Pt() > fPtMax ||
738  vect.Eta() < fEtaMin || vect.Eta() > fEtaMax ||
739  vect.Phi() < fPhiMin || vect.Phi() > fPhiMax)
740  continue;
741 
742  Int_t label = 0;
743  if (fIsAODMC) {
744  label = clus->GetLabel();
745  if (label <= 0)
746  AliDebug(1,Form("%s: Clus %d with label<=0", GetName(), i));
747  }
748  AddCluster(clus);
749  }
750  }
751  }
752 
753  if (fOutCaloCells) {
754 
755  Double_t totalEnergy = 0;
756  Int_t totalCells = 0;
757 
758  if (fCaloCells)
759  totalCells += fCaloCells->GetNumberOfCells();
760  if (fAODCaloCells)
761  totalCells += fAODCaloCells->GetNumberOfCells();
762 
763  SetNumberOfOutCells(totalCells);
764 
765  if (fCaloCells)
766  CopyCells();
767 
768  if (fAODCaloCells) {
769  for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
770  Int_t mclabel = 0;
771  Double_t efrac = 0.;
772  Double_t time = -1;
773  Short_t cellNum = -1;
774  Double_t amp = -1;
775 
776  fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
777 
778  if (fIsAODMC) {
779  if (mclabel <= 0)
780  AliDebug(1,Form("%s: Cell %d with label<=0", GetName(), i));
781  }
782  else {
783  mclabel = 0;
784  }
785 
786  AliDebug(2,Form("Adding cell with amplitude %f, absolute ID %d, time %f, mc label %d", amp, cellNum, time, mclabel));
787  AddCell(amp, cellNum, time, mclabel);
788  totalEnergy += amp;
789  }
790  }
791  AliDebug(2,Form("Added cells = %d (energy = %f), total cells = %d", fAddedCells, totalEnergy, totalCells));
792  }
793 }
794 
795 //________________________________________________________________________
796 TLorentzVector AliJetEmbeddingFromAODTask::GetLeadingJet(TClonesArray *tracks, TClonesArray *clusters)
797 {
798  TString name("kt");
799  fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
800  if (fJetAlgo == 1) {
801  name = "antikt";
802  jalgo = fastjet::antikt_algorithm;
803  }
804 
805  // setup fj wrapper
806  AliFJWrapper fjw(name, name);
807  fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
808  fjw.SetGhostArea(1); // set a very large ghost area to speed up jet finding
809  fjw.SetR(fJetRadius);
810  fjw.SetAlgorithm(jalgo);
811  fjw.SetMaxRap(1);
812  fjw.Clear();
813 
814  if (tracks) {
815  const Int_t Ntracks = tracks->GetEntries();
816  for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
817  AliVParticle *t = static_cast<AliVParticle*>(tracks->At(iTracks));
818  if (!t)
819  continue;
820 
821  AliAODMCParticle *aodmcpart = dynamic_cast<AliAODMCParticle*>(t);
822  if (aodmcpart && !aodmcpart->IsPhysicalPrimary())
823  continue;
824 
825  if ((fJetType == 1 && t->Charge() == 0) ||
826  (fJetType == 2 && t->Charge() != 0))
827  continue;
828 
829  if (t->Pt() < fJetConstituentMinPt)
830  continue;
831 
832  fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
833  }
834  }
835 
836  if (clusters && fJetType != 1) {
837  Double_t vert[3]={0};
838  if (fAODVertex)
839  ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
840 
841  const Int_t Nclus = clusters->GetEntries();
842  for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
843  AliVCluster *c = static_cast<AliVCluster*>(clusters->At(iClus));
844  if (!c)
845  continue;
846 
847  if (!c->IsEMCAL())
848  continue;
849 
850  TLorentzVector nP;
851  c->GetMomentum(nP, vert);
852 
853  if (nP.Pt() < fJetConstituentMinPt)
854  continue;
855 
856  fjw.AddInputVector(nP.Px(), nP.Py(), nP.Pz(), nP.P(), -iClus - 100);
857  }
858  }
859 
860  // run jet finder
861  fjw.Run();
862 
863  std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
864  AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
865 
866  TLorentzVector jet;
867 
868  Int_t njets = jets_incl.size();
869 
870  if (njets > 0) {
871  //std::vector<fastjet::PseudoJet> jets_incl_sorted = fastjet::sorted_by_pt(jets_incl);
872  for (Int_t i = 0; i < njets; i++) {
873  if (jet.Pt() >= jets_incl[i].perp())
874  continue;
875  if (jets_incl[i].eta() < fJetMinEta || jets_incl[i].eta() > fJetMaxEta || jets_incl[i].phi() < fJetMinPhi || jets_incl[i].phi() > fJetMaxPhi)
876  continue;
877  jet.SetPxPyPzE(jets_incl[i].px(), jets_incl[i].py(), jets_incl[i].pz(), jets_incl[i].E());
878  }
879  }
880 
881  return jet;
882 }
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:100
TH2 * fHistFileMatching
Current AOD file path being embedded.
Bool_t IsEMCAL() const
Definition: AliPicoTrack.h:52
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:101
TClonesArray * fOutTracks
! output track collection
TClonesArray * fOutMCParticles
! output MC particles collection
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
Int_t GetLabel() const
Definition: AliPicoTrack.h:40
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:95
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:42
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:99
AliNamedString * fAODFilePath
Last entry in the AOD tree.
ClassImp(AliJetEmbeddingFromAODTask) AliJetEmbeddingFromAODTask
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
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:97
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