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 <TProfile.h>
  11 #include <THnSparse.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14
  15 #include "AliVCluster.h"
  16 #include "AliVTrack.h"
  17 #include "AliEmcalJet.h"
  18 #include "AliRhoParameter.h"
  19 #include "AliLog.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 #include "AliParticleContainer.h"
  27 #include "AliEmcalTriggerPatchInfo.h"
  28 #include "AliAODHeader.h"
  29 #include "AliPicoTrack.h"
  30
  31 #include "AliAnalysisTaskEmcalJetTriggerQA.h"
  32
  33 ClassImp(AliAnalysisTaskEmcalJetTriggerQA)
  34
  35 //________________________________________________________________________
  36 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() :
  37   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
  38   fDebug(kFALSE),
  39   fTriggerClass(""),
  40   fContainerFull(0),
  41   fContainerCharged(1),
  42   fMaxPatchEnergy(0),
  43   fTriggerType(-1),
  44   fNFastOR(16),
  45   fhNEvents(0),
  46   fhTriggerbit(0),
  47   fHistRhovsCentFull(0),
  48   fHistRhovsCentCharged(0),
  49   fh3PtEtaPhiTracks(0),
  50   fh3PtEtaPhiTracksOnEmcal(0),
  51   fh3PtEtaPhiTracksToProp(0),
  52   fh3PtEtaPhiTracksProp(0),
  53   fh3PtEtaPhiTracksNoProp(0),
  54   fh2CentPtJetFull(0),
  55   fh2CentPtJetCharged(0),
  56   fh3PtEtaPhiJetFull(0),
  57   fh3PtEtaPhiJetCharged(0),
  58   fh2NJetsPtFull(0),
  59   fh2NJetsPtCharged(0),
  60   fh3PtEtaAreaJetFull(0),
  61   fh3PtEtaAreaJetCharged(0),
  62   fh2PtNConstituentsCharged(0),
  63   fh2PtNConstituents(0),
  64   fh2PtMeanPtConstituentsCharged(0),
  65   fh2PtMeanPtConstituentsNeutral(0),
  66   fh2PtNEF(0),
  67   fh3NEFEtaPhi(0),
  68   fh2NEFNConstituentsCharged(0),
  69   fh2NEFNConstituentsNeutral(0),
  70   fh2Ptz(0),
  71   fh2PtzCharged(0),
  72   fh2PtLeadJet1VsLeadJet2(0),
  73   fh3EEtaPhiCluster(0),
  74   fh3PtLeadJet1VsPatchEnergy(0),
  75   fh3PtLeadJet2VsPatchEnergy(0),
  76   fh3PatchEnergyEtaPhiCenterJ1(0),
  77   fh3PatchEnergyEtaPhiCenterJ2(0),
  78   fh3PatchEnergyEtaPhiCenterJ1J2(0),
  79   fh3PatchADCEnergyEtaPhiCenterJ1(0),
  80   fh3PatchADCEnergyEtaPhiCenterJ2(0),
  81   fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
  82   fh3PatchADCEnergyEtaPhiCenterAll(0),
  83   fh3EEtaPhiCell(0),
  84   fh2ECellVsCent(0),
  85   fh2CellEnergyVsTime(0),
  86   fh3EClusELeadingCellVsTime(0),
  87   fh3JetReacCent(0)
  88 {
  89   // Default constructor.
  90
  91   SetMakeGeneralHistograms(kTRUE);
  92 }
  93
  94 //________________________________________________________________________
  95 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
  96   AliAnalysisTaskEmcalJet(name, kTRUE),
  97   fDebug(kFALSE),
  98   fTriggerClass(""),
  99   fContainerFull(0),
 100   fContainerCharged(1),
 101   fMaxPatchEnergy(0),
 102   fTriggerType(-1),
 103   fNFastOR(16),
 104   fhNEvents(0),
 105   fhTriggerbit(0),
 106   fHistRhovsCentFull(0),
 107   fHistRhovsCentCharged(0),
 108   fh3PtEtaPhiTracks(0),
 109   fh3PtEtaPhiTracksOnEmcal(0),
 110   fh3PtEtaPhiTracksToProp(0),
 111   fh3PtEtaPhiTracksProp(0),
 112   fh3PtEtaPhiTracksNoProp(0),
 113   fh2CentPtJetFull(0),
 114   fh2CentPtJetCharged(0),
 115   fh3PtEtaPhiJetFull(0),
 116   fh3PtEtaPhiJetCharged(0),
 117   fh2NJetsPtFull(0),
 118   fh2NJetsPtCharged(0),
 119   fh3PtEtaAreaJetFull(0),
 120   fh3PtEtaAreaJetCharged(0),
 121   fh2PtNConstituentsCharged(0),
 122   fh2PtNConstituents(0),
 123   fh2PtMeanPtConstituentsCharged(0),
 124   fh2PtMeanPtConstituentsNeutral(0),
 125   fh2PtNEF(0),
 126   fh3NEFEtaPhi(0),
 127   fh2NEFNConstituentsCharged(0),
 128   fh2NEFNConstituentsNeutral(0),
 129   fh2Ptz(0),
 130   fh2PtzCharged(0),
 131   fh2PtLeadJet1VsLeadJet2(0),
 132   fh3EEtaPhiCluster(0),
 133   fh3PtLeadJet1VsPatchEnergy(0),
 134   fh3PtLeadJet2VsPatchEnergy(0),
 135   fh3PatchEnergyEtaPhiCenterJ1(0),
 136   fh3PatchEnergyEtaPhiCenterJ2(0),
 137   fh3PatchEnergyEtaPhiCenterJ1J2(0),
 138   fh3PatchADCEnergyEtaPhiCenterJ1(0),
 139   fh3PatchADCEnergyEtaPhiCenterJ2(0),
 140   fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
 141   fh3PatchADCEnergyEtaPhiCenterAll(0),
 142   fh3EEtaPhiCell(0),
 143   fh2ECellVsCent(0),
 144   fh2CellEnergyVsTime(0),
 145   fh3EClusELeadingCellVsTime(0),
 146   fh3JetReacCent(0)
 147 {
 148   // Standard constructor.
 149
 150   SetMakeGeneralHistograms(kTRUE);
 151 }
 152
 153 //________________________________________________________________________
 154 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
 155 {
 156   // Destructor.
 157
 158 }
 159
 160 //________________________________________________________________________
 161 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
 162   //
 163   // Decide if event should be selected for analysis
 164   //
 165
 166   fhNEvents->Fill(3.5);
 167
 168   if(!fTriggerClass.IsNull()) {
 169     //Check if requested trigger was fired
 170     TString trigType1 = "J1";
 171     TString trigType2 = "J2";
 172     if(fTriggerClass.Contains("G")) {
 173       trigType1 = "G1";
 174       trigType2 = "G2";
 175     }
 176
 177     TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
 178     if(fTriggerClass.Contains(trigType1.Data()) && fTriggerClass.Contains(trigType2.Data())) { //if events with J1&&J2 are requested
 179       if(!firedTrigClass.Contains(trigType1.Data()) || !firedTrigClass.Contains(trigType2.Data()) ) //check if both are fired
 180         return kFALSE;
 181     }
 182     else {
 183       if(!firedTrigClass.Contains(fTriggerClass))
 184         return kFALSE;
 185       else if(fTriggerClass.Contains(trigType1.Data()) && firedTrigClass.Contains(trigType2.Data())) //if J2 is requested also add triggers which have J1&&J2. Reject if J1 is requested and J2 is fired
 186         return kFALSE;
 187     }
 188   }
 189   fhNEvents->Fill(1.5);
 190
 191   return kTRUE;
 192 }
 193
 194 //________________________________________________________________________
 195 void AliAnalysisTaskEmcalJetTriggerQA::FillTriggerPatchHistos() {
 196
 197   //Fill trigger patch histos for main trigger
 198
 199   AliEmcalTriggerPatchInfo *patch = GetMainTriggerPatch();
 200   if(patch) {
 201     fMaxPatchEnergy = patch->GetPatchE();
 202     Double_t patchADCGeV = patch->GetADCAmpGeVRough();
 203     fh3PatchADCEnergyEtaPhiCenterAll->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
 204     if(patch->IsJetLow() && !patch->IsJetHigh()) { //main patch only fired low threshold trigger
 205       fh3PatchEnergyEtaPhiCenterJ2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 206       fh3PatchADCEnergyEtaPhiCenterJ2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
 207     }
 208     else if(patch->IsJetHigh() && !patch->IsJetLow()) { //main patch only fired high threshold trigger - should never happen
 209       fh3PatchEnergyEtaPhiCenterJ1->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 210       fh3PatchADCEnergyEtaPhiCenterJ1->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
 211     }
 212     else if(patch->IsJetHigh() && patch->IsJetLow()) { //main patch fired both triggers
 213       fh3PatchEnergyEtaPhiCenterJ1J2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
 214       fh3PatchADCEnergyEtaPhiCenterJ1J2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
 215     }
 216   }
 217 }
 218
 219 //________________________________________________________________________
 220 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
 221 {
 222   // Create user output.
 223
 224   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 225
 226   Bool_t oldStatus = TH1::AddDirectoryStatus();
 227   TH1::AddDirectory(kFALSE);
 228
 229   fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
 230   fOutput->Add(fhNEvents);
 231
 232   fhTriggerbit = new TProfile("fhTriggerbit","fhTriggerbit;;TriggerBit",1,0,1);
 233   fOutput->Add(fhTriggerbit);
 234
 235   fHistRhovsCentFull = new TH2F("fHistRhovsCentFull", "fHistRhovsCentFull", 101, -1,  100, 300, 0., 300.);
 236   fHistRhovsCentFull->GetXaxis()->SetTitle("Centrality (%)");
 237   fHistRhovsCentFull->GetYaxis()->SetTitle("s#rho_{ch} (GeV/c * rad^{-1})");
 238   fOutput->Add(fHistRhovsCentFull);
 239
 240   fHistRhovsCentCharged = new TH2F("fHistRhovsCentCharged", "fHistRhovsCentCharged", 101, -1,  100, 300, 0., 300.);
 241   fHistRhovsCentCharged->GetXaxis()->SetTitle("Centrality (%)");
 242   fHistRhovsCentCharged->GetYaxis()->SetTitle("#rho_{ch} (GeV/c * rad^{-1})");
 243   fOutput->Add(fHistRhovsCentCharged);
 244
 245   Int_t fgkNCentBins = 21;
 246   Float_t kMinCent   = 0.;
 247   Float_t kMaxCent   = 105.;
 248   Double_t *binsCent = new Double_t[fgkNCentBins+1];
 249   for(Int_t i=0; i<=fgkNCentBins; i++) binsCent[i]=(Double_t)kMinCent + (kMaxCent-kMinCent)/fgkNCentBins*(Double_t)i ;
 250   binsCent[fgkNCentBins-1] = 100.5;
 251   binsCent[fgkNCentBins] = 101.5;
 252
 253   Int_t fgkNdEPBins = 18*8;
 254   Float_t kMindEP   = 0.;
 255   Float_t kMaxdEP   = 1.*TMath::Pi()/2.;
 256   Double_t *binsdEP = new Double_t[fgkNdEPBins+1];
 257   for(Int_t i=0; i<=fgkNdEPBins; i++) binsdEP[i]=(Double_t)kMindEP + (kMaxdEP-kMindEP)/fgkNdEPBins*(Double_t)i ;
 258
 259   Int_t fgkNPtBins = 200;
 260   Float_t kMinPt   = -50.;
 261   Float_t kMaxPt   = 150.;
 262   Double_t *binsPt = new Double_t[fgkNPtBins+1];
 263   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
 264
 265   Int_t fgkNPhiBins = 18*8;
 266   Float_t kMinPhi   = 0.;
 267   Float_t kMaxPhi   = 2.*TMath::Pi();
 268   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
 269   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
 270
 271   Int_t fgkNEtaBins = 100;
 272   Float_t fgkEtaMin = -1.;
 273   Float_t fgkEtaMax =  1.;
 274   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
 275   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
 276
 277   Int_t fgkNAreaBins = 100;
 278   Float_t kMinArea   = 0.;
 279   Float_t kMaxArea   = 1.;
 280   Double_t *binsArea = new Double_t[fgkNAreaBins+1];
 281   for(Int_t i=0; i<=fgkNAreaBins; i++) binsArea[i]=(Double_t)kMinArea + (kMaxArea-kMinArea)/fgkNAreaBins*(Double_t)i ;
 282
 283   Int_t fgkNConstBins = 100;
 284   Float_t kMinConst   = 0.;
 285   Float_t kMaxConst   = 100.;
 286   Double_t *binsConst = new Double_t[fgkNConstBins+1];
 287   for(Int_t i=0; i<=fgkNConstBins; i++) binsConst[i]=(Double_t)kMinConst + (kMaxConst-kMinConst)/fgkNConstBins*(Double_t)i ;
 288
 289   Int_t fgkNMeanPtBins = 100;
 290   Float_t kMinMeanPt   = 0.;
 291   Float_t kMaxMeanPt   = 20.;
 292   Double_t *binsMeanPt = new Double_t[fgkNMeanPtBins+1];
 293   for(Int_t i=0; i<=fgkNMeanPtBins; i++) binsMeanPt[i]=(Double_t)kMinMeanPt + (kMaxMeanPt-kMinMeanPt)/fgkNMeanPtBins*(Double_t)i ;
 294
 295   Int_t fgkNNEFBins = 101;
 296   Float_t kMinNEF   = 0.;
 297   Float_t kMaxNEF   = 1.01;
 298   Double_t *binsNEF = new Double_t[fgkNNEFBins+1];
 299   for(Int_t i=0; i<=fgkNNEFBins; i++) binsNEF[i]=(Double_t)kMinNEF + (kMaxNEF-kMinNEF)/fgkNNEFBins*(Double_t)i ;
 300
 301   Int_t fgkNzBins = 101;
 302   Float_t kMinz   = 0.;
 303   Float_t kMaxz   = 1.01;
 304   Double_t *binsz = new Double_t[fgkNzBins+1];
 305   for(Int_t i=0; i<=fgkNzBins; i++) binsz[i]=(Double_t)kMinz + (kMaxz-kMinz)/fgkNzBins*(Double_t)i ;
 306
 307   Int_t fgkNJetTypeBins = 2;
 308   Float_t kMinJetType   = -0.5;
 309   Float_t kMaxJetType   = 1.5;
 310   Double_t *binsJetType = new Double_t[fgkNJetTypeBins+1];
 311   for(Int_t i=0; i<=fgkNJetTypeBins; i++) binsJetType[i]=(Double_t)kMinJetType + (kMaxJetType-kMinJetType)/fgkNJetTypeBins*(Double_t)i ;
 312
 313   Int_t fgkNTimeBins = 100;
 314   Float_t kMinTime   = -200.;
 315   Float_t kMaxTime   = 200;
 316   Double_t *binsTime = new Double_t[fgkNTimeBins+1];
 317   for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
 318
 319   Double_t enBinEdges[3][2];
 320   enBinEdges[0][0] = 1.; //10 bins
 321   enBinEdges[0][1] = 0.1;
 322   enBinEdges[1][0] = 5.; //8 bins
 323   enBinEdges[1][1] = 0.5;
 324   enBinEdges[2][0] = 100.;//95 bins
 325   enBinEdges[2][1] = 1.;
 326
 327   const Float_t enmin1 =  0;
 328   const Float_t enmax1 =  enBinEdges[0][0];
 329   const Float_t enmin2 =  enmax1 ;
 330   const Float_t enmax2 =  enBinEdges[1][0];
 331   const Float_t enmin3 =  enmax2 ;
 332   const Float_t enmax3 =  enBinEdges[2][0];//fgkEnMax;
 333   const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
 334   const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
 335   const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
 336
 337   Int_t fgkNEnBins=nbin13;
 338   Double_t *binsEn=new Double_t[fgkNEnBins+1];
 339   for(Int_t i=0; i<=fgkNEnBins; i++) {
 340     if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
 341     if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 342     if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 343   }
 344
 345   fh3PtEtaPhiTracks = new TH3F("fh3PtEtaPhiTracks","fh3PtEtaPhiTracks;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 346   fOutput->Add(fh3PtEtaPhiTracks);
 347
 348   fh3PtEtaPhiTracksOnEmcal = new TH3F("fh3PtEtaPhiTracksOnEmcal","fh3PtEtaPhiTracksOnEmcal;#it{p}_{T}^{track}_{emc};#eta_{emc};#varphi_{emc}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 349   fOutput->Add(fh3PtEtaPhiTracksOnEmcal);
 350
 351   fh3PtEtaPhiTracksToProp = new TH3F("fh3PtEtaPhiTracksToProp","fh3PtEtaPhiTracksToProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 352   fOutput->Add(fh3PtEtaPhiTracksToProp);
 353
 354   fh3PtEtaPhiTracksProp = new TH3F("fh3PtEtaPhiTracksProp","fh3PtEtaPhiTracksProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 355   fOutput->Add(fh3PtEtaPhiTracksProp);
 356
 357   fh3PtEtaPhiTracksNoProp = new TH3F("fh3PtEtaPhiTracksNoProp","fh3PtEtaPhiTracksNoProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 358   fOutput->Add(fh3PtEtaPhiTracksNoProp);
 359
 360   fh2CentPtJetFull = new TH2F("fh2CentPtJetFull","fh2CentPtJetFull;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
 361   fOutput->Add(fh2CentPtJetFull);
 362
 363   fh2CentPtJetCharged = new TH2F("fh2CentPtJetCharged","fh2CentPtJetCharged;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
 364   fOutput->Add(fh2CentPtJetCharged);
 365
 366   fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 367   fOutput->Add(fh3PtEtaPhiJetFull);
 368
 369   fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 370   fOutput->Add(fh3PtEtaPhiJetCharged);
 371
 372   fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
 373   fOutput->Add(fh2NJetsPtFull);
 374
 375   fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
 376   fOutput->Add(fh2NJetsPtCharged);
 377
 378   fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
 379   fOutput->Add(fh3PtEtaAreaJetFull);
 380
 381   fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
 382   fOutput->Add(fh3PtEtaAreaJetCharged);
 383
 384   fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
 385   fOutput->Add(fh2PtNConstituentsCharged);
 386
 387   fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
 388   fOutput->Add(fh2PtNConstituents);
 389
 390   fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
 391   fOutput->Add(fh2PtMeanPtConstituentsCharged);
 392
 393   fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
 394   fOutput->Add(fh2PtMeanPtConstituentsNeutral);
 395
 396   fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",fgkNPtBins,binsPt,fgkNNEFBins,binsNEF);
 397   fOutput->Add(fh2PtNEF);
 398
 399   fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",fgkNNEFBins,binsNEF,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 400   fOutput->Add(fh3NEFEtaPhi);
 401
 402   fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
 403   fOutput->Add(fh2NEFNConstituentsCharged);
 404
 405   fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
 406   fOutput->Add(fh2NEFNConstituentsNeutral);
 407
 408   fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
 409   fOutput->Add(fh2Ptz);
 410
 411   fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
 412   fOutput->Add(fh2PtzCharged);
 413
 414   fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",fgkNPtBins,binsPt,fgkNPtBins,binsPt);
 415   fOutput->Add(fh2PtLeadJet1VsLeadJet2);
 416
 417   fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 418   fOutput->Add(fh3EEtaPhiCluster);
 419
 420   fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
 421   fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
 422   fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
 423   fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
 424
 425   fh3PatchEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1","fh3PatchEnergyEtaPhiCenterJ1;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 426   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1);
 427
 428   fh3PatchEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ2","fh3PatchEnergyEtaPhiCenterJ2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 429   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ2);
 430
 431   fh3PatchEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1J2","fh3PatchEnergyEtaPhiCenterJ1J2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 432   fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1J2);
 433
 434   fh3PatchADCEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1","fh3PatchADCEnergyEtaPhiCenterJ1;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 435   fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1);
 436
 437   fh3PatchADCEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ2","fh3PatchADCEnergyEtaPhiCenterJ2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 438   fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ2);
 439
 440   fh3PatchADCEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1J2","fh3PatchADCEnergyEtaPhiCenterJ1J2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 441   fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1J2);
 442
 443   fh3PatchADCEnergyEtaPhiCenterAll = new TH3F("fh3PatchADCEnergyEtaPhiCenterAll","fh3PatchADCEnergyEtaPhiCenterAll;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 444   fOutput->Add(fh3PatchADCEnergyEtaPhiCenterAll);
 445
 446   fh3EEtaPhiCell = new TH3F("fh3EEtaPhiCell","fh3EEtaPhiCell;E_{cell};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 447   fOutput->Add(fh3EEtaPhiCell);
 448
 449   fh2ECellVsCent = new TH2F("fh2ECellVsCent","fh2ECellVsCent;centrality;E_{cell}",101,-1,100,500,0.,5.);
 450   fOutput->Add(fh2ECellVsCent);
 451
 452   fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
 453   fOutput->Add(fh2CellEnergyVsTime);
 454
 455   fh3EClusELeadingCellVsTime = new TH3F("fh3EClusELeadingCellVsTime","fh3EClusELeadingCellVsTime;E_{cluster};E_{leading cell};time_{leading cell}",fgkNEnBins,binsEn,fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
 456   fOutput->Add(fh3EClusELeadingCellVsTime);
 457
 458   fh3JetReacCent = new TH3F("fh3JetReacCent","fh3JetReacCent;E_{Jet};Centrality;dEP",fgkNEnBins,binsEn,fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
 459   fOutput->Add(fh3JetReacCent);
 460
 461   // =========== Switch on Sumw2 for all histos ===========
 462   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 463     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 464     if (h1){
 465       h1->Sumw2();
 466       continue;
 467     }
 468     TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
 469     if (h2){
 470       h2->Sumw2();
 471       continue;
 472     }
 473     TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
 474     if (h3){
 475       h3->Sumw2();
 476       continue;
 477     }
 478     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 479     if(hn)hn->Sumw2();
 480   }
 481
 482   TH1::AddDirectory(oldStatus);
 483
 484   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 485
 486   delete [] binsCent;
 487   if(binsdEP)               delete [] binsdEP;
 488   if(binsEn)                delete [] binsEn;
 489   if(binsPt)                delete [] binsPt;
 490   if(binsPhi)               delete [] binsPhi;
 491   if(binsEta)               delete [] binsEta;
 492   if(binsArea)              delete [] binsArea;
 493   if(binsConst)             delete [] binsConst;
 494   if(binsMeanPt)            delete [] binsMeanPt;
 495   if(binsNEF)               delete [] binsNEF;
 496   if(binsz)                 delete [] binsz;
 497   if(binsJetType)           delete [] binsJetType;
 498   if(binsTime)              delete [] binsTime;
 499
 500 }
 501
 502 //________________________________________________________________________
 503 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
 504 {
 505   // Fill histograms.
 506
 507   //Tracks
 508   AliParticleContainer *partCont = GetParticleContainer(0);
 509   if (partCont) {
 510     AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle(0));
 511     while(track) {
 512       Double_t trkphi = track->Phi()*TMath::RadToDeg();
 513       fh3PtEtaPhiTracks->Fill(track->Pt(),track->Eta(),track->Phi());
 514       //Select tracks which should be propagated
 515       if(track->Pt()>=0.350) {
 516         if (TMath::Abs(track->Eta())<=0.9 && trkphi > 10 && trkphi < 250) {
 517           fh3PtEtaPhiTracksOnEmcal->Fill(track->GetTrackPtOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPhiOnEMCal());
 518           fh3PtEtaPhiTracksToProp->Fill(track->Pt(),track->Eta(),track->Phi());
 519           if(track->GetTrackPtOnEMCal()>=0)
 520             fh3PtEtaPhiTracksProp->Fill(track->Pt(),track->Eta(),track->Phi());
 521           else
 522             fh3PtEtaPhiTracksNoProp->Fill(track->Pt(),track->Eta(),track->Phi());
 523         }
 524       }
 525       track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle());
 526     }
 527   }
 528
 529   //Clusters
 530   AliClusterContainer  *clusCont = GetClusterContainer(0);
 531   if (clusCont) {
 532     Int_t nclusters = clusCont->GetNClusters();
 533     for (Int_t ic = 0; ic < nclusters; ic++) {
 534       AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
 535       if (!cluster) {
 536         AliDebug(2,Form("Could not receive cluster %d", ic));
 537         continue;
 538       }
 539       if (!cluster->IsEMCAL()) {
 540         AliDebug(2,Form("%s: Cluster is not emcal",GetName()));
 541         continue;
 542       }
 543
 544       TLorentzVector lp;
 545       cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 546       fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
 547       if(fCaloCells) {
 548         Double_t leadCellE = GetEnergyLeadingCell(cluster);
 549         Double_t leadCellT = cluster->GetTOF();
 550         fh3EClusELeadingCellVsTime->Fill(lp.E(),leadCellE,leadCellT*1e9);
 551       }
 552     }
 553   }
 554
 555   //Cells
 556   if(fCaloCells) {
 557     const Short_t nCells   = fCaloCells->GetNumberOfCells();
 558
 559     for(Int_t iCell=0; iCell<nCells; ++iCell) {
 560       Short_t cellId = fCaloCells->GetCellNumber(iCell);
 561       Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
 562       Double_t cellT = fCaloCells->GetCellTime(cellId);
 563       TVector3 pos;
 564       fGeom->GetGlobal(cellId, pos);
 565       TLorentzVector lv(pos,cellE);
 566       Double_t cellEta = lv.Eta();
 567       Double_t cellPhi = lv.Phi();
 568       if(cellPhi<0.) cellPhi+=TMath::TwoPi();
 569       if(cellPhi>TMath::TwoPi()) cellPhi-=TMath::TwoPi();
 570
 571       AliDebug(2,Form("cell energy = %f  time = %f",cellE,cellT*1e9));
 572       fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
 573       fh3EEtaPhiCell->Fill(cellE,cellEta,cellPhi);
 574       fh2ECellVsCent->Fill(fCent,cellE);
 575     }
 576   }
 577
 578   //Jets
 579   Double_t ptLeadJet1 = 0.;
 580   Double_t ptLeadJet2 = 0.;
 581
 582   fHistRhovsCentFull->Fill(fCent,GetRhoVal(fContainerFull));
 583   fHistRhovsCentCharged->Fill(fCent,GetRhoVal(fContainerCharged));
 584
 585   TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
 586   nJetsArr->Reset(0);
 587   nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
 588
 589   if (GetJetContainer(fContainerFull)) {
 590     const Int_t njets = GetNJets(fContainerFull);
 591     for (Int_t ij = 0; ij < njets; ij++) {
 592
 593       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
 594       if (!jet)
 595         continue; //jet not selected
 596
 597       Double_t jetPt = jet->Pt() - GetRhoVal(fContainerFull)*jet->Area();
 598       if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
 599
 600       Double_t dEPJetFull = RelativeEP(jet->Phi() , fEPV0);
 601       fh3JetReacCent->Fill(jet->E(),fCent,dEPJetFull);
 602
 603       fh2CentPtJetFull->Fill(fCent,jetPt);
 604       fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
 605       fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
 606
 607       //count jets above certain pT threshold
 608       Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
 609       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
 610         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 611
 612       fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
 613       fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
 614       fh2PtNEF->Fill(jetPt,jet->NEF());
 615       fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
 616       fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
 617       fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
 618
 619       AliVParticle *vp;
 620       Double_t sumPtCh = 0.;
 621       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 622         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 623         if(!vp) continue;
 624         fh2Ptz->Fill(jetPt,GetZ(vp,jet));
 625         sumPtCh+=vp->Pt();
 626       }
 627
 628       if(jet->GetNumberOfTracks()>0)
 629         fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
 630
 631       AliVCluster *vc = 0x0;
 632       Double_t sumPtNe = 0.;
 633       if (clusCont) {
 634         for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
 635           vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
 636           if(!vc) continue;
 637           TLorentzVector lp;
 638           vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
 639           sumPtNe+=lp.Pt();
 640         }
 641
 642         if(jet->GetNumberOfClusters()>0)
 643           fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
 644       }
 645     } //full jet loop
 646
 647     for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
 648       Int_t nJetsInEvent = nJetsArr->At(i);
 649       fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
 650     }
 651   }
 652
 653   //Reset array to zero to also count charged jets
 654   nJetsArr->Reset(0);
 655
 656   //Loop over charged jets
 657   if (GetJetContainer(fContainerCharged)) {
 658     const Int_t njets = GetNJets(fContainerCharged);
 659     for (Int_t ij = 0; ij < njets; ij++) {
 660
 661       AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
 662       if (!jet)
 663         continue; //jet not selected
 664
 665       Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
 666       if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
 667       fh2CentPtJetCharged->Fill(fCent,jetPt);
 668       fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
 669       fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
 670
 671       AliVParticle *vp;
 672       for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 673         vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 674         if(!vp) continue;
 675         fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
 676       }
 677
 678       //count jets above certain pT threshold
 679       Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
 680       for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
 681         nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
 682
 683     }//ch jet loop
 684     for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
 685       Int_t nJetsInEvent = nJetsArr->At(i);
 686       fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
 687     }
 688   }
 689
 690   if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged))
 691     fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
 692
 693   fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
 694   fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
 695
 696   delete nJetsArr;
 697
 698   return kTRUE;
 699 }
 700
 701 //________________________________________________________________________
 702 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
 703 {
 704   // Run analysis code here, if needed. It will be executed before FillHistograms().
 705
 706   fhTriggerbit->Fill(0.5,GetCollisionCandidates());
 707
 708   //Check if event is selected (vertex & pile-up)
 709   if(!SelectEvent())
 710     return kFALSE;
 711
 712   if(fTriggerPatchInfo)
 713     FillTriggerPatchHistos();
 714
 715   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 716 }
 717
 718 //_______________________________________________________________________
 719 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
 720 {
 721   // Called once at the end of the analysis.
 722 }
 723 //________________________________________________________________________
 724 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
 725 {
 726   // Get Z of constituent trk
 727   return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
 728 }
 729
 730 //________________________________________________________________________
 731 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(Double_t trkPx, Double_t trkPy, Double_t trkPz, Double_t jetPx, Double_t jetPy, Double_t jetPz) const
 732 {
 733   //
 734   // Get the z of a constituent inside of a jet
 735   //
 736   Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
 737   if(pJetSq>0.)
 738     return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
 739   else {
 740     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
 741     return 0;
 742   }
 743 }
 744
 745 //________________________________________________________________________
 746 Int_t AliAnalysisTaskEmcalJetTriggerQA::GetLeadingCellId(const AliVCluster *clus) const
 747 {
 748   //Get energy of leading cell in cluster
 749
 750   if(!fCaloCells)
 751     return -1;
 752
 753   Double_t emax = -1.;
 754   Int_t iCellAbsIdMax = -1;
 755   Int_t nCells = clus->GetNCells();
 756   for(Int_t i = 0; i<nCells; i++) {
 757     Int_t absId = clus->GetCellAbsId(i);
 758     Double_t cellE = fCaloCells->GetCellAmplitude(absId);
 759     if(cellE>emax) {
 760       emax = cellE;
 761       iCellAbsIdMax = absId;
 762     }
 763   }
 764   return iCellAbsIdMax;
 765 }
 766
 767 //________________________________________________________________________
 768 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetEnergyLeadingCell(const AliVCluster *clus) const
 769 {
 770   //Get energy of leading cell in cluster
 771   if(!fCaloCells)
 772     return -1.;
 773
 774   Int_t absID = GetLeadingCellId(clus);
 775   if(absID>-1)
 776     return fCaloCells->GetCellAmplitude(absID);
 777   else
 778     return -1.;
 779 }
 780
 781 //________________________________________________________________________
 782 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetECross(Int_t absID) const {
 783
 784   //Get Ecross = sum of energy of neighbouring cells (using uncalibrated energy)
 785
 786   if(!fCaloCells)
 787     return -1.;
 788
 789   Double_t ecross = -1.;
 790
 791   Int_t absID1 = -1;
 792   Int_t absID2 = -1;
 793   Int_t absID3 = -1;
 794   Int_t absID4 = -1;
 795
 796   Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 797   fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
 798   fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 799
 800   if( iphi < AliEMCALGeoParams::fgkEMCALRows-1)
 801     absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
 802   if( iphi > 0 )
 803     absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
 804
 805   if( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) {
 806     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
 807     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta-1);
 808   }
 809   else if( ieta == 0 && imod%2 ) {
 810     absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod,   iphi, ieta+1);
 811     absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
 812   }
 813   else  {
 814     if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
 815       absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
 816     if( ieta > 0 )
 817       absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
 818   }
 819
 820   Double_t ecell1 = fCaloCells->GetCellAmplitude(absID1);
 821   Double_t ecell2 = fCaloCells->GetCellAmplitude(absID2);
 822   Double_t ecell3 = fCaloCells->GetCellAmplitude(absID3);
 823   Double_t ecell4 = fCaloCells->GetCellAmplitude(absID4);
 824
 825   ecross = ecell1+ecell2+ecell3+ecell4;
 826
 827   return ecross;
 828 }
 829
 830 //_________________________________________________________________________
 831 Float_t AliAnalysisTaskEmcalJetTriggerQA::RelativeEP(Double_t objAng, Double_t EPAng) const
 832 {
 833   // function to calculate angle between object and EP in the 1st quadrant (0,Pi/2)
 834   Double_t dphi = EPAng - objAng;
 835
 836   // ran into trouble with a few dEP<-Pi so trying this...
 837   if( dphi<-1*TMath::Pi() )
 838     dphi = dphi + 1*TMath::Pi();
 839   if( dphi>1*TMath::Pi())
 840     dphi = dphi - 1*TMath::Pi();
 841
 842   if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
 843     // Do nothing! we are in quadrant 1
 844   }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
 845     dphi = 1*TMath::Pi() - dphi;
 846   }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
 847     dphi = fabs(dphi);
 848   }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
 849     dphi = dphi + 1*TMath::Pi();
 850   }
 851
 852   return dphi;   // dphi in [0, Pi/2]
 853 }