[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliJetEmbeddingTask.cxx

// $Id$
//
// Jet embedding task.
//
// Author: Salvatore Aiola, Constantin Loizides

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

#include "AliAnalysisManager.h"
#include "AliVEvent.h"
#include "AliVCluster.h"
#include "AliEMCALDigit.h"
#include "AliEMCALRecPoint.h"
#include "AliPicoTrack.h"
#include "AliEMCALGeometry.h"
#include "AliLog.h"

#include "AliJetEmbeddingTask.h"

ClassImp(AliJetEmbeddingTask)

//________________________________________________________________________
AliJetEmbeddingTask::AliJetEmbeddingTask() : 
  AliAnalysisTaskSE("AliJetEmbeddingTask"),
  fGeomName(),
  fTracksName(),
  fOutTracksName(),
  fCaloName(),
  fOutCaloName(),
  fEtaMin(-1),
  fEtaMax(1),
  fPhiMin(0),
  fPhiMax(TMath::Pi() * 2),
  fPtMin(0),
  fPtMax(0),
  fCopyArray(kTRUE),
  fNEmbClusters(0),
  fNEmbTracks(0),
  fGeom(0),
  fClusters(0),
  fOutClusters(0),
  fTracks(0),
  fOutTracks(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliJetEmbeddingTask::AliJetEmbeddingTask(const char *name) : 
  AliAnalysisTaskSE(name),
  fGeomName(""),
  fTracksName("PicoTracks"),
  fOutTracksName("PicoTracksEmbedded"),
  fCaloName("CaloClustersCorr"),
  fOutCaloName("CaloClustersCorrEmbedded"),
  fEtaMin(-1),
  fEtaMax(1),
  fPhiMin(0),
  fPhiMax(TMath::Pi() * 2),
  fPtMin(50),
  fPtMax(60),
  fCopyArray(kTRUE),
  fNEmbClusters(0),
  fNEmbTracks(1),
  fGeom(0),
  fClusters(0),
  fOutClusters(0),
  fTracks(0),
  fOutTracks(0)
{
  // Standard constructor.

}

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

//________________________________________________________________________
void AliJetEmbeddingTask::Init()
{
  // Init task.

  if (fNEmbTracks > 0) {
    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));

    if (strcmp(fTracks->GetClass()->GetName(), "AliPicoTrack")) {
      AliError("Can only embed PicoTracks!");
      return;
    }
 
    if (!fTracks) {
      AliError(Form("Couldn't retrieve tracks with name %s!", fTracksName.Data()));
      return;
    }

    if (!fOutTracks) {
      fOutTracksName = fTracksName;
      if (fCopyArray) {
	fOutTracksName += "Embedded";
	fOutTracks = new TClonesArray(*fTracks);
	fOutTracks->SetName(fOutTracksName);
      }
      else {
	fOutTracks = fTracks;
      }
    }

    if (fCopyArray) {
      if (!(InputEvent()->FindListObject(fOutTracksName)))
	InputEvent()->AddObject(fOutTracks);
      fOutTracks->Clear();
    }
  }

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

    if (!fOutClusters) {
      fOutCaloName = fCaloName;
      if (fCopyArray) {
	fOutCaloName += "Embedded";
	fOutClusters = new TClonesArray(*fClusters);
	fOutClusters->SetName(fOutCaloName);
      }
      else {
	fOutClusters = fClusters;
      }
    }

    if (fCopyArray) {
      if (!(InputEvent()->FindListObject(fOutCaloName)))
	InputEvent()->AddObject(fOutClusters);
      fOutClusters->Clear();
    }
  }

  if (!fGeom) {
    if (fGeomName.Length()>0) {
      fGeom = AliEMCALGeometry::GetInstance(fGeomName);
      if (!fGeom)
        AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
    } else {
      fGeom = AliEMCALGeometry::GetInstance();
      if (!fGeom) 
        AliError("Could not get default geometry!");
    }
  }
}

