]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx
Changes from Salvatore
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetEmbeddingFromAODTask.cxx
index 23c424df22be478a3745c7bfcdbb162267ca2e4f..bd21277edeb57098da47641b44ec0edd0b319c73 100644 (file)
@@ -15,6 +15,7 @@
 #include <TH2C.h>
 #include <TList.h>
 #include <TStreamerInfo.h>
+#include <TRandom.h>
 
 #include "AliVEvent.h"
 #include "AliAODTrack.h"
@@ -23,6 +24,7 @@
 #include "AliVTrack.h"
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
+#include "AliAODMCParticle.h"
 #include "AliCentrality.h"
 #include "AliVHeader.h"
 #include "AliVVertex.h"
@@ -36,70 +38,104 @@ ClassImp(AliJetEmbeddingFromAODTask)
 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() : 
   AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
   fFileList(0),
+  fRandomAccess(kFALSE),
   fAODTreeName(),
   fAODHeaderName(),
   fAODVertexName(),
   fAODTrackName(),
   fAODClusName(),
   fAODCellsName(),
+  fAODMCParticlesName(),
   fMinCentrality(0),
   fMaxCentrality(10),
   fTriggerMask(AliVEvent::kAny),
   fZVertexCut(10),
   fIncludeNoITS(kTRUE),
-  fTotalFiles(2200),
+  fUseNegativeLabels(kTRUE),
+  fTrackEfficiency(1),
+  fIsMC(kFALSE),
+  fTotalFiles(2050),
+  fAttempts(5),
   fEsdTreeMode(kFALSE),
   fCurrentFileID(0),
-  fCurrentAODFileID(-1),
+  fCurrentAODFileID(0),
   fCurrentAODFile(0),
   fPicoTrackVersion(0),
+  fCurrentAODTree(0),
   fAODHeader(0),
   fAODVertex(0),
   fAODTracks(0),
   fAODClusters(0),
   fAODCaloCells(0),
-  fHistFileIDs(0)
+  fAODMCParticles(0),
+  fCurrentAODEntry(-1),
+  fHistFileMatching(0),
+  fHistAODFileError(0),
+  fHistNotEmbedded(0)
 {
   // Default constructor.
   SetSuffix("AODEmbedding");
   SetMarkMC(0);
   fAODfilterBits[0] = -1;
   fAODfilterBits[1] = -1;
+  fEtaMin = -1;
+  fEtaMax = 1;
+  fPhiMin = -10;
+  fPhiMax = 10;
+  fPtMin = 0;
+  fPtMax = 1000;
 }
 
 //________________________________________________________________________
 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) : 
   AliJetModelBaseTask(name, drawqa),
   fFileList(0),
+  fRandomAccess(kFALSE),
   fAODTreeName("aodTree"),
   fAODHeaderName("header"),
   fAODVertexName("vertices"),
   fAODTrackName("tracks"),
-  fAODClusName("caloClusters"),
+  fAODClusName(""),
   fAODCellsName("emcalCells"),
+  fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
   fMinCentrality(0),
   fMaxCentrality(10),
   fTriggerMask(AliVEvent::kAny),
   fZVertexCut(10),
   fIncludeNoITS(kTRUE),
-  fTotalFiles(2200),
+  fUseNegativeLabels(kTRUE),
+  fTrackEfficiency(1),
+  fIsMC(kFALSE),
+  fTotalFiles(2050),
+  fAttempts(5),
   fEsdTreeMode(kFALSE),
   fCurrentFileID(0),
-  fCurrentAODFileID(-1),
+  fCurrentAODFileID(0),
   fCurrentAODFile(0),
   fPicoTrackVersion(0),
+  fCurrentAODTree(0),
   fAODHeader(0),
   fAODVertex(0),
   fAODTracks(0),
   fAODClusters(0),
-  fAODCaloCells(0),
-  fHistFileIDs(0)
+  fAODCaloCells(0),  
+  fAODMCParticles(0),
+  fCurrentAODEntry(-1),
+  fHistFileMatching(0),
+  fHistAODFileError(0),
+  fHistNotEmbedded(0)
 {
   // Standard constructor.
   SetSuffix("AODEmbedding");
   SetMarkMC(0);
   fAODfilterBits[0] = -1;
   fAODfilterBits[1] = -1;
+  fEtaMin = -1;
+  fEtaMax = 1;
+  fPhiMin = -10;
+  fPhiMax = 10;
+  fPtMin = 0;
+  fPtMax = 1000;
 }
 
 //________________________________________________________________________
@@ -121,10 +157,21 @@ void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
 
   AliJetModelBaseTask::UserCreateOutputObjects();
   
