2 // Jet trigger QA analysis task.
6 #include <TClonesArray.h>
10 #include <THnSparse.h>
12 #include <TLorentzVector.h>
14 #include "AliVCluster.h"
15 #include "AliVTrack.h"
16 #include "AliEmcalJet.h"
17 #include "AliRhoParameter.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"
30 #include "AliAnalysisTaskEmcalJetTriggerQA.h"
32 ClassImp(AliAnalysisTaskEmcalJetTriggerQA)
34 //________________________________________________________________________
35 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() :
36 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
47 fHistRhovsCentFull(0),
48 fHistRhovsCentCharged(0),
50 fh3PtEtaPhiTracksOnEmcal(0),
51 fh3PtEtaPhiTracksProp(0),
52 fh3PtEtaPhiJetFull(0),
53 fh3PtEtaPhiJetCharged(0),
56 fh3PtEtaAreaJetFull(0),
57 fh3PtEtaAreaJetCharged(0),
58 fh2PtNConstituentsCharged(0),
59 fh2PtNConstituents(0),
60 fh2PtMeanPtConstituentsCharged(0),
61 fh2PtMeanPtConstituentsNeutral(0),
64 fh2NEFNConstituentsCharged(0),
65 fh2NEFNConstituentsNeutral(0),
68 fh2PtLeadJet1VsLeadJet2(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),
83 // Default constructor.
85 SetMakeGeneralHistograms(kTRUE);
88 //________________________________________________________________________
89 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
90 AliAnalysisTaskEmcalJet(name, kTRUE),
101 fHistRhovsCentFull(0),
102 fHistRhovsCentCharged(0),
103 fh3PtEtaPhiTracks(0),
104 fh3PtEtaPhiTracksOnEmcal(0),
105 fh3PtEtaPhiTracksProp(0),
106 fh3PtEtaPhiJetFull(0),
107 fh3PtEtaPhiJetCharged(0),
109 fh2NJetsPtCharged(0),
110 fh3PtEtaAreaJetFull(0),
111 fh3PtEtaAreaJetCharged(0),
112 fh2PtNConstituentsCharged(0),
113 fh2PtNConstituents(0),
114 fh2PtMeanPtConstituentsCharged(0),
115 fh2PtMeanPtConstituentsNeutral(0),
118 fh2NEFNConstituentsCharged(0),
119 fh2NEFNConstituentsNeutral(0),
122 fh2PtLeadJet1VsLeadJet2(0),
123 fh3EEtaPhiCluster(0),
124 fh3PtLeadJet1VsPatchEnergy(0),
125 fh3PtLeadJet2VsPatchEnergy(0),
126 fh3PatchEnergyEtaPhiCenterJ1(0),
127 fh3PatchEnergyEtaPhiCenterJ2(0),
128 fh3PatchEnergyEtaPhiCenterJ1J2(0),
129 fh3PatchADCEnergyEtaPhiCenterJ1(0),
130 fh3PatchADCEnergyEtaPhiCenterJ2(0),
131 fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
132 fh2CellEnergyVsTime(0),
133 fh3EClusELeadingCellVsTime(0),
137 // Standard constructor.
139 SetMakeGeneralHistograms(kTRUE);
142 //________________________________________________________________________
143 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
147 delete fOutput; // delete output object list
153 //________________________________________________________________________
154 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
156 // Decide if event should be selected for analysis
159 fhNEvents->Fill(3.5);
161 if(!fTriggerClass.IsNull()) {
162 //Check if requested trigger was fired
163 TString trigType1 = "J1";
164 TString trigType2 = "J2";
165 if(fTriggerClass.Contains("G")) {
170 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
172 if(fTriggerClass.Contains(trigType1.Data()) && fTriggerClass.Contains(trigType2.Data())) { //if events with J1&&J2 are requested
173 if(!firedTrigClass.Contains(trigType1.Data()) || !firedTrigClass.Contains(trigType2.Data()) ) //check if both are fired
177 if(!firedTrigClass.Contains(fTriggerClass))
179 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
184 fhNEvents->Fill(1.5);
190 //________________________________________________________________________
191 void AliAnalysisTaskEmcalJetTriggerQA::FindTriggerPatch() {
193 //Fill trigger patch histos for main trigger
195 AliEmcalTriggerPatchInfo *patch = GetMainTriggerPatch();
197 fMaxPatchEnergy = patch->GetPatchE();
198 Double_t patchADCGeV = patch->GetADCAmpGeVRough();
199 if(patch->IsJetLow() && !patch->IsJetHigh()) { //main patch only fired low threshold trigger
200 fh3PatchEnergyEtaPhiCenterJ2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
201 fh3PatchADCEnergyEtaPhiCenterJ2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
203 else if(patch->IsJetHigh() && !patch->IsJetLow()) { //main patch only fired high threshold trigger - should never happen
204 fh3PatchEnergyEtaPhiCenterJ1->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
205 fh3PatchADCEnergyEtaPhiCenterJ1->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
207 else if(patch->IsJetHigh() && patch->IsJetLow()) { //main patch fired both triggers
208 fh3PatchEnergyEtaPhiCenterJ1J2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
209 fh3PatchADCEnergyEtaPhiCenterJ1J2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
214 //________________________________________________________________________
215 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
217 // Create user output.
219 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
221 Bool_t oldStatus = TH1::AddDirectoryStatus();
222 TH1::AddDirectory(kFALSE);
224 fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
225 fOutput->Add(fhNEvents);
227 fHistRhovsCentFull = new TH2F("fHistRhovsCentFull", "fHistRhovsCentFull", 101, -1, 100, 300, 0., 300.);
228 fHistRhovsCentFull->GetXaxis()->SetTitle("Centrality (%)");
229 fHistRhovsCentFull->GetYaxis()->SetTitle("s#rho_{ch} (GeV/c * rad^{-1})");
230 fOutput->Add(fHistRhovsCentFull);
232 fHistRhovsCentCharged = new TH2F("fHistRhovsCentCharged", "fHistRhovsCentCharged", 101, -1, 100, 300, 0., 300.);
233 fHistRhovsCentCharged->GetXaxis()->SetTitle("Centrality (%)");
234 fHistRhovsCentCharged->GetYaxis()->SetTitle("#rho_{ch} (GeV/c * rad^{-1})");
235 fOutput->Add(fHistRhovsCentCharged);
237 Int_t fgkNCentBins = 20;
238 Float_t kMinCent = 0.;
239 Float_t kMaxCent = 100.;
240 Double_t *binsCent = new Double_t[fgkNCentBins+1];
241 for(Int_t i=0; i<=fgkNCentBins; i++) binsCent[i]=(Double_t)kMinCent + (kMaxCent-kMinCent)/fgkNCentBins*(Double_t)i ;
243 Int_t fgkNdEPBins = 18*8;
244 Float_t kMindEP = 0.;
245 Float_t kMaxdEP = 1.*TMath::Pi()/2.;
246 Double_t *binsdEP = new Double_t[fgkNdEPBins+1];
247 for(Int_t i=0; i<=fgkNdEPBins; i++) binsdEP[i]=(Double_t)kMindEP + (kMaxdEP-kMindEP)/fgkNdEPBins*(Double_t)i ;
249 Int_t fgkNPtBins = 200;
250 Float_t kMinPt = -50.;
251 Float_t kMaxPt = 150.;
252 Double_t *binsPt = new Double_t[fgkNPtBins+1];
253 for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
255 Int_t fgkNPhiBins = 18*8;
256 Float_t kMinPhi = 0.;
257 Float_t kMaxPhi = 2.*TMath::Pi();
258 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
259 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
261 Int_t fgkNEtaBins = 100;
262 Float_t fgkEtaMin = -1.;
263 Float_t fgkEtaMax = 1.;
264 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
265 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
267 Int_t fgkNAreaBins = 100;
268 Float_t kMinArea = 0.;
269 Float_t kMaxArea = 1.;
270 Double_t *binsArea = new Double_t[fgkNAreaBins+1];
271 for(Int_t i=0; i<=fgkNAreaBins; i++) binsArea[i]=(Double_t)kMinArea + (kMaxArea-kMinArea)/fgkNAreaBins*(Double_t)i ;
273 Int_t fgkNConstBins = 100;
274 Float_t kMinConst = 0.;
275 Float_t kMaxConst = 100.;
276 Double_t *binsConst = new Double_t[fgkNConstBins+1];
277 for(Int_t i=0; i<=fgkNConstBins; i++) binsConst[i]=(Double_t)kMinConst + (kMaxConst-kMinConst)/fgkNConstBins*(Double_t)i ;
279 Int_t fgkNMeanPtBins = 100;
280 Float_t kMinMeanPt = 0.;
281 Float_t kMaxMeanPt = 20.;
282 Double_t *binsMeanPt = new Double_t[fgkNMeanPtBins+1];
283 for(Int_t i=0; i<=fgkNMeanPtBins; i++) binsMeanPt[i]=(Double_t)kMinMeanPt + (kMaxMeanPt-kMinMeanPt)/fgkNMeanPtBins*(Double_t)i ;
285 Int_t fgkNNEFBins = 101;
286 Float_t kMinNEF = 0.;
287 Float_t kMaxNEF = 1.01;
288 Double_t *binsNEF = new Double_t[fgkNNEFBins+1];
289 for(Int_t i=0; i<=fgkNNEFBins; i++) binsNEF[i]=(Double_t)kMinNEF + (kMaxNEF-kMinNEF)/fgkNNEFBins*(Double_t)i ;
291 Int_t fgkNzBins = 101;
293 Float_t kMaxz = 1.01;
294 Double_t *binsz = new Double_t[fgkNzBins+1];
295 for(Int_t i=0; i<=fgkNzBins; i++) binsz[i]=(Double_t)kMinz + (kMaxz-kMinz)/fgkNzBins*(Double_t)i ;
297 Int_t fgkNJetTypeBins = 2;
298 Float_t kMinJetType = -0.5;
299 Float_t kMaxJetType = 1.5;
300 Double_t *binsJetType = new Double_t[fgkNJetTypeBins+1];
301 for(Int_t i=0; i<=fgkNJetTypeBins; i++) binsJetType[i]=(Double_t)kMinJetType + (kMaxJetType-kMinJetType)/fgkNJetTypeBins*(Double_t)i ;
303 Int_t fgkNTimeBins = 100;
304 Float_t kMinTime = -200.;
305 Float_t kMaxTime = 200;
306 Double_t *binsTime = new Double_t[fgkNTimeBins+1];
307 for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
309 Double_t enBinEdges[3][2];
310 enBinEdges[0][0] = 1.; //10 bins
311 enBinEdges[0][1] = 0.1;
312 enBinEdges[1][0] = 5.; //8 bins
313 enBinEdges[1][1] = 0.5;
314 enBinEdges[2][0] = 100.;//95 bins
315 enBinEdges[2][1] = 1.;
317 const Float_t enmin1 = 0;
318 const Float_t enmax1 = enBinEdges[0][0];
319 const Float_t enmin2 = enmax1 ;
320 const Float_t enmax2 = enBinEdges[1][0];
321 const Float_t enmin3 = enmax2 ;
322 const Float_t enmax3 = enBinEdges[2][0];//fgkEnMax;
323 const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
324 const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
325 const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
327 Int_t fgkNEnBins=nbin13;
328 Double_t *binsEn=new Double_t[fgkNEnBins+1];
329 for(Int_t i=0; i<=fgkNEnBins; i++) {
330 if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
331 if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
332 if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
335 fh3PtEtaPhiTracks = new TH3F("fh3PtEtaPhiTracks","fh3PtEtaPhiTracks;#it{p}_{T}^{track};#eta;#varphi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
336 fOutput->Add(fh3PtEtaPhiTracks);
338 fh3PtEtaPhiTracksOnEmcal = new TH3F("fh3PtEtaPhiTracksOnEmcal","fh3PtEtaPhiTracksOnEmcal;#it{p}_{T}^{track};#eta;#varphi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
339 fOutput->Add(fh3PtEtaPhiTracksOnEmcal);
341 fh3PtEtaPhiTracksProp = new TH3F("fh3PtEtaPhiTracksProp","fh3PtEtaPhiTracksProp;#it{p}_{T}^{track};#eta;#varphi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
342 fOutput->Add(fh3PtEtaPhiTracksProp);
344 fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
345 fOutput->Add(fh3PtEtaPhiJetFull);
347 fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
348 fOutput->Add(fh3PtEtaPhiJetCharged);
350 fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
351 fOutput->Add(fh2NJetsPtFull);
353 fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
354 fOutput->Add(fh2NJetsPtCharged);
356 fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
357 fOutput->Add(fh3PtEtaAreaJetFull);
359 fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
360 fOutput->Add(fh3PtEtaAreaJetCharged);
362 fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
363 fOutput->Add(fh2PtNConstituentsCharged);
365 fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
366 fOutput->Add(fh2PtNConstituents);
368 fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
369 fOutput->Add(fh2PtMeanPtConstituentsCharged);
371 fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
372 fOutput->Add(fh2PtMeanPtConstituentsNeutral);
374 fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",fgkNPtBins,binsPt,fgkNNEFBins,binsNEF);
375 fOutput->Add(fh2PtNEF);
377 fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",fgkNNEFBins,binsNEF,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
378 fOutput->Add(fh3NEFEtaPhi);
380 fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
381 fOutput->Add(fh2NEFNConstituentsCharged);
383 fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
384 fOutput->Add(fh2NEFNConstituentsNeutral);
386 fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
387 fOutput->Add(fh2Ptz);
389 fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
390 fOutput->Add(fh2PtzCharged);
392 fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",fgkNPtBins,binsPt,fgkNPtBins,binsPt);
393 fOutput->Add(fh2PtLeadJet1VsLeadJet2);
395 fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
396 fOutput->Add(fh3EEtaPhiCluster);
398 fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
399 fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
400 fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
401 fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
403 fh3PatchEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1","fh3PatchEnergyEtaPhiCenterJ1;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
404 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1);
406 fh3PatchEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ2","fh3PatchEnergyEtaPhiCenterJ2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
407 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ2);
409 fh3PatchEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1J2","fh3PatchEnergyEtaPhiCenterJ1J2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
410 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1J2);
412 fh3PatchADCEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1","fh3PatchADCEnergyEtaPhiCenterJ1;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
413 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1);
415 fh3PatchADCEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ2","fh3PatchADCEnergyEtaPhiCenterJ2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
416 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ2);
418 fh3PatchADCEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1J2","fh3PatchADCEnergyEtaPhiCenterJ1J2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
419 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1J2);
421 fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
422 fOutput->Add(fh2CellEnergyVsTime);
424 fh3EClusELeadingCellVsTime = new TH3F("fh3EClusELeadingCellVsTime","fh3EClusELeadingCellVsTime;E_{cluster};E_{leading cell};time_{leading cell}",fgkNEnBins,binsEn,fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
425 fOutput->Add(fh3EClusELeadingCellVsTime);
427 fh3JetReacCent = new TH3F("fh3JetReacCent","fh3JetReacCent;E_{Jet};Centrality;dEP",fgkNEnBins,binsEn,fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
428 fOutput->Add(fh3JetReacCent);
430 fh2FullJetCent = new TH2F("fh2FullJetCent","fh2FullJetCent;Centrality;dEP",fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
431 fOutput->Add(fh2FullJetCent);
434 // =========== Switch on Sumw2 for all histos ===========
435 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
436 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
441 TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
446 TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
451 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
455 TH1::AddDirectory(oldStatus);
457 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
459 if(binsCent) delete [] binsCent;
460 if(binsdEP) delete [] binsdEP;
461 if(binsEn) delete [] binsEn;
462 if(binsPt) delete [] binsPt;
463 if(binsPhi) delete [] binsPhi;
464 if(binsEta) delete [] binsEta;
465 if(binsArea) delete [] binsArea;
466 if(binsConst) delete [] binsConst;
467 if(binsMeanPt) delete [] binsMeanPt;
468 if(binsNEF) delete [] binsNEF;
469 if(binsz) delete [] binsz;
470 if(binsJetType) delete [] binsJetType;
471 if(binsTime) delete [] binsTime;
475 //________________________________________________________________________
476 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
481 AliParticleContainer *partCont = GetParticleContainer(0);
484 AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle(0));
486 Double_t trkphi = track->Phi()*TMath::RadToDeg();
487 fh3PtEtaPhiTracks->Fill(track->Pt(),track->Eta(),track->Phi());
488 fh3PtEtaPhiTracksOnEmcal->Fill(track->GetTrackPtOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPhiOnEMCal());
489 if(track->IsEMCAL()) {
491 if(TMath::Abs(track->Eta())<0.9 && trkphi > 10 && trkphi < 250 )
492 fh3PtEtaPhiTracksProp->Fill(track->Pt(),track->Eta(),track->Phi());
494 track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle());
499 AliClusterContainer *clusCont = GetClusterContainer(0);
501 Int_t nclusters = clusCont->GetNClusters();
502 TString arrName = clusCont->GetArrayName();
503 for (Int_t ic = 0; ic < nclusters; ic++) {
504 AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
506 AliDebug(2,Form("Could not receive cluster %d", ic));
509 if (!cluster->IsEMCAL()) {
510 AliDebug(2,Form("%s: Cluster is not emcal",GetName()));
515 cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
516 fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
518 Double_t leadCellE = GetEnergyLeadingCell(cluster);
519 Double_t leadCellT = cluster->GetTOF();
520 fh3EClusELeadingCellVsTime->Fill(lp.E(),leadCellE,leadCellT*1e9);
527 const Short_t nCells = fCaloCells->GetNumberOfCells();
529 for(Int_t iCell=0; iCell<nCells; ++iCell) {
530 Short_t cellId = fCaloCells->GetCellNumber(iCell);
531 Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
532 Double_t cellT = fCaloCells->GetCellTime(cellId);
534 AliDebug(2,Form("cell energy = %f time = %f",cellE,cellT*1e9));
535 fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
540 Double_t ptLeadJet1 = 0.;
541 Double_t ptLeadJet2 = 0.;
543 fHistRhovsCentFull->Fill(fCent,GetRhoVal(fContainerFull));
544 fHistRhovsCentCharged->Fill(fCent,GetRhoVal(fContainerCharged));
546 TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
548 nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
550 if (GetJetContainer(fContainerFull)) {
551 const Int_t njets = GetNJets(fContainerFull);
552 for (Int_t ij = 0; ij < njets; ij++) {
554 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
556 continue; //jet not selected
558 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerFull)*jet->Area();
559 Double_t jetPhi = jet->Phi();
560 Double_t dEPJetFull = -500.0;
561 dEPJetFull = RelativeEP(jetPhi , fEPV0);
563 fh3JetReacCent->Fill(jet->E(),fCent,dEPJetFull);
564 fh2FullJetCent->Fill(fCent,dEPJetFull);
566 if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
567 fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
568 fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
570 //count jets above certain pT threshold
571 Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
572 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
573 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
575 fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
576 fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
578 fh2PtNEF->Fill(jetPt,jet->NEF());
579 fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
580 fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
581 fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
584 Double_t sumPtCh = 0.;
585 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
586 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
588 fh2Ptz->Fill(jetPt,GetZ(vp,jet));
592 if(jet->GetNumberOfTracks()>0)
593 fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
595 AliVCluster *vc = 0x0;
596 Double_t sumPtNe = 0.;
598 for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
599 vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
602 vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
606 if(jet->GetNumberOfClusters()>0)
607 fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
611 for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
612 Int_t nJetsInEvent = nJetsArr->At(i);
613 fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
617 //Reset array to zero to also count charged jets
620 //Loop over charged jets
621 if (GetJetContainer(fContainerCharged)) {
622 const Int_t njets = GetNJets(fContainerCharged);
623 for (Int_t ij = 0; ij < njets; ij++) {
625 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
627 continue; //jet not selected
629 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
630 if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
631 fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
632 fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
635 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
636 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
638 fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
641 //count jets above certain pT threshold
642 Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
643 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
644 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
647 for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
648 Int_t nJetsInEvent = nJetsArr->At(i);
649 fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
653 if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged)) {
654 fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
657 fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
658 fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
660 if(nJetsArr) delete nJetsArr;
665 //________________________________________________________________________
666 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
668 // Run analysis code here, if needed. It will be executed before FillHistograms().
670 //Check if event is selected (vertex & pile-up)
674 if(fTriggerPatchInfo)
677 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
680 //_______________________________________________________________________
681 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
683 // Called once at the end of the analysis.
685 //________________________________________________________________________
686 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet) const
688 // Get Z of constituent trk
689 return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
692 //________________________________________________________________________
693 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
696 // Get the z of a constituent inside of a jet
698 Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
700 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
702 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
707 //________________________________________________________________________
708 Int_t AliAnalysisTaskEmcalJetTriggerQA::GetLeadingCellId(const AliVCluster *clus) const
710 //Get energy of leading cell in cluster
716 Int_t iCellAbsIdMax = -1;
717 Int_t nCells = clus->GetNCells();
718 for(Int_t i = 0; i<nCells; i++) {
719 Int_t absId = clus->GetCellAbsId(i);
720 Double_t cellE = fCaloCells->GetCellAmplitude(absId);
723 iCellAbsIdMax = absId;
727 return iCellAbsIdMax;
730 //________________________________________________________________________
731 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetEnergyLeadingCell(const AliVCluster *clus) const
733 //Get energy of leading cell in cluster
737 Int_t absID = GetLeadingCellId(clus);
739 return fCaloCells->GetCellAmplitude(absID);
744 //________________________________________________________________________
745 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetECross(Int_t absID) const {
747 //Get Ecross = sum of energy of neighbouring cells (using uncalibrated energy)
752 Double_t ecross = -1.;
759 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
760 fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
761 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
763 if( iphi < AliEMCALGeoParams::fgkEMCALRows-1)
764 absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
766 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
768 if( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) {
769 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
770 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
772 else if( ieta == 0 && imod%2 ) {
773 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
774 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
777 if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
778 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
780 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
783 Double_t ecell1 = fCaloCells->GetCellAmplitude(absID1);
784 Double_t ecell2 = fCaloCells->GetCellAmplitude(absID2);
785 Double_t ecell3 = fCaloCells->GetCellAmplitude(absID3);
786 Double_t ecell4 = fCaloCells->GetCellAmplitude(absID4);
788 ecross = ecell1+ecell2+ecell3+ecell4;
793 //_________________________________________________________________________
794 Float_t AliAnalysisTaskEmcalJetTriggerQA:: RelativeEP(Double_t objAng, Double_t EPAng) const
795 { // function to calculate angle between object and EP in the 1st quadrant (0,Pi/2)
796 Double_t dphi = (EPAng - objAng);
798 // ran into trouble with a few dEP<-Pi so trying this...
799 if( dphi<-1*TMath::Pi() ){
800 dphi = dphi + 1*TMath::Pi();
803 if( dphi>1*TMath::Pi()){
804 dphi = dphi - 1*TMath::Pi();
807 if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
808 // Do nothing! we are in quadrant 1
809 }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
810 dphi = 1*TMath::Pi() - dphi;
811 }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
813 }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
814 dphi = dphi + 1*TMath::Pi();
817 return dphi; // dphi in [0, Pi/2]