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", 9999);
 342         for (Int_t i = 0; i < 9999; 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",9999);
 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 = fMarkMC;
 452   else
 453     label += 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
 525   if (nlabels != 0 && labels && labels[0] >= 0) {
 526     AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
 527     if (esdClus) {
 528       TArrayI parents(nlabels, labels);
 529       esdClus->AddLabels(parents);
 530     }
 531     else {
 532       AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
 533       if (aodClus)
 534         aodClus->SetLabel(labels, nlabels);
 535     }
 536   }
 537   else if (fMarkMC != 0) {
 538     AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
 539     if (esdClus) {
 540       TArrayI parents(1, &fMarkMC);
 541       esdClus->AddLabels(parents);
 542     }
 543     else {
 544       AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
 545       if (aodClus)
 546         aodClus->SetLabel(&fMarkMC,1);
 547     }
 548   }
 549
 550   return dc;
 551 }
 552
 553 //________________________________________________________________________
 554 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi, Int_t label)
 555 {
 556   // Add a cluster to the event.
 557
 558   Int_t absId = 0;
 559   if (eta < -100 || phi < 0) {
 560     GetRandomCell(eta, phi, absId);
 561   }
 562   else {
 563     fGeom->EtaPhiFromIndex(absId, eta, phi);
 564   }
 565
 566   if (absId == -1) {
 567     AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
 568                     " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
 569                     eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 570     return 0;
 571   }
 572
 573   if (e < 0) {
 574     Double_t pt = GetRandomPt();
 575     TLorentzVector nPart;
 576     nPart.SetPtEtaPhiM(pt, eta, phi, 0);
 577     e = nPart.E();
 578   }
 579
 580   return AddCluster(e, absId, label);
 581 }
 582
 583 //________________________________________________________________________
 584 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t label)
 585 {
 586   // Add a cluster to the event.
 587
 588   const Int_t nClusters = fOutClusters->GetEntriesFast();
 589
 590   TClonesArray digits("AliEMCALDigit", 1);
 591
 592   AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
 593   digit->SetId(absId);
 594   digit->SetIndexInList(0);
 595   digit->SetType(AliEMCALDigit::kHG);
 596   digit->SetAmplitude(e);
 597
 598   AliEMCALRecPoint recPoint("");
 599   recPoint.AddDigit(*digit, e, kFALSE);
 600   recPoint.EvalGlobalPosition(0, &digits);
 601
 602   TVector3 gpos;
 603   recPoint.GetGlobalPosition(gpos);
 604   Float_t g[3];
 605   gpos.GetXYZ(g);
 606
 607   AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
 608   cluster->SetType(AliVCluster::kEMCALClusterv1);
 609   cluster->SetE(recPoint.GetEnergy());
 610   cluster->SetPosition(g);
 611   cluster->SetNCells(1);
 612   UShort_t shortAbsId = absId;
 613   cluster->SetCellsAbsId(&shortAbsId);
 614   Double32_t fract = 1;
 615   cluster->SetCellsAmplitudeFraction(&fract);
 616   cluster->SetID(nClusters);
 617   cluster->SetEmcCpvDistance(-1);
 618
 619   //MC label
 620   if (label == 0)
 621     label = fMarkMC;
 622   else
 623     label += fMCLabelShift;
 624
 625   if (fEsdMode) {
 626     AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster);
 627     TArrayI parents(1, &label);
 628     esdClus->AddLabels(parents);
 629   }
 630   else {
 631     AliAODCaloCluster *aodClus = static_cast<AliAODCaloCluster*>(cluster);
 632     aodClus->SetLabel(&label, 1);
 633   }
 634
 635   return cluster;
 636 }
 637
 638 //________________________________________________________________________
 639 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)
 640 {
 641   // Add a track to the event.
 642
 643   const Int_t nTracks = fOutTracks->GetEntriesFast();
 644
 645   if (pt < 0)
 646     pt = GetRandomPt();
 647   if (eta < -100)
 648     eta = GetRandomEta();
 649   if (phi < 0)
 650     phi = GetRandomPhi();
 651
 652   if (label == 0)
 653     label = fMarkMC;
 654   else
 655     label += fMCLabelShift;
 656
 657   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt,
 658                                                                   eta,
 659                                                                   phi,
 660                                                                   1,
 661                                                                   label,
 662                                                                   type,
 663                                                                   etaemc,
 664                                                                   phiemc,
 665                                                                   ise);
 666
 667   return track;
 668 }
 669
 670 //________________________________________________________________________
 671 AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int_t origIndex)
 672 {
 673   const Int_t nPart = fOutMCParticles->GetEntriesFast();
 674
 675   AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
 676
 677   fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift);
 678   AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)",
 679                    origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex));
 680
 681   return aodpart;
 682 }
 683
 684 //________________________________________________________________________
 685 void AliJetModelBaseTask::CopyCells()
 686 {
 687   if (!fCaloCells)
 688     return;
 689
 690   fAddedCells = 0;
 691   fCaloCells->Sort();
 692   for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
 693     Int_t mclabel = 0;
 694     Double_t efrac = 0.;
 695     Double_t time = -1;
 696     Short_t cellNum = -1;
 697     Double_t amp = -1;
 698
 699     fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
 700
 701     if (!fIsMC)
 702       mclabel = 0;
 703
 704     fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
 705     fAddedCells++;
 706   }
 707
 708   AliDebug(2, Form("%d cells from the current event", fAddedCells));
 709 }
 710
 711 //________________________________________________________________________
 712 void AliJetModelBaseTask::CopyClusters()
 713 {
 714   // Copy all the clusters in the new collection
 715   if (!fClusters)
 716     return;
 717
 718   const Int_t nClusters = fClusters->GetEntriesFast();
 719   Int_t nCopiedClusters = 0;
 720
 721   if (fEsdMode) {
 722     for (Int_t i = 0; i < nClusters; ++i) {
 723       AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
 724       if (!esdcluster || !esdcluster->IsEMCAL())
 725         continue;
 726       AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
 727       if (!fIsMC) {
 728         TArrayI *labels = clus->GetLabelsArray();
 729         if (labels)
 730           labels->Reset();
 731       }
 732       nCopiedClusters++;
 733     }
 734   }
 735   else {
 736     for (Int_t i = 0; i < nClusters; ++i) {
 737       AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
 738       if (!aodcluster || !aodcluster->IsEMCAL())
 739         continue;
 740       AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
 741       if (!fIsMC)
 742         clus->SetLabel(0,0);
 743       nCopiedClusters++;
 744     }
 745   }
 746 }
 747
 748 //________________________________________________________________________
 749 void AliJetModelBaseTask::CopyTracks()
 750 {
 751   // Copy all the tracks in the new collection
 752
 753   if (!fTracks)
 754     return;
 755
 756   const Int_t nTracks = fTracks->GetEntriesFast();
 757   Int_t nCopiedTracks = 0;
 758   for (Int_t i = 0; i < nTracks; ++i) {
 759     AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
 760     if (!picotrack)
 761       continue;
 762     AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack);
 763     if (!fIsMC && track->GetLabel() != 0)
 764       track->SetLabel(0);
 765     nCopiedTracks++;
 766   }
 767 }
 768
 769 //________________________________________________________________________
 770 void AliJetModelBaseTask::CopyMCParticles()
 771 {
 772   // Copy all the MC particles in the new collection
 773
 774   if (!fMCParticles)
 775     return;
 776
 777   const Int_t nPart = fMCParticles->GetEntriesFast();
 778   Int_t nCopiedPart = 0;
 779   for (Int_t i = 0; i < nPart; ++i) {
 780     AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fMCParticles->At(i));
 781     if (!part)
 782       continue;
 783     new ((*fOutMCParticles)[nCopiedPart]) AliAODMCParticle(*part);
 784
 785     nCopiedPart++;
 786   }
 787
 788   if (!fMCParticlesMap)
 789     return;
 790
 791   for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) {
 792     fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i);
 793     if (fMCParticlesMap->At(i) >= 0)
 794       fMCLabelShift = i;
 795   }
 796
 797   AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift));
 798 }
 799
 800 //________________________________________________________________________
 801 void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
 802 {
 803   // Get random cell.
 804
 805   Int_t repeats = 0;
 806   Double_t rndEta = eta;
 807   Double_t rndPhi = phi;
 808   do {
 809     if (eta < -100)
 810       rndEta = GetRandomEta(kTRUE);
 811     if (phi < 0)
 812       rndPhi = GetRandomPhi(kTRUE);
 813     fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
 814     repeats++;
 815   } while (absId == -1 && repeats < 100);
 816
 817   if (!(absId > -1)) {
 818     AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
 819                     "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 820   }
 821   else {
 822     eta = rndEta;
 823     phi = rndPhi;
 824   }
 825 }
 826
 827 //________________________________________________________________________
 828 Double_t AliJetModelBaseTask::GetRandomEta(Bool_t emcal)
 829 {
 830   // Get random eta.
 831
 832   Double_t etamax = fEtaMax;
 833   Double_t etamin = fEtaMin;
 834
 835   if (emcal) {
 836     const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
 837     const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
 838
 839     if (etamax > EmcalMaxEta) etamax = EmcalMaxEta;
 840     if (etamax < EmcalMinEta) etamax = EmcalMinEta;
 841     if (etamin > EmcalMaxEta) etamin = EmcalMaxEta;
 842     if (etamin < EmcalMinEta) etamin = EmcalMinEta;
 843   }
 844
 845   return gRandom->Rndm() * (etamax - etamin) + etamin;
 846 }
 847
 848 //________________________________________________________________________
 849 Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal)
 850 {
 851   // Get random phi.
 852
 853   Double_t phimax = fPhiMax;
 854   Double_t phimin = fPhiMin;
 855
 856   if (emcal) {
 857     const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 858     const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 859
 860     if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
 861     if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
 862     if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
 863     if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
 864   }
 865
 866   return gRandom->Rndm() * (phimax - phimin) + phimin;
 867 }
 868
 869 //________________________________________________________________________
 870 Double_t AliJetModelBaseTask::GetRandomPt()
 871 {
 872   // Get random pt.
 873
 874   if (fPtSpectrum)
 875     return fPtSpectrum->GetRandom();
 876   else
 877     return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 878 }
 879
 880 //________________________________________________________________________
 881 void AliJetModelBaseTask::Run()
 882 {
 883   // Run.
 884 }