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

// $Id$
//
// Jet modelling task.
//
// Author: S.Aiola, C.Loizides

#include "AliJetModelBaseTask.h"

#include <TClonesArray.h>
#include <TH1.h>
#include <TLorentzVector.h>
#include <TRandom3.h>
#include <TList.h>

#include "AliVEvent.h"
#include "AliAODCaloCluster.h"
#include "AliESDCaloCluster.h"
#include "AliVCluster.h"
#include "AliEMCALDigit.h"
#include "AliEMCALRecPoint.h"
#include "AliESDCaloCells.h"
#include "AliAODCaloCells.h"
#include "AliVCaloCells.h"
#include "AliPicoTrack.h"
#include "AliEMCALGeometry.h"
#include "AliLog.h"

ClassImp(AliJetModelBaseTask)

//________________________________________________________________________
AliJetModelBaseTask::AliJetModelBaseTask() : 
  AliAnalysisTaskSE("AliJetModelBaseTask"),
  fGeomName(),
  fTracksName(),
  fOutTracksName(),
  fCaloName(),
  fOutCaloName(),
  fCellsName(),
  fOutCellsName(),
  fSuffix(),
  fEtaMin(-1),
  fEtaMax(1),
  fPhiMin(0),
  fPhiMax(TMath::Pi() * 2),
  fPtMin(0),
  fPtMax(0),
  fCopyArray(kTRUE),
  fNClusters(0),
  fNCells(0),
  fNTracks(0),
  fMarkMC(99999),
  fPtSpectrum(0),
  fQAhistos(kFALSE),
  fIsInit(0),
  fGeom(0),
  fClusters(0),
  fOutClusters(0),
  fTracks(0),
  fOutTracks(0),
  fCaloCells(0),
  fOutCaloCells(0),
  fAddedCells(0),
  fEsdMode(kFALSE),
  fOutput(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) : 
  AliAnalysisTaskSE(name),
  fGeomName(""),
  fTracksName("PicoTracks"),
  fOutTracksName("PicoTracksEmbedded"),
  fCaloName("CaloClustersCorr"),
  fOutCaloName("CaloClustersCorrEmbedded"),
  fCellsName(""),
  fOutCellsName(""),
  fSuffix("Processed"),
  fEtaMin(-1),
  fEtaMax(1),
  fPhiMin(0),
  fPhiMax(TMath::Pi() * 2),
  fPtMin(50),
  fPtMax(60),
  fCopyArray(kTRUE),
  fNClusters(0),
  fNCells(0),
  fNTracks(1),
  fMarkMC(99999),
  fPtSpectrum(0),
  fQAhistos(drawqa),
  fIsInit(0),
  fGeom(0),
  fClusters(0),
  fOutClusters(0),
  fTracks(0),
  fOutTracks(0),
  fCaloCells(0),
  fOutCaloCells(0),
  fAddedCells(0),
  fEsdMode(kFALSE),
  fOutput(0)
{
  // Standard constructor.

  if (fQAhistos) {
    DefineOutput(1, TList::Class()); 
  }
}

//________________________________________________________________________
AliJetModelBaseTask::~AliJetModelBaseTask()
{
  // Destructor
}

//________________________________________________________________________
void AliJetModelBaseTask::UserCreateOutputObjects()
{
  // Create user output.
  if (!fQAhistos)
    return;

  OpenFile(1);
  fOutput = new TList();
  fOutput->SetOwner();

  PostData(1, fOutput);
}

//________________________________________________________________________
void AliJetModelBaseTask::UserExec(Option_t *) 
{
  // Execute per event.

  if (!fIsInit)
    fIsInit = ExecOnce();

  if (!fIsInit)
    return;

  if (fCopyArray) {
    if (fOutTracks)
      fOutTracks->Delete();
    if (fOutClusters)
      fOutClusters->Delete();
  }

  AliVCaloCells *tempCaloCells = 0;

  if (fCaloCells) {
    fAddedCells = 0;
    if (!fCopyArray) {
      tempCaloCells = fCaloCells;
      fCaloCells = static_cast<AliVCaloCells*>(tempCaloCells->Clone(Form("%s_old",fCaloCells->GetName())));
    }
  }

  Run();

  if (fCaloCells) {
    delete fCaloCells;
    fCaloCells = tempCaloCells;
  }
}

