[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetEmbeddingFromAODTask.cxx

// $Id: AliJetEmbeddingFromAODTask.cxx $
//
// Jet embedding from AOD task.
//
// Author: S.Aiola, C.Loizides

#include "AliJetEmbeddingFromAODTask.h"

#include <TFile.h>
#include <TTree.h>
#include <TClonesArray.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TGrid.h>
#include <TH2C.h>
#include <TList.h>
#include <TStreamerInfo.h>
#include <TRandom.h>

#include "AliVEvent.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"
#include "AliPicoTrack.h"
#include "AliVTrack.h"
#include "AliVCluster.h"
#include "AliVCaloCells.h"
#include "AliAODMCParticle.h"
#include "AliCentrality.h"
#include "AliVHeader.h"
#include "AliVVertex.h"
#include "AliAODHeader.h"
#include "AliLog.h"
#include "AliInputEventHandler.h"

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),
  fUseNegativeLabels(kTRUE),
  fTrackEfficiency(1),
  fIsMC(kFALSE),
  fTotalFiles(2050),
  fAttempts(5),
  fEsdTreeMode(kFALSE),
  fCurrentFileID(0),
  fCurrentAODFileID(0),
  fCurrentAODFile(0),
  fPicoTrackVersion(0),
  fCurrentAODTree(0),
  fAODHeader(0),
  fAODVertex(0),
  fAODTracks(0),
  fAODClusters(0),
  fAODCaloCells(0),
  fAODMCParticles(0),
  fCurrentAODEntry(-1),
  fHistFileMatching(0),
  fHistAODFileError(0),
  fHistNotEmbedded(0)
{
  // Default constructor.
  SetSuffix("AODEmbedding");
  SetMarkMC(0);
  fAODfilterBits[0] = -1;
  fAODfilterBits[1] = -1;
}

//________________________________________________________________________
AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) : 
  AliJetModelBaseTask(name, drawqa),
  fFileList(0),
  fRandomAccess(kFALSE),
  fAODTreeName("aodTree"),
  fAODHeaderName("header"),
  fAODVertexName("vertices"),
  fAODTrackName("tracks"),
  fAODClusName(""),
  fAODCellsName("emcalCells"),
  fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
  fMinCentrality(0),
  fMaxCentrality(10),
  fTriggerMask(AliVEvent::kAny),
  fZVertexCut(10),
  fIncludeNoITS(kTRUE),
  fUseNegativeLabels(kTRUE),
  fTrackEfficiency(1),
  fIsMC(kFALSE),
  fTotalFiles(2050),
  fAttempts(5),
  fEsdTreeMode(kFALSE),
  fCurrentFileID(0),
  fCurrentAODFileID(0),
  fCurrentAODFile(0),
  fPicoTrackVersion(0),
  fCurrentAODTree(0),
  fAODHeader(0),
  fAODVertex(0),
  fAODTracks(0),
  fAODClusters(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;
}

//________________________________________________________________________
AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
{
  // Destructor

  if (fCurrentAODFile) {
    fCurrentAODFile->Close();
    delete fCurrentAODFile;
  }
}

//________________________________________________________________________
void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
{
  if (!fQAhistos)
    return;

  AliJetModelBaseTask::UserCreateOutputObjects();
  
  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);
}

//________________________________________________________________________
Bool_t AliJetEmbeddingFromAODTask::UserNotify()
{
  TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
  if (path.EndsWith("/"))
    path.Remove(path.Length()-1);
  path.Remove(path.Last('/'));
  path.Remove(0, path.Last('/')+1);
  if (!path.IsDec()) {
    AliWarning(Form("Could not extract file number from path %s", path.Data()));
    return kTRUE;
  }
  fCurrentFileID = path.Atoi();
  if (!fRandomAccess) {
    fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles-1;
    AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
  }
  return kTRUE;
}

//________________________________________________________________________
Bool_t AliJetEmbeddingFromAODTask::ExecOnce() 
{
  if (fAODTreeName.Contains("aod"))
    fEsdTreeMode = kFALSE;
  else
    fEsdTreeMode = kTRUE;

  return AliJetModelBaseTask::ExecOnce();
}

