PWGGA/EMCALJetTasks/AliAnalysisTaskSAQA.cxx

   1 // $Id: AliAnalysisTaskSAQA.cxx 56636 2012-05-22 22:51:12Z loizides $
   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
  14 #include "AliAnalysisManager.h"
  15 #include "AliCentrality.h"
  16 #include "AliVCluster.h"
  17 #include "AliVParticle.h"
  18 #include "AliVTrack.h"
  19 #include "AliEmcalJet.h"
  20 #include "AliVEventHandler.h"
  21 #include "AliLog.h"
  22
  23 #include "AliAnalysisTaskSAQA.h"
  24
  25 ClassImp(AliAnalysisTaskSAQA)
  26
  27 //________________________________________________________________________
  28 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
  29   AliAnalysisTaskEmcal("AliAnalysisTaskSAQA"),
  30   fCellEnergyCut(0.1),
  31   fDoEoverP(kFALSE),
  32   fTrgClusName("ClustersL1GAMMAFEE"),
  33   fTrgClusters(0),
  34   fHistCentrality(0),
  35   fHistTracksCent(0),
  36   fHistClusCent(0),
  37   fHistMaxL1FastORCent(0),
  38   fHistMaxL1ClusCent(0),
  39   fHistMaxL1ThrCent(0),
  40   fHistTracksPt(0),
  41   fHistTrPhiEta(0),
  42   fHistClustersEnergy(0),
  43   fHistClusPhiEta(0),
  44   fHistJetsPt(0),
  45   fHistJetsPhiEta(0),
  46   fHistJetsPtArea(0),
  47   fHistEoverP(0),
  48   fHistCellsEnergy(0),
  49   fHistChVSneCells(0),
  50   fHistChVSneClus(0),
  51   fHistChVSneCorrCells(0)
  52 {
  53   // Default constructor.
  54
  55   for (Int_t i = 0; i < 5; i++) {
  56     fHistTrackPhi[i] = 0;
  57     fHistTrackEta[i] = 0;
  58   }
  59 }
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
  63   AliAnalysisTaskEmcal(name),
  64   fCellEnergyCut(0.1),
  65   fDoEoverP(kFALSE),
  66   fTrgClusName("ClustersL1GAMMAFEE"),
  67   fTrgClusters(0),
  68   fHistCentrality(0),
  69   fHistTracksCent(0),
  70   fHistClusCent(0),
  71   fHistMaxL1FastORCent(0),
  72   fHistMaxL1ClusCent(0),
  73   fHistMaxL1ThrCent(0),
  74   fHistTracksPt(0),
  75   fHistTrPhiEta(0),
  76   fHistClustersEnergy(0),
  77   fHistClusPhiEta(0),
  78   fHistJetsPt(0),
  79   fHistJetsPhiEta(0),
  80   fHistJetsPtArea(0),
  81   fHistEoverP(0),
  82   fHistCellsEnergy(0),
  83   fHistChVSneCells(0),
  84   fHistChVSneClus(0),
  85   fHistChVSneCorrCells(0)
  86 {
  87   // Standard constructor.
  88
  89   for (Int_t i = 0; i < 5; i++) {
  90     fHistTrackPhi[i] = 0;
  91     fHistTrackEta[i] = 0;
  92   }
  93 }
  94
  95 //________________________________________________________________________
  96 AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
  97 {
  98   // Destructor
  99 }
 100
 101 //________________________________________________________________________
 102 void AliAnalysisTaskSAQA::UserCreateOutputObjects()
 103 {
 104   // Create histograms
 105
 106   OpenFile(1);
 107   fOutput = new TList();
 108   fOutput->SetOwner();  // IMPORTANT!
 109
 110   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
 111   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 112   fHistCentrality->GetYaxis()->SetTitle("counts");
 113   fOutput->Add(fHistCentrality);
 114
 115   fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", fNbins, 0, 100, fNbins, 0, 4000);
 116   fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
 117   fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
 118   fOutput->Add(fHistTracksCent);
 119
 120   fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", fNbins, 0, 100, fNbins, 0, 2000);
 121   fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
 122   fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
 123   fOutput->Add(fHistClusCent);
 124
 125   fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
 126   fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
 127   fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
 128   fOutput->Add(fHistMaxL1FastORCent);
 129
 130   fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
 131   fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
 132   fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
 133   fOutput->Add(fHistMaxL1ClusCent);
 134
 135   fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
 136   fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
 137   fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
 138   fOutput->Add(fHistMaxL1ThrCent);
 139
 140   fHistTracksPt = new TH1F("fHistTracksPt","P_{T} spectrum of reconstructed tracks", fNbins, fMinPt, fMaxPt);
 141   fHistTracksPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 142   fHistTracksPt->GetYaxis()->SetTitle("counts");
 143   fOutput->Add(fHistTracksPt);
 144
 145   fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 20, -2, 2, 32, 0, 6.4);
 146   fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
 147   fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
 148   fOutput->Add(fHistTrPhiEta);
 149
 150   fHistClustersEnergy = new TH1F("fHistClustersEnergy","Energy spectrum of clusters", fNbins, fMinPt, fMaxPt);
 151   fHistClustersEnergy->GetXaxis()->SetTitle("E [GeV]");
 152   fHistClustersEnergy->GetYaxis()->SetTitle("counts");
 153   fOutput->Add(fHistClustersEnergy);
 154
 155   fHistClusPhiEta = new TH2F("fHistClusPhiEta","Phi-Eta distribution of clusters", 20, -2, 2, 32, 0, 6.4);
 156   fHistClusPhiEta->GetXaxis()->SetTitle("#eta");
 157   fHistClusPhiEta->GetYaxis()->SetTitle("#phi");
 158   fOutput->Add(fHistClusPhiEta);
 159
 160   fHistJetsPt = new TH1F("fHistJetsPt","P_{T} spectrum of reconstructed jets", fNbins, fMinPt, fMaxPt);
 161   fHistJetsPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 162   fHistJetsPt->GetYaxis()->SetTitle("counts");
 163   fOutput->Add(fHistJetsPt);
 164
 165   fHistJetsPhiEta = new TH2F("fHistJetsPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
 166   fHistJetsPhiEta->GetXaxis()->SetTitle("#eta");
 167   fHistJetsPhiEta->GetYaxis()->SetTitle("#phi");
 168   fOutput->Add(fHistJetsPhiEta);
 169
 170   fHistJetsPtArea = new TH2F("fHistJetsPtArea","P_{T} vs. area of reconstructed jets", fNbins, fMinPt, fMaxPt, fNbins, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
 171   fHistJetsPtArea->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 172   fHistJetsPtArea->GetYaxis()->SetTitle("area");
 173   fOutput->Add(fHistJetsPtArea);
 174
 175   if (fDoEoverP) {
 176     fHistEoverP = new TH2F("fHistEoverP","E/P vs. E", fNbins, fMinPt, fMaxPt, fNbins, 0, 10);
 177     fHistEoverP->GetXaxis()->SetTitle("E [GeV]");
 178     fHistEoverP->GetYaxis()->SetTitle("E/P [c]");
 179     fOutput->Add(fHistEoverP);
 180   }
 181
 182   fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinPt, fMaxPt);
 183   fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
 184   fHistCellsEnergy->GetYaxis()->SetTitle("counts");
 185   fOutput->Add(fHistCellsEnergy);
 186
 187   fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy", fNbins, fMinPt * 4, fMaxPt * 4, fNbins, fMinPt * 4, fMaxPt * 4);
 188   fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
 189   fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
 190   fOutput->Add(fHistChVSneCells);
 191
 192   fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy", fNbins, fMinPt * 4, fMaxPt * 4, fNbins, fMinPt * 4, fMaxPt * 4);
 193   fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
 194   fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
 195   fOutput->Add(fHistChVSneClus);
 196
 197   fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy", fNbins, fMinPt * 4, fMaxPt * 4, fNbins, fMinPt * 4, fMaxPt * 4);
 198   fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
 199   fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
 200   fOutput->Add(fHistChVSneCorrCells);
 201
 202   for (Int_t i = 0; i < 5; i++) {
 203     TString histnamephi("fHistTrackPhi_");
 204     histnamephi += i;
 205     fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 128, 0, 6.4);
 206     fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
 207     fOutput->Add(fHistTrackPhi[i]);
 208
 209     TString histnameeta("fHistTrackEta_");
 210     histnameeta += i;
 211     fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
 212     fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
 213     fOutput->Add(fHistTrackEta[i]);
 214   }
 215
 216   fHistTrackPhi[0]->SetLineColor(kRed);
 217   fHistTrackEta[0]->SetLineColor(kRed);
 218   fHistTrackPhi[1]->SetLineColor(kBlue);
 219   fHistTrackEta[1]->SetLineColor(kBlue);
 220   fHistTrackPhi[2]->SetLineColor(kGreen);
 221   fHistTrackEta[2]->SetLineColor(kGreen);
 222   fHistTrackPhi[3]->SetLineColor(kOrange);
 223   fHistTrackEta[3]->SetLineColor(kOrange);
 224   fHistTrackPhi[4]->SetLineColor(kBlack);
 225   fHistTrackEta[4]->SetLineColor(kBlack);
 226
 227   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 228 }
 229
 230 //________________________________________________________________________
 231 void AliAnalysisTaskSAQA::RetrieveEventObjects()
 232 {
 233   AliAnalysisTaskEmcal::RetrieveEventObjects();
 234
 235   if (strcmp(fTrgClusName,"")) {
 236     fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
 237     if (!fTrgClusters) {
 238       AliWarning(Form("Could not retrieve trigger clusters!"));
 239     }
 240   }
 241 }
 242
 243 //________________________________________________________________________
 244 void AliAnalysisTaskSAQA::FillHistograms()
 245 {
 246   fHistCentrality->Fill(fCent);
 247   if (fTracks)
 248     fHistTracksCent->Fill(fCent, fTracks->GetEntriesFast());
 249   if (fCaloClusters)
 250     fHistClusCent->Fill(fCent, fCaloClusters->GetEntriesFast());
 251
 252   Float_t clusSum = DoClusterLoop();
 253
 254   Float_t trackSum = DoTrackLoop();
 255
 256   Float_t cellSum = 0, cellCutSum = 0;
 257   DoCellLoop(cellSum, cellCutSum);
 258
 259   fHistChVSneCells->Fill(cellSum, trackSum);
 260   fHistChVSneClus->Fill(clusSum, trackSum);
 261   fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
 262
 263   DoJetLoop();
 264
 265   Float_t maxTrgClus = DoTriggerClusLoop();
 266   fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
 267
 268   Int_t maxL1amp = -1;
 269   Int_t maxL1thr = -1;
 270
 271   DoTriggerPrimitives(maxL1amp, maxL1thr);
 272
 273   if (maxL1amp > -1)
 274     fHistMaxL1FastORCent->Fill(fCent, maxL1amp);
 275
 276   if (maxL1thr > -1)
 277     fHistMaxL1ThrCent->Fill(fCent, maxL1thr);
 278 }
 279
 280 //________________________________________________________________________
 281 void AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
 282 {
 283   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 284
 285   if (!cells)
 286     return;
 287
 288   Int_t ncells = cells->GetNumberOfCells();
 289
 290   for (Int_t pos = 0; pos < ncells; pos++) {
 291
 292     Float_t amp = cells->GetAmplitude(pos);
 293
 294     fHistCellsEnergy->Fill(amp);
 295
 296     sum += amp;
 297
 298     if (amp < fCellEnergyCut)
 299       continue;
 300
 301     sum_cut += amp;
 302
 303   }
 304 }
 305
 306 //________________________________________________________________________
 307 Float_t AliAnalysisTaskSAQA::DoClusterLoop()
 308 {
 309   if (!fCaloClusters)
 310     return 0;
 311
 312   Float_t sum = 0;
 313
 314   // Cluster loop
 315   Int_t nclusters =  fCaloClusters->GetEntriesFast();
 316
 317   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 318     AliVCluster* cluster = dynamic_cast<AliVCluster*>(fCaloClusters->At(iClusters));
 319     if (!cluster) {
 320       AliError(Form("Could not receive cluster %d", iClusters));
 321       continue;
 322     }
 323
 324     if (!(cluster->IsEMCAL())) continue;
 325
 326     if (!AcceptCluster(cluster, kTRUE)) continue;
 327
 328     fHistClustersEnergy->Fill(cluster->E());
 329
 330     sum += cluster->E();
 331
 332     Float_t pos[3];
 333     cluster->GetPosition(pos);
 334     TVector3 clusVec(pos);
 335     fHistClusPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
 336
 337   } //cluster loop
 338
 339   return sum;
 340 }
 341
 342 //________________________________________________________________________
 343 Float_t AliAnalysisTaskSAQA::DoTrackLoop()
 344 {
 345   if (!fTracks)
 346     return 0;
 347
 348   Float_t sum = 0;
 349
 350   // Track loop
 351   Int_t ntracks = fTracks->GetEntriesFast();
 352   Int_t nclusters = 0;
 353   if (fCaloClusters)
 354     nclusters = fCaloClusters->GetEntriesFast();
 355
 356   for(Int_t i = 0; i < ntracks; i++) {
 357
 358     AliVTrack* track = dynamic_cast<AliVTrack*>(fTracks->At(i)); // pointer to reconstructed to track
 359     if(!track) {
 360       AliError(Form("Could not retrieve track %d",i));
 361       continue;
 362     }
 363
 364     if (!AcceptTrack(track, kTRUE)) continue;
 365
 366     fHistTracksPt->Fill(track->Pt());
 367
 368     sum += track->P();
 369
 370     Int_t label = track->GetLabel();
 371
 372     fHistTrPhiEta->Fill(track->Eta(), track->Phi());
 373
 374     fHistTrackEta[4]->Fill(track->Eta());
 375     fHistTrackPhi[4]->Fill(track->Phi());
 376
 377     if (label >= 0 && label < 4) {
 378       fHistTrackEta[label]->Fill(track->Eta());
 379       fHistTrackPhi[label]->Fill(track->Phi());
 380     }
 381
 382     if (!fDoEoverP)
 383       continue;
 384
 385     Int_t clId = track->GetEMCALcluster();
 386     if (clId > -1 && clId < nclusters && fCaloClusters) {
 387       AliVCluster* cluster = dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
 388       if (cluster) {
 389         fHistEoverP->Fill(cluster->E(), cluster->E() / track->P());
 390       }
 391     }
 392   }
 393
 394   return sum;
 395 }
 396
 397 //________________________________________________________________________
 398 void AliAnalysisTaskSAQA::DoJetLoop()
 399 {
 400   if (!fJets)
 401     return;
 402
 403   Int_t njets = fJets->GetEntriesFast();
 404
 405   for (Int_t ij = 0; ij < njets; ij++) {
 406
 407     AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(ij));
 408
 409     if (!jet) {
 410       AliError(Form("Could not receive jet %d", ij));
 411       continue;
 412     }
 413
 414     if (!AcceptJet(jet))
 415       continue;
 416
 417     fHistJetsPt->Fill(jet->Pt());
 418     fHistJetsPtArea->Fill(jet->Pt(), jet->Area());
 419     fHistJetsPhiEta->Fill(jet->Eta(), jet->Phi());
 420   }
 421 }
 422
 423 //________________________________________________________________________
 424 Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
 425 {
 426   if (!fTrgClusters)
 427     return 0;
 428
 429   Int_t ntrgclusters = fTrgClusters->GetEntriesFast();
 430   Float_t maxTrgClus = 0;
 431
 432   for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
 433     AliVCluster* cluster = dynamic_cast<AliVCluster*>(fTrgClusters->At(iClusters));
 434     if (!cluster) {
 435       AliError(Form("Could not receive cluster %d", iClusters));
 436       continue;
 437     }
 438
 439     if (!(cluster->IsEMCAL())) continue;
 440
 441     if (cluster->E() > maxTrgClus)
 442       maxTrgClus = cluster->E();
 443
 444   }
 445   return maxTrgClus;
 446 }
 447
 448 //________________________________________________________________________
 449 void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
 450 {
 451   AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
 452
 453   if (!triggers || triggers->GetEntries() == 0)
 454     return;
 455
 456   triggers->Reset();
 457   Int_t L1amp = 0;
 458   Int_t L1thr = 0;
 459   maxL1amp = -1;
 460   maxL1thr = -1;
 461
 462   while (triggers->Next()) {
 463     triggers->GetL1TimeSum(L1amp);
 464     if (maxL1amp < L1amp)
 465       maxL1amp = L1amp;
 466
 467     triggers->GetL1Threshold(L1thr);
 468     if (maxL1thr < L1thr)
 469       maxL1thr = L1thr;
 470   }
 471 }
 472
 473 //________________________________________________________________________
 474 void AliAnalysisTaskSAQA::Terminate(Option_t *)
 475 {
 476   // Called once at the end of the analysis.
 477 }