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