-  fHistFileIDs = new TH2C("fHistFileIDs", "fHistFileIDs", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
-  fHistFileIDs->GetXaxis()->SetTitle("File no. (PYTHIA)");
-  fHistFileIDs->GetYaxis()->SetTitle("File no. (Embedded AOD)");
-  fOutput->Add(fHistFileIDs);
+  fHistFileMatching = new TH2C("fHistFileMatching", "fHistFileMatching", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
+  fHistFileMatching->GetXaxis()->SetTitle("File no. (PYTHIA)");
+  fHistFileMatching->GetYaxis()->SetTitle("File no. (Embedded AOD)");
+  fHistFileMatching->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistFileMatching);
+
+  fHistAODFileError = new TH1C("fHistAODFileError", "fHistAODFileError", fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
+  fHistAODFileError->GetXaxis()->SetTitle("File no. (Embedded AOD)");
+  fHistAODFileError->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistAODFileError);
+
+  fHistNotEmbedded = new TH1C("fHistNotEmbedded", "fHistNotEmbedded", fTotalFiles, -0.5, fTotalFiles - 0.5);
+  fHistNotEmbedded->GetXaxis()->SetTitle("File no. (PYTHIA)");
+  fHistNotEmbedded->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistNotEmbedded);
 
   PostData(1, fOutput);
 }
@@ -142,8 +189,10 @@ Bool_t AliJetEmbeddingFromAODTask::UserNotify()
     return kTRUE;
   }
   fCurrentFileID = path.Atoi();
-  fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles;
-  AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
+  if (!fRandomAccess) {
+    fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles-1;
+    AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
+  }
   return kTRUE;
 }
 
@@ -167,22 +216,35 @@ Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
     fCurrentAODFile = 0;
   }
 
-  if (fCurrentAODFileID >= fFileList->GetEntriesFast())
-    return kFALSE;
+  Int_t i = 0;
+  TString fileName;
 
-  TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
-  TString fileName = objFileName->GetString();
-  
-  if (fileName.BeginsWith("alien://") && !gGrid) {
+  do {
+    if (i>0) {
+      AliDebug(3,Form("Failed to open file %s...", fileName.Data()));
+      if (fHistAODFileError)
+       fHistAODFileError->Fill(fCurrentAODFileID);
+    }
+
+    fileName = GetNextFileName();
+
+    if (fileName.IsNull())
+      break;
+    
+    if (fileName.BeginsWith("alien://") && !gGrid) {
       AliInfo("Trying to connect to AliEn ...");
       TGrid::Connect("alien://");
-  }
+    }
 
-  fCurrentAODFile = TFile::Open(fileName);
+    AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
+    fCurrentAODFile = TFile::Open(fileName);
 
-  if (!fCurrentAODFile || fCurrentAODFile->IsZombie()) {
+    i++;
+
+  } while ((!fCurrentAODFile || fCurrentAODFile->IsZombie()) && i < fAttempts);
+
+  if (!fCurrentAODFile || fCurrentAODFile->IsZombie())
     return kFALSE;
-  }
 
   const TList *clist = fCurrentAODFile->GetStreamerInfoCache();
   if(clist) {
@@ -193,59 +255,77 @@ Bool_t AliJetEmbeddingFromAODTask::OpenNextFile()
       fPicoTrackVersion = 0;
   }
 
-  if (fQAhistos)
-    fHistFileIDs->Fill(fCurrentFileID, fCurrentAODFileID);
+  fCurrentAODTree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
+  if (!fCurrentAODTree)
+    return kFALSE;
+
+  if (!fAODHeaderName.IsNull()) 
+    fCurrentAODTree->SetBranchAddress(fAODHeaderName, &fAODHeader);
+  
+  if (!fAODVertexName.IsNull()) 
+    fCurrentAODTree->SetBranchAddress(fAODVertexName, &fAODVertex);
+      
+  if (!fAODTrackName.IsNull()) 
+    fCurrentAODTree->SetBranchAddress(fAODTrackName, &fAODTracks);
+  
+  if (!fAODClusName.IsNull()) 
+    fCurrentAODTree->SetBranchAddress(fAODClusName, &fAODClusters);
+  
+  if (!fAODCellsName.IsNull()) 
+    fCurrentAODTree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
+  
+  if (!fAODMCParticlesName.IsNull()) 
+    fCurrentAODTree->SetBranchAddress(fAODMCParticlesName, &fAODMCParticles);
+  
+  if (fRandomAccess)
+    fCurrentAODEntry = TMath::Nint(gRandom->Rndm()*fCurrentAODTree->GetEntries());
+  else
+    fCurrentAODEntry = 0;
+
+  AliDebug(2,Form("Will start embedding from entry %d", fCurrentAODEntry));
+  
+  if (fHistFileMatching)
+    fHistFileMatching->Fill(fCurrentFileID, fCurrentAODFileID-1);
   
