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