PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.cxx

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