-  fCurrentAODFileID++;
   return kTRUE;
 }
 
 //________________________________________________________________________
-Bool_t AliJetEmbeddingFromAODTask::GetNextEntry() 
+TString AliJetEmbeddingFromAODTask::GetNextFileName()
 {
-  static TTree *tree = 0;
-  static Int_t entry = 0;
+    if (fRandomAccess) 
+      fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*fFileList->GetEntriesFast());
+    else
+      fCurrentAODFileID++;
 
+    if (fCurrentAODFileID >= fFileList->GetEntriesFast())
+      return "";
+    
+    TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
+    return objFileName->GetString();
+}
+
+//________________________________________________________________________
+Bool_t AliJetEmbeddingFromAODTask::GetNextEntry() 
+{
   Int_t attempts = 0;
 
   do {
-    
-    if (!fCurrentAODFile || !tree || entry >= tree->GetEntries()) {
+    if (!fCurrentAODFile || !fCurrentAODTree || fCurrentAODEntry >= fCurrentAODTree->GetEntries()) {
       if (!OpenNextFile())
        return kFALSE;
-      
-      tree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
-      if (!tree)
-       return kFALSE;
-
-      if (!fAODHeaderName.IsNull()) 
-       tree->SetBranchAddress(fAODHeaderName, &fAODHeader);
-
-      if (!fAODVertexName.IsNull()) 
-       tree->SetBranchAddress(fAODVertexName, &fAODVertex);
-      
-      if (!fAODTrackName.IsNull()) 
-       tree->SetBranchAddress(fAODTrackName, &fAODTracks);
-      
-      if (!fAODClusName.IsNull()) 
-       tree->SetBranchAddress(fAODClusName, &fAODClusters);
-      
-      if (!fAODCellsName.IsNull()) 
-       tree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
-      
-      entry = 0;
     }
     
-    tree->GetEntry(entry);
-    entry++;
+    fCurrentAODTree->GetEntry(fCurrentAODEntry);
+    fCurrentAODEntry++;
     attempts++;
 
     if (attempts == 1000) 
       AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
 
-  } while (!IsAODEventSelected() && tree);
+  } while (!IsAODEventSelected() && fCurrentAODTree);
 
-  return (tree!=0);
+  return (fCurrentAODTree!=0);
 }
 
 //________________________________________________________________________
@@ -254,18 +334,22 @@ Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
   if (!fEsdTreeMode && fAODHeader) {
     AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
 
-    UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
-  
-    if ((offlineTrigger & fTriggerMask) == 0) {
-      AliDebug(2, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask));
-      return kFALSE;
+    if (fTriggerMask != AliVEvent::kAny) {
+      UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
+      
+      if ((offlineTrigger & fTriggerMask) == 0) {
+       AliDebug(2,Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask));
+       return kFALSE;
+      }
     }
     
-    AliCentrality *cent = aodHeader->GetCentralityP();
-    Float_t centVal = cent->GetCentralityPercentile("V0M");
-    if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
-      AliDebug(2, Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality));
-      return kFALSE;
+    if (fMinCentrality >= 0) {
+      AliCentrality *cent = aodHeader->GetCentralityP();
+      Float_t centVal = cent->GetCentralityPercentile("V0M");
+      if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
+       AliDebug(2,Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality));
+       return kFALSE;
+      }
     }
   }
 
