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