//________________________________________________________________________
void AliJetEmbeddingTask::Embed() 
{
  // Embed particles.

  if (fNEmbClusters > 0 && fOutClusters) {
    const Int_t nClusters = fOutClusters->GetEntriesFast();
    TClonesArray digits("AliEMCALDigit", 1);
    for (Int_t i = 0; i < fNEmbClusters; ++i) {
      Double_t pt  = gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
      Double_t eta = gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
      Double_t phi = gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
    
      TLorentzVector nPart;
      nPart.SetPtEtaPhiM(pt, eta, phi, 0);
      Double_t e = nPart.E();
      
      Int_t absId = 0;
      fGeom->GetAbsCellIdFromEtaPhi(eta, phi, absId);

      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));
	continue;
      }
      
      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 + i));
      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 + i);
      cluster->SetEmcCpvDistance(-1);
      cluster->SetChi2(100); // MC flag!
    }
  }
 
  if (fNEmbTracks > 0 && fOutTracks) {
    Int_t nTracks = fOutTracks->GetEntriesFast();
    for (Int_t i = 0; i < fNEmbTracks; ++i) {
      Double_t pt  = gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
      Double_t eta = gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
      Double_t phi = gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
      
      new ((*fOutTracks)[nTracks + i]) AliPicoTrack(pt, 
						    eta, 
						    phi, 
						    1, 
						    100,    // MC flag!      
						    0, 
						    0, 
						    kFALSE);
    }
  }
}

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

  Init();

  Embed();
}
Commit	Line	Data
44629b14	1	// $Id$
	2	//
	3	// Jet embedding task.
	4	//
	5	// Author: Salvatore Aiola, Constantin Loizides
	6
	7	#include <TClonesArray.h>
	8	#include <TLorentzVector.h>
	9	#include <TRandom3.h>
	10
	11	#include "AliAnalysisManager.h"
	12	#include "AliVEvent.h"
	13	#include "AliVCluster.h"
	14	#include "AliEMCALDigit.h"
	15	#include "AliEMCALRecPoint.h"
	16	#include "AliPicoTrack.h"
	17	#include "AliEMCALGeometry.h"
	18	#include "AliLog.h"
	19
	20	#include "AliJetEmbeddingTask.h"
	21
	22	ClassImp(AliJetEmbeddingTask)
	23
	24	//________________________________________________________________________
	25	AliJetEmbeddingTask::AliJetEmbeddingTask() :
	26	AliAnalysisTaskSE("AliJetEmbeddingTask"),
	27	fGeomName(),
	28	fTracksName(),
	29	fOutTracksName(),
	30	fCaloName(),
	31	fOutCaloName(),
	32	fEtaMin(-1),
	33	fEtaMax(1),
	34	fPhiMin(0),
	35	fPhiMax(TMath::Pi() * 2),
	36	fPtMin(0),
	37	fPtMax(0),
	38	fCopyArray(kTRUE),
	39	fNEmbClusters(0),
	40	fNEmbTracks(0),
	41	fGeom(0),
	42	fClusters(0),
	43	fOutClusters(0),
	44	fTracks(0),
	45	fOutTracks(0)
	46	{
	47	// Default constructor.
	48	}
	49
	50	//________________________________________________________________________
	51	AliJetEmbeddingTask::AliJetEmbeddingTask(const char *name) :
	52	AliAnalysisTaskSE(name),
	53	fGeomName(""),
	54	fTracksName("PicoTracks"),
	55	fOutTracksName("PicoTracksEmbedded"),
	56	fCaloName("CaloClustersCorr"),
	57	fOutCaloName("CaloClustersCorrEmbedded"),
	58	fEtaMin(-1),
	59	fEtaMax(1),
	60	fPhiMin(0),
	61	fPhiMax(TMath::Pi() * 2),
	62	fPtMin(50),
	63	fPtMax(60),
	64	fCopyArray(kTRUE),
