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