//________________________________________________________________________
Bool_t AliJetModelBaseTask::ExecOnce()
{
  // Init task.

  delete gRandom;
  gRandom = new TRandom3(0);

  fEsdMode = InputEvent()->InheritsFrom("AliESDEvent");

  if (fPtMax < fPtMin) {
    AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
    fPtMax = fPtMin;
  }

  if (fEtaMax < fEtaMin) {
    AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
    fEtaMax = fEtaMin;
  }

  if (fPhiMax < fPhiMin) {
    AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
    fPhiMax = fPhiMin;
  }

  if (!fCaloCells && !fCellsName.IsNull()) {
    fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCellsName));
    if (!fCaloCells) {
      AliError(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data()));
      return kFALSE;
    }
    
    if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
      AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data())); 
      fCaloCells = 0;
      return kFALSE;
    }

    if (!fOutCaloCells) {
      fOutCellsName = fCellsName;
      if (fCopyArray) {
	fOutCellsName += fSuffix;

	if (fEsdMode) 
	  fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName);
	else
	  fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName);

	if (InputEvent()->FindListObject(fOutCellsName)) {
	  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data())); 
	  return kFALSE;
	}
	else {
	  InputEvent()->AddObject(fOutCaloCells);
	}
      }
      else {
	fOutCaloCells = fCaloCells;
      }
    }
  }

  if (fNTracks > 0 && !fTracks && !fTracksName.IsNull()) {
    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
    if (!fTracks) {
      AliError(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
      return kFALSE;
    }
    
    if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
      AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data())); 
      fTracks = 0;
      return kFALSE;
    }

    if (!fOutTracks) {
      fOutTracksName = fTracksName;
      if (fCopyArray) {
	fOutTracksName += fSuffix;
	fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize());
	fOutTracks->SetName(fOutTracksName);
	if (InputEvent()->FindListObject(fOutTracksName)) {
	  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data())); 
	  return kFALSE;
	}
	else {
	  InputEvent()->AddObject(fOutTracks);
	}
      }
      else {
	fOutTracks = fTracks;
      }
    }
  }

  if (fNClusters > 0 && !fClusters && !fCaloName.IsNull()) {
    fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 
    if (!fClusters) {
      AliError(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
      return kFALSE;
    }

    if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
      AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data())); 
      fClusters = 0;
      return kFALSE;
    }

    if (!fOutClusters) {
      fOutCaloName = fCaloName;
      if (fCopyArray) {
	fOutCaloName += fSuffix;
	fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
	fOutClusters->SetName(fOutCaloName);
	if (InputEvent()->FindListObject(fOutCaloName)) {
	  AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); 
	  return kFALSE;
	}
	else {
	  InputEvent()->AddObject(fOutClusters);
	}
      }
      else {
	fOutClusters = fClusters;
      }
    }
  }

  if (!fCaloName.IsNull() || !fCellsName.IsNull()) {
    if (!fGeom) {
      if (fGeomName.Length() > 0) {
	fGeom = AliEMCALGeometry::GetInstance(fGeomName);
	if (!fGeom) {
	  AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data()));
	  return kFALSE;
	}
      } else {
	fGeom = AliEMCALGeometry::GetInstance();
	if (!fGeom) {
	  AliFatal("Could not get default geometry!");
	  return kFALSE;
	}
      }
    }
  
    const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
    const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
    const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
    const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
    
    if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta;
    if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta;
    if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta;
    if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta;
    
    if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi;
    if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi;
    if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi;
    if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi;
  }

  return kTRUE;
}

//________________________________________________________________________
Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
{
  if (fOutCaloCells->GetNumberOfCells() < n) {
    fOutCaloCells->DeleteContainer();
    fOutCaloCells->CreateContainer(n);
  }
  else {
    fOutCaloCells->SetNumberOfCells(n);
  }

  fAddedCells = 0;

  return n;
}

//________________________________________________________________________
void AliJetModelBaseTask::CopyCells()
{
  for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
    Short_t mclabel = 0;
    Double_t efrac = 0.;
    Double_t time = -1;
    Short_t cellNum = -1;
    Double_t amp = -1;

    fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
    fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
  }

  fAddedCells = fCaloCells->GetNumberOfCells();

  AliDebug(2, Form("%d cells from the PYTHIA event", fAddedCells));
}

