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 "AliEmcalParticle.h"
  20 #include "AliAODCaloTrigger.h"
  21 #include "AliEMCALGeometry.h"
  22 #include "AliVCaloCells.h"
  23 #include "AliJetContainer.h"
  24 #include "AliClusterContainer.h"
  25 #include "AliEmcalTriggerPatchInfo.h"
  26 #include "AliAODHeader.h"
  27
  28 #include "AliAnalysisTaskEmcalJetTriggerQA.h"
  29
  30 ClassImp(AliAnalysisTaskEmcalJetTriggerQA)
  31
  32 //________________________________________________________________________
  33 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() :
  34   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
  35   fDebug(kFALSE),
  36   fTriggerClass(""),
  37   fBitJ1((1<<8)),
  38   fBitJ2((1<<11)),
  39   fContainerFull(0),
  40   fContainerCharged(1),
  41   fMaxPatchEnergy(0),
  42   fTriggerType(-1),
  43   fNFastOR(16),
  44   fhNEvents(0),
  45   fh3PtEtaPhiJetFull(0),
  46   fh3PtEtaPhiJetCharged(0),
  47   fh2NJetsPtFull(0),
  48   fh2NJetsPtCharged(0),
  49   fh3PtEtaAreaJetFull(0),
  50   fh3PtEtaAreaJetCharged(0),
  51   fh2PtNConstituentsCharged(0),
  52   fh2PtNConstituents(0),
  53   fh2PtMeanPtConstituentsCharged(0),
  54   fh2PtMeanPtConstituentsNeutral(0),
  55   fh2PtNEF(0),
  56   fh3NEFEtaPhi(0),
  57   fh2NEFNConstituentsCharged(0),
  58   fh2NEFNConstituentsNeutral(0),
  59   fh2Ptz(0),
  60   fh2PtzCharged(0),
  61   fh2PtLeadJet1VsLeadJet2(0),
  62   fh3EEtaPhiCluster(0),
  63   fh3PtLeadJet1VsPatchEnergy(0),
  64   fh3PtLeadJet2VsPatchEnergy(0),
  65   fh3PatchEnergyEtaPhiCenterJ1(0),
  66   fh3PatchEnergyEtaPhiCenterJ2(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   fTriggerClass(""),
  80   fBitJ1((1<<8)),
  81   fBitJ2((1<<11)),
  82   fContainerFull(0),
  83   fContainerCharged(1),
  84   fMaxPatchEnergy(0),
  85   fTriggerType(-1),
  86   fNFastOR(16),
  87   fhNEvents(0),
  88   fh3PtEtaPhiJetFull(0),
  89   fh3PtEtaPhiJetCharged(0),
  90   fh2NJetsPtFull(0),
  91   fh2NJetsPtCharged(0),
  92   fh3PtEtaAreaJetFull(0),
  93   fh3PtEtaAreaJetCharged(0),
  94   fh2PtNConstituentsCharged(0),
  95   fh2PtNConstituents(0),
  96   fh2PtMeanPtConstituentsCharged(0),
  97   fh2PtMeanPtConstituentsNeutral(0),
  98   fh2PtNEF(0),
  99   fh3NEFEtaPhi(0),
 100   fh2NEFNConstituentsCharged(0),
 101   fh2NEFNConstituentsNeutral(0),
 102   fh2Ptz(0),
 103   fh2PtzCharged(0),
 104   fh2PtLeadJet1VsLeadJet2(0),
 105   fh3EEtaPhiCluster(0),
 106   fh3PtLeadJet1VsPatchEnergy(0),
 107   fh3PtLeadJet2VsPatchEnergy(0),
 108   fh3PatchEnergyEtaPhiCenterJ1(0),
 109   fh3PatchEnergyEtaPhiCenterJ2(0),
 110   fh2CellEnergyVsTime(0),
 111   fh3EClusELeadingCellVsTime(0)
 112 {
 113   // Standard constructor.
 114
 115   SetMakeGeneralHistograms(kTRUE);
 116 }
 117
 118 //________________________________________________________________________
 119 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
 120 {
 121   // Destructor.
 122 }
 123
 124 //________________________________________________________________________
 125 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
 126   //
 127   // Decide if event should be selected for analysis
 128   //
 129
 130   fhNEvents->Fill(3.5);
 131
 132   if(!fTriggerClass.IsNull()) {
 133     //Check if requested trigger was fired
 134     TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
 135
 136     if(fTriggerClass.Contains("J1") && fTriggerClass.Contains("J2")) {
 137       if(!firedTrigClass.Contains("J1") || !firedTrigClass.Contains("J2") )
 138         return kFALSE;
 139     }
 140     else {
 141
 142       if(!firedTrigClass.Contains(fTriggerClass))
 143         return kFALSE;
 144       if(fTriggerClass.Contains("J2") && firedTrigClass.Contains("J1")) //only accept J2 triggers which were not fired by J1 as well
 145         return kFALSE;
 146     }
 147   }
 148
 149   fhNEvents->Fill(1.5);
 150
 151   return kTRUE;
 152
 153 }
 154
 155 //________________________________________________________________________
 156 void AliAnalysisTaskEmcalJetTriggerQA::FindTriggerPatch() {
 157
 158   //Code to get position of trigger
 159
 160   AliEmcalTriggerPatchInfo *patch = GetMainTriggerPatch();
 161   if(patch) {
 162     fMaxPatchEnergy = patch->GetPatchE();
 163     if(patch->IsJetLow() && !patch->IsJetHigh())
 164       fh3PatchEnergyEtaPhiCenterJ2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 165     if(patch->IsJetHigh())
 166       fh3PatchEnergyEtaPhiCenterJ1->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 167   }
 168 }
 169
 170 //________________________________________________________________________
 171 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
 172 {
 173   // Create user output.
 174
 175   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 176
 177   Bool_t oldStatus = TH1::AddDirectoryStatus();
 178   TH1::AddDirectory(kFALSE);
 179
 180   fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
 181   fOutput->Add(fhNEvents);
 182
 183   Int_t fgkNPtBins = 200;
 184   Float_t kMinPt   = -50.;
 185   Float_t kMaxPt   = 150.;
 186   Double_t *binsPt = new Double_t[fgkNPtBins+1];
 187   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
 188
 189   Int_t fgkNPhiBins = 18*8;
 190   Float_t kMinPhi   = 0.;
 191   Float_t kMaxPhi   = 2.*TMath::Pi();
 192   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
 193   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
 194
 195   Int_t fgkNEtaBins = 100;
 196   Float_t fgkEtaMin = -1.;
 197   Float_t fgkEtaMax =  1.;
 198   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
 199   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
 200
 201   Int_t fgkNAreaBins = 100;
 202   Float_t kMinArea   = 0.;
 203   Float_t kMaxArea   = 1.;
 204   Double_t *binsArea = new Double_t[fgkNAreaBins+1];
 205   for(Int_t i=0; i<=fgkNAreaBins; i++) binsArea[i]=(Double_t)kMinArea + (kMaxArea-kMinArea)/fgkNAreaBins*(Double_t)i ;
 206
 207   Int_t fgkNConstBins = 100;
 208   Float_t kMinConst   = 0.;
 209   Float_t kMaxConst   = 100.;
 210   Double_t *binsConst = new Double_t[fgkNConstBins+1];
 211   for(Int_t i=0; i<=fgkNConstBins; i++) binsConst[i]=(Double_t)kMinConst + (kMaxConst-kMinConst)/fgkNConstBins*(Double_t)i ;
 212
 213   Int_t fgkNMeanPtBins = 200;
 214   Float_t kMinMeanPt   = 0.;
 215   Float_t kMaxMeanPt   = 20.;
 216   Double_t *binsMeanPt = new Double_t[fgkNMeanPtBins+1];
 217   for(Int_t i=0; i<=fgkNMeanPtBins; i++) binsMeanPt[i]=(Double_t)kMinMeanPt + (kMaxMeanPt-kMinMeanPt)/fgkNMeanPtBins*(Double_t)i ;
 218
 219   Int_t fgkNNEFBins = 101;
 220   Float_t kMinNEF   = 0.;
 221   Float_t kMaxNEF   = 1.01;
 222   Double_t *binsNEF = new Double_t[fgkNNEFBins+1];
 223   for(Int_t i=0; i<=fgkNNEFBins; i++) binsNEF[i]=(Double_t)kMinNEF + (kMaxNEF-kMinNEF)/fgkNNEFBins*(Double_t)i ;
 224
 225   Int_t fgkNzBins = 101;
 226   Float_t kMinz   = 0.;
 227   Float_t kMaxz   = 1.01;
 228   Double_t *binsz = new Double_t[fgkNzBins+1];
 229   for(Int_t i=0; i<=fgkNzBins; i++) binsz[i]=(Double_t)kMinz + (kMaxz-kMinz)/fgkNzBins*(Double_t)i ;
 230
 231   Int_t fgkNJetTypeBins = 2;
 232   Float_t kMinJetType   = -0.5;
 233   Float_t kMaxJetType   = 1.5;
 234   Double_t *binsJetType = new Double_t[fgkNJetTypeBins+1];
 235   for(Int_t i=0; i<=fgkNJetTypeBins; i++) binsJetType[i]=(Double_t)kMinJetType + (kMaxJetType-kMinJetType)/fgkNJetTypeBins*(Double_t)i ;
 236
 237   Int_t fgkNTimeBins = 700;
 238   Float_t kMinTime   = -400.;
 239   Float_t kMaxTime   = 1000;
 240   Double_t *binsTime = new Double_t[fgkNTimeBins+1];
 241   for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
 242
 243   Double_t enBinEdges[3][2];
 244   enBinEdges[0][0] = 1.; //10 bins
 245   enBinEdges[0][1] = 0.1;
 246   enBinEdges[1][0] = 5.; //8 bins
 247   enBinEdges[1][1] = 0.5;
 248   enBinEdges[2][0] = 100.;//95 bins
 249   enBinEdges[2][1] = 1.;
 250
 251   const Float_t enmin1 =  0;
 252   const Float_t enmax1 =  enBinEdges[0][0];
 253   const Float_t enmin2 =  enmax1 ;
 254   const Float_t enmax2 =  enBinEdges[1][0];
 255   const Float_t enmin3 =  enmax2 ;
 256   const Float_t enmax3 =  enBinEdges[2][0];//fgkEnMax;
 257   const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
 258   const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
 259   const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
 260
 261   Int_t fgkNEnBins=nbin13;
 262   Double_t *binsEn=new Double_t[fgkNEnBins+1];
 263   for(Int_t i=0; i<=fgkNEnBins; i++) {
 264     if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
 265     if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 266     if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 267   }
 268
 269
 270   fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 271   fOutput->Add(fh3PtEtaPhiJetFull);
 272
 273   fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 274   fOutput->Add(fh3PtEtaPhiJetCharged);
 275
 276   fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
 277   fOutput->Add(fh2NJetsPtFull);
 278
 279   fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
 280   fOutput->Add(fh2NJetsPtCharged);
 281
 282   fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
 283   fOutput->Add(fh3PtEtaAreaJetFull);
 284
 285   fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
 286   fOutput->Add(fh3PtEtaAreaJetCharged);
 287
 288   fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
 289   fOutput->Add(fh2PtNConstituentsCharged);
 290
 291   fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
 292   fOutput->Add(fh2PtNConstituents);
 293
 294   fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
 295   fOutput->Add(fh2PtMeanPtConstituentsCharged);
 296
 297   fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
 298   fOutput->Add(fh2PtMeanPtConstituentsNeutral);
 299
 300   fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",fgkNPtBins,binsPt,fgkNNEFBins,binsNEF);
 301   fOutput->Add(fh2PtNEF);
 302
 303   fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",fgkNNEFBins,binsNEF,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 304   fOutput->Add(fh3NEFEtaPhi);
 305
 306   fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
 307   fOutput->Add(fh2NEFNConstituentsCharged);
 308
 309   fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
 310   fOutput->Add(fh2NEFNConstituentsNeutral);
 311
 312   fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
 313   fOutput->Add(fh2Ptz);
 314
 315   fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
 316   fOutput->Add(fh2PtzCharged);
 317
 318   fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",fgkNPtBins,binsPt,fgkNPtBins,binsPt);
 319   fOutput->Add(fh2PtLeadJet1VsLeadJet2);
 320
 321   fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 322   fOutput->Add(fh3EEtaPhiCluster);
 323
 324   fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
 325   fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
 326   fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
 327   fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
 328
 329   fh3PatchEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1","fh3PatchEnergyEtaPhiCenterJ1;E_{patch};#eta;#phi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 330   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1);
 331
 332   fh3PatchEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ2","fh3PatchEnergyEtaPhiCenterJ2;E_{patch};#eta;#phi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 333   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ2);
 334
 335   fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
 336   fOutput->Add(fh2CellEnergyVsTime);
 337
 338   fh3EClusELeadingCellVsTime = new TH3F("fh3EClusELeadingCellVsTime","fh3EClusELeadingCellVsTime;E_{cluster};E_{leading cell};time_{leading cell}",fgkNEnBins,binsEn,fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
 339   fOutput->Add(fh3EClusELeadingCellVsTime);
 340
 341
 342   // =========== Switch on Sumw2 for all histos ===========
 343   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 344     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 345     if (h1){
 346       h1->Sumw2();
 347       continue;
 348     }
 349     TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
 350     if (h2){
 351       h2->Sumw2();
 352       continue;
 353     }
 354     TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
 355     if (h3){
 356       h3->Sumw2();
 357       continue;
 358     }
 359     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 360     if(hn)hn->Sumw2();
 361   }
 362
 363   TH1::AddDirectory(oldStatus);
 364
 365   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 366
 367   if(binsEn)                delete [] binsEn;
 368   if(binsPt)                delete [] binsPt;
 369   if(binsPhi)               delete [] binsPhi;
 370   if(binsEta)               delete [] binsEta;
 371   if(binsArea)              delete [] binsArea;
 372   if(binsConst)             delete [] binsConst;
 373   if(binsMeanPt)            delete [] binsMeanPt;
 374   if(binsNEF)               delete [] binsNEF;
 375   if(binsz)                 delete [] binsz;
 376   if(binsJetType)           delete [] binsJetType;
 377   if(binsTime)              delete [] binsTime;
 378
 379 }
 380
 381 //________________________________________________________________________
 382 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
 383 {
 384   // Fill histograms.
 385
 386   AliClusterContainer  *clusCont = GetClusterContainer(0);
 387   if (clusCont) {
 388     Int_t nclusters = clusCont->GetNClusters();
 389     TString arrName = clusCont->GetArrayName();
 390     for (Int_t ic = 0; ic < nclusters; ic++) {
 391       AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
 392       if (!cluster) {
 393         AliDebug(2,Form("Could not receive cluster %d", ic));
 394         continue;
 395       }
 396       if (!cluster->IsEMCAL()) {
 397         AliDebug(2,Form("%s: Cluster is not emcal",GetName()));
 398         continue;
 399       }
 400
 401       TLorentzVector lp;
 402       cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 403
 404       //Fill eta,phi,E of clusters here
 405       fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
 406
 407       if(fCaloCells) {
 408         Double_t leadCellE = GetEnergyLeadingCell(cluster);
 409         Double_t leadCellT = cluster->GetTOF();
 410         fh3EClusELeadingCellVsTime->Fill(lp.E(),leadCellE,leadCellT*1e9);
 411       }
 412     }
 413   }
 414
 415   if(fCaloCells) {
 416     const Short_t nCells   = fCaloCells->GetNumberOfCells();
 417
 418     for(Int_t iCell=0; iCell<nCells; ++iCell) {
 419       Short_t cellId = fCaloCells->GetCellNumber(iCell);
 420       Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
 421       Double_t cellT = fCaloCells->GetCellTime(cellId);
 422
 423       AliDebug(2,Form("cell energy = %f  time = %f",cellE,cellT*1e9));
 424       fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
 425     }
 426   }
 427
 428   Double_t ptLeadJet1 = 0.;
 429   Double_t ptLeadJet2 = 0.;
 430
 431   TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
 432   nJetsArr->Reset(0);
 433   nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
 434
 435   if (GetJetContainer(fContainerFull)) {
 436     const Int_t njets = GetNJets(fContainerFull);
 437     for (Int_t ij = 0; ij < njets; ij++) {
 438
 439       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
 440       if (!jet)
 441         continue; //jet not selected
 442
 443       Double_t jetPt = jet->Pt() - GetRhoVal(fContainerFull)*jet->Area();
 444       if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
 445       fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
 446       fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
 447
 448       //count jets above certain pT threshold
 449       Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
 450       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
 451         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 452
 453       fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
 454       fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
 455
 456       fh2PtNEF->Fill(jetPt,jet->NEF());
 457       fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
 458       fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
 459       fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
 460
 461       AliVParticle *vp;
 462       Double_t sumPtCh = 0.;
 463       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 464         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 465         if(!vp) continue;
 466         fh2Ptz->Fill(jetPt,GetZ(vp,jet));
 467         sumPtCh+=vp->Pt();
 468       }
 469
 470       if(jet->GetNumberOfTracks()>0)
 471         fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
 472
 473
 474       AliVCluster *vc = 0x0;
 475       Double_t sumPtNe = 0.;
 476       if (clusCont) {
 477         for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
 478           vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
 479           if(!vc) continue;
 480
 481           TLorentzVector lp;
 482           vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 483           sumPtNe+=lp.Pt();
 484
 485         }
 486
 487         if(jet->GetNumberOfClusters()>0)
 488           fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
 489       }
 490     } //full jet loop
 491
 492     for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
 493       Int_t nJetsInEvent = nJetsArr->At(i);
 494       fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
 495     }
 496
 497   }
 498
 499   //Reset array to zero to also count charged jets
 500   nJetsArr->Reset(0);
 501
 502   //Loop over charged jets
 503   if (GetJetContainer(fContainerCharged)) {
 504     const Int_t njets = GetNJets(fContainerCharged);
 505     for (Int_t ij = 0; ij < njets; ij++) {
 506
 507       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
 508       if (!jet)
 509         continue; //jet not selected
 510
 511       Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
 512       if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
 513       fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
 514       fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
 515
 516       AliVParticle *vp;
 517       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 518         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 519         if(!vp) continue;
 520         fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
 521       }
 522
 523       //count jets above certain pT threshold
 524       Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
 525       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
 526         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 527
 528     }//ch jet loop
 529     for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
 530       Int_t nJetsInEvent = nJetsArr->At(i);
 531       fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
 532     }
 533   }
 534
 535   if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged)) {
 536     fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
 537   }
 538
 539   fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
 540   fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
 541
 542   if(nJetsArr) delete nJetsArr;
 543
 544   return kTRUE;
 545 }
 546
 547 //________________________________________________________________________
 548 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
 549 {
 550   // Run analysis code here, if needed. It will be executed before FillHistograms().
 551
 552   //Check if event is selected (vertex & pile-up)
 553   if(!SelectEvent())
 554     return kFALSE;
 555
 556   if(!fTriggerClass.IsNull())
 557     FindTriggerPatch();
 558
 559   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 560 }
 561
 562 //_______________________________________________________________________
 563 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
 564 {
 565   // Called once at the end of the analysis.
 566 }
 567 //________________________________________________________________________
 568 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
 569 {
 570   // Get Z of constituent trk
 571
 572   return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
 573 }
 574
 575 //________________________________________________________________________
 576 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
 577 {
 578   //
 579   // Get the z of a constituent inside of a jet
 580   //
 581
 582   Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
 583
 584   if(pJetSq>0.)
 585     return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
 586   else {
 587     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
 588     return 0;
 589   }
 590 }
 591
 592 //________________________________________________________________________
 593 Int_t AliAnalysisTaskEmcalJetTriggerQA::GetLeadingCellId(const AliVCluster *clus) const
 594 {
 595   //Get energy of leading cell in cluster
 596
 597   if(!fCaloCells)
 598     return -1;
 599
 600   Double_t emax = -1.;
 601   Int_t iCellAbsIdMax = -1;
 602   Int_t nCells = clus->GetNCells();
 603   for(Int_t i = 0; i<nCells; i++) {
 604     Int_t absId = clus->GetCellAbsId(i);
 605     Double_t cellE = fCaloCells->GetCellAmplitude(absId);
 606     if(cellE>emax) {
 607       emax = cellE;
 608       iCellAbsIdMax = absId;
 609     }
 610   }
 611
 612   return iCellAbsIdMax;
 613 }
 614
 615 //________________________________________________________________________
 616 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetEnergyLeadingCell(const AliVCluster *clus) const
 617 {
 618   //Get energy of leading cell in cluster
 619   if(!fCaloCells)
 620     return -1.;
 621
 622   Int_t absID = GetLeadingCellId(clus);
 623   if(absID>-1)
 624     return fCaloCells->GetCellAmplitude(absID);
 625   else
 626     return -1.;
 627
 628 }