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