//________________________________________________________________________
Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi)
{
  // Add a cell to the event.

  Int_t absId = 0;
  if (eta < -100 || phi < 0) {
    GetRandomCell(eta, phi, absId);
  }
  else {
    fGeom->EtaPhiFromIndex(absId, eta, phi);
  }

  if (absId == -1) {
    AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!"
		    " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])", 
		    eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
    return 0;
  } 

  if (e < 0) {
    Double_t pt = GetRandomPt();
    TLorentzVector nPart;
    nPart.SetPtEtaPhiM(pt, eta, phi, 0);
    e = nPart.E();
  }

  return AddCell(e, absId);
}

//________________________________________________________________________
Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time)
{
  // Add a cell to the event.

  Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, fMarkMC, 0);

  if (r) {
    fAddedCells++;
    return fAddedCells;
  }
  else {
    return 0;
  }
}


//________________________________________________________________________
AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi)
{
  // Add a cluster to the event.

  Int_t absId = 0;
  if (eta < -100 || phi < 0) {
    GetRandomCell(eta, phi, absId);
  }
  else {
    fGeom->EtaPhiFromIndex(absId, eta, phi);
  }

  if (absId == -1) {
    AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
		    " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])", 
		    eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
    return 0;
  } 

  if (e < 0) {
    Double_t pt = GetRandomPt();
    TLorentzVector nPart;
    nPart.SetPtEtaPhiM(pt, eta, phi, 0);
    e = nPart.E();
  }

  return AddCluster(e, absId);
}
      
//________________________________________________________________________
AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
{
  // Add a cluster to the event.

  const Int_t nClusters = fOutClusters->GetEntriesFast();

  TClonesArray digits("AliEMCALDigit", 1);

  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
  digit->SetId(absId);
  digit->SetIndexInList(0);
  digit->SetType(AliEMCALDigit::kHG);
  digit->SetAmplitude(e);
      
  AliEMCALRecPoint recPoint("");
  recPoint.AddDigit(*digit, e, kFALSE);
  recPoint.EvalGlobalPosition(0, &digits);

  TVector3 gpos;
  recPoint.GetGlobalPosition(gpos);
  Float_t g[3];
  gpos.GetXYZ(g);
      
  AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
  cluster->SetType(AliVCluster::kEMCALClusterv1);
  cluster->SetE(recPoint.GetEnergy());
  cluster->SetPosition(g);
  cluster->SetNCells(1);
  UShort_t shortAbsId = absId;
  cluster->SetCellsAbsId(&shortAbsId);
  Double32_t fract = 1;
  cluster->SetCellsAmplitudeFraction(&fract);
  cluster->SetID(nClusters);
  cluster->SetEmcCpvDistance(-1);

  //MC label
  if (fEsdMode) {
    AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster);
    TArrayI parents(1, &fMarkMC);
    esdClus->AddLabels(parents); 
  }
  else {
    AliAODCaloCluster *aodClus = static_cast<AliAODCaloCluster*>(cluster);
    aodClus->SetLabel(&fMarkMC, 1); 
  }
  
  return cluster;
}

//________________________________________________________________________
AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Bool_t ise)
{
  // Add a track to the event.

  const Int_t nTracks = fOutTracks->GetEntriesFast();
  
  if (pt < 0) 
    pt = GetRandomPt();
  if (eta < -100) 
    eta = GetRandomEta();
  if (phi < 0) 
    phi = GetRandomPhi();

  AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, 
								  eta, 
								  phi, 
								  1,
								  fMarkMC,    // MC flag!
								  type, 
								  etaemc, 
								  phiemc, 
								  ise);

  return track;
}

//________________________________________________________________________
void AliJetModelBaseTask::CopyClusters()
{
  // Copy all the clusters in the new collection

  const Int_t nClusters = fClusters->GetEntriesFast();
  Int_t nCopiedClusters = 0;
  
  if (fEsdMode) {
    for (Int_t i = 0; i < nClusters; ++i) {
      AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
      if (!esdcluster || !esdcluster->IsEMCAL())
	continue;
      new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
      nCopiedClusters++;
    }
  }
  else {
    for (Int_t i = 0; i < nClusters; ++i) {
      AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
      if (!aodcluster || !aodcluster->IsEMCAL())
	continue;
      new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
      nCopiedClusters++;
    }
  }
}

