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