]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetModelBaseTask.cxx
CommitLineData
762e8424 1// $Id$
2//
3// Jet modelling task.
4//
297edd60 5// Author: S.Aiola, C.Loizides
762e8424 6
ffe32451 7#include "AliJetModelBaseTask.h"
8
762e8424 9#include <TClonesArray.h>
787a3c4f 10#include <TH1I.h>
762e8424 11#include <TLorentzVector.h>
12#include <TRandom3.h>
b16bb001 13#include <TList.h>
43032ce2 14#include <TF2.h>
762e8424 15
b16bb001 16#include "AliVEvent.h"
e44e8726 17#include "AliAODCaloCluster.h"
b16bb001 18#include "AliESDCaloCluster.h"
19#include "AliVCluster.h"
762e8424 20#include "AliEMCALDigit.h"
ffe32451 21#include "AliEMCALRecPoint.h"
b16bb001 22#include "AliESDCaloCells.h"
23#include "AliAODCaloCells.h"
787a3c4f 24#include "AliAODMCParticle.h"
b16bb001 25#include "AliVCaloCells.h"
ffe32451 26#include "AliPicoTrack.h"
b16bb001 27#include "AliEMCALGeometry.h"
28#include "AliLog.h"
5be3857d 29#include "AliNamedArrayI.h"
762e8424 30
31ClassImp(AliJetModelBaseTask)
32
33//________________________________________________________________________
34AliJetModelBaseTask::AliJetModelBaseTask() :
35 AliAnalysisTaskSE("AliJetModelBaseTask"),
36 fGeomName(),
37 fTracksName(),
38 fOutTracksName(),
39 fCaloName(),
40 fOutCaloName(),
b16bb001 41 fCellsName(),
42 fOutCellsName(),
787a3c4f 43 fMCParticlesName(),
44 fOutMCParticlesName(),
4358e58a 45 fIsMC(kFALSE),
762e8424 46 fSuffix(),
47 fEtaMin(-1),
48 fEtaMax(1),
49 fPhiMin(0),
50 fPhiMax(TMath::Pi() * 2),
51 fPtMin(0),
52 fPtMax(0),
53 fCopyArray(kTRUE),
54 fNClusters(0),
b16bb001 55 fNCells(0),
762e8424 56 fNTracks(0),
7f76e479 57 fMarkMC(99999),
e44e8726 58 fPtSpectrum(0),
43032ce2 59 fPtPhiEvPlDistribution(0),
27f77648 60 fDensitySpectrum(0),
b16bb001 61 fQAhistos(kFALSE),
43032ce2 62 fPsi(0),
ffe32451 63 fIsInit(0),
762e8424 64 fGeom(0),
65 fClusters(0),
66 fOutClusters(0),
67 fTracks(0),
b16bb001 68 fOutTracks(0),
69 fCaloCells(0),
70 fOutCaloCells(0),
71 fAddedCells(0),
787a3c4f 72 fMCParticles(0),
73 fMCParticlesMap(0),
74 fOutMCParticles(0),
75 fOutMCParticlesMap(0),
76 fMCLabelShift(0),
b16bb001 77 fEsdMode(kFALSE),
78 fOutput(0)
762e8424 79{
80 // Default constructor.
8e48f7f5 81
82 fVertex[0] = 0;
83 fVertex[1] = 0;
84 fVertex[2] = 0;
762e8424 85}
86
87//________________________________________________________________________
b16bb001 88AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
762e8424 89 AliAnalysisTaskSE(name),
90 fGeomName(""),
91 fTracksName("PicoTracks"),
92 fOutTracksName("PicoTracksEmbedded"),
93 fCaloName("CaloClustersCorr"),
94 fOutCaloName("CaloClustersCorrEmbedded"),
b16bb001 95 fCellsName(""),
96 fOutCellsName(""),
787a3c4f 97 fMCParticlesName(""),
98 fOutMCParticlesName(""),
4358e58a 99 fIsMC(kFALSE),
762e8424 100 fSuffix("Processed"),
101 fEtaMin(-1),
102 fEtaMax(1),
103 fPhiMin(0),
104 fPhiMax(TMath::Pi() * 2),
105 fPtMin(50),
106 fPtMax(60),
107 fCopyArray(kTRUE),
108 fNClusters(0),
b16bb001 109 fNCells(0),
762e8424 110 fNTracks(1),
7f76e479 111 fMarkMC(99999),
e44e8726 112 fPtSpectrum(0),
43032ce2 113 fPtPhiEvPlDistribution(0),
27f77648 114 fDensitySpectrum(0),
b16bb001 115 fQAhistos(drawqa),
43032ce2 116 fPsi(0),
ffe32451 117 fIsInit(0),
762e8424 118 fGeom(0),
119 fClusters(0),
120 fOutClusters(0),
121 fTracks(0),
b16bb001 122 fOutTracks(0),
123 fCaloCells(0),
124 fOutCaloCells(0),
125 fAddedCells(0),
787a3c4f 126 fMCParticles(0),
127 fMCParticlesMap(0),
128 fOutMCParticles(0),
129 fOutMCParticlesMap(0),
130 fMCLabelShift(0),
b16bb001 131 fEsdMode(kFALSE),
132 fOutput(0)
762e8424 133{
134 // Standard constructor.
b16bb001 135
136 if (fQAhistos) {
137 DefineOutput(1, TList::Class());
138 }
8e48f7f5 139
140 fVertex[0] = 0;
141 fVertex[1] = 0;
142 fVertex[2] = 0;
762e8424 143}
144
145//________________________________________________________________________
146AliJetModelBaseTask::~AliJetModelBaseTask()
147{
148 // Destructor
149}
150
b16bb001 151//________________________________________________________________________
152void AliJetModelBaseTask::UserCreateOutputObjects()
153{
154 // Create user output.
155 if (!fQAhistos)
156 return;
157
158 OpenFile(1);
159 fOutput = new TList();
160 fOutput->SetOwner();
161
162 PostData(1, fOutput);
163}
164
762e8424 165//________________________________________________________________________
ffe32451 166void AliJetModelBaseTask::UserExec(Option_t *)
167{
168 // Execute per event.
169
d63c8c07 170 if (!fIsInit)
171 fIsInit = ExecOnce();
172
173 if (!fIsInit)
174 return;
175
6fd2f1a9 176 fVertex[0] = 0;
177 fVertex[1] = 0;
178 fVertex[2] = 0;
179
180 const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
181 if (vert)
182 vert->GetXYZ(fVertex);
183
d63c8c07 184 if (fCopyArray) {
185 if (fOutTracks)
186 fOutTracks->Delete();
187 if (fOutClusters)
188 fOutClusters->Delete();
787a3c4f 189 if (fOutMCParticles)
190 fOutMCParticles->Delete();
b4ffb4d2 191 }
27f77648 192
193 if (fDensitySpectrum) {
43032ce2 194 fNTracks = TMath::Nint(fDensitySpectrum->GetRandom());
195 fNCells = TMath::Nint(fDensitySpectrum->GetRandom());
196 fNClusters = TMath::Nint(fDensitySpectrum->GetRandom());
27f77648 197 }
63fac07f 198
5be3857d 199 // Clear map
200 if (fOutMCParticlesMap)
4358e58a 201 fOutMCParticlesMap->Clear();
d63c8c07 202
cd6431de 203 AliVCaloCells *tempCaloCells = 0;
204
b16bb001 205 if (fCaloCells) {
206 fAddedCells = 0;
cd6431de 207 if (!fCopyArray) {
208 tempCaloCells = fCaloCells;
209 fCaloCells = static_cast<AliVCaloCells*>(tempCaloCells->Clone(Form("%s_old",fCaloCells->GetName())));
210 }
b16bb001 211 }
212
43032ce2 213 if (fPtPhiEvPlDistribution)
214 fPsi = gRandom->Rndm() * TMath::Pi();
215
ffe32451 216 Run();
cd6431de 217
fde82e42 218 if (fCaloCells && !fCopyArray) {
cd6431de 219 delete fCaloCells;
220 fCaloCells = tempCaloCells;
221 }
ffe32451 222}
223
b16bb001 224//________________________________________________________________________
225Bool_t AliJetModelBaseTask::ExecOnce()
226{
227 // Init task.
228
229 delete gRandom;
230 gRandom = new TRandom3(0);
231
232 fEsdMode = InputEvent()->InheritsFrom("AliESDEvent");
233
234 if (fPtMax < fPtMin) {
235 AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
236 fPtMax = fPtMin;
237 }
238
239 if (fEtaMax < fEtaMin) {
240 AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
241 fEtaMax = fEtaMin;
242 }
243
244 if (fPhiMax < fPhiMin) {
245 AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
246 fPhiMax = fPhiMin;
247 }
248
787a3c4f 249 if (!fCellsName.IsNull()) {
b16bb001 250 fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCellsName));
251 if (!fCaloCells) {
787a3c4f 252 AliWarning(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data()));
b16bb001 253 }
787a3c4f 254 else if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
b16bb001 255 AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data()));
256 fCaloCells = 0;
257 return kFALSE;
258 }
259
260 if (!fOutCaloCells) {
261 fOutCellsName = fCellsName;
787a3c4f 262 if (fCopyArray)
b16bb001 263 fOutCellsName += fSuffix;
787a3c4f 264 if (fCopyArray || !fCaloCells) {
b16bb001 265 if (fEsdMode)
266 fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName);
267 else
268 fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName);
269
270 if (InputEvent()->FindListObject(fOutCellsName)) {
271 AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data()));
272 return kFALSE;
273 }
274 else {
275 InputEvent()->AddObject(fOutCaloCells);
276 }
277 }
278 else {
279 fOutCaloCells = fCaloCells;
280 }
281 }
282 }
283
27f77648 284 if (!fTracksName.IsNull()) {
b16bb001 285 fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
286 if (!fTracks) {
787a3c4f 287 AliWarning(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
b16bb001 288 }
787a3c4f 289 else if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
b16bb001 290 AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data()));
291 fTracks = 0;
292 return kFALSE;
293 }
294
295 if (!fOutTracks) {
296 fOutTracksName = fTracksName;
787a3c4f 297 if (fCopyArray)
b16bb001 298 fOutTracksName += fSuffix;
787a3c4f 299 if (fCopyArray || !fTracks) {
300 fOutTracks = new TClonesArray("AliPicoTrack");
b16bb001 301 fOutTracks->SetName(fOutTracksName);
302 if (InputEvent()->FindListObject(fOutTracksName)) {
303 AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data()));
304 return kFALSE;
305 }
306 else {
307 InputEvent()->AddObject(fOutTracks);
308 }
309 }
310 else {
311 fOutTracks = fTracks;
312 }
313 }
314 }
315
27f77648 316 if (!fCaloName.IsNull()) {
b16bb001 317 fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
318
319 if (!fClusters) {
787a3c4f 320 AliWarning(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
b16bb001 321 }
787a3c4f 322 else if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
b16bb001 323 AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data()));
324 fClusters = 0;
325 return kFALSE;
326 }
327
328 if (!fOutClusters) {
329 fOutCaloName = fCaloName;
787a3c4f 330 if (fCopyArray)
b16bb001 331 fOutCaloName += fSuffix;
507f74bc 332 TString className;
333 if (fClusters)
334 className = fClusters->GetClass()->GetName();
335 else if (fEsdMode)
336 className = "AliESDCaloCluster";
337 else
338 className = "AliAODCaloCluster";
339
787a3c4f 340 if (fCopyArray || !fClusters) {
507f74bc 341 fOutClusters = new TClonesArray(className.Data());
b16bb001 342 fOutClusters->SetName(fOutCaloName);
343 if (InputEvent()->FindListObject(fOutCaloName)) {
344 AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data()));
345 return kFALSE;
346 }
347 else {
348 InputEvent()->AddObject(fOutClusters);
349 }
350 }
351 else {
352 fOutClusters = fClusters;
353 }
354 }
355 }
356
787a3c4f 357 if (!fMCParticlesName.IsNull()) {
358 fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
359 if (!fMCParticles) {
360 AliWarning(Form("%s: Couldn't retrieve MC particles with name %s!", GetName(), fMCParticlesName.Data()));
361 }
362 else {
363 if (!fMCParticles->GetClass()->GetBaseClass("AliAODMCParticle")) {
364 AliError(Form("%s: Collection %s does not contain AliAODMCParticle objects!", GetName(), fMCParticlesName.Data()));
365 fMCParticles = 0;
366 return kFALSE;
367 }
368
5be3857d 369 fMCParticlesMap = dynamic_cast<AliNamedArrayI*>(InputEvent()->FindListObject(fMCParticlesName + "_Map"));
787a3c4f 370
371 if (!fMCParticlesMap) {
372 AliWarning(Form("%s: Could not retrieve map for MC particles %s! Will assume MC labels consistent with indexes...", GetName(), fMCParticlesName.Data()));
a7477843 373 fMCParticlesMap = new AliNamedArrayI(fMCParticlesName + "_Map", 99999);
374 for (Int_t i = 0; i < 99999; i++) {
5be3857d 375 fMCParticlesMap->AddAt(i,i);
787a3c4f 376 }
377 }
378 }
379
380 if (!fOutMCParticles) {
381 fOutMCParticlesName = fMCParticlesName;
382 if (fCopyArray)
383 fOutMCParticlesName += fSuffix;
384 if (fCopyArray || !fMCParticles) {
385
386 fOutMCParticles = new TClonesArray("AliAODMCParticle");
387 fOutMCParticles->SetName(fOutMCParticlesName);
388 if (InputEvent()->FindListObject(fOutMCParticlesName)) {
389 AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutMCParticlesName.Data()));
390 return kFALSE;
391 }
392 else {
393 InputEvent()->AddObject(fOutMCParticles);
394 }
395
a7477843 396 fOutMCParticlesMap = new AliNamedArrayI(fOutMCParticlesName + "_Map",99999);
787a3c4f 397 if (InputEvent()->FindListObject(fOutMCParticlesName + "_Map")) {
398 AliFatal(Form("%s: Map %s_Map is already present in the event!", GetName(), fOutMCParticlesName.Data()));
399 return kFALSE;
400 }
401 else {
402 InputEvent()->AddObject(fOutMCParticlesMap);
403 }
404 }
405 else {
406 fOutMCParticles = fMCParticles;
407 fOutMCParticlesMap = fMCParticlesMap;
408 }
409 }
410 }
411
2103dc6a 412 if (!fGeom) {
413 if (fGeomName.Length() > 0) {
414 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
415 if (!fGeom) {
416 AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data()));
417 return kFALSE;
418 }
419 } else {
4358e58a 420 fGeom = AliEMCALGeometry::GetInstance();
421 if (!fGeom) {
422 AliFatal("Could not get default geometry!");
423 return kFALSE;
424 }
b16bb001 425 }
cfc2ac24 426 }
b16bb001 427
428 return kTRUE;
429}
430
431//________________________________________________________________________
432Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
433{
434 if (fOutCaloCells->GetNumberOfCells() < n) {
435 fOutCaloCells->DeleteContainer();
436 fOutCaloCells->CreateContainer(n);
437 }
438 else {
439 fOutCaloCells->SetNumberOfCells(n);
440 }
441
442 fAddedCells = 0;
443
444 return n;
445}
446
b16bb001 447//________________________________________________________________________
448Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi)
449{
450 // Add a cell to the event.
451
452 Int_t absId = 0;
453 if (eta < -100 || phi < 0) {
454 GetRandomCell(eta, phi, absId);
455 }
456 else {
457 fGeom->EtaPhiFromIndex(absId, eta, phi);
458 }
459
460 if (absId == -1) {
461 AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!"
462 " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
463 eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
464 return 0;
465 }
466
467 if (e < 0) {
468 Double_t pt = GetRandomPt();
469 TLorentzVector nPart;
470 nPart.SetPtEtaPhiM(pt, eta, phi, 0);
471 e = nPart.E();
472 }
473
474 return AddCell(e, absId);
475}
476
477//________________________________________________________________________
787a3c4f 478Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time, Int_t label)
b16bb001 479{
480 // Add a cell to the event.
481
a7477843 482 if (label < 0)
483 label = 0;
484
485 label += fMarkMC + fMCLabelShift;
787a3c4f 486
fde82e42 487 Short_t pos = -1;
488 if (fCaloCells)
489 pos = fCaloCells->GetCellPosition(absId);
490
4358e58a 491 Double_t efrac = 1;
492 Bool_t increaseOnSuccess = kFALSE;
493
494 if (pos < 0) {
495 increaseOnSuccess = kTRUE;
496 pos = fAddedCells;
497 }
498 else {
499 Short_t cellNumber = -1;
500 Double_t old_e = 0;
501 Double_t old_time = 0;
502 Int_t old_label = 0;
503 Double_t old_efrac = 0;
504 fOutCaloCells->GetCell(pos, cellNumber, old_e, old_time, old_label, old_efrac);
505
fde82e42 506 efrac = e / (old_e + e);
507
508 if (old_label > 0 && e < old_efrac * old_e) {
4358e58a 509 label = old_label;
fde82e42 510 efrac = old_efrac;
4358e58a 511 time = old_time;
512 }
513
4358e58a 514 e += old_e;
515 }
516
517 Bool_t r = fOutCaloCells->SetCell(pos, absId, e, time, label, efrac);
b16bb001 518
519 if (r) {
4358e58a 520 if (increaseOnSuccess)
521 fAddedCells++;
b16bb001 522 return fAddedCells;
523 }
4358e58a 524 else
b16bb001 525 return 0;
b16bb001 526}
527
507f74bc 528//________________________________________________________________________
529AliVCluster* AliJetModelBaseTask::AddCluster(AliVCluster *oc)
530{
531 // Add a cluster to the event.
532
533 const Int_t nClusters = fOutClusters->GetEntriesFast();
534 AliVCluster *dc = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
535 dc->SetType(AliVCluster::kEMCALClusterv1);
536 dc->SetE(oc->E());
537 Float_t pos[3] = {0};
538 oc->GetPosition(pos);
539 dc->SetPosition(pos);
540 dc->SetNCells(oc->GetNCells());
541 dc->SetCellsAbsId(oc->GetCellsAbsId());
542 dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
543 dc->SetID(oc->GetID());
544 dc->SetDispersion(oc->GetDispersion());
545 dc->SetEmcCpvDistance(-1);
546 dc->SetChi2(-1);
547 dc->SetTOF(oc->GetTOF()); //time-of-flight
548 dc->SetNExMax(oc->GetNExMax()); //number of local maxima
549 dc->SetM02(oc->GetM02());
550 dc->SetM20(oc->GetM20());
551 dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel());
552
553 //MC
554 UInt_t nlabels = oc->GetNLabels();
555 Int_t *labels = oc->GetLabels();
a7477843 556 TArrayI parents;
557 if (nlabels != 0 && labels)
558 parents.Set(nlabels, labels);
559 else {
560 nlabels = 1;
561 parents.Set(1);
562 parents[0] = 0;
563 }
507f74bc 564
a7477843 565 if (fMarkMC+fMCLabelShift != 0) {
566 for (UInt_t i = 0; i < nlabels; i++) {
567 parents[i] += fMarkMC+fMCLabelShift;
507f74bc 568 }
569 }
a7477843 570
571 AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
572 if (esdClus) {
573 esdClus->AddLabels(parents);
574 }
575 else {
576 AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
577 if (aodClus)
578 aodClus->SetLabel(parents.GetArray(), nlabels);
fde82e42 579 }
507f74bc 580
581 return dc;
582}
b16bb001 583
ffe32451 584//________________________________________________________________________
787a3c4f 585AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi, Int_t label)
ffe32451 586{
587 // Add a cluster to the event.
588
589 Int_t absId = 0;
590 if (eta < -100 || phi < 0) {
591 GetRandomCell(eta, phi, absId);
592 }
593 else {
594 fGeom->EtaPhiFromIndex(absId, eta, phi);
595 }
596
597 if (absId == -1) {
598 AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
599 " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
600 eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
601 return 0;
602 }
603
604 if (e < 0) {
605 Double_t pt = GetRandomPt();
606 TLorentzVector nPart;
607 nPart.SetPtEtaPhiM(pt, eta, phi, 0);
608 e = nPart.E();
609 }
610
787a3c4f 611 return AddCluster(e, absId, label);
ffe32451 612}
613
614//________________________________________________________________________
787a3c4f 615AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t label)
ffe32451 616{
617 // Add a cluster to the event.
618
619 const Int_t nClusters = fOutClusters->GetEntriesFast();
620
621 TClonesArray digits("AliEMCALDigit", 1);
622
623 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
624 digit->SetId(absId);
625 digit->SetIndexInList(0);
626 digit->SetType(AliEMCALDigit::kHG);
627 digit->SetAmplitude(e);
628
629 AliEMCALRecPoint recPoint("");
630 recPoint.AddDigit(*digit, e, kFALSE);
631 recPoint.EvalGlobalPosition(0, &digits);
632
633 TVector3 gpos;
634 recPoint.GetGlobalPosition(gpos);
635 Float_t g[3];
636 gpos.GetXYZ(g);
637
638 AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
639 cluster->SetType(AliVCluster::kEMCALClusterv1);
640 cluster->SetE(recPoint.GetEnergy());
641 cluster->SetPosition(g);
642 cluster->SetNCells(1);
643 UShort_t shortAbsId = absId;
644 cluster->SetCellsAbsId(&shortAbsId);
645 Double32_t fract = 1;
646 cluster->SetCellsAmplitudeFraction(&fract);
647 cluster->SetID(nClusters);
648 cluster->SetEmcCpvDistance(-1);
ffe32451 649
f660c2d6 650 //MC label
a7477843 651 if (label < 0)
652 label = 0;
653
654 label += fMarkMC+fMCLabelShift;
787a3c4f 655
f660c2d6 656 if (fEsdMode) {
657 AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster);
787a3c4f 658 TArrayI parents(1, &label);
f660c2d6 659 esdClus->AddLabels(parents);
660 }
661 else {
662 AliAODCaloCluster *aodClus = static_cast<AliAODCaloCluster*>(cluster);
787a3c4f 663 aodClus->SetLabel(&label, 1);
f660c2d6 664 }
665
ffe32451 666 return cluster;
667}
668
669//________________________________________________________________________
56bd3193 670AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Double_t ptemc, Bool_t ise, Int_t label, Short_t charge)
ffe32451 671{
672 // Add a track to the event.
ffe32451 673
43032ce2 674 if (pt < 0 && eta < -100 && phi < 0) {
675 GetRandomParticle(pt,eta,phi);
676 }
677 else {
678 if (pt < 0)
679 pt = GetRandomPt();
680 if (eta < -100)
681 eta = GetRandomEta();
682 if (phi < 0)
683 phi = GetRandomPhi(pt);
684 }
ffe32451 685
a7477843 686 if (label >= 0)
687 label += fMarkMC+fMCLabelShift;
688 else if (label < 0)
689 label -= fMarkMC+fMCLabelShift;
43032ce2 690
691 const Int_t nTracks = fOutTracks->GetEntriesFast();
787a3c4f 692
ffe32451 693 AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt,
694 eta,
695 phi,
6fd2f1a9 696 charge,
787a3c4f 697 label,
b16bb001 698 type,
699 etaemc,
56bd3193 700 phiemc,
701 ptemc,
b16bb001 702 ise);
703
ffe32451 704 return track;
705}
706
787a3c4f 707//________________________________________________________________________
708AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int_t origIndex)
709{
710 const Int_t nPart = fOutMCParticles->GetEntriesFast();
711
712 AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
713
a7477843 714 if (origIndex + fMCLabelShift >= fOutMCParticlesMap->GetSize())
715 fOutMCParticlesMap->Set((origIndex + fMCLabelShift)*2);
716
5be3857d 717 fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift);
5ce8ae64 718 AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)",
719 origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex));
787a3c4f 720
721 return aodpart;
722}
723
4358e58a 724//________________________________________________________________________
725void AliJetModelBaseTask::CopyCells()
726{
727 if (!fCaloCells)
728 return;
729
fde82e42 730 fAddedCells = 0;
731 fCaloCells->Sort();
4358e58a 732 for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
733 Int_t mclabel = 0;
734 Double_t efrac = 0.;
735 Double_t time = -1;
736 Short_t cellNum = -1;
737 Double_t amp = -1;
738
739 fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
740
741 if (!fIsMC)
742 mclabel = 0;
743
744 fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
fde82e42 745 fAddedCells++;
4358e58a 746 }
747
4358e58a 748 AliDebug(2, Form("%d cells from the current event", fAddedCells));
749}
750
ffe32451 751//________________________________________________________________________
752void AliJetModelBaseTask::CopyClusters()
753{
754 // Copy all the clusters in the new collection
787a3c4f 755 if (!fClusters)
756 return;
ffe32451 757
ffe32451 758 const Int_t nClusters = fClusters->GetEntriesFast();
56cc432c 759 Int_t nCopiedClusters = 0;
ffe32451 760
b16bb001 761 if (fEsdMode) {
ffe32451 762 for (Int_t i = 0; i < nClusters; ++i) {
763 AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
56cc432c 764 if (!esdcluster || !esdcluster->IsEMCAL())
ffe32451 765 continue;
4358e58a 766 AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
767 if (!fIsMC) {
768 TArrayI *labels = clus->GetLabelsArray();
769 if (labels)
770 labels->Reset();
771 }
56cc432c 772 nCopiedClusters++;
ffe32451 773 }
774 }
775 else {
776 for (Int_t i = 0; i < nClusters; ++i) {
777 AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
56cc432c 778 if (!aodcluster || !aodcluster->IsEMCAL())
ffe32451 779 continue;
4358e58a 780 AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
781 if (!fIsMC)
782 clus->SetLabel(0,0);
56cc432c 783 nCopiedClusters++;
ffe32451 784 }
785 }
786}
787
788//________________________________________________________________________
789void AliJetModelBaseTask::CopyTracks()
790{
787a3c4f 791 // Copy all the tracks in the new collection
792
793 if (!fTracks)
794 return;
795
ffe32451 796 const Int_t nTracks = fTracks->GetEntriesFast();
56cc432c 797 Int_t nCopiedTracks = 0;
ffe32451 798 for (Int_t i = 0; i < nTracks; ++i) {
4358e58a 799 AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
800 if (!picotrack)
ffe32451 801 continue;
4358e58a 802 AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack);
803 if (!fIsMC && track->GetLabel() != 0)
804 track->SetLabel(0);
56cc432c 805 nCopiedTracks++;
ffe32451 806 }
807}
808
787a3c4f 809//________________________________________________________________________
810void AliJetModelBaseTask::CopyMCParticles()
811{
812 // Copy all the MC particles in the new collection
813
814 if (!fMCParticles)
815 return;
816
817 const Int_t nPart = fMCParticles->GetEntriesFast();
818 Int_t nCopiedPart = 0;
819 for (Int_t i = 0; i < nPart; ++i) {
820 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fMCParticles->At(i));
821 if (!part)
822 continue;
823 new ((*fOutMCParticles)[nCopiedPart]) AliAODMCParticle(*part);
824
825 nCopiedPart++;
826 }
827
a7477843 828 if (!fMCParticlesMap || !fOutMCParticlesMap)
787a3c4f 829 return;
830
a7477843 831 if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize())
832 fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2);
833
5be3857d 834 for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) {
835 fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i);
836 if (fMCParticlesMap->At(i) >= 0)
5ce8ae64 837 fMCLabelShift = i;
787a3c4f 838 }
839
5ce8ae64 840 AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift));
787a3c4f 841}
842
762e8424 843//________________________________________________________________________
844void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
845{
acd859a2 846 // Get random cell.
847
762e8424 848 Int_t repeats = 0;
05debff6 849 Double_t rndEta = eta;
850 Double_t rndPhi = phi;
762e8424 851 do {
05debff6 852 if (eta < -100)
2103dc6a 853 rndEta = GetRandomEta(kTRUE);
05debff6 854 if (phi < 0)
2103dc6a 855 rndPhi = GetRandomPhi(kTRUE);
05debff6 856 fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
762e8424 857 repeats++;
858 } while (absId == -1 && repeats < 100);
859
860 if (!(absId > -1)) {
861 AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
862 "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
863 }
05debff6 864 else {
865 eta = rndEta;
866 phi = rndPhi;
867 }
762e8424 868}
869
870//________________________________________________________________________
2103dc6a 871Double_t AliJetModelBaseTask::GetRandomEta(Bool_t emcal)
762e8424 872{
acd859a2 873 // Get random eta.
874
2103dc6a 875 Double_t etamax = fEtaMax;
876 Double_t etamin = fEtaMin;
877
878 if (emcal) {
879 const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
880 const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
881
882 if (etamax > EmcalMaxEta) etamax = EmcalMaxEta;
883 if (etamax < EmcalMinEta) etamax = EmcalMinEta;
884 if (etamin > EmcalMaxEta) etamin = EmcalMaxEta;
885 if (etamin < EmcalMinEta) etamin = EmcalMinEta;
886 }
887
888 return gRandom->Rndm() * (etamax - etamin) + etamin;
762e8424 889}
890
891//________________________________________________________________________
2103dc6a 892Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal)
762e8424 893{
acd859a2 894 // Get random phi.
2103dc6a 895
896 Double_t phimax = fPhiMax;
897 Double_t phimin = fPhiMin;
898
899 if (emcal) {
900 const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
901 const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
902
903 if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
904 if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
905 if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
906 if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
907 }
acd859a2 908
43032ce2 909 Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
910
911 return result;
762e8424 912}
913
914//________________________________________________________________________
915Double_t AliJetModelBaseTask::GetRandomPt()
916{
acd859a2 917 // Get random pt.
918
e44e8726 919 if (fPtSpectrum)
920 return fPtSpectrum->GetRandom();
921 else
922 return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
762e8424 923}
924
43032ce2 925//________________________________________________________________________
926void AliJetModelBaseTask::GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal)
927{
928 // Get a random particle.
929
930 eta = GetRandomEta(emcal);
931
932 if (fPtPhiEvPlDistribution) {
933 Double_t phimax = fPhiMax;
934 Double_t phimin = fPhiMin;
935
936 if (emcal) {
937 const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
938 const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
939
940 if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
941 if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
942 if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
943 if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
944 }
945
946 if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) {
947 AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName()));
948 pt = GetRandomPt();
949 phi = GetRandomPhi(emcal);
950 }
951 else {
952 do {
953 fPtPhiEvPlDistribution->GetRandom2(pt,phi);
954 phi += fPsi;
955 if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2;
956 } while (phi > phimax || phi < phimin);
957 }
958 }
959 else {
960 pt = GetRandomPt();
961 phi = GetRandomPhi(emcal);
962 }
963}
964
762e8424 965//________________________________________________________________________
966void AliJetModelBaseTask::Run()
967{
acd859a2 968 // Run.
762e8424 969}