//________________________________________________________________________
void AliJetModelBaseTask::CopyTracks()
{
  const Int_t nTracks = fTracks->GetEntriesFast();
  Int_t nCopiedTracks = 0;
  for (Int_t i = 0; i < nTracks; ++i) {
    AliPicoTrack *track = static_cast<AliPicoTrack*>(fTracks->At(i));
    if (!track)
      continue;
    new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*track);
    nCopiedTracks++;
  }
}

//________________________________________________________________________
void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
{
  // Get random cell.

  Int_t repeats = 0;
  Double_t rndEta = eta;
  Double_t rndPhi = phi;
  do {
    if (eta < -100)
      rndEta = GetRandomEta();
    if (phi < 0)
      rndPhi = GetRandomPhi();
    fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);  
    repeats++;
  } while (absId == -1 && repeats < 100);
  
  if (!(absId > -1)) {
    AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
		    "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
  }
  else {
    eta = rndEta;
    phi = rndPhi;
  }
}

//________________________________________________________________________
Double_t AliJetModelBaseTask::GetRandomEta()
{
  // Get random eta.

  return gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
}

//________________________________________________________________________
Double_t AliJetModelBaseTask::GetRandomPhi()
{
  // Get random phi.

  return gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
}

//________________________________________________________________________
Double_t AliJetModelBaseTask::GetRandomPt()
{
  // Get random pt.

  if (fPtSpectrum)
    return fPtSpectrum->GetRandom();
  else
    return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
}

//________________________________________________________________________
void AliJetModelBaseTask::Run() 
{
  // Run.
}
Commit	Line	Data
762e8424	1	// $Id$
	2	//
	3	// Jet modelling task.
	4	//
297edd60	5	// Author: S.Aiola, C.Loizides
762e8424	6
ffe32451	7	#include "AliJetModelBaseTask.h"
ffe32451	8
762e8424	9	#include <TClonesArray.h>
65bb5510	10	#include <TH1.h>
762e8424	11	#include <TLorentzVector.h>
762e8424	12	#include <TRandom3.h>
b16bb001	13	#include <TList.h>
762e8424	14
b16bb001	15	#include "AliVEvent.h"
e44e8726	16	#include "AliAODCaloCluster.h"
b16bb001	17	#include "AliESDCaloCluster.h"
b16bb001	18	#include "AliVCluster.h"
762e8424	19	#include "AliEMCALDigit.h"
ffe32451	20	#include "AliEMCALRecPoint.h"
b16bb001	21	#include "AliESDCaloCells.h"
	22	#include "AliAODCaloCells.h"
	23	#include "AliVCaloCells.h"
ffe32451	24	#include "AliPicoTrack.h"
b16bb001	25	#include "AliEMCALGeometry.h"
b16bb001	26	#include "AliLog.h"
762e8424	27
	28	ClassImp(AliJetModelBaseTask)
	29
	30	//________________________________________________________________________
	31	AliJetModelBaseTask::AliJetModelBaseTask() :
	32	AliAnalysisTaskSE("AliJetModelBaseTask"),
	33	fGeomName(),
	34	fTracksName(),
	35	fOutTracksName(),
	36	fCaloName(),
	37	fOutCaloName(),
b16bb001	38	fCellsName(),
b16bb001	39	fOutCellsName(),
762e8424	40	fSuffix(),
	41	fEtaMin(-1),
	42	fEtaMax(1),
	43	fPhiMin(0),
	44	fPhiMax(TMath::Pi() * 2),
	45	fPtMin(0),
	46	fPtMax(0),
	47	fCopyArray(kTRUE),
	48	fNClusters(0),
b16bb001	49	fNCells(0),
762e8424	50	fNTracks(0),
7f76e479	51	fMarkMC(99999),
e44e8726	52	fPtSpectrum(0),
b16bb001	53	fQAhistos(kFALSE),
ffe32451	54	fIsInit(0),
762e8424	55	fGeom(0),
	56	fClusters(0),
	57	fOutClusters(0),
	58	fTracks(0),
b16bb001	59	fOutTracks(0),
	60	fCaloCells(0),
	61	fOutCaloCells(0),
	62	fAddedCells(0),
	63	fEsdMode(kFALSE),
	64	fOutput(0)
762e8424	65	{
	66	// Default constructor.
	67	}
	68
	69	//________________________________________________________________________
b16bb001	70	AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
762e8424	71	AliAnalysisTaskSE(name),
	72	fGeomName(""),
	73	fTracksName("PicoTracks"),
	74	fOutTracksName("PicoTracksEmbedded"),
	75	fCaloName("CaloClustersCorr"),
	76	fOutCaloName("CaloClustersCorrEmbedded"),
