PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTriggerQA.cxx

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