65	fNEmbClusters(0),
66	fNEmbTracks(1),
67	fGeom(0),
68	fClusters(0),
69	fOutClusters(0),
70	fTracks(0),
71	fOutTracks(0)
72	{
73	// Standard constructor.
74
75	}
76
77	//________________________________________________________________________
78	AliJetEmbeddingTask::~AliJetEmbeddingTask()
79	{
80	// Destructor
81	}
82
83	//________________________________________________________________________
84	void AliJetEmbeddingTask::Init()
85	{
86	// Init task.
87
88	if (fNEmbTracks > 0) {
89	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
90
91	if (strcmp(fTracks->GetClass()->GetName(), "AliPicoTrack")) {
92	AliError("Can only embed PicoTracks!");
93	return;
94	}
95
96	if (!fTracks) {
97	AliError(Form("Couldn't retrieve tracks with name %s!", fTracksName.Data()));
98	return;
99	}
100
101	if (!fOutTracks) {
102	fOutTracksName = fTracksName;
103	if (fCopyArray) {
104	fOutTracksName += "Embedded";
105	fOutTracks = new TClonesArray(*fTracks);
106	fOutTracks->SetName(fOutTracksName);
107	}
108	else {
109	fOutTracks = fTracks;
110	}
111	}
112
113	if (fCopyArray) {
114	if (!(InputEvent()->FindListObject(fOutTracksName)))
115	InputEvent()->AddObject(fOutTracks);
116	fOutTracks->Clear();
117	}
118	}
119
120	if (fNEmbClusters > 0) {
121	fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
122
123	if (!fClusters) {
124	AliError(Form("Couldn't retrieve clusters with name %s!", fCaloName.Data()));
125	return;
126	}
127
128	if (!fOutClusters) {
129	fOutCaloName = fCaloName;
130	if (fCopyArray) {
131	fOutCaloName += "Embedded";
132	fOutClusters = new TClonesArray(*fClusters);
133	fOutClusters->SetName(fOutCaloName);
134	}
135	else {
136	fOutClusters = fClusters;
137	}
138	}
139
140	if (fCopyArray) {
141	if (!(InputEvent()->FindListObject(fOutCaloName)))
142	InputEvent()->AddObject(fOutClusters);
143	fOutClusters->Clear();
144	}
145	}
146
808ba0db	147	if (!fGeom) {
	148	if (fGeomName.Length()>0) {
	149	fGeom = AliEMCALGeometry::GetInstance(fGeomName);
	150	if (!fGeom)
	151	AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
	152	} else {
	153	fGeom = AliEMCALGeometry::GetInstance();
	154	if (!fGeom)
	155	AliError("Could not get default geometry!");
44629b14	156	}
808ba0db	157	}
44629b14	158	}
	159
	160	//________________________________________________________________________
	161	void AliJetEmbeddingTask::Embed()
	162	{
	163	// Embed particles.
	164
	165	if (fNEmbClusters > 0 && fOutClusters) {
44629b14	166	const Int_t nClusters = fOutClusters->GetEntriesFast();
44629b14	167	TClonesArray digits("AliEMCALDigit", 1);
d69437c0	168	for (Int_t i = 0; i < fNEmbClusters; ++i) {
44629b14	169	Double_t pt = gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
	170	Double_t eta = gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
	171	Double_t phi = gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
	172
	173	TLorentzVector nPart;
	174	nPart.SetPtEtaPhiM(pt, eta, phi, 0);
	175	Double_t e = nPart.E();
	176
	177	Int_t absId = 0;
	178	fGeom->GetAbsCellIdFromEtaPhi(eta, phi, absId);
808ba0db	179
	180	if (absId == -1) {
	181	AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
	182	" Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
	183	eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
	184	continue;
	185	}
44629b14	186
	187	AliEMCALDigit digit = static_cast<AliEMCALDigit>(digits.New(0));
	188	digit->SetId(absId);
	189	digit->SetIndexInList(0);
	190	digit->SetType(AliEMCALDigit::kHG);
	191	digit->SetAmplitude(e);
	192
808ba0db	193	AliEMCALRecPoint recPoint("");
d69437c0	194	recPoint.AddDigit(*digit, e, kFALSE);
d69437c0	195	recPoint.EvalGlobalPosition(0, &digits);
44629b14	196
44629b14	197	TVector3 gpos;
d69437c0	198	recPoint.GetGlobalPosition(gpos);
44629b14	199	Float_t g[3];
	200	gpos.GetXYZ(g);
	201
	202	AliVCluster cluster = static_cast<AliVCluster>(fOutClusters->New(nClusters + i));
	203	cluster->SetType(AliVCluster::kEMCALClusterv1);
d69437c0	204	cluster->SetE(recPoint.GetEnergy());
44629b14	205	cluster->SetPosition(g);
	206	cluster->SetNCells(1);
	207	UShort_t shortAbsId = absId;
	208	cluster->SetCellsAbsId(&shortAbsId);
808ba0db	209	Double32_t fract = 1;
808ba0db	210	cluster->SetCellsAmplitudeFraction(&fract);
44629b14	211	cluster->SetID(nClusters + i);
	212	cluster->SetEmcCpvDistance(-1);
	213	cluster->SetChi2(100); // MC flag!
	214	}
	215	}
	216
	217	if (fNEmbTracks > 0 && fOutTracks) {
	218	Int_t nTracks = fOutTracks->GetEntriesFast();
d69437c0	219	for (Int_t i = 0; i < fNEmbTracks; ++i) {
44629b14	220	Double_t pt = gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
	221	Double_t eta = gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
	222	Double_t phi = gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
	223
	224	new ((*fOutTracks)[nTracks + i]) AliPicoTrack(pt,
	225	eta,
	226	phi,
	227	1,
	228	100, // MC flag!
	229	0,
	230	0,
	231	kFALSE);
	232	}
	233	}
	234	}
	235
	236	//________________________________________________________________________
	237	void AliJetEmbeddingTask::UserExec(Option_t *)
	238	{
	239	// Execute per event.
	240
	241	Init();
	242
	243	Embed();
	244	}