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