b16bb001	77	fCellsName(""),
b16bb001	78	fOutCellsName(""),
762e8424	79	fSuffix("Processed"),
	80	fEtaMin(-1),
	81	fEtaMax(1),
	82	fPhiMin(0),
	83	fPhiMax(TMath::Pi() * 2),
	84	fPtMin(50),
	85	fPtMax(60),
	86	fCopyArray(kTRUE),
	87	fNClusters(0),
b16bb001	88	fNCells(0),
762e8424	89	fNTracks(1),
7f76e479	90	fMarkMC(99999),
e44e8726	91	fPtSpectrum(0),
b16bb001	92	fQAhistos(drawqa),
ffe32451	93	fIsInit(0),
762e8424	94	fGeom(0),
	95	fClusters(0),
	96	fOutClusters(0),
	97	fTracks(0),
b16bb001	98	fOutTracks(0),
	99	fCaloCells(0),
	100	fOutCaloCells(0),
	101	fAddedCells(0),
	102	fEsdMode(kFALSE),
	103	fOutput(0)
762e8424	104	{
762e8424	105	// Standard constructor.
b16bb001	106
	107	if (fQAhistos) {
	108	DefineOutput(1, TList::Class());
	109	}
762e8424	110	}
	111
	112	//________________________________________________________________________
	113	AliJetModelBaseTask::~AliJetModelBaseTask()
	114	{
	115	// Destructor
	116	}
	117
b16bb001	118	//________________________________________________________________________
	119	void AliJetModelBaseTask::UserCreateOutputObjects()
	120	{
	121	// Create user output.
	122	if (!fQAhistos)
	123	return;
	124
	125	OpenFile(1);
	126	fOutput = new TList();
	127	fOutput->SetOwner();
	128
	129	PostData(1, fOutput);
	130	}
	131
