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