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