PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJet.cxx

   1 // $Id$
   2 //
   3 // Emcal jet analysis base task.
   4 //
   5 // Author: S.Aiola, M. Verweij
   6
   7 #include "AliAnalysisTaskEmcalJet.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TList.h>
  11
  12 #include "AliAnalysisManager.h"
  13 #include "AliCentrality.h"
  14 #include "AliEMCALGeometry.h"
  15 #include "AliESDEvent.h"
  16 #include "AliEmcalJet.h"
  17 #include "AliLog.h"
  18 #include "AliRhoParameter.h"
  19 #include "AliLocalRhoParameter.h"
  20 #include "AliVCluster.h"
  21 #include "AliVEventHandler.h"
  22 #include "AliVParticle.h"
  23 #include "AliJetContainer.h"
  24
  25 ClassImp(AliAnalysisTaskEmcalJet)
  26
  27 //________________________________________________________________________
  28 AliAnalysisTaskEmcalJet::AliAnalysisTaskEmcalJet() :
  29   AliAnalysisTaskEmcal("AliAnalysisTaskEmcalJet"),
  30   fRhoName(),
  31   fLocalRhoName(),
  32   fJetCollArray(),
  33   fJets(0),
  34   fRho(0),
  35   fLocalRho(0),
  36   fRhoVal(0)
  37 {
  38   // Default constructor.
  39
  40   fJetCollArray.SetOwner(kTRUE);
  41 }
  42
  43 //________________________________________________________________________
  44 AliAnalysisTaskEmcalJet::AliAnalysisTaskEmcalJet(const char *name, Bool_t histo) :
  45   AliAnalysisTaskEmcal(name, histo),
  46   fRhoName(),
  47   fLocalRhoName(),
  48   fJetCollArray(),
  49   fJets(0),
  50   fRho(0),
  51   fLocalRho(0),
  52   fRhoVal(0)
  53 {
  54   // Standard constructor.
  55
  56   fJetCollArray.SetOwner(kTRUE);
  57 }
  58
  59 //________________________________________________________________________
  60 AliAnalysisTaskEmcalJet::~AliAnalysisTaskEmcalJet()
  61 {
  62   // Destructor
  63 }
  64
  65 //________________________________________________________________________
  66 Bool_t AliAnalysisTaskEmcalJet::AcceptBiasJet(AliEmcalJet *jet, Int_t c)
  67 {
  68   // Accept jet with a bias.
  69
  70   AliJetContainer *cont = GetJetContainer(c);
  71   if(!cont) {
  72     AliError(Form("%s:Container %d not found",GetName(),c));
  73     return 0;
  74   }
  75
  76   return cont->AcceptBiasJet(jet);
  77 }
  78
  79 //________________________________________________________________________
  80 Float_t* AliAnalysisTaskEmcalJet::GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const
  81 {
  82   Float_t *bins = new Float_t[n+1];
  83
  84   Float_t binWidth = (max-min)/n;
  85   bins[0] = min;
  86   for (Int_t i = 1; i <= n; i++) {
  87     bins[i] = bins[i-1]+binWidth;
  88   }
  89
  90   return bins;
  91 }
  92
  93 //________________________________________________________________________
  94 Double_t AliAnalysisTaskEmcalJet::GetLeadingHadronPt(AliEmcalJet *jet, Int_t c)
  95 {
  96
  97   AliJetContainer *cont = GetJetContainer(c);
  98   if(!cont) {
  99     AliError(Form("%s:Container %d not found",GetName(),c));
 100     return 0;
 101   }
 102
 103   return cont->GetLeadingHadronPt(jet);
 104 }
 105
 106 //________________________________________________________________________
 107 Bool_t AliAnalysisTaskEmcalJet::AcceptJet(AliEmcalJet *jet, Int_t c)
 108 {
 109   // Return true if jet is accepted.
 110   if (!jet)
 111     return kFALSE;
 112
 113   AliJetContainer *cont = GetJetContainer(c);
 114   if(!cont) {
 115     AliError(Form("%s:Container %d not found",GetName(),c));
 116     return 0;
 117   }
 118
 119   return cont->AcceptJet(jet);
 120 }
 121
 122 //________________________________________________________________________
 123 AliRhoParameter *AliAnalysisTaskEmcalJet::GetRhoFromEvent(const char *name)
 124 {
 125   // Get rho from event.
 126
 127   AliRhoParameter *rho = 0;
 128   TString sname(name);
 129   if (!sname.IsNull()) {
 130     rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
 131     if (!rho) {
 132       AliWarning(Form("%s: Could not retrieve rho with name %s!", GetName(), name));
 133       return 0;
 134     }
 135   }
 136   return rho;
 137 }
 138
 139 //________________________________________________________________________
 140 AliLocalRhoParameter *AliAnalysisTaskEmcalJet::GetLocalRhoFromEvent(const char *name)
 141 {
 142   // Get local rho from event.
 143   AliLocalRhoParameter *rho = 0;
 144   TString sname(name);
 145   if (!sname.IsNull()) {
 146     rho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(sname));
 147     if (!rho) {
 148       AliWarning(Form("%s: Could not retrieve local rho with name %s!", GetName(), name));
 149       return 0;
 150     }
 151   }
 152   return rho;
 153 }
 154
 155
 156 //________________________________________________________________________
 157 void AliAnalysisTaskEmcalJet::ExecOnce()
 158 {
 159   // Init the analysis.
 160
 161   AliAnalysisTaskEmcal::ExecOnce();
 162
 163   if (!fRhoName.IsNull() && !fRho) { // get rho from the event
 164     fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
 165     if (!fRho) {
 166       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
 167       fInitialized = kFALSE;
 168       return;
 169     }
 170   }
 171
 172   if (!fLocalRhoName.IsNull() && !fLocalRho) {
 173     fLocalRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fLocalRhoName));
 174     if (!fLocalRho) {
 175       AliError(Form("%s: Could not retrieve local rho %s!", GetName(), fLocalRhoName.Data()));
 176       fInitialized = kFALSE;
 177       return;
 178     }
 179   }
 180
 181   //Load all requested jet branches - each container knows name already
 182   if(fJetCollArray.GetEntriesFast()==0) {
 183     AliWarning("There are no jet collections");
 184     return;
 185   }
 186
 187   for(Int_t i =0; i<fJetCollArray.GetEntriesFast(); i++) {
 188     AliJetContainer *cont = static_cast<AliJetContainer*>(fJetCollArray.At(i));
 189     cont->SetRunNumber(InputEvent()->GetRunNumber());
 190     cont->SetEMCALGeometry();
 191     cont->SetArray(InputEvent());
 192     cont->LoadRho(InputEvent());
 193   }
 194
 195   //Get Jets, cuts and rho for first jet container
 196   AliJetContainer *cont = GetJetContainer(0);
 197
 198   if (!cont->GetArrayName().IsNull()) {
 199     fJets = cont->GetArray();
 200     if(!fJets && fJetCollArray.GetEntriesFast()>0) {
 201       AliError(Form("%s: Could not retrieve first jet branch!", GetName()));
 202       fInitialized = kFALSE;
 203       return;
 204     }
 205   }
 206
 207   if (!fRho) { // if rho name is not provided, tries to use the rho object of the first jet branch
 208     fRhoName = cont->GetRhoName();
 209     fRho = cont->GetRhoParameter();
 210   }
 211 }
 212
 213 //________________________________________________________________________
 214 Bool_t AliAnalysisTaskEmcalJet::IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted) const
 215 {
 216   // Return true if cluster is in jet.
 217
 218   for (Int_t i = 0; i < jet->GetNumberOfClusters(); ++i) {
 219     Int_t ijetclus = jet->ClusterAt(i);
 220     if (sorted && ijetclus > iclus)
 221       return kFALSE;
 222     if (ijetclus == iclus)
 223       return kTRUE;
 224   }
 225   return kFALSE;
 226 }
 227
 228 //________________________________________________________________________
 229 Bool_t AliAnalysisTaskEmcalJet::IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted) const
 230 {
 231   // Return true if track is in jet.
 232
 233   for (Int_t i = 0; i < jet->GetNumberOfTracks(); ++i) {
 234     Int_t ijettrack = jet->TrackAt(i);
 235     if (sorted && ijettrack > itrack)
 236       return kFALSE;
 237     if (ijettrack == itrack)
 238       return kTRUE;
 239   }
 240   return kFALSE;
 241 }
 242
 243 //________________________________________________________________________
 244 Bool_t AliAnalysisTaskEmcalJet::RetrieveEventObjects()
 245 {
 246   // Retrieve objects from event.
 247
 248   if (!AliAnalysisTaskEmcal::RetrieveEventObjects())
 249     return kFALSE;
 250
 251   if (fRho)
 252     fRhoVal = fRho->GetVal();
 253
 254   return kTRUE;
 255 }
 256
 257 //________________________________________________________________________
 258 AliJetContainer* AliAnalysisTaskEmcalJet::AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius) {
 259
 260   // Add particle container
 261   // will be called in AddTask macro
 262
 263   TString tmp = TString(n);
 264   if(tmp.IsNull()) return 0;
 265
 266   AliJetContainer *cont = 0x0;
 267   cont = new AliJetContainer();
 268   cont->SetArrayName(n);
 269   cont->SetJetRadius(jetRadius);
 270
 271   defaultCutType.ToUpper();
 272
 273   if(!defaultCutType.IsNull() && !defaultCutType.EqualTo("USER")) {
 274     if(defaultCutType.EqualTo("TPC"))
 275       cont->SetJetAcceptanceType(AliJetContainer::kTPC);
 276     else if(defaultCutType.EqualTo("EMCAL"))
 277       cont->SetJetAcceptanceType(AliJetContainer::kEMCAL);
 278     else
 279       AliWarning(Form("%s: default cut type %s not recognized. Not setting cuts.",GetName(),defaultCutType.Data()));
 280   } else
 281     cont->SetJetAcceptanceType(AliJetContainer::kUser);
 282
 283   fJetCollArray.Add(cont);
 284
 285   return cont;
 286 }
 287
 288 //________________________________________________________________________
 289 AliJetContainer* AliAnalysisTaskEmcalJet::GetJetContainer(Int_t i) const{
 290   // Get i^th jet container
 291
 292   if(i<0 || i>=fJetCollArray.GetEntriesFast()) return 0;
 293   AliJetContainer *cont = static_cast<AliJetContainer*>(fJetCollArray.At(i));
 294   return cont;
 295 }
 296
 297 //________________________________________________________________________
 298 AliJetContainer* AliAnalysisTaskEmcalJet::GetJetContainer(const char* name) const{
 299   // Get the jet container with name
 300
 301   AliJetContainer *cont = static_cast<AliJetContainer*>(fJetCollArray.FindObject(name));
 302   return cont;
 303 }
 304
 305 //________________________________________________________________________
 306 void AliAnalysisTaskEmcalJet::SetJetAcceptanceType(UInt_t t, Int_t c)
 307 {
 308   // Set acceptance cuts
 309   AliJetContainer *cont = GetJetContainer(c);
 310   if (cont) {
 311     cont->SetJetAcceptanceType((AliJetContainer::JetAcceptanceType)t);
 312   }
 313   else {
 314     AliError(Form("%s in SetJetAcceptanceType(...): container %d not found!",GetName(),c));
 315   }
 316 }
 317
 318 //________________________________________________________________________
 319 void AliAnalysisTaskEmcalJet::SetJetAcceptanceType(TString cutType, Int_t c) {
 320   //set acceptance cuts
 321
 322   AliJetContainer *cont = GetJetContainer(c);
 323   if (!cont) {
 324     AliError(Form("%s in SetJetAcceptanceType(...): container %d not found",GetName(),c));
 325     return;
 326   }
 327
 328   cutType.ToUpper();
 329
 330   if(!cutType.IsNull() && !cutType.EqualTo("USER")) {
 331     if(cutType.EqualTo("TPC"))
 332      cont->SetJetAcceptanceType(AliJetContainer::kTPC);
 333     else if(cutType.EqualTo("EMCAL"))
 334       cont->SetJetAcceptanceType(AliJetContainer::kEMCAL);
 335     else
 336       AliWarning(Form("%s: default cut type %s not recognized. Not setting cuts.",GetName(),cutType.Data()));
 337   } else
 338     cont->SetJetAcceptanceType(AliJetContainer::kUser);
 339 }
 340
 341 //________________________________________________________________________
 342 void AliAnalysisTaskEmcalJet::SetRhoName(const char *n, Int_t c)
 343 {
 344   if (c >= 0) {
 345     AliJetContainer *cont = GetJetContainer(c);
 346     if (cont) cont->SetRhoName(n);
 347     else AliError(Form("%s in SetRhoName(...): container %d not found",GetName(),c));
 348   }
 349   else {
 350     fRhoName = n;
 351   }
 352 }
 353
 354 //________________________________________________________________________
 355 void AliAnalysisTaskEmcalJet::SetJetEtaLimits(Float_t min, Float_t max, Int_t c)
 356 {
 357   AliJetContainer *cont = GetJetContainer(c);
 358   if (cont) cont->SetJetEtaLimits(min,max);
 359   else AliError(Form("%s in SetJetEtaLimits(...): container %d not found",GetName(),c));
 360 }
 361
 362 //________________________________________________________________________
 363 void AliAnalysisTaskEmcalJet::SetJetPhiLimits(Float_t min, Float_t max, Int_t c)
 364 {
 365   AliJetContainer *cont = GetJetContainer(c);
 366   if (cont) cont->SetJetPhiLimits(min,max);
 367   else AliError(Form("%s in SetJetPhiLimits(...): container %d not found",GetName(),c));
 368 }
 369
 370 //________________________________________________________________________
 371 void AliAnalysisTaskEmcalJet::SetJetAreaCut(Float_t cut, Int_t c)
 372 {
 373   AliJetContainer *cont = GetJetContainer(c);
 374   if (cont) cont->SetJetAreaCut(cut);
 375   else AliError(Form("%s in SetJetAreaCut(...): container %d not found",GetName(),c));
 376 }
 377
 378 //________________________________________________________________________
 379 void AliAnalysisTaskEmcalJet::SetPercAreaCut(Float_t p, Int_t c)
 380 {
 381   AliJetContainer *cont = GetJetContainer(c);
 382   if (cont) cont->SetPercAreaCut(p);
 383   else AliError(Form("%s in SetPercAreaCut(...): container %d not found",GetName(),c));
 384 }
 385
 386 //________________________________________________________________________
 387 void AliAnalysisTaskEmcalJet::SetZLeadingCut(Float_t zemc, Float_t zch, Int_t c)
 388 {
 389   AliJetContainer *cont = GetJetContainer(c);
 390   if (cont) cont->SetZLeadingCut(zemc,zch);
 391   else AliError(Form("%s in SetZLeadingCut(...): container %d not found",GetName(),c));
 392 }
 393
 394 //________________________________________________________________________
 395 void AliAnalysisTaskEmcalJet::SetNEFCut(Float_t min, Float_t max, Int_t c)
 396 {
 397   AliJetContainer *cont = GetJetContainer(c);
 398   if (cont) cont->SetNEFCut(min,max);
 399   else AliError(Form("%s in SetNEFCut(...): container %d not found",GetName(),c));
 400 }
 401
 402 //________________________________________________________________________
 403 void AliAnalysisTaskEmcalJet::SetAreaEmcCut(Double_t a, Int_t c)
 404 {
 405   AliJetContainer *cont = GetJetContainer(c);
 406   if (cont) cont->SetAreaEmcCut(a);
 407   else AliError(Form("%s in SetAreaEmcCut(...): container %d not found",GetName(),c));
 408 }
 409
 410 //________________________________________________________________________
 411 void AliAnalysisTaskEmcalJet::SetJetPtCut(Float_t cut, Int_t c)
 412 {
 413   AliJetContainer *cont = GetJetContainer(c);
 414   if (cont) cont->SetJetPtCut(cut);
 415   else AliError(Form("%s in SetJetPtCut(...): container %d not found",GetName(),c));
 416 }
 417
 418 //________________________________________________________________________
 419 void AliAnalysisTaskEmcalJet::SetJetRadius(Float_t r, Int_t c)
 420 {
 421   AliJetContainer *cont = GetJetContainer(c);
 422   if (cont) cont->SetJetRadius(r);
 423   else AliError(Form("%s in SetJetRadius(...): container %d not found",GetName(),c));
 424 }
 425
 426 //________________________________________________________________________
 427 void AliAnalysisTaskEmcalJet::SetMaxClusterPt(Float_t cut, Int_t c)
 428 {
 429   AliJetContainer *cont = GetJetContainer(c);
 430   if (cont) cont->SetMaxClusterPt(cut);
 431   else AliError(Form("%s in SetMaxClusterPt(...): container %d not found",GetName(),c));
 432 }
 433
 434 //________________________________________________________________________
 435 void AliAnalysisTaskEmcalJet::SetMaxTrackPt(Float_t cut, Int_t c)
 436 {
 437   AliJetContainer *cont = GetJetContainer(c);
 438   if (cont) cont->SetMaxTrackPt(cut);
 439   else AliError(Form("%s in SetMaxTrackPt(...): container %d not found",GetName(),c));
 440 }
 441
 442 //________________________________________________________________________
 443 void AliAnalysisTaskEmcalJet::SetPtBiasJetClus(Float_t cut, Int_t c)
 444 {
 445   AliJetContainer *cont = GetJetContainer(c);
 446   if (cont) cont->SetPtBiasJetClus(cut);
 447   else AliError(Form("%s in SetPtBiasJetClus(...): container %d not found",GetName(),c));
 448 }
 449
 450 //________________________________________________________________________
 451 void AliAnalysisTaskEmcalJet::SetPtBiasJetTrack(Float_t cut, Int_t c)
 452 {
 453   AliJetContainer *cont = GetJetContainer(c);
 454   if (cont) cont->SetPtBiasJetTrack(cut);
 455   else AliError(Form("%s in SetPtBiasJetTrack(...): container %d not found",GetName(),c));
 456 }
 457
 458 //________________________________________________________________________
 459 void AliAnalysisTaskEmcalJet::SetLeadingHadronType(Int_t t, Int_t c)
 460 {
 461   AliJetContainer *cont = GetJetContainer(c);
 462   if (cont) cont->SetLeadingHadronType(t);
 463   else AliError(Form("%s in SetLeadingHadronType(...): container %d not found",GetName(),c));
 464 }
 465
 466 //________________________________________________________________________
 467 void AliAnalysisTaskEmcalJet::SetNLeadingJets(Int_t t, Int_t c)
 468 {
 469   AliJetContainer *cont = GetJetContainer(c);
 470   if (cont) cont->SetNLeadingJets(t);
 471   else AliError(Form("%s in SetNLeadingJets(...): container %d not found",GetName(),c));
 472 }
 473
 474 //________________________________________________________________________
 475 void AliAnalysisTaskEmcalJet::SetJetBitMap(UInt_t m, Int_t c)
 476 {
 477   AliJetContainer *cont = GetJetContainer(c);
 478   if (cont) cont->SetJetBitMap(m);
 479   else AliError(Form("%s in SetJetBitMap(...): container %d not found",GetName(),c));
 480 }
 481
 482 //________________________________________________________________________
 483 void AliAnalysisTaskEmcalJet::SetIsParticleLevel(Bool_t b, Int_t c)
 484 {
 485   AliJetContainer *cont = GetJetContainer(c);
 486   if (cont) cont->SetIsParticleLevel(b);
 487   else AliError(Form("%s in SetIsParticleLevel(...): container %d not found",GetName(),c));
 488 }
 489
 490 //________________________________________________________________________
 491 const TString& AliAnalysisTaskEmcalJet::GetRhoName(Int_t c) const
 492 {
 493   if (c >= 0) {
 494     AliJetContainer *cont = GetJetContainer(c);
 495     if (cont) return cont->GetRhoName();
 496     else { AliError(Form("%s in GetRhoName(...): container %d not found. Returning fRhoName...",GetName(),c)); return fRhoName; }
 497   }
 498   else {
 499     return fRhoName;
 500   }
 501 }
 502
 503 //________________________________________________________________________
 504 TClonesArray* AliAnalysisTaskEmcalJet::GetJetArray(Int_t i) const {
 505   // Get i^th TClonesArray with AliEmcalJet
 506
 507   AliJetContainer *cont = GetJetContainer(i);
 508   if(!cont) {
 509     AliError(Form("%s:Container %d not found",GetName(),i));
 510     return 0;
 511   }
 512   return cont->GetArray();
 513 }
 514
 515 //________________________________________________________________________
 516 Double_t AliAnalysisTaskEmcalJet::GetJetRadius(Int_t i) const {
 517   // Get jet radius from jet container i
 518
 519   AliJetContainer *cont = GetJetContainer(i);
 520   if(!cont) {
 521     AliError(Form("%s:Container %d not found",GetName(),i));
 522     return 0;
 523   }
 524
 525   return cont->GetJetRadius();
 526 }
 527
 528 //________________________________________________________________________
 529 AliEmcalJet* AliAnalysisTaskEmcalJet::GetJetFromArray(Int_t j, Int_t c) const {
 530   // Get jet j if accepted from  container c
 531   // If jet not accepted return 0
 532
 533   AliJetContainer *cont = GetJetContainer(c);
 534   if(!cont) {
 535     AliError(Form("%s:Container %d not found",GetName(),c));
 536     return 0;
 537   }
 538   AliEmcalJet *jet = cont->GetJet(j);
 539
 540   return jet;
 541 }
 542
 543 //________________________________________________________________________
 544 AliEmcalJet* AliAnalysisTaskEmcalJet::GetAcceptJetFromArray(Int_t j, Int_t c) const {
 545   // Get jet j if accepted from  container c
 546   // If jet not accepted return 0
 547
 548   AliJetContainer *cont = GetJetContainer(c);
 549   if(!cont) {
 550     AliError(Form("%s:Container %d not found",GetName(),c));
 551     return 0;
 552   }
 553   AliEmcalJet *jet = cont->GetAcceptJet(j);
 554
 555   return jet;
 556 }
 557
 558 //________________________________________________________________________
 559 Int_t AliAnalysisTaskEmcalJet::GetNJets(Int_t i) const {
 560   // Get number of entries in jet array i
 561
 562   AliJetContainer *cont = GetJetContainer(i);
 563   if(!cont) {
 564     AliError(Form("%s:Container %d not found",GetName(),i));
 565     return 0;
 566   }
 567   return cont->GetNJets();
 568
 569 }
 570
 571 //________________________________________________________________________
 572 Double_t AliAnalysisTaskEmcalJet::GetRhoVal(Int_t i) const {
 573   // Get rho corresponding to jet array i
 574
 575   AliJetContainer *cont = GetJetContainer(i);
 576   if(!cont) {
 577     AliError(Form("%s:Container %d not found",GetName(),i));
 578     return 0;
 579   }
 580   return cont->GetRhoVal();
 581 }