762e8424	132	//________________________________________________________________________
ffe32451	133	void AliJetModelBaseTask::UserExec(Option_t *)
	134	{
	135	// Execute per event.
	136
d63c8c07	137	if (!fIsInit)
	138	fIsInit = ExecOnce();
	139
	140	if (!fIsInit)
	141	return;
	142
	143	if (fCopyArray) {
	144	if (fOutTracks)
	145	fOutTracks->Delete();
	146	if (fOutClusters)
	147	fOutClusters->Delete();
b4ffb4d2	148	}
d63c8c07	149
cd6431de	150	AliVCaloCells *tempCaloCells = 0;
cd6431de	151
b16bb001	152	if (fCaloCells) {
b16bb001	153	fAddedCells = 0;
cd6431de	154	if (!fCopyArray) {
	155	tempCaloCells = fCaloCells;
	156	fCaloCells = static_cast<AliVCaloCells*>(tempCaloCells->Clone(Form("%s_old",fCaloCells->GetName())));
	157	}
b16bb001	158	}
b16bb001	159
ffe32451	160	Run();
cd6431de	161
	162	if (fCaloCells) {
	163	delete fCaloCells;
	164	fCaloCells = tempCaloCells;
	165	}
ffe32451	166	}
ffe32451	167
b16bb001	168	//________________________________________________________________________
	169	Bool_t AliJetModelBaseTask::ExecOnce()
	170	{
	171	// Init task.
	172
	173	delete gRandom;
	174	gRandom = new TRandom3(0);
	175
	176	fEsdMode = InputEvent()->InheritsFrom("AliESDEvent");
	177
	178	if (fPtMax < fPtMin) {
	179	AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
	180	fPtMax = fPtMin;
	181	}
	182
	183	if (fEtaMax < fEtaMin) {
	184	AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
	185	fEtaMax = fEtaMin;
	186	}
	187
	188	if (fPhiMax < fPhiMin) {
	189	AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
	190	fPhiMax = fPhiMin;
	191	}
	192
	193	if (!fCaloCells && !fCellsName.IsNull()) {
	194	fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCellsName));
	195	if (!fCaloCells) {
	196	AliError(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data()));
	197	return kFALSE;
	198	}
	199
	200	if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
	201	AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data()));
	202	fCaloCells = 0;
	203	return kFALSE;
	204	}
	205
	206	if (!fOutCaloCells) {
	207	fOutCellsName = fCellsName;
	208	if (fCopyArray) {
	209	fOutCellsName += fSuffix;
	210
	211	if (fEsdMode)
	212	fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName);
	213	else
	214	fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName);
	215
	216	if (InputEvent()->FindListObject(fOutCellsName)) {
	217	AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data()));
	218	return kFALSE;
	219	}
	220	else {
	221	InputEvent()->AddObject(fOutCaloCells);
	222	}
	223	}
	224	else {
	225	fOutCaloCells = fCaloCells;
	226	}
	227	}
	228	}
	229
	230	if (fNTracks > 0 && !fTracks && !fTracksName.IsNull()) {
	231	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
232	if (!fTracks) {
233	AliError(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
234	return kFALSE;
235	}
236
237	if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
238	AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data()));
239	fTracks = 0;
240	return kFALSE;
241	}
242
243	if (!fOutTracks) {
244	fOutTracksName = fTracksName;
245	if (fCopyArray) {
246	fOutTracksName += fSuffix;
247	fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize());
248	fOutTracks->SetName(fOutTracksName);
249	if (InputEvent()->FindListObject(fOutTracksName)) {
250	AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data()));
251	return kFALSE;
252	}
253	else {
254	InputEvent()->AddObject(fOutTracks);
255	}
256	}
257	else {
258	fOutTracks = fTracks;
259	}
260	}
261	}
262
263	if (fNClusters > 0 && !fClusters && !fCaloName.IsNull()) {
264	fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
265
266	if (!fClusters) {
267	AliError(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
268	return kFALSE;
269	}
270
271	if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
272	AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data()));
273	fClusters = 0;
274	return kFALSE;
275	}
276
277	if (!fOutClusters) {
278	fOutCaloName = fCaloName;
279	if (fCopyArray) {
280	fOutCaloName += fSuffix;
281	fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
282	fOutClusters->SetName(fOutCaloName);
283	if (InputEvent()->FindListObject(fOutCaloName)) {
284	AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data()));
285	return kFALSE;
286	}
287	else {
288	InputEvent()->AddObject(fOutClusters);
289	}
290	}
291	else {
292	fOutClusters = fClusters;
293	}
294	}
295	}
296
cfc2ac24	297	if (!fCaloName.IsNull() \|\| !fCellsName.IsNull()) {
	298	if (!fGeom) {
	299	if (fGeomName.Length() > 0) {
	300	fGeom = AliEMCALGeometry::GetInstance(fGeomName);
7f76e479	301	if (!fGeom) {
	302	AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data()));
	303	return kFALSE;
	304	}
cfc2ac24	305	} else {
cfc2ac24	306	fGeom = AliEMCALGeometry::GetInstance();
7f76e479	307	if (!fGeom) {
	308	AliFatal("Could not get default geometry!");
	309	return kFALSE;
	310	}
cfc2ac24	311	}
b16bb001	312	}
b16bb001	313
cfc2ac24	314	const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
	315	const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
	316	const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
	317	const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
	318
	319	if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta;
	320	if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta;
	321	if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta;
	322	if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta;
	323
	324	if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi;
	325	if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi;
	326	if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi;
	327	if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi;
	328	}
b16bb001	329
	330	return kTRUE;
	331	}
	332
	333	//________________________________________________________________________
	334	Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
	335	{
	336	if (fOutCaloCells->GetNumberOfCells() < n) {
	337	fOutCaloCells->DeleteContainer();
	338	fOutCaloCells->CreateContainer(n);
	339	}
	340	else {
	341	fOutCaloCells->SetNumberOfCells(n);
	342	}
	343
	344	fAddedCells = 0;
	345
	346	return n;
	347	}
	348
	349	//________________________________________________________________________
	350	void AliJetModelBaseTask::CopyCells()
	351	{
	352	for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
7f76e479	353	Short_t mclabel = 0;
b16bb001	354	Double_t efrac = 0.;
	355	Double_t time = -1;
	356	Short_t cellNum = -1;
	357	Double_t amp = -1;
	358
	359	fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
	360	fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
	361	}
	362
	363	fAddedCells = fCaloCells->GetNumberOfCells();