//________________________________________________________________________
Bool_t AliJetEmbeddingFromAODTask::OpenNextFile() 
{
  if (fCurrentAODFile) {
    fCurrentAODFile->Close();
    delete fCurrentAODFile;
    fCurrentAODFile = 0;
  }

  Int_t i = 0;
  TString fileName;

  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++;

  } while ((!fCurrentAODFile || fCurrentAODFile->IsZombie()) && i < fAttempts);

  if (!fCurrentAODFile || fCurrentAODFile->IsZombie())
    return kFALSE;

  const TList *clist = fCurrentAODFile->GetStreamerInfoCache();
  if(clist) {
    TStreamerInfo *cinfo = static_cast<TStreamerInfo*>(clist->FindObject("AliPicoTrack"));
    if(cinfo) 
      fPicoTrackVersion = cinfo->GetClassVersion();
    else
      fPicoTrackVersion = 0;
  }

  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);
  
  return kTRUE;
}

//________________________________________________________________________
TString AliJetEmbeddingFromAODTask::GetNextFileName()
{
    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 || !fCurrentAODTree || fCurrentAODEntry >= fCurrentAODTree->GetEntries()) {
      if (!OpenNextFile())
	return kFALSE;
    }
    
    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() && fCurrentAODTree);

  return (fCurrentAODTree!=0);
}

//________________________________________________________________________
Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
{
  if (!fEsdTreeMode && fAODHeader) {
    AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);

    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;
      }
    }
    
    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;
      }
    }
  }

  if (fAODVertex) {
    Double_t vert[3]={0};
    ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
    if (TMath::Abs(vert[2]) > fZVertexCut) {
      AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f", vert[2], fZVertexCut));
      return kFALSE;
    }
  }

  return kTRUE;
}

//________________________________________________________________________
void AliJetEmbeddingFromAODTask::Run() 
{
  if (!GetNextEntry()) {
    if (fHistNotEmbedded)
      fHistNotEmbedded->Fill(fCurrentFileID);
    AliError("Unable to get the AOD event to embed. Nothing will be embedded.");
    return;
  }

  if (fOutMCParticles) {

    if (fCopyArray && fMCParticles)
      CopyMCParticles();

    if (fAODMCParticles) {
      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;
	}
	
	if (!part->IsPhysicalPrimary()) 
	  continue;

	AliDebug(3, Form("Embedding MC particle with pT = %f, eta = %f, phi = %f", part->Pt(), part->Eta(), part->Phi()));
	AddMCParticle(part, i);
      }
    }
  }

  if (fOutTracks) {

    if (fCopyArray && fTracks)
      CopyTracks();

    if (fAODTracks) {
      AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
      for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
	AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
	if (!track) {
	  AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
	  continue;
	}
	
	Int_t type = 0;
	Bool_t isEmc = kFALSE;
	if (!fEsdTreeMode) {
	  AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
	  if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
	    type = 0;
	  }
	  else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
	    if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
	      if (fIncludeNoITS)
		type = 2;
	      else
		continue;
	    }
	    else {
	      type = 1;
	    }
	  }
	  else { /*not a good track*/
	    continue;
	  }

	  if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 && 
	      aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() && 
	      aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
	    isEmc = kTRUE;
	}
	else if (fPicoTrackVersion > 0) { /*not AOD mode, let's see if it is PicoTrack*/
	  AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
	  if (ptrack) {
	    if (fPicoTrackVersion >= 3)
	      type = ptrack->GetTrackType();
	    else
	      type = ptrack->GetLabel();
	    isEmc = ptrack->IsEMCAL();
	  }
	}
	
	if (fTrackEfficiency < 1) {
	  Double_t r = gRandom->Rndm();
	  if (fTrackEfficiency < r) 
	    continue;
	}
	
	Int_t label = 0;
	if (fIsMC) {
	  if (fUseNegativeLabels)
	    label = track->GetLabel();
	  else 
	    label = TMath::Abs(track->GetLabel());
	  
	  if (label == 0) {
	    AliWarning(Form("%s: Track %d with label==0", GetName(), i));
	    label = 99999;
	  }
	}

	AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f, label = %d", track->Pt(), track->Eta(), track->Phi(), label));
	AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, label);
      }
    }
  }

  if (fOutClusters) {

    if (fCopyArray && fClusters)
      CopyClusters();

    if (fAODClusters) {
      for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
	AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
	if (!clus) {
	  AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
	  continue;
	}
	TLorentzVector vect;
	Double_t vert[3] = {0,0,0};
	clus->GetMomentum(vect,vert);
	AddCluster(clus->E(), vect.Eta(), vect.Phi(), clus->GetLabel());
      }
    }
  }

  if (fOutCaloCells) {

    Int_t totalCells = 0;

    if (fCaloCells)
      totalCells += fCaloCells->GetNumberOfCells();
    if (fAODCaloCells)
      totalCells += fAODCaloCells->GetNumberOfCells();

    SetNumberOfOutCells(totalCells);

    if (fCaloCells)
      CopyCells();

    if (fAODCaloCells) {
      for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
	Int_t mclabel = 0;
	Double_t efrac = 0.;
	Double_t time = -1;
	Short_t cellNum = -1;
	Double_t amp = -1;
	
	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, mclabel);
      }
    }

    AliDebug(2,Form("Added cells = %d, total cells = %d", fAddedCells, totalCells));
  }
}
Commit	Line	Data
b16bb001	1	// $Id: AliJetEmbeddingFromAODTask.cxx $
	2	//
	3	// Jet embedding from AOD task.
	4	//
	5	// Author: S.Aiola, C.Loizides
	6
	7	#include "AliJetEmbeddingFromAODTask.h"
	8
