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