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