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