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(kTRUE),
  72   fUseNegativeLabels(kTRUE),
  73   fTrackEfficiency(1),
  74   fIsAODMC(kFALSE),
  75   fTotalFiles(2050),
  76   fAttempts(5),
  77   fEsdTreeMode(kFALSE),
  78   fCurrentFileID(0),
  79   fCurrentAODFileID(0),
  80   fCurrentAODFile(0),
  81   fPicoTrackVersion(0),
  82   fCurrentAODTree(0),
  83   fAODHeader(0),
  84   fAODVertex(0),
  85   fAODTracks(0),
  86   fAODClusters(0),
  87   fAODCaloCells(0),
  88   fAODMCParticles(0),
  89   fCurrentAODEntry(-1),
  90   fHistFileMatching(0),
  91   fHistAODFileError(0),
  92   fHistNotEmbedded(0),
  93   fHistEmbeddingQA(0)
  94 {
  95   // Default constructor.
  96   SetSuffix("AODEmbedding");
  97   SetMarkMC(0);
  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(kTRUE),
 135   fUseNegativeLabels(kTRUE),
 136   fTrackEfficiency(1),
 137   fIsAODMC(kFALSE),
 138   fTotalFiles(2050),
 139   fAttempts(5),
 140   fEsdTreeMode(kFALSE),
 141   fCurrentFileID(0),
 142   fCurrentAODFileID(0),
 143   fCurrentAODFile(0),
 144   fPicoTrackVersion(0),
 145   fCurrentAODTree(0),
 146   fAODHeader(0),
 147   fAODVertex(0),
 148   fAODTracks(0),
 149   fAODClusters(0),
 150   fAODCaloCells(0),
 151   fAODMCParticles(0),
 152   fCurrentAODEntry(-1),
 153   fHistFileMatching(0),
 154   fHistAODFileError(0),
 155   fHistNotEmbedded(0),
 156   fHistEmbeddingQA(0)
 157 {
 158   // Standard constructor.
 159   SetSuffix("AODEmbedding");
 160   SetMarkMC(0);
 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
 537           if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
 538               aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
 539               aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
 540             isEmc = kTRUE;
 541         }
 542         else if (fPicoTrackVersion > 0) { /*not AOD mode, let's see if it is PicoTrack*/
 543           AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
 544           if (ptrack) {
 545             if (fPicoTrackVersion >= 3)
 546               type = ptrack->GetTrackType();
 547             else
 548               type = ptrack->GetLabel();
 549             isEmc = ptrack->IsEMCAL();
 550           }
 551         }
 552
 553         if (track->Pt() < fPtMin || track->Pt() > fPtMax ||
 554             track->Eta() < fEtaMin || track->Eta() > fEtaMax ||
 555             track->Phi() < fPhiMin || track->Phi() > fPhiMax) {
 556           AliDebug(3, "Track not embedded because out of limits.");
 557           continue;
 558         }
 559
 560         if (fTrackEfficiency < 1) {
 561           Double_t r = gRandom->Rndm();
 562           if (fTrackEfficiency < r) {
 563             AliDebug(3, "Track not embedded because of artificial inefiiciency.");
 564             continue;
 565           }
 566         }
 567
 568         Int_t label = 0;
 569         if (fIsAODMC) {
 570           if (fUseNegativeLabels)
 571             label = track->GetLabel();
 572           else
 573             label = TMath::Abs(track->GetLabel());
 574
 575           if (label == 0)
 576             AliWarning(Form("%s: Track %d with label==0", GetName(), i));
 577         }
 578
 579         AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, label);
 580         AliDebug(3, "Track embedded!");
 581       }
 582     }
 583   }
 584
 585   if (fOutClusters) {
 586
 587     if (fCopyArray && fClusters)
 588       CopyClusters();
 589
 590     if (fAODClusters) {
 591       for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
 592         AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
 593         if (!clus) {
 594           AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
 595           continue;
 596         }
 597         TLorentzVector vect;
 598         Double_t vert[3] = {0,0,0};
 599         clus->GetMomentum(vect,vert);
 600
 601         if (vect.Pt() < fPtMin || vect.Pt() > fPtMax ||
 602             vect.Eta() < fEtaMin || vect.Eta() > fEtaMax ||
 603             vect.Phi() < fPhiMin || vect.Phi() > fPhiMax)
 604           continue;
 605
 606         Int_t label = 0;
 607         if (fIsAODMC) {
 608           label = clus->GetLabel();
 609           if (label <= 0)
 610             AliWarning(Form("%s: Clus %d with label<=0", GetName(), i));
 611         }
 612
 613         AddCluster(clus);
 614       }
 615     }
 616   }
 617
 618   if (fOutCaloCells) {
 619
 620     Double_t totalEnergy = 0;
 621     Int_t totalCells = 0;
 622
 623     if (fCaloCells)
 624       totalCells += fCaloCells->GetNumberOfCells();
 625     if (fAODCaloCells)
 626       totalCells += fAODCaloCells->GetNumberOfCells();
 627
 628     SetNumberOfOutCells(totalCells);
 629
 630     if (fCaloCells)
 631       CopyCells();
 632
 633     if (fAODCaloCells) {
 634       for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
 635         Int_t mclabel = 0;
 636         Double_t efrac = 0.;
 637         Double_t time = -1;
 638         Short_t cellNum = -1;
 639         Double_t amp = -1;
 640
 641         fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
 642
 643         if (fIsAODMC) {
 644           if (mclabel <= 0)
 645             AliWarning(Form("%s: Cell %d with label<=0", GetName(), i));
 646         }
 647         else {
 648           mclabel = 0;
 649         }
 650
 651         AliDebug(2,Form("Adding cell with amplitude %f, absolute ID %d, time %f, mc label %d", amp, cellNum, time, mclabel));
 652         AddCell(amp, cellNum, time, mclabel);
 653         totalEnergy += amp;
 654       }
 655     }
 656
 657     AliDebug(2,Form("Added cells = %d (energy = %f), total cells = %d", fAddedCells, totalEnergy, totalCells));
 658   }
 659 }
 660
 661 //________________________________________________________________________
 662 TLorentzVector AliJetEmbeddingFromAODTask::GetLeadingJet(TClonesArray *tracks, TClonesArray *clusters)
 663 {
 664   TString name("kt");
 665   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
 666   if (fJetAlgo == 1) {
 667     name  = "antikt";
 668     jalgo = fastjet::antikt_algorithm;
 669   }
 670
 671   // setup fj wrapper
 672   AliFJWrapper fjw(name, name);
 673   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
 674   fjw.SetGhostArea(1);  // set a very large ghost area to speed up jet finding
 675   fjw.SetR(fJetRadius);
 676   fjw.SetAlgorithm(jalgo);
 677   fjw.SetMaxRap(1);
 678   fjw.Clear();
 679
 680   if (tracks) {
 681     const Int_t Ntracks = tracks->GetEntries();
 682     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 683       AliVParticle *t = static_cast<AliVParticle*>(tracks->At(iTracks));
 684       if (!t)
 685         continue;
 686
 687       AliAODMCParticle *aodmcpart = dynamic_cast<AliAODMCParticle*>(t);
 688       if (aodmcpart && !aodmcpart->IsPhysicalPrimary())
 689         continue;
 690
 691       if ((fJetType == 1 && t->Charge() == 0) ||
 692           (fJetType == 2 && t->Charge() != 0))
 693         continue;
 694
 695       if (t->Pt() < fJetConstituentMinPt)
 696         continue;
 697
 698       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
 699     }
 700   }
 701
 702   if (clusters && fJetType != 1) {
 703     Double_t vert[3]={0};
 704     if (fAODVertex)
 705       ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
 706
 707     const Int_t Nclus = clusters->GetEntries();
 708     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 709       AliVCluster *c = static_cast<AliVCluster*>(clusters->At(iClus));
 710       if (!c)
 711         continue;
 712
 713       if (!c->IsEMCAL())
 714         continue;
 715
 716       TLorentzVector nP;
 717       c->GetMomentum(nP, vert);
 718
 719       if (nP.Pt() < fJetConstituentMinPt)
 720         continue;
 721
 722       fjw.AddInputVector(nP.Px(), nP.Py(), nP.Pz(), nP.P(), -iClus - 100);
 723     }
 724   }
 725
 726   // run jet finder
 727   fjw.Run();
 728
 729   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
 730   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
 731
 732   TLorentzVector jet;
 733
 734   Int_t njets = jets_incl.size();
 735
 736   if (njets > 0) {
 737     //std::vector<fastjet::PseudoJet> jets_incl_sorted = fastjet::sorted_by_pt(jets_incl);
 738     for (Int_t i = 0; i < njets; i++) {
 739       if (jet.Pt() >= jets_incl[i].perp())
 740         continue;
 741       if (jets_incl[i].eta() < fJetMinEta || jets_incl[i].eta() > fJetMaxEta || jets_incl[i].phi() < fJetMinPhi || jets_incl[i].phi() > fJetMaxPhi)
 742         continue;
 743       jet.SetPxPyPzE(jets_incl[i].px(), jets_incl[i].py(), jets_incl[i].pz(), jets_incl[i].E());
 744     }
 745   }
 746
 747   return jet;
 748 }