PWGGA/EMCALTasks/AliAnalysisTaskSAQA.cxx

   1 // $Id$
   2 //
   3 // General QA task (S.Aiola).
   4 //
   5 //
   6
   7 #include <TChain.h>
   8 #include <TClonesArray.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TParticle.h>
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliCentrality.h"
  17 #include "AliVCluster.h"
  18 #include "AliESDtrack.h"
  19 #include "AliEmcalJet.h"
  20 #include "AliAODTrack.h"
  21 #include "AliESDtrack.h"
  22 #include "AliVEventHandler.h"
  23 #include "AliPicoTrack.h"
  24
  25 #include "AliAnalysisTaskSAQA.h"
  26
  27 ClassImp(AliAnalysisTaskSAQA)
  28
  29 //________________________________________________________________________
  30 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
  31   AliAnalysisTaskSE("AliAnalysisTaskSAQA"),
  32   fOutput(0),
  33   fCellEnergyCut(0.1),
  34   fTracksName("Tracks"),
  35   fCaloName("CaloClusters"),
  36   fTrgClusName("ClustersL1GAMMAFEE"),
  37   fTracks(0),
  38   fCaloClusters(0),
  39   fJets(0),
  40   fTrgClusters(0),
  41   fCent(0),
  42   fHistCentrality(0),
  43   fHistTracksCent(0),
  44   fHistClusCent(0),
  45   fHistMaxL1FastORCent(0),
  46   fHistMaxL1ClusCent(0),
  47   fHistMaxL1ThrCent(0),
  48   fHistTracksPt(0),
  49   fHistCellsEnergy(0),
  50   fHistClustersEnergy(0),
  51   fHistEoverP(0),
  52   fHistTrPhiEta(0),
  53   fHistClusPhiEta(0),
  54   fHistChVSneCells(0),
  55   fHistChVSneClus(0),
  56   fHistChVSneCorrCells(0),
  57   fNbins(250),
  58   fMinPt(0),
  59   fMaxPt(25)
  60 {
  61   // Default constructor.
  62
  63   for (Int_t i = 0; i < 5; i++) {
  64     fHistTrackPhi[i] = 0;
  65     fHistTrackEta[i] = 0;
  66   }
  67
  68   // Output slot #1 writes into a TH1 container
  69   DefineOutput(1, TList::Class());
  70 }
  71
  72 //________________________________________________________________________
  73 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
  74   AliAnalysisTaskSE(name),
  75   fOutput(0),
  76   fCellEnergyCut(0.1),
  77   fTracksName("Tracks"),
  78   fCaloName("CaloClusters"),
  79   fTrgClusName("ClustersL1GAMMAFEE"),
  80   fTracks(0),
  81   fCaloClusters(0),
  82   fJets(0),
  83   fTrgClusters(0),
  84   fCent(0),
  85   fHistCentrality(0),
  86   fHistTracksCent(0),
  87   fHistClusCent(0),
  88   fHistMaxL1FastORCent(0),
  89   fHistMaxL1ClusCent(0),
  90   fHistMaxL1ThrCent(0),
  91   fHistTracksPt(0),
  92   fHistCellsEnergy(0),
  93   fHistClustersEnergy(0),
  94   fHistEoverP(0),
  95   fHistTrPhiEta(0),
  96   fHistClusPhiEta(0),
  97   fHistChVSneCells(0),
  98   fHistChVSneClus(0),
  99   fHistChVSneCorrCells(0),
 100   fNbins(250),
 101   fMinPt(0),
 102   fMaxPt(25)
 103 {
 104   // Standard constructor.
 105
 106   for (Int_t i = 0; i < 5; i++) {
 107     fHistTrackPhi[i] = 0;
 108     fHistTrackEta[i] = 0;
 109   }
 110
 111   // Output slot #1 writes into a TH1 container
 112   DefineOutput(1, TList::Class());
 113 }
 114
 115 //________________________________________________________________________
 116 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
 117 {
 118   // Destructor
 119 }
 120
 121 //________________________________________________________________________
 122 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 123 {
 124   // Create histograms
 125
 126   fOutput = new TList();
 127   fOutput->SetOwner();  // IMPORTANT!
 128
 129   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
 130   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 131   fHistCentrality->GetYaxis()->SetTitle("counts");
 132   fOutput->Add(fHistCentrality);
 133
 134   fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", fNbins, 0, 100, fNbins, 0, 4000);
 135   fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
 136   fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
 137   fOutput->Add(fHistTracksCent);
 138
 139   fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", fNbins, 0, 100, fNbins, 0, 2000);
 140   fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
 141   fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
 142   fOutput->Add(fHistClusCent);
 143
 144   fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
 145   fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
 146   fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
 147   fOutput->Add(fHistMaxL1FastORCent);
 148
 149   fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
 150   fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
 151   fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
 152   fOutput->Add(fHistMaxL1ClusCent);
 153
 154   fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
 155   fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
 156   fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
 157   fOutput->Add(fHistMaxL1ThrCent);
 158
 159   fHistTracksPt = new TH1F("fHistTracksPt","P_{T} spectrum of reconstructed tracks", fNbins, fMinPt, fMaxPt);
 160   fHistTracksPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 161   fHistTracksPt->GetYaxis()->SetTitle("counts");
 162   fOutput->Add(fHistTracksPt);
 163
 164   fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinPt, fMaxPt);
 165   fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
 166   fHistCellsEnergy->GetYaxis()->SetTitle("counts");
 167   fOutput->Add(fHistCellsEnergy);
 168
 169   fHistClustersEnergy = new TH1F("fHistClustersEnergy","Energy spectrum of clusters", fNbins, fMinPt, fMaxPt);
 170   fHistClustersEnergy->GetXaxis()->SetTitle("E [GeV]");
 171   fHistClustersEnergy->GetYaxis()->SetTitle("counts");
 172   fOutput->Add(fHistClustersEnergy);
 173
 174   fHistEoverP = new TH2F("fHistEoverP","E/P vs. E", fNbins, fMinPt, fMaxPt, fNbins, 0, 10);
 175   fHistEoverP->GetXaxis()->SetTitle("E [GeV]");
 176   fHistEoverP->GetYaxis()->SetTitle("E/P [c]");
 177   fOutput->Add(fHistEoverP);
 178
 179   fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 20, -2, 2, 32, 0, 6.4);
 180   fHistTrPhiEta->GetXaxis()->SetTitle("Eta");
 181   fHistTrPhiEta->GetYaxis()->SetTitle("Phi");
 182   fOutput->Add(fHistTrPhiEta);
 183
 184   fHistClusPhiEta = new TH2F("fHistClusPhiEta","Phi-Eta distribution of clusters", 20, -2, 2, 32, 0, 6.4);
 185   fHistClusPhiEta->GetXaxis()->SetTitle("Eta");
 186   fHistClusPhiEta->GetYaxis()->SetTitle("Phi");
 187   fOutput->Add(fHistClusPhiEta);
 188
 189   fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy", fNbins, fMinPt * 10, fMaxPt * 10, fNbins, fMinPt * 10, fMaxPt * 10);
 190   fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
 191   fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
 192   fOutput->Add(fHistChVSneCells);
 193
 194   fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy", fNbins, fMinPt * 10, fMaxPt * 10, fNbins, fMinPt * 10, fMaxPt * 10);
 195   fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
 196   fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
 197   fOutput->Add(fHistChVSneClus);
 198
 199   fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy", fNbins, fMinPt * 10, fMaxPt * 10, fNbins, fMinPt * 10, fMaxPt * 10);
 200   fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
 201   fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
 202   fOutput->Add(fHistChVSneCorrCells);
 203
 204   for (Int_t i = 0; i < 5; i++) {
 205     TString histnamephi("fHistTrackPhi_");
 206     histnamephi += i;
 207     fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 128, 0, 6.4);
 208     fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
 209     fOutput->Add(fHistTrackPhi[i]);
 210
 211     TString histnameeta("fHistTrackEta_");
 212     histnameeta += i;
 213     fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
 214     fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
 215     fOutput->Add(fHistTrackEta[i]);
 216   }
 217
 218   fHistTrackPhi[0]->SetLineColor(kRed);
 219   fHistTrackEta[0]->SetLineColor(kRed);
 220   fHistTrackPhi[1]->SetLineColor(kBlue);
 221   fHistTrackEta[1]->SetLineColor(kBlue);
 222   fHistTrackPhi[2]->SetLineColor(kGreen);
 223   fHistTrackEta[2]->SetLineColor(kGreen);
 224   fHistTrackPhi[3]->SetLineColor(kOrange);
 225   fHistTrackEta[3]->SetLineColor(kOrange);
 226   fHistTrackPhi[4]->SetLineColor(kBlack);
 227   fHistTrackEta[4]->SetLineColor(kBlack);
 228
 229   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 230 }
 231
 232 void AliAnalysisTaskSAQA::RetrieveEventObjects()
 233 {
 234   fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 235   if (!fCaloClusters) {
 236     AliWarning(Form("Could not retrieve clusters!"));
 237   }
 238
 239   fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 240   if (!fTracks) {
 241     AliWarning(Form("Could not retrieve tracks!"));
 242   }
 243
 244   if (strcmp(fTrgClusName,"")) {
 245     fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
 246     if (!fTrgClusters) {
 247       AliWarning(Form("Could not retrieve trigger clusters!"));
 248     }
 249   }
 250
 251   fCent = InputEvent()->GetCentrality();
 252 }
 253
 254 AliVTrack* AliAnalysisTaskSAQA::GetTrack(const Int_t i) const
 255 {
 256   if (fTracks)
 257     return dynamic_cast<AliVTrack*>(fTracks->At(i));
 258   else
 259     return 0;
 260 }
 261
 262 Int_t AliAnalysisTaskSAQA::GetNumberOfTracks() const
 263 {
 264   if (fTracks)
 265     return fTracks->GetEntriesFast();
 266   else
 267     return 0;
 268 }
 269
 270 AliVCluster* AliAnalysisTaskSAQA::GetCaloCluster(const Int_t i) const
 271 {
 272   if (fCaloClusters)
 273     return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
 274   else
 275     return 0;
 276 }
 277
 278 Int_t AliAnalysisTaskSAQA::GetNumberOfCaloClusters() const
 279 {
 280   if (fCaloClusters)
 281     return fCaloClusters->GetEntriesFast();
 282   else
 283     return 0;
 284 }
 285
 286 AliVCluster* AliAnalysisTaskSAQA::GetTrgCluster(const Int_t i) const
 287 {
 288   if (fTrgClusters)
 289     return dynamic_cast<AliVCluster*>(fTrgClusters->At(i));
 290   else
 291     return 0;
 292 }
 293
 294 Int_t AliAnalysisTaskSAQA::GetNumberOfTrgClusters() const
 295 {
 296   if (fTrgClusters)
 297     return fTrgClusters->GetEntriesFast();
 298   else
 299     return 0;
 300 }
 301
 302 //________________________________________________________________________
 303 Float_t AliAnalysisTaskSAQA::GetAcceptanceNormFactor() const
 304 {
 305   const Float_t EmcalEtaMin = -0.7;
 306   const Float_t EmcalEtaMax = 0.7;
 307   const Float_t EmcalPhiMin = 80 * TMath::DegToRad();
 308   const Float_t EmcalPhiMax = 180 * TMath::DegToRad();
 309
 310   const Float_t TpcEtaMin = -0.9;
 311   const Float_t TpcEtaMax = 0.9;
 312   const Float_t TpcPhiMin = 0;
 313   const Float_t TpcPhiMax = 2 * TMath::Pi();
 314
 315   Float_t emcalArea = (EmcalEtaMax - EmcalEtaMin) * (EmcalPhiMax - EmcalPhiMin);
 316   Float_t tpcArea   = (TpcEtaMax - TpcEtaMin) * (TpcPhiMax - TpcPhiMin);
 317
 318   return emcalArea / tpcArea;
 319 }
 320
 321 void AliAnalysisTaskSAQA::FillHistograms()
 322 {
 323   Float_t cent = 100;
 324
 325   if (fCent)
 326     cent = fCent->GetCentralityPercentile("V0M");
 327   else
 328     AliWarning("Centrality not available!");
 329
 330   fHistCentrality->Fill(cent);
 331   fHistTracksCent->Fill(cent, GetNumberOfTracks());
 332   fHistClusCent->Fill(cent, GetNumberOfCaloClusters());
 333
 334   Float_t clusSum = DoClusterLoop();
 335
 336   Float_t trackSum = DoTrackLoop();
 337
 338   //Normalization to EMCal acceptance
 339   trackSum *= GetAcceptanceNormFactor();
 340
 341   Float_t cellSum = 0, cellCutSum = 0;
 342   DoCellLoop(cellSum, cellCutSum);
 343
 344   fHistChVSneCells->Fill(cellSum, trackSum);
 345   fHistChVSneClus->Fill(clusSum, trackSum);
 346   fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
 347
 348   Float_t maxTrgClus = DoTriggerClusLoop();
 349   fHistMaxL1ClusCent->Fill(cent, maxTrgClus);
 350
 351   Int_t maxL1amp = -1;
 352   Int_t maxL1thr = -1;
 353
 354   DoTriggerPrimitives(maxL1amp, maxL1thr);
 355
 356   if (maxL1amp > -1)
 357     fHistMaxL1FastORCent->Fill(cent, maxL1amp);
 358
 359   if (maxL1thr > -1)
 360     fHistMaxL1ThrCent->Fill(cent, maxL1thr);
 361 }
 362
 363 //________________________________________________________________________
 364 void AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 365 {
 366   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 367
 368   if (!cells)
 369     return;
 370
 371   Int_t ncells = cells->GetNumberOfCells();
 372
 373   for (Int_t pos = 0; pos < ncells; pos++) {
 374
 375     Float_t amp = cells->GetAmplitude(pos);
 376
 377     fHistCellsEnergy->Fill(amp);
 378
 379     sum += amp;
 380
 381     if (amp < fCellEnergyCut)
 382       continue;
 383
 384     sum_cut += amp;
 385
 386   }
 387 }
 388
 389 //________________________________________________________________________
 390 Float_t AliAnalysisTaskSAQA::DoClusterLoop()
 391 {
 392   Float_t sum = 0;
 393
 394   // get primary vertex
 395   //Double_t vertex[3] = {0, 0, 0};
 396   //InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 397
 398   // Cluster loop
 399   Int_t nclusters =  GetNumberOfCaloClusters();
 400
 401   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 402     AliVCluster* cluster = GetCaloCluster(iClusters);
 403     if (!cluster) {
 404       AliError(Form("Could not receive cluster %d", iClusters));
 405       continue;
 406     }
 407
 408     if (!(cluster->IsEMCAL())) continue;
 409
 410     fHistClustersEnergy->Fill(cluster->E());
 411
 412     //TLorentzVector nPart;
 413     //cluster->GetMomentum(nPart, vertex);
 414     sum += cluster->E();
 415
 416     Float_t pos[3];
 417     cluster->GetPosition(pos);
 418     TVector3 clusVec(pos);
 419     fHistClusPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
 420
 421   } //cluster loop
 422
 423   return sum;
 424 }
 425
 426 //________________________________________________________________________
 427 Float_t AliAnalysisTaskSAQA::DoTrackLoop()
 428 {
 429   Float_t sum = 0;
 430
 431   // Track loop
 432   Int_t nclusters =  GetNumberOfCaloClusters();
 433   Int_t ntracks = GetNumberOfTracks();
 434
 435   for(Int_t i = 0; i < ntracks; i++) {
 436
 437     AliVTrack* track = GetTrack(i); // pointer to reconstructed to track
 438     if(!track) {
 439       AliError(Form("Could not retrieve esdtrack %d",i));
 440       continue;
 441     }
 442
 443     if (!AcceptTrack(track)) continue;
 444
 445     fHistTracksPt->Fill(track->Pt());
 446
 447     sum += track->P();
 448
 449     Int_t clId = track->GetEMCALcluster();
 450     if (clId > -1 && clId < nclusters) {
 451       AliVCluster* cluster = GetCaloCluster(clId);
 452       if (cluster) {
 453         fHistEoverP->Fill(cluster->E(), cluster->E() / track->P());
 454       }
 455     }
 456
 457     Int_t label = track->GetLabel();
 458
 459     fHistTrPhiEta->Fill(track->Eta(), track->Phi());
 460
 461     fHistTrackEta[4]->Fill(track->Eta());
 462     fHistTrackPhi[4]->Fill(track->Phi());
 463
 464     if (label >= 0 && label < 4) {
 465       fHistTrackEta[label]->Fill(track->Eta());
 466       fHistTrackPhi[label]->Fill(track->Phi());
 467     }
 468   }
 469
 470   return sum;
 471 }
 472
 473 //________________________________________________________________________
 474 Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
 475 {
 476   Int_t ntrgclusters =  GetNumberOfTrgClusters();
 477   Float_t maxTrgClus = 0;
 478
 479   for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
 480     AliVCluster* cluster = GetTrgCluster(iClusters);
 481     if (!cluster) {
 482       AliError(Form("Could not receive cluster %d", iClusters));
 483       continue;
 484     }
 485
 486     if (!(cluster->IsEMCAL())) continue;
 487
 488     if (cluster->E() > maxTrgClus)
 489       maxTrgClus = cluster->E();
 490
 491   }
 492   return maxTrgClus;
 493 }
 494
 495 //________________________________________________________________________
 496 void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
 497 {
 498   AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
 499
 500   if (!triggers || triggers->GetEntries() == 0)
 501     return;
 502
 503   triggers->Reset();
 504   //Float_t L0FastORamp = 0;
 505   Int_t L1amp = 0;
 506   Int_t L1thr = 0;
 507   maxL1amp = -1;
 508   maxL1thr = -1;
 509
 510   while (triggers->Next()) {
 511     /*
 512     triggers->GetAmplitude(L0FastORamp);
 513
 514     if (L0FastORamp < 0)
 515       continue;
 516
 517     Int_t ntimes = 0;
 518
 519     triggers->GetNL0Times(ntimes);
 520
 521     if (ntimes > 0) {
 522       Int_t trgtimes[25];
 523       triggers->GetL0Times(trgtimes);
 524
 525       Int_t mintime = trgtimes[0];
 526       Int_t maxtime = trgtimes[0];
 527
 528       for (Int_t i = 0; i < ntimes; ++i) {
 529         if (trgtimes[i] < mintime)
 530           mintime = trgtimes[i];
 531         if (maxtime < trgtimes[i])
 532           maxtime = trgtimes[i];
 533       }
 534     }
 535
 536     Int_t gCol = 0, gRow = 0;
 537     triggers->GetPosition(gCol, gRow);
 538
 539     Int_t find = -1;
 540     fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol, gRow, find);
 541
 542     if (find < 0)
 543       continue;
 544
 545     Int_t cidx[4] = {-1};
 546     Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
 547
 548     if (!ret)
 549       continue;
 550     */
 551
 552     triggers->GetL1TimeSum(L1amp);
 553     if (maxL1amp < L1amp)
 554       maxL1amp = L1amp;
 555
 556     triggers->GetL1Threshold(L1thr);
 557     if (maxL1thr < L1thr)
 558       maxL1thr = L1thr;
 559   }
 560 }
 561
 562 //________________________________________________________________________
 563 Bool_t AliAnalysisTaskSAQA::AcceptTrack(AliVTrack* /*track*/) const
 564 {
 565   return kTRUE;
 566 }
 567
 568 //________________________________________________________________________
 569 void AliAnalysisTaskSAQA::UserExec(Option_t *)
 570 {
 571   // Main loop, called for each event.
 572   // Add jets to event if not yet there
 573
 574   RetrieveEventObjects();
 575
 576   FillHistograms();
 577
 578   // information for this iteration of the UserExec in the container
 579   PostData(1, fOutput);
 580 }
 581
 582 //________________________________________________________________________
 583 void AliAnalysisTaskSAQA::Terminate(Option_t *)
 584 {
 585   // Called once at the end of the analysis.
 586 }