]> 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 c36ea78b7f312f8ff648b08ca18f89927153fba1..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"
@@ -37,6 +38,7 @@ ClassImp(AliJetEmbeddingFromAODTask)
 AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() : 
   AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
   fFileList(0),
+  fRandomAccess(kFALSE),
   fAODTreeName(),
   fAODHeaderName(),
   fAODVertexName(),
@@ -49,31 +51,46 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
   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),
   fAODMCParticles(0),
-  fHistFileIDs(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"),
@@ -86,25 +103,39 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t
   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),  
   fAODMCParticles(0),
-  fHistFileIDs(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;
 }
 
 //________________________________________________________________________
@@ -126,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);
 }
@@ -147,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;
 }
 
@@ -172,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://");
-  }
+    }
+
+    AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
+    fCurrentAODFile = TFile::Open(fileName);
+
+    i++;
 
-  fCurrentAODFile = TFile::Open(fileName);
+  } while ((!fCurrentAODFile || fCurrentAODFile->IsZombie()) && i < fAttempts);
 
-  if (!fCurrentAODFile || fCurrentAODFile->IsZombie()) {
+  if (!fCurrentAODFile || fCurrentAODFile->IsZombie())
     return kFALSE;
-  }
 
   const TList *clist = fCurrentAODFile->GetStreamerInfoCache();
   if(clist) {
@@ -198,62 +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);
-
-      if (!fAODMCParticlesName.IsNull()) 
-       tree->SetBranchAddress(fAODMCParticlesName, &fAODMCParticles);
-      
-      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);
 }
 
 //________________________________________________________________________
@@ -262,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;
+      }
     }
   }
 
@@ -293,7 +369,9 @@ 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;
   }
 
@@ -303,6 +381,7 @@ void AliJetEmbeddingFromAODTask::Run()
       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) {
@@ -310,10 +389,18 @@ void AliJetEmbeddingFromAODTask::Run()
          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!");
       }
     }
   }
@@ -323,6 +410,8 @@ void AliJetEmbeddingFromAODTask::Run()
     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++) {
@@ -332,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) {
@@ -343,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;
          }
 
@@ -370,8 +464,36 @@ 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, track->GetLabel());
+       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!");
       }
     }
   }
@@ -391,6 +513,12 @@ void AliJetEmbeddingFromAODTask::Run()
        TLorentzVector vect;
        Double_t vert[3] = {0,0,0};
        clus->GetMomentum(vect,vert);
+
+       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());
       }
     }
@@ -412,7 +540,7 @@ void AliJetEmbeddingFromAODTask::Run()
 
     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;