]> 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
set phys sel for jet finders
[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
507f6e75 184 if (fCopyArray) {
185 if (fOutTracks)
186 fOutTracks->Delete();
187 if (fOutClusters)
188 fOutClusters->Delete();
189 if (fOutMCParticles)
190 fOutMCParticles->Delete();
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//________________________________________________________________________
396ffc82 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, Double_t mass)
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,
396ffc82
MV		 702 ise,
703 mass);
b16bb001 704
ffe32451 705 return track;
706}
707
787a3c4f 708//________________________________________________________________________
709AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int_t origIndex)
710{
711 const Int_t nPart = fOutMCParticles->GetEntriesFast();
712
713 AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
714
a7477843 715 if (origIndex + fMCLabelShift >= fOutMCParticlesMap->GetSize())
716 fOutMCParticlesMap->Set((origIndex + fMCLabelShift)*2);
717
5be3857d 718 fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift);
5ce8ae64 719 AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)",
720 origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex));
787a3c4f 721
722 return aodpart;
723}
724
4358e58a 725//________________________________________________________________________
726void AliJetModelBaseTask::CopyCells()
727{
728 if (!fCaloCells)
729 return;
730
fde82e42 731 fAddedCells = 0;
732 fCaloCells->Sort();
4358e58a 733 for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
734 Int_t mclabel = 0;
735 Double_t efrac = 0.;
736 Double_t time = -1;
737 Short_t cellNum = -1;
738 Double_t amp = -1;
739
740 fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
741
742 if (!fIsMC)
743 mclabel = 0;
744
745 fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
fde82e42 746 fAddedCells++;
4358e58a 747 }
748
4358e58a 749 AliDebug(2, Form("%d cells from the current event", fAddedCells));
750}
751
ffe32451 752//________________________________________________________________________
753void AliJetModelBaseTask::CopyClusters()
754{
755 // Copy all the clusters in the new collection
787a3c4f 756 if (!fClusters)
757 return;
ffe32451 758
ffe32451 759 const Int_t nClusters = fClusters->GetEntriesFast();
56cc432c 760 Int_t nCopiedClusters = 0;
ffe32451 761
b16bb001 762 if (fEsdMode) {
ffe32451 763 for (Int_t i = 0; i < nClusters; ++i) {
764 AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
56cc432c 765 if (!esdcluster || !esdcluster->IsEMCAL())
ffe32451 766 continue;
4358e58a 767 AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
768 if (!fIsMC) {
769 TArrayI *labels = clus->GetLabelsArray();
770 if (labels)
771 labels->Reset();
772 }
56cc432c 773 nCopiedClusters++;
ffe32451 774 }
775 }
776 else {
777 for (Int_t i = 0; i < nClusters; ++i) {
778 AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
56cc432c 779 if (!aodcluster || !aodcluster->IsEMCAL())
ffe32451 780 continue;
4358e58a 781 AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
782 if (!fIsMC)
783 clus->SetLabel(0,0);
56cc432c 784 nCopiedClusters++;
ffe32451 785 }
786 }
787}
788
789//________________________________________________________________________
790void AliJetModelBaseTask::CopyTracks()
791{
787a3c4f 792 // Copy all the tracks in the new collection
793
794 if (!fTracks)
795 return;
796
ffe32451 797 const Int_t nTracks = fTracks->GetEntriesFast();
56cc432c 798 Int_t nCopiedTracks = 0;
ffe32451 799 for (Int_t i = 0; i < nTracks; ++i) {
4358e58a 800 AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
801 if (!picotrack)
ffe32451 802 continue;
4358e58a 803 AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack);
804 if (!fIsMC && track->GetLabel() != 0)
805 track->SetLabel(0);
56cc432c 806 nCopiedTracks++;
ffe32451 807 }
808}
809
787a3c4f 810//________________________________________________________________________
811void AliJetModelBaseTask::CopyMCParticles()
812{
813 // Copy all the MC particles in the new collection
814
815 if (!fMCParticles)
816 return;
817
818 const Int_t nPart = fMCParticles->GetEntriesFast();
819 Int_t nCopiedPart = 0;
820 for (Int_t i = 0; i < nPart; ++i) {
821 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fMCParticles->At(i));
822 if (!part)
823 continue;
824 new ((*fOutMCParticles)[nCopiedPart]) AliAODMCParticle(*part);
825
826 nCopiedPart++;
827 }
828
a7477843 829 if (!fMCParticlesMap || !fOutMCParticlesMap)
787a3c4f 830 return;
831
a7477843 832 if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize())
833 fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2);
834
5be3857d 835 for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) {
836 fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i);
837 if (fMCParticlesMap->At(i) >= 0)
5ce8ae64 838 fMCLabelShift = i;
787a3c4f 839 }
840
5ce8ae64 841 AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift));
787a3c4f 842}
843
762e8424 844//________________________________________________________________________
845void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
846{
acd859a2 847 // Get random cell.
848
762e8424 849 Int_t repeats = 0;
05debff6 850 Double_t rndEta = eta;
851 Double_t rndPhi = phi;
762e8424 852 do {
05debff6 853 if (eta < -100)
2103dc6a 854 rndEta = GetRandomEta(kTRUE);
05debff6 855 if (phi < 0)
2103dc6a 856 rndPhi = GetRandomPhi(kTRUE);
05debff6 857 fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
762e8424 858 repeats++;
859 } while (absId == -1 && repeats < 100);
860
861 if (!(absId > -1)) {
862 AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
863 "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
864 }
05debff6 865 else {
866 eta = rndEta;
867 phi = rndPhi;
868 }
762e8424 869}
870
871//________________________________________________________________________
2103dc6a 872Double_t AliJetModelBaseTask::GetRandomEta(Bool_t emcal)
762e8424 873{
acd859a2 874 // Get random eta.
875
2103dc6a 876 Double_t etamax = fEtaMax;
877 Double_t etamin = fEtaMin;
878
879 if (emcal) {
880 const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
881 const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
882
883 if (etamax > EmcalMaxEta) etamax = EmcalMaxEta;
884 if (etamax < EmcalMinEta) etamax = EmcalMinEta;
885 if (etamin > EmcalMaxEta) etamin = EmcalMaxEta;
886 if (etamin < EmcalMinEta) etamin = EmcalMinEta;
887 }
888
889 return gRandom->Rndm() * (etamax - etamin) + etamin;
762e8424 890}
891
892//________________________________________________________________________
2103dc6a 893Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal)
762e8424 894{
acd859a2 895 // Get random phi.
2103dc6a 896
897 Double_t phimax = fPhiMax;
898 Double_t phimin = fPhiMin;
899
900 if (emcal) {
901 const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
902 const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
903
904 if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
905 if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
906 if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
907 if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
908 }
acd859a2 909
43032ce2 910 Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
911
912 return result;
762e8424 913}
914
915//________________________________________________________________________
916Double_t AliJetModelBaseTask::GetRandomPt()
917{
acd859a2 918 // Get random pt.
919
e44e8726 920 if (fPtSpectrum)
921 return fPtSpectrum->GetRandom();
922 else
923 return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
762e8424 924}
925
43032ce2 926//________________________________________________________________________
927void AliJetModelBaseTask::GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal)
928{
929 // Get a random particle.
930
931 eta = GetRandomEta(emcal);
932
933 if (fPtPhiEvPlDistribution) {
934 Double_t phimax = fPhiMax;
935 Double_t phimin = fPhiMin;
936
937 if (emcal) {
938 const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
939 const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
940
941 if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
942 if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
943 if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
944 if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
945 }
946
947 if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) {
948 AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName()));
949 pt = GetRandomPt();
950 phi = GetRandomPhi(emcal);
951 }
952 else {
953 do {
954 fPtPhiEvPlDistribution->GetRandom2(pt,phi);
955 phi += fPsi;
956 if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2;
957 } while (phi > phimax || phi < phimin);
958 }
959 }
960 else {
961 pt = GetRandomPt();
962 phi = GetRandomPhi(emcal);
963 }
964}
965
762e8424 966//________________________________________________________________________
967void AliJetModelBaseTask::Run()
968{
acd859a2 969 // Run.
762e8424 970}