PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx

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