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, Double_t mass)
 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                                                                   mass);
 704
 705   return track;
 706 }
 707
 708 //________________________________________________________________________
 709 AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int_t origIndex)
 710 {
 711   const Int_t nPart = fOutMCParticles->GetEntriesFast();
 712
 713   AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
 714
 715   if (origIndex + fMCLabelShift >= fOutMCParticlesMap->GetSize())
 716     fOutMCParticlesMap->Set((origIndex + fMCLabelShift)*2);
 717
 718   fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift);
 719   AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)",
 720                    origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex));
 721
 722   return aodpart;
 723 }
 724
 725 //________________________________________________________________________
 726 void AliJetModelBaseTask::CopyCells()
 727 {
 728   if (!fCaloCells)
 729     return;
 730
 731   fAddedCells = 0;
 732   fCaloCells->Sort();
 733   for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
 734     Int_t mclabel = 0;
 735     Double_t efrac = 0.;
 736     Double_t time = -1;
 737     Short_t cellNum = -1;
 738     Double_t amp = -1;
 739
 740     fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
 741
 742     if (!fIsMC)
 743       mclabel = 0;
 744
 745     fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
 746     fAddedCells++;
 747   }
 748
 749   AliDebug(2, Form("%d cells from the current event", fAddedCells));
 750 }
 751
 752 //________________________________________________________________________
 753 void AliJetModelBaseTask::CopyClusters()
 754 {
 755   // Copy all the clusters in the new collection
 756   if (!fClusters)
 757     return;
 758
 759   const Int_t nClusters = fClusters->GetEntriesFast();
 760   Int_t nCopiedClusters = 0;
 761
 762   if (fEsdMode) {
 763     for (Int_t i = 0; i < nClusters; ++i) {
 764       AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
 765       if (!esdcluster || !esdcluster->IsEMCAL())
 766         continue;
 767       AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
 768       if (!fIsMC) {
 769         TArrayI *labels = clus->GetLabelsArray();
 770         if (labels)
 771           labels->Reset();
 772       }
 773       nCopiedClusters++;
 774     }
 775   }
 776   else {
 777     for (Int_t i = 0; i < nClusters; ++i) {
 778       AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
 779       if (!aodcluster || !aodcluster->IsEMCAL())
 780         continue;
 781       AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
 782       if (!fIsMC)
 783         clus->SetLabel(0,0);
 784       nCopiedClusters++;
 785     }
 786   }
 787 }
 788
 789 //________________________________________________________________________
 790 void AliJetModelBaseTask::CopyTracks()
 791 {
 792   // Copy all the tracks in the new collection
 793
 794   if (!fTracks)
 795     return;
 796
 797   const Int_t nTracks = fTracks->GetEntriesFast();
 798   Int_t nCopiedTracks = 0;
 799   for (Int_t i = 0; i < nTracks; ++i) {
 800     AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
 801     if (!picotrack)
 802       continue;
 803     AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack);
 804     if (!fIsMC && track->GetLabel() != 0)
 805       track->SetLabel(0);
 806     nCopiedTracks++;
 807   }
 808 }
 809
 810 //________________________________________________________________________
 811 void AliJetModelBaseTask::CopyMCParticles()
 812 {
 813   // Copy all the MC particles in the new collection
 814
 815   if (!fMCParticles)
 816     return;
 817
 818   const Int_t nPart = fMCParticles->GetEntriesFast();
 819   Int_t nCopiedPart = 0;
 820   for (Int_t i = 0; i < nPart; ++i) {
 821     AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fMCParticles->At(i));
 822     if (!part)
 823       continue;
 824     new ((*fOutMCParticles)[nCopiedPart]) AliAODMCParticle(*part);
 825
 826     nCopiedPart++;
 827   }
 828
 829   if (!fMCParticlesMap || !fOutMCParticlesMap)
 830     return;
 831
 832   if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize())
 833     fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2);
 834
 835   for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) {
 836     fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i);
 837     if (fMCParticlesMap->At(i) >= 0)
 838       fMCLabelShift = i;
 839   }
 840
 841   AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift));
 842 }
 843
 844 //________________________________________________________________________
 845 void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
 846 {
 847   // Get random cell.
 848
 849   Int_t repeats = 0;
 850   Double_t rndEta = eta;
 851   Double_t rndPhi = phi;
 852   do {
 853     if (eta < -100)
 854       rndEta = GetRandomEta(kTRUE);
 855     if (phi < 0)
 856       rndPhi = GetRandomPhi(kTRUE);
 857     fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
 858     repeats++;
 859   } while (absId == -1 && repeats < 100);
 860
 861   if (!(absId > -1)) {
 862     AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
 863                     "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 864   }
 865   else {
 866     eta = rndEta;
 867     phi = rndPhi;
 868   }
 869 }
 870
 871 //________________________________________________________________________
 872 Double_t AliJetModelBaseTask::GetRandomEta(Bool_t emcal)
 873 {
 874   // Get random eta.
 875
 876   Double_t etamax = fEtaMax;
 877   Double_t etamin = fEtaMin;
 878
 879   if (emcal) {
 880     const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
 881     const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
 882
 883     if (etamax > EmcalMaxEta) etamax = EmcalMaxEta;
 884     if (etamax < EmcalMinEta) etamax = EmcalMinEta;
 885     if (etamin > EmcalMaxEta) etamin = EmcalMaxEta;
 886     if (etamin < EmcalMinEta) etamin = EmcalMinEta;
 887   }
 888
 889   return gRandom->Rndm() * (etamax - etamin) + etamin;
 890 }
 891
 892 //________________________________________________________________________
 893 Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal)
 894 {
 895   // Get random phi.
 896
 897   Double_t phimax = fPhiMax;
 898   Double_t phimin = fPhiMin;
 899
 900   if (emcal) {
 901     const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 902     const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 903
 904     if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
 905     if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
 906     if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
 907     if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
 908   }
 909
 910   Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
 911
 912   return result;
 913 }
 914
 915 //________________________________________________________________________
 916 Double_t AliJetModelBaseTask::GetRandomPt()
 917 {
 918   // Get random pt.
 919
 920   if (fPtSpectrum)
 921     return fPtSpectrum->GetRandom();
 922   else
 923     return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 924 }
 925
 926 //________________________________________________________________________
 927 void AliJetModelBaseTask::GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal)
 928 {
 929   // Get a random particle.
 930
 931   eta = GetRandomEta(emcal);
 932
 933   if (fPtPhiEvPlDistribution) {
 934     Double_t phimax = fPhiMax;
 935     Double_t phimin = fPhiMin;
 936
 937     if (emcal) {
 938       const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 939       const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 940
 941       if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
 942       if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
 943       if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
 944       if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
 945     }
 946
 947     if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) {
 948       AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName()));
 949       pt = GetRandomPt();
 950       phi = GetRandomPhi(emcal);
 951     }
 952     else {
 953       do {
 954         fPtPhiEvPlDistribution->GetRandom2(pt,phi);
 955         phi += fPsi;
 956         if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2;
 957       } while (phi > phimax || phi < phimin);
 958     }
 959   }
 960   else {
 961     pt = GetRandomPt();
 962     phi = GetRandomPhi(emcal);
 963   }
 964 }
 965
 966 //________________________________________________________________________
 967 void AliJetModelBaseTask::Run()
 968 {
 969   // Run.
 970 }