PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx

   1 // $Id$
   2 //
   3 // General QA task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TChain.h>
   8 #include <TClonesArray.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TH3F.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliCentrality.h"
  17 #include "AliVCluster.h"
  18 #include "AliVParticle.h"
  19 #include "AliVTrack.h"
  20 #include "AliEmcalJet.h"
  21 #include "AliVEventHandler.h"
  22 #include "AliAODEvent.h"
  23 #include "AliExternalTrackParam.h"
  24 #include "AliTrackerBase.h"
  25 #include "AliLog.h"
  26
  27 #include "AliAnalysisTaskSAQA.h"
  28
  29 ClassImp(AliAnalysisTaskSAQA)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
  33   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
  34   fCellEnergyCut(0.1),
  35   fNclusters(0),
  36   fNtracks(0),
  37   fNjets(0),
  38   fHistTracksCent(0),
  39   fHistClusCent(0),
  40   fHistJetsCent(0),
  41   fHistClusTracks(0),
  42   fHistJetsParts(0),
  43   fHistCellsCent(0),
  44   fHistCellsTracks(0),
  45   fHistTrEmcPhiEta(0),
  46   fHistTrPhiEtaNonProp(0),
  47   fHistDeltaEtaPt(0),
  48   fHistDeltaPhiPt(0),
  49   fHistNCellsEnergy(0),
  50   fHistClusTimeEnergy(0),
  51   fHistCellsEnergy(0),
  52   fHistChVSneCells(0),
  53   fHistChVSneClus(0),
  54   fHistChVSneCorrCells(0)
  55 {
  56   // Default constructor.
  57
  58   for (Int_t i = 0; i < 5; i++) {
  59     fHistTrackPhi[i] = 0;
  60     fHistTrackEta[i] = 0;
  61   }
  62
  63   for (Int_t i = 0; i < 4; i++) {
  64     fHistTrPhiEtaPt[i] = 0;
  65     fHistClusPhiEtaEnergy[i] = 0;
  66     fHistJetsPhiEtaPt[i] = 0;
  67     fHistJetsPtArea[i] = 0;
  68   }
  69
  70   SetMakeGeneralHistograms(kTRUE);
  71 }
  72
  73 //________________________________________________________________________
  74 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
  75   AliAnalysisTaskEmcalJet(name, kTRUE),
  76   fCellEnergyCut(0.1),
  77   fNclusters(0),
  78   fNtracks(0),
  79   fNjets(0),
  80   fHistTracksCent(0),
  81   fHistClusCent(0),
  82   fHistJetsCent(0),
  83   fHistClusTracks(0),
  84   fHistJetsParts(0),
  85   fHistCellsCent(0),
  86   fHistCellsTracks(0),
  87   fHistTrEmcPhiEta(0),
  88   fHistTrPhiEtaNonProp(0),
  89   fHistDeltaEtaPt(0),
  90   fHistDeltaPhiPt(0),
  91   fHistNCellsEnergy(0),
  92   fHistClusTimeEnergy(0),
  93   fHistCellsEnergy(0),
  94   fHistChVSneCells(0),
  95   fHistChVSneClus(0),
  96   fHistChVSneCorrCells(0)
  97 {
  98   // Standard constructor.
  99
 100   for (Int_t i = 0; i < 5; i++) {
 101     fHistTrackPhi[i] = 0;
 102     fHistTrackEta[i] = 0;
 103   }
 104
 105   for (Int_t i = 0; i < 4; i++) {
 106     fHistTrPhiEtaPt[i] = 0;
 107     fHistClusPhiEtaEnergy[i] = 0;
 108     fHistJetsPhiEtaPt[i] = 0;
 109     fHistJetsPtArea[i] = 0;
 110   }
 111
 112   SetMakeGeneralHistograms(kTRUE);
 113 }
 114
 115 //________________________________________________________________________
 116 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
 117 {
 118   // Destructor
 119 }
 120
 121 //________________________________________________________________________
 122 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 123 {
 124   // Create histograms
 125
 126   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 127
 128   if (!fTracksName.IsNull()) {
 129     fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
 130     fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
 131     fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
 132     fOutput->Add(fHistTracksCent);
 133
 134     TString histname;
 135
 136     for (Int_t i = 0; i < 4; i++) {
 137       histname = "fHistTrPhiEtaPt_";
 138       histname += i;
 139       fHistTrPhiEtaPt[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
 140       fHistTrPhiEtaPt[i] ->GetXaxis()->SetTitle("#eta");
 141       fHistTrPhiEtaPt[i] ->GetYaxis()->SetTitle("#phi");
 142       fHistTrPhiEtaPt[i] ->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 143       fOutput->Add(fHistTrPhiEtaPt[i]);
 144     }
 145
 146     fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 147     fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
 148     fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
 149     fOutput->Add(fHistTrEmcPhiEta);
 150
 151     fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 152     fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
 153     fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
 154     fOutput->Add(fHistTrPhiEtaNonProp);
 155
 156     fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
 157     fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 158     fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
 159     fOutput->Add(fHistDeltaEtaPt);
 160
 161     fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
 162     fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 163     fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
 164     fOutput->Add(fHistDeltaPhiPt);
 165
 166     for (Int_t i = 0; i < 5; i++) {
 167       TString histnamephi("fHistTrackPhi_");
 168       histnamephi += i;
 169       fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 201, 0, TMath::Pi() * 2.01);
 170       fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
 171       fOutput->Add(fHistTrackPhi[i]);
 172
 173       TString histnameeta("fHistTrackEta_");
 174       histnameeta += i;
 175       fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -1, 1);
 176       fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
 177       fOutput->Add(fHistTrackEta[i]);
 178     }
 179
 180     fHistTrackPhi[0]->SetLineColor(kRed);
 181     fHistTrackEta[0]->SetLineColor(kRed);
 182     fHistTrackPhi[1]->SetLineColor(kBlue);
 183     fHistTrackEta[1]->SetLineColor(kBlue);
 184     fHistTrackPhi[2]->SetLineColor(kGreen);
 185     fHistTrackEta[2]->SetLineColor(kGreen);
 186     fHistTrackPhi[3]->SetLineColor(kOrange);
 187     fHistTrackEta[3]->SetLineColor(kOrange);
 188     fHistTrackPhi[4]->SetLineColor(kBlack);
 189     fHistTrackEta[4]->SetLineColor(kBlack);
 190   }
 191
 192   if (!fCaloName.IsNull()) {
 193     fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
 194     fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
 195     fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
 196     fOutput->Add(fHistClusCent);
 197
 198     fHistClusTracks = new TH2F("fHistClusTracks","Clusters vs. tracks", fNbins, 0, 4000, fNbins, 0, 2000);
 199     fHistClusTracks->GetXaxis()->SetTitle("No. of tracks");
 200     fHistClusTracks->GetYaxis()->SetTitle("No. of clusters");
 201     fOutput->Add(fHistClusTracks);
 202
 203     fHistCellsCent = new TH2F("fHistCellsCent","Cells vs. centrality", 100, 0, 100, fNbins, 0, 6000);
 204     fHistCellsCent->GetXaxis()->SetTitle("Centrality (%)");
 205     fHistCellsCent->GetYaxis()->SetTitle("No. of EMCal cells");
 206     fOutput->Add(fHistCellsCent);
 207
 208     fHistCellsTracks = new TH2F("fHistCellsTracks","Cells vs. tracks", fNbins, 0, 4000, fNbins, 0, 6000);
 209     fHistCellsTracks->GetXaxis()->SetTitle("No. of tracks");
 210     fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
 211     fOutput->Add(fHistCellsTracks);
 212
 213     TString histname;
 214
 215     for (Int_t i = 0; i < 4; i++) {
 216       histname = "fHistClusPhiEtaEnergy_";
 217       histname += i;
 218       fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
 219       fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
 220       fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
 221       fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("Energy (GeV)");
 222       fOutput->Add(fHistClusPhiEtaEnergy[i]);
 223     }
 224
 225     fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins,  -1e-6, 1e-6);
 226     fHistClusTimeEnergy->GetXaxis()->SetTitle("Energy (GeV)");
 227     fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
 228     fOutput->Add(fHistClusTimeEnergy);
 229
 230     fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
 231     fHistNCellsEnergy->GetXaxis()->SetTitle("Energy (GeV)");
 232     fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
 233     fOutput->Add(fHistNCellsEnergy);
 234
 235     fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
 236     fHistCellsEnergy->GetXaxis()->SetTitle("Energy (GeV)");
 237     fHistCellsEnergy->GetYaxis()->SetTitle("counts");
 238     fOutput->Add(fHistCellsEnergy);
 239
 240     fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy",
 241                                 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 242     fHistChVSneCells->GetXaxis()->SetTitle("Energy (GeV)");
 243     fHistChVSneCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
 244     fOutput->Add(fHistChVSneCells);
 245
 246     fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy",
 247                                (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 248     fHistChVSneClus->GetXaxis()->SetTitle("Energy (GeV)");
 249     fHistChVSneClus->GetYaxis()->SetTitle("Momentum (GeV/c)");
 250     fOutput->Add(fHistChVSneClus);
 251
 252     fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy",
 253                                     (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
 254     fHistChVSneCorrCells->GetXaxis()->SetTitle("Energy (GeV)");
 255     fHistChVSneCorrCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
 256     fOutput->Add(fHistChVSneCorrCells);
 257   }
 258
 259   if (!fJetsName.IsNull()) {
 260     fHistJetsCent = new TH2F("fHistJetsCent","Jets vs. centrality", 100, 0, 100, 150, -0.5, 149.5);
 261     fHistJetsCent->GetXaxis()->SetTitle("Centrality (%)");
 262     fHistJetsCent->GetYaxis()->SetTitle("No. of jets");
 263     fOutput->Add(fHistJetsCent);
 264
 265     fHistJetsParts = new TH2F("fHistJetsParts","Jets vs. centrality", fNbins, 0, 6000, 150, -0.5, 149.5);
 266     fHistJetsParts->GetXaxis()->SetTitle("No. of particles");
 267     fHistJetsParts->GetYaxis()->SetTitle("No. of jets");
 268     fOutput->Add(fHistJetsParts);
 269
 270     TString histname;
 271
 272     for (Int_t i = 0; i < 4; i++) {
 273       histname = "fHistJetsPhiEtaPt_";
 274       histname += i;
 275       fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01,  (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
 276       fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
 277       fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
 278       fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
 279       fOutput->Add(fHistJetsPhiEtaPt[i]);
 280
 281       histname = "fHistJetsPtArea_";
 282       histname += i;
 283       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
 284       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 285       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 286       fOutput->Add(fHistJetsPtArea[i]);
 287     }
 288   }
 289
 290   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 291 }
 292
 293 //________________________________________________________________________
 294 Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
 295 {
 296   // Retrieve event objects.
 297
 298   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 299     return kFALSE;
 300
 301   fNclusters = 0;
 302   fNtracks = 0;
 303   fNjets = 0;
 304
 305   return kTRUE;
 306 }
 307
 308
 309 //________________________________________________________________________
 310 Bool_t AliAnalysisTaskSAQA::FillHistograms()
 311 {
 312   // Fill histograms.
 313
 314   Float_t clusSum = 0;
 315   Float_t trackSum = 0;
 316
 317   if (fTracks) {
 318     trackSum = DoTrackLoop();
 319
 320     fHistTracksCent->Fill(fCent, fNtracks);
 321   }
 322
 323   if (fCaloClusters) {
 324     clusSum = DoClusterLoop();
 325
 326     fHistClusCent->Fill(fCent, fNclusters);
 327     fHistClusTracks->Fill(fNtracks, fNclusters);
 328
 329     Float_t cellSum = 0, cellCutSum = 0;
 330
 331     Int_t ncells = DoCellLoop(cellSum, cellCutSum);
 332
 333     if (fTracks)
 334       fHistCellsTracks->Fill(fNtracks, ncells);
 335
 336     fHistCellsCent->Fill(fCent, ncells);
 337
 338     fHistChVSneCells->Fill(cellSum, trackSum);
 339     fHistChVSneClus->Fill(clusSum, trackSum);
 340     fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
 341   }
 342
 343   if (fJets) {
 344     DoJetLoop();
 345
 346     fHistJetsCent->Fill(fCent, fNjets);
 347     fHistJetsParts->Fill(fNtracks + fNclusters, fNjets);
 348   }
 349
 350   return kTRUE;
 351 }
 352
 353 //________________________________________________________________________
 354 Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 355 {
 356   // Do cell loop.
 357
 358   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 359
 360   if (!cells)
 361     return 0;
 362
 363   const Int_t ncells = cells->GetNumberOfCells();
 364
 365   for (Int_t pos = 0; pos < ncells; pos++) {
 366     Float_t amp = cells->GetAmplitude(pos);
 367     fHistCellsEnergy->Fill(amp);
 368     sum += amp;
 369     if (amp < fCellEnergyCut)
 370       continue;
 371     sum_cut += amp;
 372   }
 373
 374   return ncells;
 375 }
 376
 377 //________________________________________________________________________
 378 Float_t AliAnalysisTaskSAQA::DoClusterLoop()
 379 {
 380   // Do cluster loop.
 381
 382   if (!fCaloClusters)
 383     return 0;
 384
 385   Float_t sum = 0;
 386
 387   // Cluster loop
 388   Int_t nclusters =  fCaloClusters->GetEntriesFast();
 389
 390   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 391     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
 392     if (!cluster) {
 393       AliError(Form("Could not receive cluster %d", iClusters));
 394       continue;
 395     }
 396
 397     if (!AcceptCluster(cluster, kTRUE))
 398       continue;
 399
 400     sum += cluster->E();
 401
 402     TLorentzVector nPart;
 403     cluster->GetMomentum(nPart, fVertex);
 404
 405     fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
 406     fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
 407
 408     fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
 409
 410     fNclusters++;
 411   }
 412
 413   return sum;
 414 }
 415
 416 //________________________________________________________________________
 417 Float_t AliAnalysisTaskSAQA::DoTrackLoop()
 418 {
 419   // Do track loop.
 420
 421   if (!fTracks)
 422     return 0;
 423
 424   Float_t sum = 0;
 425
 426   Int_t ntracks = fTracks->GetEntriesFast();
 427   Int_t nclusters = 0;
 428   if (fCaloClusters)
 429     nclusters = fCaloClusters->GetEntriesFast();
 430
 431   for (Int_t i = 0; i < ntracks; i++) {
 432
 433     AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
 434
 435     if (!track) {
 436       AliError(Form("Could not retrieve track %d",i));
 437       continue;
 438     }
 439
 440     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 441
 442     if (vtrack && !AcceptTrack(vtrack, kTRUE))
 443       continue;
 444
 445     fNtracks++;
 446
 447     sum += track->P();
 448
 449     fHistTrPhiEtaPt[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
 450
 451     Int_t label = track->GetLabel();
 452
 453     if (label >= 0 && label < 4) {
 454       fHistTrackEta[label]->Fill(track->Eta());
 455       fHistTrackPhi[label]->Fill(track->Phi());
 456     }
 457
 458     fHistTrackEta[4]->Fill(track->Eta());
 459     fHistTrackPhi[4]->Fill(track->Phi());
 460
 461     if (!vtrack)
 462       continue;
 463
 464     if (vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999)
 465       fHistTrPhiEtaNonProp->Fill(vtrack->Eta(), vtrack->Phi());
 466
 467     fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
 468     fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
 469     fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
 470   }
 471
 472   return sum;
 473 }
 474
 475 //________________________________________________________________________
 476 void AliAnalysisTaskSAQA::DoJetLoop()
 477 {
 478   // Do jet loop.
 479
 480   if (!fJets)
 481     return;
 482
 483   Int_t njets = fJets->GetEntriesFast();
 484
 485   for (Int_t ij = 0; ij < njets; ij++) {
 486
 487     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 488
 489     if (!jet) {
 490       AliError(Form("Could not receive jet %d", ij));
 491       continue;
 492     }
 493
 494     if (!AcceptJet(jet))
 495       continue;
 496
 497     fNjets++;
 498
 499     fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
 500     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 501   }
 502 }
 503
 504
 505 //________________________________________________________________________
 506 void AliAnalysisTaskSAQA::Terminate(Option_t *)
 507 {
 508   // Called once at the end of the analysis.
 509 }