b16bb001	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>
cd6431de	16	#include <TList.h>
cd6431de	17	#include <TStreamerInfo.h>
6a20534a	18	#include <TRandom.h>
b16bb001	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"
787a3c4f	27	#include "AliAODMCParticle.h"
b16bb001	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),
6a20534a	41	fRandomAccess(kFALSE),
b16bb001	42	fAODTreeName(),
	43	fAODHeaderName(),
	44	fAODVertexName(),
	45	fAODTrackName(),
	46	fAODClusName(),
	47	fAODCellsName(),
787a3c4f	48	fAODMCParticlesName(),
b16bb001	49	fMinCentrality(0),
	50	fMaxCentrality(10),
	51	fTriggerMask(AliVEvent::kAny),
	52	fZVertexCut(10),
	53	fIncludeNoITS(kTRUE),
5ce8ae64	54	fUseNegativeLabels(kTRUE),
	55	fTrackEfficiency(1),
	56	fIsMC(kFALSE),
6a20534a	57	fTotalFiles(2050),
6a20534a	58	fAttempts(5),
b16bb001	59	fEsdTreeMode(kFALSE),
b16bb001	60	fCurrentFileID(0),
6a20534a	61	fCurrentAODFileID(0),
b16bb001	62	fCurrentAODFile(0),
cd6431de	63	fPicoTrackVersion(0),
6a20534a	64	fCurrentAODTree(0),
b16bb001	65	fAODHeader(0),
	66	fAODVertex(0),
	67	fAODTracks(0),
	68	fAODClusters(0),
	69	fAODCaloCells(0),
787a3c4f	70	fAODMCParticles(0),
6a20534a	71	fCurrentAODEntry(-1),
	72	fHistFileMatching(0),
	73	fHistAODFileError(0),
	74	fHistNotEmbedded(0)
b16bb001	75	{
	76	// Default constructor.
	77	SetSuffix("AODEmbedding");
7f76e479	78	SetMarkMC(0);
b16bb001	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),
6a20534a	87	fRandomAccess(kFALSE),
b16bb001	88	fAODTreeName("aodTree"),
	89	fAODHeaderName("header"),
	90	fAODVertexName("vertices"),
	91	fAODTrackName("tracks"),
