PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx

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