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