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