PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTriggerQA.cxx

   1 //
   2 // Jet trigger QA analysis task.
   3 //
   4 // Author: M.Verweij
   5
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <THnSparse.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13
  14 #include "AliVCluster.h"
  15 #include "AliVTrack.h"
  16 #include "AliEmcalJet.h"
  17 #include "AliRhoParameter.h"
  18 #include "AliLog.h"
  19 #include "AliAnalysisUtils.h"
  20 #include "AliEmcalParticle.h"
  21 #include "AliAODCaloTrigger.h"
  22 #include "AliEMCALGeometry.h"
  23 #include "AliVCaloCells.h"
  24 #include "AliJetContainer.h"
  25 #include "AliClusterContainer.h"
  26
  27 #include "AliAnalysisTaskEmcalJetTriggerQA.h"
  28
  29 ClassImp(AliAnalysisTaskEmcalJetTriggerQA)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() :
  33   AliAnalysisTaskEmcalJetDev("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
  34   fDebug(kFALSE),
  35   fUseAnaUtils(kTRUE),
  36   fAnalysisUtils(0),
  37   fTriggerClass(""),
  38   fBitJ1((1<<8)),
  39   fBitJ2((1<<11)),
  40   fContainerFull(0),
  41   fContainerCharged(1),
  42   fMaxPatchEnergy(0),
  43   fTriggerType(-1),
  44   fNFastOR(16),
  45   fhNEvents(0),
  46   fh3PtEtaPhiJetFull(0),
  47   fh3PtEtaPhiJetCharged(0),
  48   fh2NJetsPtFull(0),
  49   fh2NJetsPtCharged(0),
  50   fh3PtEtaAreaJetFull(0),
  51   fh3PtEtaAreaJetCharged(0),
  52   fh2PtNConstituentsCharged(0),
  53   fh2PtNConstituents(0),
  54   fh2PtMeanPtConstituentsCharged(0),
  55   fh2PtMeanPtConstituentsNeutral(0),
  56   fh2PtNEF(0),
  57   fh3NEFEtaPhi(0),
  58   fh2NEFNConstituentsCharged(0),
  59   fh2NEFNConstituentsNeutral(0),
  60   fh2Ptz(0),
  61   fh2PtzCharged(0),
  62   fh2PtLeadJet1VsLeadJet2(0),
  63   fh3EEtaPhiCluster(0),
  64   fh3PtLeadJet1VsPatchEnergy(0),
  65   fh3PtLeadJet2VsPatchEnergy(0),
  66   fh3PatchEnergyEtaPhiCenter(0),
  67   fh2CellEnergyVsTime(0)
  68 {
  69   // Default constructor.
  70
  71
  72   SetMakeGeneralHistograms(kTRUE);
  73 }
  74
  75 //________________________________________________________________________
  76 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
  77   AliAnalysisTaskEmcalJetDev(name, kTRUE),
  78   fDebug(kFALSE),
  79   fUseAnaUtils(kTRUE),
  80   fAnalysisUtils(0),
  81   fTriggerClass(""),
  82   fBitJ1((1<<8)),
  83   fBitJ2((1<<11)),
  84   fContainerFull(0),
  85   fContainerCharged(1),
  86   fMaxPatchEnergy(0),
  87   fTriggerType(-1),
  88   fNFastOR(16),
  89   fhNEvents(0),
  90   fh3PtEtaPhiJetFull(0),
  91   fh3PtEtaPhiJetCharged(0),
  92   fh2NJetsPtFull(0),
  93   fh2NJetsPtCharged(0),
  94   fh3PtEtaAreaJetFull(0),
  95   fh3PtEtaAreaJetCharged(0),
  96   fh2PtNConstituentsCharged(0),
  97   fh2PtNConstituents(0),
  98   fh2PtMeanPtConstituentsCharged(0),
  99   fh2PtMeanPtConstituentsNeutral(0),
 100   fh2PtNEF(0),
 101   fh3NEFEtaPhi(0),
 102   fh2NEFNConstituentsCharged(0),
 103   fh2NEFNConstituentsNeutral(0),
 104   fh2Ptz(0),
 105   fh2PtzCharged(0),
 106   fh2PtLeadJet1VsLeadJet2(0),
 107   fh3EEtaPhiCluster(0),
 108   fh3PtLeadJet1VsPatchEnergy(0),
 109   fh3PtLeadJet2VsPatchEnergy(0),
 110   fh3PatchEnergyEtaPhiCenter(0),
 111   fh2CellEnergyVsTime(0)
 112 {
 113   // Standard constructor.
 114
 115   SetMakeGeneralHistograms(kTRUE);
 116 }
 117
 118 //________________________________________________________________________
 119 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
 120 {
 121   // Destructor.
 122 }
 123
 124 //----------------------------------------------------------------------------------------------
 125 void AliAnalysisTaskEmcalJetTriggerQA::InitOnce() {
 126   //
 127   // only initialize once
 128   //
 129
 130   // Initialize analysis util class for vertex selection
 131   if(fUseAnaUtils) {
 132     fAnalysisUtils = new AliAnalysisUtils();
 133     fAnalysisUtils->SetMinVtxContr(2);
 134     fAnalysisUtils->SetMaxVtxZ(10.);
 135   }
 136 }
 137
 138 //________________________________________________________________________
 139 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
 140   //
 141   // Decide if event should be selected for analysis
 142   //
 143
 144   fhNEvents->Fill(0.5);
 145
 146   if(fAnalysisUtils) {
 147     if(!fAnalysisUtils->IsVertexSelected2013pA(InputEvent()))
 148       return kFALSE;
 149
 150     fhNEvents->Fill(2.5);
 151
 152     if(fAnalysisUtils->IsPileUpEvent(InputEvent()))
 153       return kFALSE;
 154   }
 155   else{
 156     if(fUseAnaUtils)
 157       AliError("fAnalysisUtils not initialized. Call AliAnalysisTaskEmcalJetTriggerQA::InitOnce()");
 158   }
 159
 160   fhNEvents->Fill(3.5);
 161
 162   if(!fTriggerClass.IsNull()) {
 163     //Check if requested trigger was fired
 164     TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
 165     if(!firedTrigClass.Contains(fTriggerClass))
 166       return kFALSE;
 167   }
 168
 169   fhNEvents->Fill(1.5);
 170
 171   return kTRUE;
 172
 173 }
 174
 175 //________________________________________________________________________
 176 void AliAnalysisTaskEmcalJetTriggerQA::FindTriggerPatch() {
 177
 178   //Code to get position of trigger
 179
 180   if(!fGeom)
 181     fGeom = AliEMCALGeometry::GetInstance();
 182
 183
 184   TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
 185
 186   AliAODCaloTrigger *trg = dynamic_cast<AliAODCaloTrigger*>(InputEvent()->GetCaloTrigger("EMCAL"));
 187   trg->Reset();
 188
 189   int col, row; //FASTOR position
 190
 191   Int_t nPatchNotEmpty = 0; //counter number of patches which are not empty
 192   fMaxPatchEnergy = 0.;
 193   while (trg->Next())
 194     {
 195
 196       trg->GetPosition(col, row); //col (0 to 63), row (0 to 47)
 197
 198       if (col > -1 && row > -1)
 199         {
 200           Int_t id = -1; //FASTOR index
 201           Int_t cellIndex[4] = {-1};
 202           fGeom->GetAbsFastORIndexFromPositionInEMCAL(col,row,id); //phi is column, eta is row
 203           if(id<0) continue;
 204
 205           fGeom->GetCellIndexFromFastORIndex(id,cellIndex);
 206
 207           Int_t ts = 0;
 208           trg->GetL1TimeSum(ts);
 209
 210           Float_t ampTrg = 0.;
 211           trg->GetAmplitude(ampTrg);
 212
 213           //L1 analysis
 214           Int_t bit = 0;
 215           trg->GetTriggerBits(bit);
 216
 217           Bool_t bTrigJ = TestFilterBit(bit,fBitJ1|fBitJ2);
 218           if(!bTrigJ)
 219             continue;
 220
 221           Bool_t bTrigJ1 = TestFilterBit(bit,fBitJ1);
 222           Bool_t bTrigJ2 = TestFilterBit(bit,fBitJ2);
 223
 224           if(bTrigJ1) fTriggerType = 0;
 225           if(bTrigJ2) fTriggerType = 1;
 226           if(!bTrigJ1 && !bTrigJ2) fTriggerType = -1;
 227
 228           Double_t minPhiPatch =  10.;
 229           Double_t maxPhiPatch = -10.;
 230           Double_t minEtaPatch =  10.;
 231           Double_t maxEtaPatch = -10.;
 232
 233           //Get energy in trigger patch 8x8 FASTOR
 234           Double_t patchEnergy = 0.;
 235           Double_t sumAmp = 0.;
 236           //      const Int_t nFastOR = 8;//16;
 237           for(Int_t fastrow = 0; fastrow<fNFastOR; fastrow++) {
 238             for(Int_t fastcol = 0; fastcol<fNFastOR; fastcol++) {
 239               Int_t nrow = row+fastrow;
 240               Int_t ncol = col+fastcol;
 241               fGeom->GetAbsFastORIndexFromPositionInEMCAL(ncol,nrow,id);
 242
 243               if(id<0) {
 244                 AliWarning(Form("%s: id smaller than 0 %d",GetName(),id));
 245                 continue;
 246               }
 247
 248               fGeom->GetCellIndexFromFastORIndex(id,cellIndex);
 249               for(Int_t icell=0; icell<4; icell++) {
 250
 251                 if(!fGeom->CheckAbsCellId(cellIndex[icell])) continue;
 252
 253                 Double_t amp =0., time = 0., efrac = 0;
 254                 Int_t mclabel = -1;
 255                 Short_t absId = -1;
 256                 Int_t nSupMod = -1, nModule = -1, nIphi = -1, nIeta = -1;
 257
 258                 fGeom->GetCellIndex(cellIndex[icell], nSupMod, nModule, nIphi, nIeta );
 259
 260                 fCaloCells->GetCell(cellIndex[icell], absId, amp, time,mclabel,efrac);
 261
 262                 Double_t eta, phi;
 263                 fGeom->EtaPhiFromIndex(cellIndex[icell], eta, phi);
 264
 265                 if(phi<minPhiPatch) minPhiPatch = phi;
 266                 if(phi>maxPhiPatch) maxPhiPatch = phi;
 267                 if(eta<minEtaPatch) minEtaPatch = eta;
 268                 if(eta>maxEtaPatch) maxEtaPatch = eta;
 269
 270                 sumAmp+=fCaloCells->GetAmplitude(cellIndex[icell]);
 271                 patchEnergy+=amp;
 272
 273               }//cells in fastor loop
 274             }//fastor col loop
 275           }//fastor row loop
 276
 277           Double_t etaCent = minEtaPatch + 0.5*(maxEtaPatch-minEtaPatch);
 278           Double_t phiCent = minPhiPatch + 0.5*(maxPhiPatch-minPhiPatch);
 279           fh3PatchEnergyEtaPhiCenter->Fill(patchEnergy,etaCent,phiCent);
 280
 281           if(patchEnergy>0)
 282             if(fDebug>2) AliInfo(Form("%s: patch edges eta: %f - %f   phi: %f - %f ampTrg: %f patchEnergy: %f sumAmp: %f",GetName(),minEtaPatch,maxEtaPatch,minPhiPatch,maxPhiPatch,ampTrg,patchEnergy,sumAmp));
 283
 284           if(patchEnergy>0) nPatchNotEmpty++;
 285
 286           if(patchEnergy>fMaxPatchEnergy) fMaxPatchEnergy=patchEnergy;
 287
 288         }
 289     }
 290
 291 }
 292
 293 //________________________________________________________________________
 294 void AliAnalysisTaskEmcalJetTriggerQA::LoadExtraBranches() {
 295   //
 296   // load extra brances
 297   //
 298
 299 }
 300
 301 //________________________________________________________________________
 302 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
 303 {
 304   // Create user output.
 305
 306   InitOnce();
 307
 308   AliAnalysisTaskEmcalJetDev::UserCreateOutputObjects();
 309
 310   Bool_t oldStatus = TH1::AddDirectoryStatus();
 311   TH1::AddDirectory(kFALSE);
 312
 313   fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
 314   fOutput->Add(fhNEvents);
 315
 316   Int_t nBinsPt = 120;
 317   Double_t minPt = -20.;
 318   Double_t maxPt = 100.;
 319   Int_t nBinsEta = 40;
 320   Double_t minEta = -1.;
 321   Double_t maxEta = 1.;
 322   Int_t nBinsPhi = 18*6;
 323   Double_t minPhi = 0.;
 324   Double_t maxPhi = TMath::TwoPi();
 325   Int_t nBinsArea = 100;
 326   Double_t minArea = 0.;
 327   Double_t maxArea = 1.;
 328
 329   fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 330   fOutput->Add(fh3PtEtaPhiJetFull);
 331
 332   fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 333   fOutput->Add(fh3PtEtaPhiJetCharged);
 334
 335   fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,nBinsPt,minPt,maxPt);
 336   fOutput->Add(fh2NJetsPtFull);
 337
 338   fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,nBinsPt,minPt,maxPt);
 339   fOutput->Add(fh2NJetsPtCharged);
 340
 341   fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsArea,minArea,maxArea);
 342   fOutput->Add(fh3PtEtaAreaJetFull);
 343
 344   fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsArea,minArea,maxArea);
 345   fOutput->Add(fh3PtEtaAreaJetCharged);
 346
 347   Int_t nBinsConst =100;
 348   Double_t minConst = 0.;
 349   Double_t maxConst = 100.;
 350
 351   Int_t nBinsMeanPt = 200;
 352   Double_t minMeanPt = 0.;
 353   Double_t maxMeanPt = 20.;
 354
 355   Int_t nBinsNEF = 101;
 356   Double_t minNEF = 0.;
 357   Double_t maxNEF = 1.01;
 358
 359   Int_t nBinsz = 101;
 360   Double_t minz = 0.;
 361   Double_t maxz = 1.01;
 362
 363   Int_t nBinsECluster =100;
 364   Double_t minECluster = 0.;
 365   Double_t maxECluster = 100.;
 366
 367
 368   fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",nBinsPt,minPt,maxPt,nBinsConst,minConst,maxConst);
 369   fOutput->Add(fh2PtNConstituentsCharged);
 370
 371   fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",nBinsPt,minPt,maxPt,nBinsConst,minConst,maxConst);
 372   fOutput->Add(fh2PtNConstituents);
 373
 374   fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",nBinsPt,minPt,maxPt,nBinsMeanPt,minMeanPt,maxMeanPt);
 375   fOutput->Add(fh2PtMeanPtConstituentsCharged);
 376
 377   fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",nBinsPt,minPt,maxPt,nBinsMeanPt,minMeanPt,maxMeanPt);
 378   fOutput->Add(fh2PtMeanPtConstituentsNeutral);
 379
 380   fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",nBinsPt,minPt,maxPt,nBinsNEF,minNEF,maxNEF);
 381   fOutput->Add(fh2PtNEF);
 382
 383   fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",nBinsNEF,minNEF,maxNEF,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 384   fOutput->Add(fh3NEFEtaPhi);
 385
 386   fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",nBinsNEF,minNEF,maxNEF,nBinsConst,minConst,maxConst);
 387   fOutput->Add(fh2NEFNConstituentsCharged);
 388
 389   fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",nBinsNEF,minNEF,maxNEF,nBinsConst,minConst,maxConst);
 390   fOutput->Add(fh2NEFNConstituentsNeutral);
 391
 392   fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",nBinsPt,minPt,maxPt,nBinsz,minz,maxz);
 393   fOutput->Add(fh2Ptz);
 394
 395   fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",nBinsPt,minPt,maxPt,nBinsz,minz,maxz);
 396   fOutput->Add(fh2PtzCharged);
 397
 398   fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 399   fOutput->Add(fh2PtLeadJet1VsLeadJet2);
 400
 401   fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",nBinsECluster,minECluster,maxECluster,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 402   fOutput->Add(fh3EEtaPhiCluster);
 403
 404   fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,2,-0.5,1.5);
 405   fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
 406   fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,2,-0.5,1.5);
 407   fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
 408
 409   fh3PatchEnergyEtaPhiCenter = new TH3F("fh3PatchEnergyEtaPhiCenter","fh3PatchEnergyEtaPhiCenter;E_{patch};#eta;#phi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 410   fOutput->Add(fh3PatchEnergyEtaPhiCenter);
 411
 412   fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",100,0.,100.,700,-400,1000);
 413   fOutput->Add(fh2CellEnergyVsTime);
 414
 415
 416   // =========== Switch on Sumw2 for all histos ===========
 417   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 418     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 419     if (h1){
 420       h1->Sumw2();
 421       continue;
 422     }
 423     TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
 424     if (h2){
 425       h2->Sumw2();
 426       continue;
 427     }
 428     TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
 429     if (h3){
 430       h3->Sumw2();
 431       continue;
 432     }
 433     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 434     if(hn)hn->Sumw2();
 435   }
 436
 437   TH1::AddDirectory(oldStatus);
 438
 439   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 440 }
 441
 442 //________________________________________________________________________
 443 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
 444 {
 445   // Fill histograms.
 446
 447   AliClusterContainer  *clusCont = GetClusterContainer(0);
 448
 449   if (clusCont) {
 450     Int_t nclusters = clusCont->GetNClusters();
 451     for (Int_t ic = 0; ic < nclusters; ic++) {
 452       AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
 453       if (!cluster) {
 454         AliError(Form("Could not receive cluster %d", ic));
 455         continue;
 456       }
 457       if (!cluster->IsEMCAL()) {
 458         AliDebug(11,Form("%s: Cluster is not emcal",GetName()));
 459         continue;
 460       }
 461
 462       TLorentzVector lp;
 463       cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 464
 465       //Fill eta,phi,E of clusters here
 466       fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
 467     }
 468   }
 469
 470   if(fCaloCells) {
 471     const Short_t nCells   = fCaloCells->GetNumberOfCells();
 472
 473     for(Int_t iCell=0; iCell<nCells; ++iCell) {
 474       Short_t cellId = fCaloCells->GetCellNumber(iCell);
 475       Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
 476       Double_t cellT = fCaloCells->GetCellTime(cellId);
 477
 478       AliDebug(2,Form("cell energy = %f  time = %f",cellE,cellT*1e9));
 479       fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
 480
 481     }
 482
 483   }
 484
 485   Double_t ptLeadJet1 = 0.;
 486   Double_t ptLeadJet2 = 0.;
 487
 488   TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
 489   nJetsArr->Reset(0);
 490   nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
 491
 492   if (GetJetContainer(fContainerFull)) {
 493     const Int_t njets = GetNJets(fContainerFull);
 494     for (Int_t ij = 0; ij < njets; ij++) {
 495
 496       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
 497       if (!jet)
 498         continue; //jet not selected
 499
 500       Double_t jetPt = jet->Pt();
 501       if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
 502       fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
 503       fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
 504
 505       //count jets above certain pT threshold
 506       Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
 507       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
 508         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 509
 510       fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
 511       fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
 512
 513       fh2PtNEF->Fill(jetPt,jet->NEF());
 514       fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
 515       fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
 516       fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
 517
 518       AliVParticle *vp;
 519       Double_t sumPtCh = 0.;
 520       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 521         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 522         if(!vp) continue;
 523         fh2Ptz->Fill(jetPt,GetZ(vp,jet));
 524         sumPtCh+=vp->Pt();
 525       }
 526
 527       if(jet->GetNumberOfTracks()>0)
 528         fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
 529
 530
 531       AliVCluster *vc = 0x0;
 532       Double_t sumPtNe = 0.;
 533       if (clusCont) {
 534         for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
 535           vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
 536           if(!vc) continue;
 537
 538           TLorentzVector lp;
 539           vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 540           sumPtNe+=lp.Pt();
 541
 542         }
 543
 544         if(jet->GetNumberOfClusters()>0)
 545           fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
 546       }
 547     } //full jet loop
 548
 549     for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
 550       Int_t nJetsInEvent = nJetsArr->At(i);
 551       fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
 552     }
 553
 554   }
 555
 556   //Reset array to zero to also count charged jets
 557   nJetsArr->Reset(0);
 558
 559   //Loop over charged jets
 560   if (GetJetContainer(fContainerCharged)) {
 561     const Int_t njets = GetNJets(fContainerCharged);
 562     for (Int_t ij = 0; ij < njets; ij++) {
 563
 564       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
 565       if (!jet)
 566         continue; //jet not selected
 567
 568       Double_t jetPt = jet->Pt();
 569       if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
 570       fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
 571       fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
 572
 573       AliVParticle *vp;
 574       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 575         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 576         if(!vp) continue;
 577         fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
 578       }
 579
 580       //count jets above certain pT threshold
 581       Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
 582       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
 583         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 584
 585     }//ch jet loop
 586     for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
 587       Int_t nJetsInEvent = nJetsArr->At(i);
 588       fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
 589     }
 590   }
 591
 592   if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged)) {
 593     fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
 594   }
 595
 596   fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
 597   fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
 598
 599   if(nJetsArr) delete nJetsArr;
 600
 601   return kTRUE;
 602 }
 603
 604 //________________________________________________________________________
 605 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
 606 {
 607   // Run analysis code here, if needed. It will be executed before FillHistograms().
 608
 609   //Check if event is selected (vertex & pile-up)
 610   if(!SelectEvent())
 611     return kFALSE;
 612
 613   LoadExtraBranches();
 614
 615   if(!fTriggerClass.IsNull())
 616     FindTriggerPatch();
 617
 618   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 619 }
 620
 621 //_______________________________________________________________________
 622 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
 623 {
 624   // Called once at the end of the analysis.
 625 }
 626 //________________________________________________________________________
 627 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
 628 {
 629   // Get Z of constituent trk
 630
 631   return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
 632 }
 633
 634 //________________________________________________________________________
 635 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
 636 {
 637   //
 638   // Get the z of a constituent inside of a jet
 639   //
 640
 641   Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
 642
 643   if(pJetSq>0.)
 644     return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
 645   else {
 646     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
 647     return 0;
 648   }
 649 }