@@ -285,15 +369,49 @@ Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
 void AliJetEmbeddingFromAODTask::Run() 
 {
   if (!GetNextEntry()) {
-    AliError("Unable to get AOD event to embed. Nothing will be embedded.");
+    if (fHistNotEmbedded)
+      fHistNotEmbedded->Fill(fCurrentFileID);
+    AliError("Unable to get the AOD event to embed. Nothing will be embedded.");
     return;
   }
 
-  if (fTracks && fOutTracks) {
+  if (fOutMCParticles) {
+
+    if (fCopyArray && fMCParticles)
+      CopyMCParticles();
+
+    if (fAODMCParticles) {
+      AliDebug(2, Form("%d MC particles will be processed for embedding.", fAODMCParticles->GetEntriesFast()));
+      for (Int_t i = 0; i < fAODMCParticles->GetEntriesFast(); i++) {
+       AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fAODMCParticles->At(i));
+       if (!part) {
+         AliError(Form("Could not find MC particle %d in branch %s of tree %s!", i, fAODMCParticlesName.Data(), fAODTreeName.Data()));
+         continue;
+       }
+       
+       AliDebug(3, Form("Processing MC particle %d with pT = %f, eta = %f, phi = %f", i, part->Pt(), part->Eta(), part->Phi()));
+       
+       if (!part->IsPhysicalPrimary()) 
+         continue;
+
+       if (part->Pt() < fPtMin || part->Pt() > fPtMax ||
+           part->Eta() < fEtaMin || part->Eta() > fEtaMax ||
+           part->Phi() < fPhiMin || part->Phi() > fPhiMax)
+         continue;
+       
+       AddMCParticle(part, i);
+       AliDebug(3, "Embedded!");
+      }
+    }
+  }
+
+  if (fOutTracks) {
 
-    if (fCopyArray)
+    if (fCopyArray && fTracks)
       CopyTracks();
 
+    AliDebug(2, Form("Start embedding with %d tracks.", fOutTracks->GetEntriesFast()));
+
     if (fAODTracks) {
       AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
       for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
@@ -303,6 +421,8 @@ void AliJetEmbeddingFromAODTask::Run()
          continue;
        }
        
+       AliDebug(3, Form("Processing track %d with pT = %f, eta = %f, phi = %f, label = %d", i, track->Pt(), track->Eta(), track->Phi(), track->GetLabel()));
+       
        Int_t type = 0;
        Bool_t isEmc = kFALSE;
        if (!fEsdTreeMode) {
@@ -314,14 +434,17 @@ void AliJetEmbeddingFromAODTask::Run()
            if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
              if (fIncludeNoITS)
                type = 2;
-             else
+             else {
+               AliDebug(3, "Track not embedded because ITS refit failed.");
                continue;
            }
+           }
            else {
              type = 1;
            }
          }
          else { /*not a good track*/
+           AliDebug(3, "Track not embedded because not an hybrid track.");
            continue;
          }
 
@@ -341,15 +464,43 @@ void AliJetEmbeddingFromAODTask::Run()
          }
        }
        
-       AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f", track->Pt(), track->Eta(), track->Phi()));
-       AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc);
+       if (track->Pt() < fPtMin || track->Pt() > fPtMax ||
+           track->Eta() < fEtaMin || track->Eta() > fEtaMax ||
+           track->Phi() < fPhiMin || track->Phi() > fPhiMax) {
+         AliDebug(3, "Track not embedded because out of limits.");
+         continue;
+       }
+       
+       if (fTrackEfficiency < 1) {
+         Double_t r = gRandom->Rndm();
+         if (fTrackEfficiency < r) {
+           AliDebug(3, "Track not embedded because of artificial inefiiciency.");
+           continue;
+       }
+       }
+       
+       Int_t label = 0;
+       if (fIsMC) {
+         if (fUseNegativeLabels)
+           label = track->GetLabel();
+         else 
+           label = TMath::Abs(track->GetLabel());
+         
+         if (label == 0) {
+           AliDebug(2,Form("%s: Track %d with label==0", GetName(), i));
+           label = 99999;
+         }
+       }
+
+       AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, label);
+       AliDebug(3, "Track embedded!");
       }
     }
   }
 
-  if (fClusters && fOutClusters) {
+  if (fOutClusters) {
 
-    if (fCopyArray)
+    if (fCopyArray && fClusters)
       CopyClusters();
 
     if (fAODClusters) {
@@ -362,24 +513,34 @@ void AliJetEmbeddingFromAODTask::Run()
        TLorentzVector vect;
        Double_t vert[3] = {0,0,0};
        clus->GetMomentum(vect,vert);
-       AddCluster(clus->E(), vect.Eta(), vect.Phi());
+
+       if (vect.Pt() < fPtMin || vect.Pt() > fPtMax ||
+           vect.Eta() < fEtaMin || vect.Eta() > fEtaMax ||
+           vect.Phi() < fPhiMin || vect.Phi() > fPhiMax)
+         continue;
+
+       AddCluster(clus->E(), vect.Eta(), vect.Phi(), clus->GetLabel());
       }
     }
   }
 
-  if (fCaloCells && fOutCaloCells) {
+  if (fOutCaloCells) {
+
+    Int_t totalCells = 0;
 
-    Int_t totalCells = fCaloCells->GetNumberOfCells();
+    if (fCaloCells)
+      totalCells += fCaloCells->GetNumberOfCells();
     if (fAODCaloCells)
       totalCells += fAODCaloCells->GetNumberOfCells();
 
     SetNumberOfOutCells(totalCells);
 
-    CopyCells();
+    if (fCaloCells)
+      CopyCells();
 
     if (fAODCaloCells) {
       for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
-       Short_t mclabel = 0;
+       Int_t mclabel = 0;
        Double_t efrac = 0.;
        Double_t time = -1;
        Short_t cellNum = -1;
@@ -387,7 +548,7 @@ void AliJetEmbeddingFromAODTask::Run()
        
        fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
        AliDebug(3,Form("Adding cell with amplitude %f, absolute ID %d, time %f", amp, cellNum, time));
-       AddCell(amp, cellNum, time);
+       AddCell(amp, cellNum, time, mclabel);
       }
     }