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 #include <TFile.h>
  10 #include <TTree.h>
  11 #include <TClonesArray.h>
  12 #include <TObjArray.h>
  13 #include <TObjString.h>
  14 #include <TGrid.h>
  15 #include <TH2C.h>
  16 #include <TList.h>
  17 #include <TStreamerInfo.h>
  18 #include <TRandom.h>
  19
  20 #include "AliVEvent.h"
  21 #include "AliAODTrack.h"
  22 #include "AliESDtrack.h"
  23 #include "AliPicoTrack.h"
  24 #include "AliVTrack.h"
  25 #include "AliVCluster.h"
  26 #include "AliVCaloCells.h"
  27 #include "AliAODMCParticle.h"
  28 #include "AliCentrality.h"
  29 #include "AliVHeader.h"
  30 #include "AliVVertex.h"
  31 #include "AliAODHeader.h"
  32 #include "AliLog.h"
  33 #include "AliInputEventHandler.h"
  34
  35 ClassImp(AliJetEmbeddingFromAODTask)
  36
  37 //________________________________________________________________________
  38 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
  39   AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
  40   fFileList(0),
  41   fRandomAccess(kFALSE),
  42   fAODTreeName(),
  43   fAODHeaderName(),
  44   fAODVertexName(),
  45   fAODTrackName(),
  46   fAODClusName(),
  47   fAODCellsName(),
  48   fAODMCParticlesName(),
  49   fMinCentrality(0),
  50   fMaxCentrality(10),
  51   fTriggerMask(AliVEvent::kAny),
  52   fZVertexCut(10),
  53   fIncludeNoITS(kTRUE),
  54   fUseNegativeLabels(kTRUE),
  55   fTrackEfficiency(1),
  56   fIsMC(kFALSE),
  57   fTotalFiles(2050),
  58   fAttempts(5),
  59   fEsdTreeMode(kFALSE),
  60   fCurrentFileID(0),
  61   fCurrentAODFileID(0),
  62   fCurrentAODFile(0),
  63   fPicoTrackVersion(0),
  64   fCurrentAODTree(0),
  65   fAODHeader(0),
  66   fAODVertex(0),
  67   fAODTracks(0),
  68   fAODClusters(0),
  69   fAODCaloCells(0),
  70   fAODMCParticles(0),
  71   fCurrentAODEntry(-1),
  72   fHistFileMatching(0),
  73   fHistAODFileError(0),
  74   fHistNotEmbedded(0)
  75 {
  76   // Default constructor.
  77   SetSuffix("AODEmbedding");
  78   SetMarkMC(0);
  79   fAODfilterBits[0] = -1;
  80   fAODfilterBits[1] = -1;
  81 }
  82
  83 //________________________________________________________________________
  84 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) :
  85   AliJetModelBaseTask(name, drawqa),
  86   fFileList(0),
  87   fRandomAccess(kFALSE),
  88   fAODTreeName("aodTree"),
  89   fAODHeaderName("header"),
  90   fAODVertexName("vertices"),
  91   fAODTrackName("tracks"),
  92   fAODClusName(""),
  93   fAODCellsName("emcalCells"),
  94   fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
  95   fMinCentrality(0),
  96   fMaxCentrality(10),
  97   fTriggerMask(AliVEvent::kAny),
  98   fZVertexCut(10),
  99   fIncludeNoITS(kTRUE),
 100   fUseNegativeLabels(kTRUE),
 101   fTrackEfficiency(1),
 102   fIsMC(kFALSE),
 103   fTotalFiles(2050),
 104   fAttempts(5),
 105   fEsdTreeMode(kFALSE),
 106   fCurrentFileID(0),
 107   fCurrentAODFileID(0),
 108   fCurrentAODFile(0),
 109   fPicoTrackVersion(0),
 110   fCurrentAODTree(0),
 111   fAODHeader(0),
 112   fAODVertex(0),
 113   fAODTracks(0),
 114   fAODClusters(0),
 115   fAODCaloCells(0),
 116   fAODMCParticles(0),
 117   fCurrentAODEntry(-1),
 118   fHistFileMatching(0),
 119   fHistAODFileError(0),
 120   fHistNotEmbedded(0)
 121 {
 122   // Standard constructor.
 123   SetSuffix("AODEmbedding");
 124   SetMarkMC(0);
 125   fAODfilterBits[0] = -1;
 126   fAODfilterBits[1] = -1;
 127 }
 128
 129 //________________________________________________________________________
 130 AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
 131 {
 132   // Destructor
 133
 134   if (fCurrentAODFile) {
 135     fCurrentAODFile->Close();
 136     delete fCurrentAODFile;
 137   }
 138 }
 139
 140 //________________________________________________________________________
 141 void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
 142 {
 143   if (!fQAhistos)
 144     return;
 145
 146   AliJetModelBaseTask::UserCreateOutputObjects();
 147
 148   fHistFileMatching = new TH2C("fHistFileMatching", "fHistFileMatching", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
 149   fHistFileMatching->GetXaxis()->SetTitle("File no. (PYTHIA)");
 150   fHistFileMatching->GetYaxis()->SetTitle("File no. (Embedded AOD)");
 151   fHistFileMatching->GetZaxis()->SetTitle("counts");
 152   fOutput->Add(fHistFileMatching);
 153
 154   fHistAODFileError = new TH1C("fHistAODFileError", "fHistAODFileError", fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
 155   fHistAODFileError->GetXaxis()->SetTitle("File no. (Embedded AOD)");
 156   fHistAODFileError->GetYaxis()->SetTitle("counts");
 157   fOutput->Add(fHistAODFileError);
 158
 159   fHistNotEmbedded = new TH1C("fHistNotEmbedded", "fHistNotEmbedded", fTotalFiles, -0.5, fTotalFiles - 0.5);
 160   fHistNotEmbedded->GetXaxis()->SetTitle("File no. (PYTHIA)");
 161   fHistNotEmbedded->GetYaxis()->SetTitle("counts");
 162   fOutput->Add(fHistNotEmbedded);
 163
 164   PostData(1, fOutput);
 165 }
 166
 167 //________________________________________________________________________
 168 Bool_t AliJetEmbeddingFromAODTask::UserNotify()
 169 {
 170   TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
 171   if (path.EndsWith("/"))
 172     path.Remove(path.Length()-1);
 173   path.Remove(path.Last('/'));
 174   path.Remove(0, path.Last('/')+1);
 175   if (!path.IsDec()) {
 176     AliWarning(Form("Could not extract file number from path %s", path.Data()));
 177     return kTRUE;
 178   }
 179   fCurrentFileID = path.Atoi();
 180   if (!fRandomAccess) {
 181     fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles-1;
 182     AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
 183   }
 184   return kTRUE;
 185 }
 186
 187 //________________________________________________________________________
 188 Bool_t AliJetEmbeddingFromAODTask::ExecOnce()
 189 {
 190   if (fAODTreeName.Contains("aod"))
 191     fEsdTreeMode = kFALSE;
 192   else
 193     fEsdTreeMode = kTRUE;
 194
 195   return AliJetModelBaseTask::ExecOnce();
 196 }
 197
 198 //________________________________________________________________________
 199 Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
 200 {
 201   if (fCurrentAODFile) {
 202     fCurrentAODFile->Close();
 203     delete fCurrentAODFile;
 204     fCurrentAODFile = 0;
 205   }
 206
 207   Int_t i = 0;
 208   TString fileName;
 209
 210   do {
 211     if (i>0) {
 212       AliDebug(3,Form("Failed to open file %s...", fileName.Data()));
 213       if (fHistAODFileError)
 214         fHistAODFileError->Fill(fCurrentAODFileID);
 215     }
 216
 217     fileName = GetNextFileName();
 218
 219     if (fileName.IsNull())
 220       break;
 221
 222     if (fileName.BeginsWith("alien://") && !gGrid) {
 223       AliInfo("Trying to connect to AliEn ...");
 224       TGrid::Connect("alien://");
 225     }
 226
 227     AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
 228     fCurrentAODFile = TFile::Open(fileName);
 229
 230     i++;
 231
 232   } while ((!fCurrentAODFile || fCurrentAODFile->IsZombie()) && i < fAttempts);
 233
 234   if (!fCurrentAODFile || fCurrentAODFile->IsZombie())
 235     return kFALSE;
 236
 237   const TList *clist = fCurrentAODFile->GetStreamerInfoCache();
 238   if(clist) {
 239     TStreamerInfo *cinfo = static_cast<TStreamerInfo*>(clist->FindObject("AliPicoTrack"));
 240     if(cinfo)
 241       fPicoTrackVersion = cinfo->GetClassVersion();
 242     else
 243       fPicoTrackVersion = 0;
 244   }
 245
 246   fCurrentAODTree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
 247   if (!fCurrentAODTree)
 248     return kFALSE;
 249
 250   if (!fAODHeaderName.IsNull())
 251     fCurrentAODTree->SetBranchAddress(fAODHeaderName, &fAODHeader);
 252
 253   if (!fAODVertexName.IsNull())
 254     fCurrentAODTree->SetBranchAddress(fAODVertexName, &fAODVertex);
 255
 256   if (!fAODTrackName.IsNull())
 257     fCurrentAODTree->SetBranchAddress(fAODTrackName, &fAODTracks);
 258
 259   if (!fAODClusName.IsNull())
 260     fCurrentAODTree->SetBranchAddress(fAODClusName, &fAODClusters);
 261
 262   if (!fAODCellsName.IsNull())
 263     fCurrentAODTree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
 264
 265   if (!fAODMCParticlesName.IsNull())
 266     fCurrentAODTree->SetBranchAddress(fAODMCParticlesName, &fAODMCParticles);
 267
 268   if (fRandomAccess)
 269     fCurrentAODEntry = TMath::Nint(gRandom->Rndm()*fCurrentAODTree->GetEntries());
 270   else
 271     fCurrentAODEntry = 0;
 272
 273   AliDebug(2,Form("Will start embedding from entry %d", fCurrentAODEntry));
 274
 275   if (fHistFileMatching)
 276     fHistFileMatching->Fill(fCurrentFileID, fCurrentAODFileID-1);
 277
 278   return kTRUE;
 279 }
 280
 281 //________________________________________________________________________
 282 TString AliJetEmbeddingFromAODTask::GetNextFileName()
 283 {
 284     if (fRandomAccess)
 285       fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*fFileList->GetEntriesFast());
 286     else
 287       fCurrentAODFileID++;
 288
 289     if (fCurrentAODFileID >= fFileList->GetEntriesFast())
 290       return "";
 291
 292     TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
 293     return objFileName->GetString();
 294 }
 295
 296 //________________________________________________________________________
 297 Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
 298 {
 299   Int_t attempts = 0;
 300
 301   do {
 302     if (!fCurrentAODFile || !fCurrentAODTree || fCurrentAODEntry >= fCurrentAODTree->GetEntries()) {
 303       if (!OpenNextFile())
 304         return kFALSE;
 305     }
 306
 307     fCurrentAODTree->GetEntry(fCurrentAODEntry);
 308     fCurrentAODEntry++;
 309     attempts++;
 310
 311     if (attempts == 1000)
 312       AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
 313
 314   } while (!IsAODEventSelected() && fCurrentAODTree);
 315
 316   return (fCurrentAODTree!=0);
 317 }
 318
 319 //________________________________________________________________________
 320 Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
 321 {
 322   if (!fEsdTreeMode && fAODHeader) {
 323     AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
 324
 325     if (fTriggerMask != AliVEvent::kAny) {
 326       UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
 327
 328       if ((offlineTrigger & fTriggerMask) == 0) {
 329         AliDebug(2,Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask));
 330         return kFALSE;
 331       }
 332     }
 333
 334     if (fMinCentrality >= 0) {
 335       AliCentrality *cent = aodHeader->GetCentralityP();
 336       Float_t centVal = cent->GetCentralityPercentile("V0M");
 337       if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
 338         AliDebug(2,Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality));
 339         return kFALSE;
 340       }
 341     }
 342   }
 343
 344   if (fAODVertex) {
 345     Double_t vert[3]={0};
 346     ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
 347     if (TMath::Abs(vert[2]) > fZVertexCut) {
 348       AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f", vert[2], fZVertexCut));
 349       return kFALSE;
 350     }
 351   }
 352
 353   return kTRUE;
 354 }
 355
 356 //________________________________________________________________________
 357 void AliJetEmbeddingFromAODTask::Run()
 358 {
 359   if (!GetNextEntry()) {
 360     if (fHistNotEmbedded)
 361       fHistNotEmbedded->Fill(fCurrentFileID);
 362     AliError("Unable to get the AOD event to embed. Nothing will be embedded.");
 363     return;
 364   }
 365
 366   if (fOutMCParticles) {
 367
 368     if (fCopyArray && fMCParticles)
 369       CopyMCParticles();
 370
 371     if (fAODMCParticles) {
 372       for (Int_t i = 0; i < fAODMCParticles->GetEntriesFast(); i++) {
 373         AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fAODMCParticles->At(i));
 374         if (!part) {
 375           AliError(Form("Could not find MC particle %d in branch %s of tree %s!", i, fAODMCParticlesName.Data(), fAODTreeName.Data()));
 376           continue;
 377         }
 378
 379         if (!part->IsPhysicalPrimary())
 380           continue;
 381
 382         AliDebug(3, Form("Embedding MC particle with pT = %f, eta = %f, phi = %f", part->Pt(), part->Eta(), part->Phi()));
 383         AddMCParticle(part, i);
 384       }
 385     }
 386   }
 387
 388   if (fOutTracks) {
 389
 390     if (fCopyArray && fTracks)
 391       CopyTracks();
 392
 393     if (fAODTracks) {
 394       AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
 395       for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
 396         AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
 397         if (!track) {
 398           AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
 399           continue;
 400         }
 401
 402         Int_t type = 0;
 403         Bool_t isEmc = kFALSE;
 404         if (!fEsdTreeMode) {
 405           AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
 406           if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
 407             type = 0;
 408           }
 409           else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
 410             if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
 411               if (fIncludeNoITS)
 412                 type = 2;
 413               else
 414                 continue;
 415             }
 416             else {
 417               type = 1;
 418             }
 419           }
 420           else { /*not a good track*/
 421             continue;
 422           }
 423
 424           if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 &&
 425               aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() &&
 426               aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
 427             isEmc = kTRUE;
 428         }
 429         else if (fPicoTrackVersion > 0) { /*not AOD mode, let's see if it is PicoTrack*/
 430           AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
 431           if (ptrack) {
 432             if (fPicoTrackVersion >= 3)
 433               type = ptrack->GetTrackType();
 434             else
 435               type = ptrack->GetLabel();
 436             isEmc = ptrack->IsEMCAL();
 437           }
 438         }
 439
 440         if (fTrackEfficiency < 1) {
 441           Double_t r = gRandom->Rndm();
 442           if (fTrackEfficiency < r)
 443             continue;
 444         }
 445
 446         Int_t label = 0;
 447         if (fIsMC) {
 448           if (fUseNegativeLabels)
 449             label = track->GetLabel();
 450           else
 451             label = TMath::Abs(track->GetLabel());
 452
 453           if (label == 0) {
 454             AliWarning(Form("%s: Track %d with label==0", GetName(), i));
 455             label = 99999;
 456           }
 457         }
 458
 459         AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f, label = %d", track->Pt(), track->Eta(), track->Phi(), label));
 460         AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, label);
 461       }
 462     }
 463   }
 464
 465   if (fOutClusters) {
 466
 467     if (fCopyArray && fClusters)
 468       CopyClusters();
 469
 470     if (fAODClusters) {
 471       for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
 472         AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
 473         if (!clus) {
 474           AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
 475           continue;
 476         }
 477         TLorentzVector vect;
 478         Double_t vert[3] = {0,0,0};
 479         clus->GetMomentum(vect,vert);
 480         AddCluster(clus->E(), vect.Eta(), vect.Phi(), clus->GetLabel());
 481       }
 482     }
 483   }
 484
 485   if (fOutCaloCells) {
 486
 487     Int_t totalCells = 0;
 488
 489     if (fCaloCells)
 490       totalCells += fCaloCells->GetNumberOfCells();
 491     if (fAODCaloCells)
 492       totalCells += fAODCaloCells->GetNumberOfCells();
 493
 494     SetNumberOfOutCells(totalCells);
 495
 496     if (fCaloCells)
 497       CopyCells();
 498
 499     if (fAODCaloCells) {
 500       for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
 501         Int_t mclabel = 0;
 502         Double_t efrac = 0.;
 503         Double_t time = -1;
 504         Short_t cellNum = -1;
 505         Double_t amp = -1;
 506
 507         fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
 508         AliDebug(3,Form("Adding cell with amplitude %f, absolute ID %d, time %f", amp, cellNum, time));
 509         AddCell(amp, cellNum, time, mclabel);
 510       }
 511     }
 512
 513     AliDebug(2,Form("Added cells = %d, total cells = %d", fAddedCells, totalCells));
 514   }
 515 }