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