787a3c4f	92	fAODClusName(""),
b16bb001	93	fAODCellsName("emcalCells"),
787a3c4f	94	fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
b16bb001	95	fMinCentrality(0),
	96	fMaxCentrality(10),
	97	fTriggerMask(AliVEvent::kAny),
	98	fZVertexCut(10),
	99	fIncludeNoITS(kTRUE),
5ce8ae64	100	fUseNegativeLabels(kTRUE),
	101	fTrackEfficiency(1),
	102	fIsMC(kFALSE),
6a20534a	103	fTotalFiles(2050),
6a20534a	104	fAttempts(5),
b16bb001	105	fEsdTreeMode(kFALSE),
b16bb001	106	fCurrentFileID(0),
6a20534a	107	fCurrentAODFileID(0),
b16bb001	108	fCurrentAODFile(0),
cd6431de	109	fPicoTrackVersion(0),
6a20534a	110	fCurrentAODTree(0),
b16bb001	111	fAODHeader(0),
	112	fAODVertex(0),
	113	fAODTracks(0),
	114	fAODClusters(0),
787a3c4f	115	fAODCaloCells(0),
787a3c4f	116	fAODMCParticles(0),
6a20534a	117	fCurrentAODEntry(-1),
	118	fHistFileMatching(0),
	119	fHistAODFileError(0),
	120	fHistNotEmbedded(0)
b16bb001	121	{
	122	// Standard constructor.
	123	SetSuffix("AODEmbedding");
7f76e479	124	SetMarkMC(0);
b16bb001	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
6a20534a	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);
b16bb001	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();
6a20534a	180	if (!fRandomAccess) {
	181	fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles-1;
	182	AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
	183	}
b16bb001	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
6a20534a	207	Int_t i = 0;
6a20534a	208	TString fileName;
b16bb001	209
6a20534a	210	do {
6a20534a	211	if (i>0) {
7030f36f	212	AliDebug(3,Form("Failed to open file %s...", fileName.Data()));
7030f36f	213	if (fHistAODFileError)
6a20534a	214	fHistAODFileError->Fill(fCurrentAODFileID);
	215	}
	216
	217	fileName = GetNextFileName();
	218
	219	if (fileName.IsNull())
	220	break;
	221
	222	if (fileName.BeginsWith("alien://") && !gGrid) {
b16bb001	223	AliInfo("Trying to connect to AliEn ...");
b16bb001	224	TGrid::Connect("alien://");
6a20534a	225	}
6a20534a	226
7030f36f	227	AliDebug(3,Form("Trying to open file %s...", fileName.Data()));
6a20534a	228	fCurrentAODFile = TFile::Open(fileName);
b16bb001	229
7030f36f	230	i++;
7030f36f	231
6a20534a	232	} while ((!fCurrentAODFile \|\| fCurrentAODFile->IsZombie()) && i < fAttempts);
b16bb001	233
6a20534a	234	if (!fCurrentAODFile \|\| fCurrentAODFile->IsZombie())
b16bb001	235	return kFALSE;
b16bb001	236
cd6431de	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
6a20534a	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;
7030f36f	272
7030f36f	273	AliDebug(2,Form("Will start embedding from entry %d", fCurrentAODEntry));
6a20534a	274
7030f36f	275	if (fHistFileMatching)
6a20534a	276	fHistFileMatching->Fill(fCurrentFileID, fCurrentAODFileID-1);
b16bb001	277
b16bb001	278	return kTRUE;
	279	}
	280
	281	//________________________________________________________________________
