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