PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx

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