6a20534a	282	TString AliJetEmbeddingFromAODTask::GetNextFileName()
b16bb001	283	{
6a20534a	284	if (fRandomAccess)
	285	fCurrentAODFileID = TMath::Nint(gRandom->Rndm()*fFileList->GetEntriesFast());
	286	else
	287	fCurrentAODFileID++;
b16bb001	288
6a20534a	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	{
b16bb001	299	Int_t attempts = 0;
	300
	301	do {
6a20534a	302	if (!fCurrentAODFile \|\| !fCurrentAODTree \|\| fCurrentAODEntry >= fCurrentAODTree->GetEntries()) {
b16bb001	303	if (!OpenNextFile())
b16bb001	304	return kFALSE;
b16bb001	305	}
b16bb001	306
6a20534a	307	fCurrentAODTree->GetEntry(fCurrentAODEntry);
6a20534a	308	fCurrentAODEntry++;
b16bb001	309	attempts++;
	310
	311	if (attempts == 1000)
	312	AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
	313
6a20534a	314	} while (!IsAODEventSelected() && fCurrentAODTree);
b16bb001	315
6a20534a	316	return (fCurrentAODTree!=0);
b16bb001	317	}
	318
	319	//________________________________________________________________________
	320	Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
	321	{
	322	if (!fEsdTreeMode && fAODHeader) {
	323	AliAODHeader aodHeader = static_cast<AliAODHeader>(fAODHeader);
	324
6a20534a	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	}
b16bb001	332	}
b16bb001	333
6a20534a	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	}
b16bb001	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()) {
7030f36f	360	if (fHistNotEmbedded)
7030f36f	361	fHistNotEmbedded->Fill(fCurrentFileID);
6a20534a	362	AliError("Unable to get the AOD event to embed. Nothing will be embedded.");
b16bb001	363	return;
	364	}
	365
787a3c4f	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
ca5c29fa	382	AliDebug(3, Form("Embedding MC particle with pT = %f, eta = %f, phi = %f", part->Pt(), part->Eta(), part->Phi()));
787a3c4f	383	AddMCParticle(part, i);
	384	}
	385	}
	386	}
	387
	388	if (fOutTracks) {
b16bb001	389
787a3c4f	390	if (fCopyArray && fTracks)
b16bb001	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	}
cd6431de	429	else if (fPicoTrackVersion > 0) { /not AOD mode, let's see if it is PicoTrack/
b16bb001	430	AliPicoTrack ptrack = dynamic_cast<AliPicoTrack>(track);
b16bb001	431	if (ptrack) {
cd6431de	432	if (fPicoTrackVersion >= 3)
	433	type = ptrack->GetTrackType();
	434	else
	435	type = ptrack->GetLabel();
b16bb001	436	isEmc = ptrack->IsEMCAL();
	437	}
	438	}
	439
5ce8ae64	440	if (fTrackEfficiency < 1) {
	441	Double_t r = gRandom->Rndm();
	442	if (fTrackEfficiency < r)
	443	continue;
	444	}
	445
7030f36f	446	Int_t label = 0;
5ce8ae64	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));
7030f36f	455	label = 99999;
5ce8ae64	456	}
7030f36f	457	}
7030f36f	458
ca5c29fa	459	AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f, label = %d", track->Pt(), track->Eta(), track->Phi(), label));
7030f36f	460	AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, label);
b16bb001	461	}
	462	}
	463	}
	464
787a3c4f	465	if (fOutClusters) {
b16bb001	466
787a3c4f	467	if (fCopyArray && fClusters)
b16bb001	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);
787a3c4f	480	AddCluster(clus->E(), vect.Eta(), vect.Phi(), clus->GetLabel());
b16bb001	481	}
	482	}
	483	}
	484
787a3c4f	485	if (fOutCaloCells) {
	486
	487	Int_t totalCells = 0;
b16bb001	488
787a3c4f	489	if (fCaloCells)
787a3c4f	490	totalCells += fCaloCells->GetNumberOfCells();
b16bb001	491	if (fAODCaloCells)
	492	totalCells += fAODCaloCells->GetNumberOfCells();
	493
	494	SetNumberOfOutCells(totalCells);
	495
787a3c4f	496	if (fCaloCells)
787a3c4f	497	CopyCells();
b16bb001	498
	499	if (fAODCaloCells) {
	500	for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
60d77596	501	Int_t mclabel = 0;
b16bb001	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));
787a3c4f	509	AddCell(amp, cellNum, time, mclabel);
b16bb001	510	}
b16bb001	511	}
cd6431de	512
cd6431de	513	AliDebug(2,Form("Added cells = %d, total cells = %d", fAddedCells, totalCells));
b16bb001	514	}
b16bb001	515	}