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