cd6431de	364
cd6431de	365	AliDebug(2, Form("%d cells from the PYTHIA event", fAddedCells));
b16bb001	366	}
	367
	368	//________________________________________________________________________
	369	Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi)
	370	{
	371	// Add a cell to the event.
	372
	373	Int_t absId = 0;
	374	if (eta < -100 \|\| phi < 0) {
	375	GetRandomCell(eta, phi, absId);
	376	}
	377	else {
	378	fGeom->EtaPhiFromIndex(absId, eta, phi);
	379	}
	380
	381	if (absId == -1) {
	382	AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!"
	383	" Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
	384	eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
	385	return 0;
	386	}
	387
	388	if (e < 0) {
	389	Double_t pt = GetRandomPt();
	390	TLorentzVector nPart;
	391	nPart.SetPtEtaPhiM(pt, eta, phi, 0);
	392	e = nPart.E();
	393	}
	394
	395	return AddCell(e, absId);
	396	}
	397
	398	//________________________________________________________________________
	399	Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time)
	400	{
	401	// Add a cell to the event.
	402
f660c2d6	403	Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, fMarkMC, 0);
b16bb001	404
	405	if (r) {
	406	fAddedCells++;
	407	return fAddedCells;
	408	}
	409	else {
	410	return 0;
	411	}
	412	}
	413
	414
ffe32451	415	//________________________________________________________________________
	416	AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi)
	417	{
	418	// Add a cluster to the event.
	419
	420	Int_t absId = 0;
	421	if (eta < -100 \|\| phi < 0) {
	422	GetRandomCell(eta, phi, absId);
	423	}
	424	else {
	425	fGeom->EtaPhiFromIndex(absId, eta, phi);
	426	}
	427
	428	if (absId == -1) {
	429	AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
	430	" Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
	431	eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
	432	return 0;
	433	}
	434
	435	if (e < 0) {
	436	Double_t pt = GetRandomPt();
	437	TLorentzVector nPart;
	438	nPart.SetPtEtaPhiM(pt, eta, phi, 0);
	439	e = nPart.E();
	440	}
	441
	442	return AddCluster(e, absId);
	443	}
	444
	445	//________________________________________________________________________
	446	AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
	447	{
	448	// Add a cluster to the event.
	449
	450	const Int_t nClusters = fOutClusters->GetEntriesFast();
	451
	452	TClonesArray digits("AliEMCALDigit", 1);
	453
	454	AliEMCALDigit digit = static_cast<AliEMCALDigit>(digits.New(0));
	455	digit->SetId(absId);
	456	digit->SetIndexInList(0);
	457	digit->SetType(AliEMCALDigit::kHG);
	458	digit->SetAmplitude(e);
	459
	460	AliEMCALRecPoint recPoint("");
	461	recPoint.AddDigit(*digit, e, kFALSE);
	462	recPoint.EvalGlobalPosition(0, &digits);
	463
	464	TVector3 gpos;
	465	recPoint.GetGlobalPosition(gpos);
	466	Float_t g[3];
	467	gpos.GetXYZ(g);
	468
	469	AliVCluster cluster = static_cast<AliVCluster>(fOutClusters->New(nClusters));
	470	cluster->SetType(AliVCluster::kEMCALClusterv1);
	471	cluster->SetE(recPoint.GetEnergy());
	472	cluster->SetPosition(g);
	473	cluster->SetNCells(1);
	474	UShort_t shortAbsId = absId;
	475	cluster->SetCellsAbsId(&shortAbsId);
	476	Double32_t fract = 1;
	477	cluster->SetCellsAmplitudeFraction(&fract);
	478	cluster->SetID(nClusters);
479	cluster->SetEmcCpvDistance(-1);
ffe32451	480
f660c2d6	481	//MC label
	482	if (fEsdMode) {
	483	AliESDCaloCluster esdClus = static_cast<AliESDCaloCluster>(cluster);
	484	TArrayI parents(1, &fMarkMC);
	485	esdClus->AddLabels(parents);
	486	}
	487	else {
	488	AliAODCaloCluster aodClus = static_cast<AliAODCaloCluster>(cluster);
	489	aodClus->SetLabel(&fMarkMC, 1);
	490	}
	491
ffe32451	492	return cluster;
	493	}
	494
	495	//________________________________________________________________________
