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