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 <TH1.h>
  11 #include <TLorentzVector.h>
  12 #include <TRandom3.h>
  13 #include <TList.h>
  14
  15 #include "AliVEvent.h"
  16 #include "AliAODCaloCluster.h"
  17 #include "AliESDCaloCluster.h"
  18 #include "AliVCluster.h"
  19 #include "AliEMCALDigit.h"
  20 #include "AliEMCALRecPoint.h"
  21 #include "AliESDCaloCells.h"
  22 #include "AliAODCaloCells.h"
  23 #include "AliVCaloCells.h"
  24 #include "AliPicoTrack.h"
  25 #include "AliEMCALGeometry.h"
  26 #include "AliLog.h"
  27
  28 ClassImp(AliJetModelBaseTask)
  29
  30 //________________________________________________________________________
  31 AliJetModelBaseTask::AliJetModelBaseTask() :
  32   AliAnalysisTaskSE("AliJetModelBaseTask"),
  33   fGeomName(),
  34   fTracksName(),
  35   fOutTracksName(),
  36   fCaloName(),
  37   fOutCaloName(),
  38   fCellsName(),
  39   fOutCellsName(),
  40   fSuffix(),
  41   fEtaMin(-1),
  42   fEtaMax(1),
  43   fPhiMin(0),
  44   fPhiMax(TMath::Pi() * 2),
  45   fPtMin(0),
  46   fPtMax(0),
  47   fCopyArray(kTRUE),
  48   fNClusters(0),
  49   fNCells(0),
  50   fNTracks(0),
  51   fMarkMC(100),
  52   fPtSpectrum(0),
  53   fQAhistos(kFALSE),
  54   fIsInit(0),
  55   fGeom(0),
  56   fClusters(0),
  57   fOutClusters(0),
  58   fTracks(0),
  59   fOutTracks(0),
  60   fCaloCells(0),
  61   fOutCaloCells(0),
  62   fAddedCells(0),
  63   fEsdMode(kFALSE),
  64   fOutput(0)
  65 {
  66   // Default constructor.
  67 }
  68
  69 //________________________________________________________________________
  70 AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
  71   AliAnalysisTaskSE(name),
  72   fGeomName(""),
  73   fTracksName("PicoTracks"),
  74   fOutTracksName("PicoTracksEmbedded"),
  75   fCaloName("CaloClustersCorr"),
  76   fOutCaloName("CaloClustersCorrEmbedded"),
  77   fCellsName(""),
  78   fOutCellsName(""),
  79   fSuffix("Processed"),
  80   fEtaMin(-1),
  81   fEtaMax(1),
  82   fPhiMin(0),
  83   fPhiMax(TMath::Pi() * 2),
  84   fPtMin(50),
  85   fPtMax(60),
  86   fCopyArray(kTRUE),
  87   fNClusters(0),
  88   fNCells(0),
  89   fNTracks(1),
  90   fMarkMC(100),
  91   fPtSpectrum(0),
  92   fQAhistos(drawqa),
  93   fIsInit(0),
  94   fGeom(0),
  95   fClusters(0),
  96   fOutClusters(0),
  97   fTracks(0),
  98   fOutTracks(0),
  99   fCaloCells(0),
 100   fOutCaloCells(0),
 101   fAddedCells(0),
 102   fEsdMode(kFALSE),
 103   fOutput(0)
 104 {
 105   // Standard constructor.
 106
 107   if (fQAhistos) {
 108     DefineOutput(1, TList::Class());
 109   }
 110 }
 111
 112 //________________________________________________________________________
 113 AliJetModelBaseTask::~AliJetModelBaseTask()
 114 {
 115   // Destructor
 116 }
 117
 118 //________________________________________________________________________
 119 void AliJetModelBaseTask::UserCreateOutputObjects()
 120 {
 121   // Create user output.
 122   if (!fQAhistos)
 123     return;
 124
 125   OpenFile(1);
 126   fOutput = new TList();
 127   fOutput->SetOwner();
 128
 129   PostData(1, fOutput);
 130 }
 131
 132 //________________________________________________________________________
 133 void AliJetModelBaseTask::UserExec(Option_t *)
 134 {
 135   // Execute per event.
 136
 137   if (!fIsInit)
 138     fIsInit = ExecOnce();
 139
 140   if (!fIsInit)
 141     return;
 142
 143   if (fCopyArray) {
 144     if (fOutTracks)
 145       fOutTracks->Delete();
 146     if (fOutClusters)
 147       fOutClusters->Delete();
 148   }
 149
 150   AliVCaloCells *tempCaloCells = 0;
 151
 152   if (fCaloCells) {
 153     fAddedCells = 0;
 154     if (!fCopyArray) {
 155       tempCaloCells = fCaloCells;
 156       fCaloCells = static_cast<AliVCaloCells*>(tempCaloCells->Clone(Form("%s_old",fCaloCells->GetName())));
 157     }
 158   }
 159
 160   Run();
 161
 162   if (fCaloCells) {
 163     delete fCaloCells;
 164     fCaloCells = tempCaloCells;
 165   }
 166 }
 167
 168 //________________________________________________________________________
 169 Bool_t AliJetModelBaseTask::ExecOnce()
 170 {
 171   // Init task.
 172
 173   delete gRandom;
 174   gRandom = new TRandom3(0);
 175
 176   fEsdMode = InputEvent()->InheritsFrom("AliESDEvent");
 177
 178   if (fPtMax < fPtMin) {
 179     AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
 180     fPtMax = fPtMin;
 181   }
 182
 183   if (fEtaMax < fEtaMin) {
 184     AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
 185     fEtaMax = fEtaMin;
 186   }
 187
 188   if (fPhiMax < fPhiMin) {
 189     AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
 190     fPhiMax = fPhiMin;
 191   }
 192
 193   if (!fCaloCells && !fCellsName.IsNull()) {
 194     fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCellsName));
 195     if (!fCaloCells) {
 196       AliError(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data()));
 197       return kFALSE;
 198     }
 199
 200     if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
 201       AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data()));
 202       fCaloCells = 0;
 203       return kFALSE;
 204     }
 205
 206     if (!fOutCaloCells) {
 207       fOutCellsName = fCellsName;
 208       if (fCopyArray) {
 209         fOutCellsName += fSuffix;
 210
 211         if (fEsdMode)
 212           fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName);
 213         else
 214           fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName);
 215
 216         if (InputEvent()->FindListObject(fOutCellsName)) {
 217           AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data()));
 218           return kFALSE;
 219         }
 220         else {
 221           InputEvent()->AddObject(fOutCaloCells);
 222         }
 223       }
 224       else {
 225         fOutCaloCells = fCaloCells;
 226       }
 227     }
 228   }
 229
 230   if (fNTracks > 0 && !fTracks && !fTracksName.IsNull()) {
 231     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 232     if (!fTracks) {
 233       AliError(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
 234       return kFALSE;
 235     }
 236
 237     if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
 238       AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data()));
 239       fTracks = 0;
 240       return kFALSE;
 241     }
 242
 243     if (!fOutTracks) {
 244       fOutTracksName = fTracksName;
 245       if (fCopyArray) {
 246         fOutTracksName += fSuffix;
 247         fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize());
 248         fOutTracks->SetName(fOutTracksName);
 249         if (InputEvent()->FindListObject(fOutTracksName)) {
 250           AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data()));
 251           return kFALSE;
 252         }
 253         else {
 254           InputEvent()->AddObject(fOutTracks);
 255         }
 256       }
 257       else {
 258         fOutTracks = fTracks;
 259       }
 260     }
 261   }
 262
 263   if (fNClusters > 0 && !fClusters && !fCaloName.IsNull()) {
 264     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 265
 266     if (!fClusters) {
 267       AliError(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
 268       return kFALSE;
 269     }
 270
 271     if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
 272       AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data()));
 273       fClusters = 0;
 274       return kFALSE;
 275     }
 276
 277     if (!fOutClusters) {
 278       fOutCaloName = fCaloName;
 279       if (fCopyArray) {
 280         fOutCaloName += fSuffix;
 281         fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
 282         fOutClusters->SetName(fOutCaloName);
 283         if (InputEvent()->FindListObject(fOutCaloName)) {
 284           AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data()));
 285           return kFALSE;
 286         }
 287         else {
 288           InputEvent()->AddObject(fOutClusters);
 289         }
 290       }
 291       else {
 292         fOutClusters = fClusters;
 293       }
 294     }
 295   }
 296
 297   if (!fCaloName.IsNull() || !fCellsName.IsNull()) {
 298     if (!fGeom) {
 299       if (fGeomName.Length() > 0) {
 300         fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 301         if (!fGeom)
 302           AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
 303       } else {
 304         fGeom = AliEMCALGeometry::GetInstance();
 305         if (!fGeom)
 306           AliError("Could not get default geometry!");
 307       }
 308     }
 309
 310     const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
 311     const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
 312     const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
 313     const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
 314
 315     if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta;
 316     if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta;
 317     if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta;
 318     if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta;
 319
 320     if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi;
 321     if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi;
 322     if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi;
 323     if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi;
 324   }
 325
 326   return kTRUE;
 327 }
 328
 329 //________________________________________________________________________
 330 Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
 331 {
 332   if (fOutCaloCells->GetNumberOfCells() < n) {
 333     fOutCaloCells->DeleteContainer();
 334     fOutCaloCells->CreateContainer(n);
 335   }
 336   else {
 337     fOutCaloCells->SetNumberOfCells(n);
 338   }
 339
 340   fAddedCells = 0;
 341
 342   return n;
 343 }
 344
 345 //________________________________________________________________________
 346 void AliJetModelBaseTask::CopyCells()
 347 {
 348   for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
 349     Short_t mclabel = -1;
 350     Double_t efrac = 0.;
 351     Double_t time = -1;
 352     Short_t cellNum = -1;
 353     Double_t amp = -1;
 354
 355     fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
 356     fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
 357   }
 358
 359   fAddedCells = fCaloCells->GetNumberOfCells();
 360
 361   AliDebug(2, Form("%d cells from the PYTHIA event", fAddedCells));
 362 }
 363
 364 //________________________________________________________________________
 365 Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi)
 366 {
 367   // Add a cell to the event.
 368
 369   Int_t absId = 0;
 370   if (eta < -100 || phi < 0) {
 371     GetRandomCell(eta, phi, absId);
 372   }
 373   else {
 374     fGeom->EtaPhiFromIndex(absId, eta, phi);
 375   }
 376
 377   if (absId == -1) {
 378     AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!"
 379                     " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
 380                     eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 381     return 0;
 382   }
 383
 384   if (e < 0) {
 385     Double_t pt = GetRandomPt();
 386     TLorentzVector nPart;
 387     nPart.SetPtEtaPhiM(pt, eta, phi, 0);
 388     e = nPart.E();
 389   }
 390
 391   return AddCell(e, absId);
 392 }
 393
 394 //________________________________________________________________________
 395 Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time)
 396 {
 397   // Add a cell to the event.
 398
 399   Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, fMarkMC, 0);
 400
 401   if (r) {
 402     fAddedCells++;
 403     return fAddedCells;
 404   }
 405   else {
 406     return 0;
 407   }
 408 }
 409
 410
 411 //________________________________________________________________________
 412 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi)
 413 {
 414   // Add a cluster to the event.
 415
 416   Int_t absId = 0;
 417   if (eta < -100 || phi < 0) {
 418     GetRandomCell(eta, phi, absId);
 419   }
 420   else {
 421     fGeom->EtaPhiFromIndex(absId, eta, phi);
 422   }
 423
 424   if (absId == -1) {
 425     AliWarning(Form("Unable to embed cluster in eta = %f, phi = %f!"
 426                     " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])",
 427                     eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 428     return 0;
 429   }
 430
 431   if (e < 0) {
 432     Double_t pt = GetRandomPt();
 433     TLorentzVector nPart;
 434     nPart.SetPtEtaPhiM(pt, eta, phi, 0);
 435     e = nPart.E();
 436   }
 437
 438   return AddCluster(e, absId);
 439 }
 440
 441 //________________________________________________________________________
 442 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
 443 {
 444   // Add a cluster to the event.
 445
 446   const Int_t nClusters = fOutClusters->GetEntriesFast();
 447
 448   TClonesArray digits("AliEMCALDigit", 1);
 449
 450   AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(digits.New(0));
 451   digit->SetId(absId);
 452   digit->SetIndexInList(0);
 453   digit->SetType(AliEMCALDigit::kHG);
 454   digit->SetAmplitude(e);
 455
 456   AliEMCALRecPoint recPoint("");
 457   recPoint.AddDigit(*digit, e, kFALSE);
 458   recPoint.EvalGlobalPosition(0, &digits);
 459
 460   TVector3 gpos;
 461   recPoint.GetGlobalPosition(gpos);
 462   Float_t g[3];
 463   gpos.GetXYZ(g);
 464
 465   AliVCluster *cluster = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
 466   cluster->SetType(AliVCluster::kEMCALClusterv1);
 467   cluster->SetE(recPoint.GetEnergy());
 468   cluster->SetPosition(g);
 469   cluster->SetNCells(1);
 470   UShort_t shortAbsId = absId;
 471   cluster->SetCellsAbsId(&shortAbsId);
 472   Double32_t fract = 1;
 473   cluster->SetCellsAmplitudeFraction(&fract);
 474   cluster->SetID(nClusters);
 475   cluster->SetEmcCpvDistance(-1);
 476
 477   //MC label
 478   if (fEsdMode) {
 479     AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster);
 480     TArrayI parents(1, &fMarkMC);
 481     esdClus->AddLabels(parents);
 482   }
 483   else {
 484     AliAODCaloCluster *aodClus = static_cast<AliAODCaloCluster*>(cluster);
 485     aodClus->SetLabel(&fMarkMC, 1);
 486   }
 487
 488   return cluster;
 489 }
 490
 491 //________________________________________________________________________
 492 AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Bool_t ise)
 493 {
 494   // Add a track to the event.
 495
 496   const Int_t nTracks = fOutTracks->GetEntriesFast();
 497
 498   if (pt < 0)
 499     pt = GetRandomPt();
 500   if (eta < -100)
 501     eta = GetRandomEta();
 502   if (phi < 0)
 503     phi = GetRandomPhi();
 504
 505   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt,
 506                                                                   eta,
 507                                                                   phi,
 508                                                                   1,
 509                                                                   fMarkMC,    // MC flag!
 510                                                                   type,
 511                                                                   etaemc,
 512                                                                   phiemc,
 513                                                                   ise);
 514
 515   return track;
 516 }
 517
 518 //________________________________________________________________________
 519 void AliJetModelBaseTask::CopyClusters()
 520 {
 521   // Copy all the clusters in the new collection
 522
 523   const Int_t nClusters = fClusters->GetEntriesFast();
 524   Int_t nCopiedClusters = 0;
 525
 526   if (fEsdMode) {
 527     for (Int_t i = 0; i < nClusters; ++i) {
 528       AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
 529       if (!esdcluster || !esdcluster->IsEMCAL())
 530         continue;
 531       new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
 532       nCopiedClusters++;
 533     }
 534   }
 535   else {
 536     for (Int_t i = 0; i < nClusters; ++i) {
 537       AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
 538       if (!aodcluster || !aodcluster->IsEMCAL())
 539         continue;
 540       new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
 541       nCopiedClusters++;
 542     }
 543   }
 544 }
 545
 546 //________________________________________________________________________
 547 void AliJetModelBaseTask::CopyTracks()
 548 {
 549   const Int_t nTracks = fTracks->GetEntriesFast();
 550   Int_t nCopiedTracks = 0;
 551   for (Int_t i = 0; i < nTracks; ++i) {
 552     AliPicoTrack *track = static_cast<AliPicoTrack*>(fTracks->At(i));
 553     if (!track)
 554       continue;
 555     new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*track);
 556     nCopiedTracks++;
 557   }
 558 }
 559
 560 //________________________________________________________________________
 561 void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
 562 {
 563   // Get random cell.
 564
 565   Int_t repeats = 0;
 566   Double_t rndEta = eta;
 567   Double_t rndPhi = phi;
 568   do {
 569     if (eta < -100)
 570       rndEta = GetRandomEta();
 571     if (phi < 0)
 572       rndPhi = GetRandomPhi();
 573     fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);
 574     repeats++;
 575   } while (absId == -1 && repeats < 100);
 576
 577   if (!(absId > -1)) {
 578     AliWarning(Form("Could not extract random cluster! Random eta-phi extracted more than 100 times!\n"
 579                     "eta [%f, %f], phi [%f, %f]\n", fEtaMin, fEtaMax, fPhiMin, fPhiMax));
 580   }
 581   else {
 582     eta = rndEta;
 583     phi = rndPhi;
 584   }
 585 }
 586
 587 //________________________________________________________________________
 588 Double_t AliJetModelBaseTask::GetRandomEta()
 589 {
 590   // Get random eta.
 591
 592   return gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
 593 }
 594
 595 //________________________________________________________________________
 596 Double_t AliJetModelBaseTask::GetRandomPhi()
 597 {
 598   // Get random phi.
 599
 600   return gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
 601 }
 602
 603 //________________________________________________________________________
 604 Double_t AliJetModelBaseTask::GetRandomPt()
 605 {
 606   // Get random pt.
 607
 608   if (fPtSpectrum)
 609     return fPtSpectrum->GetRandom();
 610   else
 611     return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 612 }
 613
 614 //________________________________________________________________________
 615 void AliJetModelBaseTask::Run()
 616 {
 617   // Run.
 618 }