b16bb001	496	AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Bool_t ise)
ffe32451	497	{
	498	// Add a track to the event.
	499
	500	const Int_t nTracks = fOutTracks->GetEntriesFast();
	501
	502	if (pt < 0)
	503	pt = GetRandomPt();
	504	if (eta < -100)
	505	eta = GetRandomEta();
	506	if (phi < 0)
	507	phi = GetRandomPhi();
	508
	509	AliPicoTrack track = new ((fOutTracks)[nTracks]) AliPicoTrack(pt,
	510	eta,
	511	phi,
b16bb001	512	1,
f660c2d6	513	fMarkMC, // MC flag!
b16bb001	514	type,
	515	etaemc,
	516	phiemc,
	517	ise);
	518
ffe32451	519	return track;
	520	}
	521
	522	//________________________________________________________________________
	523	void AliJetModelBaseTask::CopyClusters()
	524	{
	525	// Copy all the clusters in the new collection
	526
ffe32451	527	const Int_t nClusters = fClusters->GetEntriesFast();
56cc432c	528	Int_t nCopiedClusters = 0;
ffe32451	529
b16bb001	530	if (fEsdMode) {
ffe32451	531	for (Int_t i = 0; i < nClusters; ++i) {
ffe32451	532	AliESDCaloCluster esdcluster = static_cast<AliESDCaloCluster>(fClusters->At(i));
56cc432c	533	if (!esdcluster \|\| !esdcluster->IsEMCAL())
ffe32451	534	continue;
56cc432c	535	new ((fOutClusters)[nCopiedClusters]) AliESDCaloCluster(esdcluster);
56cc432c	536	nCopiedClusters++;
ffe32451	537	}
	538	}
	539	else {
	540	for (Int_t i = 0; i < nClusters; ++i) {
	541	AliAODCaloCluster aodcluster = static_cast<AliAODCaloCluster>(fClusters->At(i));
56cc432c	542	if (!aodcluster \|\| !aodcluster->IsEMCAL())
ffe32451	543	continue;
56cc432c	544	new ((fOutClusters)[nCopiedClusters]) AliAODCaloCluster(aodcluster);
56cc432c	545	nCopiedClusters++;
ffe32451	546	}
	547	}
	548	}
	549
	550	//________________________________________________________________________
	551	void AliJetModelBaseTask::CopyTracks()
	552	{
	553	const Int_t nTracks = fTracks->GetEntriesFast();
56cc432c	554	Int_t nCopiedTracks = 0;
ffe32451	555	for (Int_t i = 0; i < nTracks; ++i) {
	556	AliPicoTrack track = static_cast<AliPicoTrack>(fTracks->At(i));
	557	if (!track)
	558	continue;
56cc432c	559	new ((fOutTracks)[nCopiedTracks]) AliPicoTrack(track);
56cc432c	560	nCopiedTracks++;
ffe32451	561	}
	562	}
	563
762e8424	564	//________________________________________________________________________
	565	void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
	566	{
acd859a2	567	// Get random cell.
acd859a2	568
762e8424	569	Int_t repeats = 0;
05debff6	570	Double_t rndEta = eta;
05debff6	571	Double_t rndPhi = phi;
762e8424	572	do {
05debff6	573	if (eta < -100)
	574	rndEta = GetRandomEta();
	575	if (phi < 0)
	576	rndPhi = GetRandomPhi();
	577	fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
762e8424	578	repeats++;
	579	} while (absId == -1 && repeats < 100);
	580
	581	if (!(absId > -1)) {
	582	AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
	583	"eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
	584	}
05debff6	585	else {
	586	eta = rndEta;
	587	phi = rndPhi;
	588	}
762e8424	589	}
	590
	591	//________________________________________________________________________
	592	Double_t AliJetModelBaseTask::GetRandomEta()
	593	{
acd859a2	594	// Get random eta.
acd859a2	595
762e8424	596	return gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
	597	}
	598
	599	//________________________________________________________________________
	600	Double_t AliJetModelBaseTask::GetRandomPhi()
	601	{
acd859a2	602	// Get random phi.
acd859a2	603
762e8424	604	return gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
	605	}
	606
	607	//________________________________________________________________________
	608	Double_t AliJetModelBaseTask::GetRandomPt()
	609	{
acd859a2	610	// Get random pt.
acd859a2	611
e44e8726	612	if (fPtSpectrum)
	613	return fPtSpectrum->GetRandom();
	614	else
	615	return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
762e8424	616	}
762e8424	617
762e8424	618	//________________________________________________________________________
	619	void AliJetModelBaseTask::Run()
	620	{
acd859a2	621	// Run.
762e8424	622	}