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