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