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