2 // Jet trigger QA analysis task.
6 #include <TClonesArray.h>
11 #include <THnSparse.h>
13 #include <TLorentzVector.h>
15 #include "AliVCluster.h"
16 #include "AliVTrack.h"
17 #include "AliEmcalJet.h"
18 #include "AliRhoParameter.h"
20 #include "AliEmcalParticle.h"
21 #include "AliAODCaloTrigger.h"
22 #include "AliEMCALGeometry.h"
23 #include "AliVCaloCells.h"
24 #include "AliJetContainer.h"
25 #include "AliClusterContainer.h"
26 #include "AliParticleContainer.h"
27 #include "AliEmcalTriggerPatchInfo.h"
28 #include "AliAODHeader.h"
29 #include "AliPicoTrack.h"
31 #include "AliAnalysisTaskEmcalJetTriggerQA.h"
33 ClassImp(AliAnalysisTaskEmcalJetTriggerQA)
35 //________________________________________________________________________
36 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() :
37 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
47 fHistRhovsCentFull(0),
48 fHistRhovsCentCharged(0),
50 fh3PtEtaPhiTracksOnEmcal(0),
51 fh3PtEtaPhiTracksToProp(0),
52 fh3PtEtaPhiTracksProp(0),
53 fh3PtEtaPhiTracksNoProp(0),
55 fh2CentPtJetCharged(0),
56 fh3PtEtaPhiJetFull(0),
57 fh3PtEtaPhiJetCharged(0),
60 fh3PtEtaAreaJetFull(0),
61 fh3PtEtaAreaJetCharged(0),
62 fh2PtNConstituentsCharged(0),
63 fh2PtNConstituents(0),
64 fh2PtMeanPtConstituentsCharged(0),
65 fh2PtMeanPtConstituentsNeutral(0),
68 fh2NEFNConstituentsCharged(0),
69 fh2NEFNConstituentsNeutral(0),
72 fh2PtLeadJet1VsLeadJet2(0),
74 fh3PtLeadJet1VsPatchEnergy(0),
75 fh3PtLeadJet2VsPatchEnergy(0),
76 fh3PatchEnergyEtaPhiCenterJ1(0),
77 fh3PatchEnergyEtaPhiCenterJ2(0),
78 fh3PatchEnergyEtaPhiCenterJ1J2(0),
79 fh3PatchADCEnergyEtaPhiCenterJ1(0),
80 fh3PatchADCEnergyEtaPhiCenterJ2(0),
81 fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
82 fh3PatchADCEnergyEtaPhiCenterAll(0),
84 fh2CellEnergyVsTime(0),
85 fh3EClusELeadingCellVsTime(0),
88 // Default constructor.
90 SetMakeGeneralHistograms(kTRUE);
93 //________________________________________________________________________
94 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
95 AliAnalysisTaskEmcalJet(name, kTRUE),
105 fHistRhovsCentFull(0),
106 fHistRhovsCentCharged(0),
107 fh3PtEtaPhiTracks(0),
108 fh3PtEtaPhiTracksOnEmcal(0),
109 fh3PtEtaPhiTracksToProp(0),
110 fh3PtEtaPhiTracksProp(0),
111 fh3PtEtaPhiTracksNoProp(0),
113 fh2CentPtJetCharged(0),
114 fh3PtEtaPhiJetFull(0),
115 fh3PtEtaPhiJetCharged(0),
117 fh2NJetsPtCharged(0),
118 fh3PtEtaAreaJetFull(0),
119 fh3PtEtaAreaJetCharged(0),
120 fh2PtNConstituentsCharged(0),
121 fh2PtNConstituents(0),
122 fh2PtMeanPtConstituentsCharged(0),
123 fh2PtMeanPtConstituentsNeutral(0),
126 fh2NEFNConstituentsCharged(0),
127 fh2NEFNConstituentsNeutral(0),
130 fh2PtLeadJet1VsLeadJet2(0),
131 fh3EEtaPhiCluster(0),
132 fh3PtLeadJet1VsPatchEnergy(0),
133 fh3PtLeadJet2VsPatchEnergy(0),
134 fh3PatchEnergyEtaPhiCenterJ1(0),
135 fh3PatchEnergyEtaPhiCenterJ2(0),
136 fh3PatchEnergyEtaPhiCenterJ1J2(0),
137 fh3PatchADCEnergyEtaPhiCenterJ1(0),
138 fh3PatchADCEnergyEtaPhiCenterJ2(0),
139 fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
140 fh3PatchADCEnergyEtaPhiCenterAll(0),
142 fh2CellEnergyVsTime(0),
143 fh3EClusELeadingCellVsTime(0),
146 // Standard constructor.
148 SetMakeGeneralHistograms(kTRUE);
151 //________________________________________________________________________
152 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
158 //________________________________________________________________________
159 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
161 // Decide if event should be selected for analysis
164 fhNEvents->Fill(3.5);
166 if(!fTriggerClass.IsNull()) {
167 //Check if requested trigger was fired
168 TString trigType1 = "J1";
169 TString trigType2 = "J2";
170 if(fTriggerClass.Contains("G")) {
175 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
176 if(fTriggerClass.Contains(trigType1.Data()) && fTriggerClass.Contains(trigType2.Data())) { //if events with J1&&J2 are requested
177 if(!firedTrigClass.Contains(trigType1.Data()) || !firedTrigClass.Contains(trigType2.Data()) ) //check if both are fired
181 if(!firedTrigClass.Contains(fTriggerClass))
183 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
187 fhNEvents->Fill(1.5);
192 //________________________________________________________________________
193 void AliAnalysisTaskEmcalJetTriggerQA::FillTriggerPatchHistos() {
195 //Fill trigger patch histos for main trigger
197 AliEmcalTriggerPatchInfo *patch = GetMainTriggerPatch();
199 fMaxPatchEnergy = patch->GetPatchE();
200 Double_t patchADCGeV = patch->GetADCAmpGeVRough();
201 fh3PatchADCEnergyEtaPhiCenterAll->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
202 if(patch->IsJetLow() && !patch->IsJetHigh()) { //main patch only fired low threshold trigger
203 fh3PatchEnergyEtaPhiCenterJ2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
204 fh3PatchADCEnergyEtaPhiCenterJ2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
206 else if(patch->IsJetHigh() && !patch->IsJetLow()) { //main patch only fired high threshold trigger - should never happen
207 fh3PatchEnergyEtaPhiCenterJ1->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
208 fh3PatchADCEnergyEtaPhiCenterJ1->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
210 else if(patch->IsJetHigh() && patch->IsJetLow()) { //main patch fired both triggers
211 fh3PatchEnergyEtaPhiCenterJ1J2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
212 fh3PatchADCEnergyEtaPhiCenterJ1J2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
217 //________________________________________________________________________
218 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
220 // Create user output.
222 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
224 Bool_t oldStatus = TH1::AddDirectoryStatus();
225 TH1::AddDirectory(kFALSE);
227 fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
228 fOutput->Add(fhNEvents);
230 fhTriggerbit = new TProfile("fhTriggerbit","fhTriggerbit;;TriggerBit",1,0,1);
231 fOutput->Add(fhTriggerbit);
233 fHistRhovsCentFull = new TH2F("fHistRhovsCentFull", "fHistRhovsCentFull", 101, -1, 100, 300, 0., 300.);
234 fHistRhovsCentFull->GetXaxis()->SetTitle("Centrality (%)");
235 fHistRhovsCentFull->GetYaxis()->SetTitle("s#rho_{ch} (GeV/c * rad^{-1})");
236 fOutput->Add(fHistRhovsCentFull);
238 fHistRhovsCentCharged = new TH2F("fHistRhovsCentCharged", "fHistRhovsCentCharged", 101, -1, 100, 300, 0., 300.);
239 fHistRhovsCentCharged->GetXaxis()->SetTitle("Centrality (%)");
240 fHistRhovsCentCharged->GetYaxis()->SetTitle("#rho_{ch} (GeV/c * rad^{-1})");
241 fOutput->Add(fHistRhovsCentCharged);
243 Int_t fgkNCentBins = 21;
244 Float_t kMinCent = 0.;
245 Float_t kMaxCent = 105.;
246 Double_t *binsCent = new Double_t[fgkNCentBins+1];
247 for(Int_t i=0; i<=fgkNCentBins; i++) binsCent[i]=(Double_t)kMinCent + (kMaxCent-kMinCent)/fgkNCentBins*(Double_t)i ;
248 binsCent[fgkNCentBins-1] = 100.5;
249 binsCent[fgkNCentBins] = 101.5;
251 Int_t fgkNdEPBins = 18*8;
252 Float_t kMindEP = 0.;
253 Float_t kMaxdEP = 1.*TMath::Pi()/2.;
254 Double_t *binsdEP = new Double_t[fgkNdEPBins+1];
255 for(Int_t i=0; i<=fgkNdEPBins; i++) binsdEP[i]=(Double_t)kMindEP + (kMaxdEP-kMindEP)/fgkNdEPBins*(Double_t)i ;
257 Int_t fgkNPtBins = 200;
258 Float_t kMinPt = -50.;
259 Float_t kMaxPt = 150.;
260 Double_t *binsPt = new Double_t[fgkNPtBins+1];
261 for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
263 Int_t fgkNPhiBins = 18*8;
264 Float_t kMinPhi = 0.;
265 Float_t kMaxPhi = 2.*TMath::Pi();
266 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
267 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
269 Int_t fgkNEtaBins = 100;
270 Float_t fgkEtaMin = -1.;
271 Float_t fgkEtaMax = 1.;
272 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
273 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
275 Int_t fgkNAreaBins = 100;
276 Float_t kMinArea = 0.;
277 Float_t kMaxArea = 1.;
278 Double_t *binsArea = new Double_t[fgkNAreaBins+1];
279 for(Int_t i=0; i<=fgkNAreaBins; i++) binsArea[i]=(Double_t)kMinArea + (kMaxArea-kMinArea)/fgkNAreaBins*(Double_t)i ;
281 Int_t fgkNConstBins = 100;
282 Float_t kMinConst = 0.;
283 Float_t kMaxConst = 100.;
284 Double_t *binsConst = new Double_t[fgkNConstBins+1];
285 for(Int_t i=0; i<=fgkNConstBins; i++) binsConst[i]=(Double_t)kMinConst + (kMaxConst-kMinConst)/fgkNConstBins*(Double_t)i ;
287 Int_t fgkNMeanPtBins = 100;
288 Float_t kMinMeanPt = 0.;
289 Float_t kMaxMeanPt = 20.;
290 Double_t *binsMeanPt = new Double_t[fgkNMeanPtBins+1];
291 for(Int_t i=0; i<=fgkNMeanPtBins; i++) binsMeanPt[i]=(Double_t)kMinMeanPt + (kMaxMeanPt-kMinMeanPt)/fgkNMeanPtBins*(Double_t)i ;
293 Int_t fgkNNEFBins = 101;
294 Float_t kMinNEF = 0.;
295 Float_t kMaxNEF = 1.01;
296 Double_t *binsNEF = new Double_t[fgkNNEFBins+1];
297 for(Int_t i=0; i<=fgkNNEFBins; i++) binsNEF[i]=(Double_t)kMinNEF + (kMaxNEF-kMinNEF)/fgkNNEFBins*(Double_t)i ;
299 Int_t fgkNzBins = 101;
301 Float_t kMaxz = 1.01;
302 Double_t *binsz = new Double_t[fgkNzBins+1];
303 for(Int_t i=0; i<=fgkNzBins; i++) binsz[i]=(Double_t)kMinz + (kMaxz-kMinz)/fgkNzBins*(Double_t)i ;
305 Int_t fgkNJetTypeBins = 2;
306 Float_t kMinJetType = -0.5;
307 Float_t kMaxJetType = 1.5;
308 Double_t *binsJetType = new Double_t[fgkNJetTypeBins+1];
309 for(Int_t i=0; i<=fgkNJetTypeBins; i++) binsJetType[i]=(Double_t)kMinJetType + (kMaxJetType-kMinJetType)/fgkNJetTypeBins*(Double_t)i ;
311 Int_t fgkNTimeBins = 100;
312 Float_t kMinTime = -200.;
313 Float_t kMaxTime = 200;
314 Double_t *binsTime = new Double_t[fgkNTimeBins+1];
315 for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
317 Double_t enBinEdges[3][2];
318 enBinEdges[0][0] = 1.; //10 bins
319 enBinEdges[0][1] = 0.1;
320 enBinEdges[1][0] = 5.; //8 bins
321 enBinEdges[1][1] = 0.5;
322 enBinEdges[2][0] = 100.;//95 bins
323 enBinEdges[2][1] = 1.;
325 const Float_t enmin1 = 0;
326 const Float_t enmax1 = enBinEdges[0][0];
327 const Float_t enmin2 = enmax1 ;
328 const Float_t enmax2 = enBinEdges[1][0];
329 const Float_t enmin3 = enmax2 ;
330 const Float_t enmax3 = enBinEdges[2][0];//fgkEnMax;
331 const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
332 const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
333 const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
335 Int_t fgkNEnBins=nbin13;
336 Double_t *binsEn=new Double_t[fgkNEnBins+1];
337 for(Int_t i=0; i<=fgkNEnBins; i++) {
338 if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
339 if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
340 if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
343 fh3PtEtaPhiTracks = new TH3F("fh3PtEtaPhiTracks","fh3PtEtaPhiTracks;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
344 fOutput->Add(fh3PtEtaPhiTracks);
346 fh3PtEtaPhiTracksOnEmcal = new TH3F("fh3PtEtaPhiTracksOnEmcal","fh3PtEtaPhiTracksOnEmcal;#it{p}_{T}^{track}_{emc};#eta_{emc};#varphi_{emc}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
347 fOutput->Add(fh3PtEtaPhiTracksOnEmcal);
349 fh3PtEtaPhiTracksToProp = new TH3F("fh3PtEtaPhiTracksToProp","fh3PtEtaPhiTracksToProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
350 fOutput->Add(fh3PtEtaPhiTracksToProp);
352 fh3PtEtaPhiTracksProp = new TH3F("fh3PtEtaPhiTracksProp","fh3PtEtaPhiTracksProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
353 fOutput->Add(fh3PtEtaPhiTracksProp);
355 fh3PtEtaPhiTracksNoProp = new TH3F("fh3PtEtaPhiTracksNoProp","fh3PtEtaPhiTracksNoProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
356 fOutput->Add(fh3PtEtaPhiTracksNoProp);
358 fh2CentPtJetFull = new TH2F("fh2CentPtJetFull","fh2CentPtJetFull;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
359 fOutput->Add(fh2CentPtJetFull);
361 fh2CentPtJetCharged = new TH2F("fh2CentPtJetCharged","fh2CentPtJetCharged;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
362 fOutput->Add(fh2CentPtJetCharged);
364 fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
365 fOutput->Add(fh3PtEtaPhiJetFull);
367 fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
368 fOutput->Add(fh3PtEtaPhiJetCharged);
370 fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
371 fOutput->Add(fh2NJetsPtFull);
373 fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
374 fOutput->Add(fh2NJetsPtCharged);
376 fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
377 fOutput->Add(fh3PtEtaAreaJetFull);
379 fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
380 fOutput->Add(fh3PtEtaAreaJetCharged);
382 fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
383 fOutput->Add(fh2PtNConstituentsCharged);
385 fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
386 fOutput->Add(fh2PtNConstituents);
388 fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
389 fOutput->Add(fh2PtMeanPtConstituentsCharged);
391 fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
392 fOutput->Add(fh2PtMeanPtConstituentsNeutral);
394 fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",fgkNPtBins,binsPt,fgkNNEFBins,binsNEF);
395 fOutput->Add(fh2PtNEF);
397 fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",fgkNNEFBins,binsNEF,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
398 fOutput->Add(fh3NEFEtaPhi);
400 fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
401 fOutput->Add(fh2NEFNConstituentsCharged);
403 fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
404 fOutput->Add(fh2NEFNConstituentsNeutral);
406 fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
407 fOutput->Add(fh2Ptz);
409 fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
410 fOutput->Add(fh2PtzCharged);
412 fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",fgkNPtBins,binsPt,fgkNPtBins,binsPt);
413 fOutput->Add(fh2PtLeadJet1VsLeadJet2);
415 fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
416 fOutput->Add(fh3EEtaPhiCluster);
418 fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
419 fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
420 fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
421 fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
423 fh3PatchEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1","fh3PatchEnergyEtaPhiCenterJ1;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
424 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1);
426 fh3PatchEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ2","fh3PatchEnergyEtaPhiCenterJ2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
427 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ2);
429 fh3PatchEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1J2","fh3PatchEnergyEtaPhiCenterJ1J2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
430 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1J2);
432 fh3PatchADCEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1","fh3PatchADCEnergyEtaPhiCenterJ1;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
433 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1);
435 fh3PatchADCEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ2","fh3PatchADCEnergyEtaPhiCenterJ2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
436 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ2);
438 fh3PatchADCEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1J2","fh3PatchADCEnergyEtaPhiCenterJ1J2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
439 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1J2);
441 fh3PatchADCEnergyEtaPhiCenterAll = new TH3F("fh3PatchADCEnergyEtaPhiCenterAll","fh3PatchADCEnergyEtaPhiCenterAll;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
442 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterAll);
444 fh3EEtaPhiCell = new TH3F("fh3EEtaPhiCell","fh3EEtaPhiCell;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
445 fOutput->Add(fh3EEtaPhiCell);
447 fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
448 fOutput->Add(fh2CellEnergyVsTime);
450 fh3EClusELeadingCellVsTime = new TH3F("fh3EClusELeadingCellVsTime","fh3EClusELeadingCellVsTime;E_{cluster};E_{leading cell};time_{leading cell}",fgkNEnBins,binsEn,fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
451 fOutput->Add(fh3EClusELeadingCellVsTime);
453 fh3JetReacCent = new TH3F("fh3JetReacCent","fh3JetReacCent;E_{Jet};Centrality;dEP",fgkNEnBins,binsEn,fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
454 fOutput->Add(fh3JetReacCent);
456 // =========== Switch on Sumw2 for all histos ===========
457 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
458 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
463 TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
468 TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
473 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
477 TH1::AddDirectory(oldStatus);
479 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
482 if(binsdEP) delete [] binsdEP;
483 if(binsEn) delete [] binsEn;
484 if(binsPt) delete [] binsPt;
485 if(binsPhi) delete [] binsPhi;
486 if(binsEta) delete [] binsEta;
487 if(binsArea) delete [] binsArea;
488 if(binsConst) delete [] binsConst;
489 if(binsMeanPt) delete [] binsMeanPt;
490 if(binsNEF) delete [] binsNEF;
491 if(binsz) delete [] binsz;
492 if(binsJetType) delete [] binsJetType;
493 if(binsTime) delete [] binsTime;
497 //________________________________________________________________________
498 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
503 AliParticleContainer *partCont = GetParticleContainer(0);
505 AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle(0));
507 Double_t trkphi = track->Phi()*TMath::RadToDeg();
508 fh3PtEtaPhiTracks->Fill(track->Pt(),track->Eta(),track->Phi());
509 //Select tracks which should be propagated
510 if(track->Pt()>=0.350) {
511 if (TMath::Abs(track->Eta())<=0.9 && trkphi > 10 && trkphi < 250) {
512 fh3PtEtaPhiTracksOnEmcal->Fill(track->GetTrackPtOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPhiOnEMCal());
513 fh3PtEtaPhiTracksToProp->Fill(track->Pt(),track->Eta(),track->Phi());
514 if(track->GetTrackPtOnEMCal()>=0)
515 fh3PtEtaPhiTracksProp->Fill(track->Pt(),track->Eta(),track->Phi());
517 fh3PtEtaPhiTracksNoProp->Fill(track->Pt(),track->Eta(),track->Phi());
520 track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle());
525 AliClusterContainer *clusCont = GetClusterContainer(0);
527 Int_t nclusters = clusCont->GetNClusters();
528 for (Int_t ic = 0; ic < nclusters; ic++) {
529 AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
531 AliDebug(2,Form("Could not receive cluster %d", ic));
534 if (!cluster->IsEMCAL()) {
535 AliDebug(2,Form("%s: Cluster is not emcal",GetName()));
540 cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
541 fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
543 Double_t leadCellE = GetEnergyLeadingCell(cluster);
544 Double_t leadCellT = cluster->GetTOF();
545 fh3EClusELeadingCellVsTime->Fill(lp.E(),leadCellE,leadCellT*1e9);
552 const Short_t nCells = fCaloCells->GetNumberOfCells();
554 for(Int_t iCell=0; iCell<nCells; ++iCell) {
555 Short_t cellId = fCaloCells->GetCellNumber(iCell);
556 Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
557 Double_t cellT = fCaloCells->GetCellTime(cellId);
559 fGeom->GetGlobal(cellId, pos);
560 TLorentzVector lv(pos,cellE);
561 Double_t cellEta = lv.Eta();
562 Double_t cellPhi = lv.Phi();
563 if(cellPhi<0.) cellPhi+=TMath::TwoPi();
564 if(cellPhi>TMath::TwoPi()) cellPhi-=TMath::TwoPi();
566 AliDebug(2,Form("cell energy = %f time = %f",cellE,cellT*1e9));
567 fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
568 fh3EEtaPhiCell->Fill(cellE,cellEta,cellPhi);
573 Double_t ptLeadJet1 = 0.;
574 Double_t ptLeadJet2 = 0.;
576 fHistRhovsCentFull->Fill(fCent,GetRhoVal(fContainerFull));
577 fHistRhovsCentCharged->Fill(fCent,GetRhoVal(fContainerCharged));
579 TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
581 nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
583 if (GetJetContainer(fContainerFull)) {
584 const Int_t njets = GetNJets(fContainerFull);
585 for (Int_t ij = 0; ij < njets; ij++) {
587 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
589 continue; //jet not selected
591 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerFull)*jet->Area();
592 if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
594 Double_t dEPJetFull = RelativeEP(jet->Phi() , fEPV0);
595 fh3JetReacCent->Fill(jet->E(),fCent,dEPJetFull);
597 fh2CentPtJetFull->Fill(fCent,jetPt);
598 fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
599 fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
601 //count jets above certain pT threshold
602 Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
603 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
604 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
606 fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
607 fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
608 fh2PtNEF->Fill(jetPt,jet->NEF());
609 fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
610 fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
611 fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
614 Double_t sumPtCh = 0.;
615 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
616 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
618 fh2Ptz->Fill(jetPt,GetZ(vp,jet));
622 if(jet->GetNumberOfTracks()>0)
623 fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
625 AliVCluster *vc = 0x0;
626 Double_t sumPtNe = 0.;
628 for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
629 vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
632 vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
636 if(jet->GetNumberOfClusters()>0)
637 fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
641 for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
642 Int_t nJetsInEvent = nJetsArr->At(i);
643 fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
647 //Reset array to zero to also count charged jets
650 //Loop over charged jets
651 if (GetJetContainer(fContainerCharged)) {
652 const Int_t njets = GetNJets(fContainerCharged);
653 for (Int_t ij = 0; ij < njets; ij++) {
655 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
657 continue; //jet not selected
659 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
660 if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
661 fh2CentPtJetCharged->Fill(fCent,jetPt);
662 fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
663 fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
666 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
667 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
669 fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
672 //count jets above certain pT threshold
673 Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
674 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
675 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
678 for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
679 Int_t nJetsInEvent = nJetsArr->At(i);
680 fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
684 if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged))
685 fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
687 fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
688 fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
695 //________________________________________________________________________
696 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
698 // Run analysis code here, if needed. It will be executed before FillHistograms().
700 fhTriggerbit->Fill(0.5,GetCollisionCandidates());
702 //Check if event is selected (vertex & pile-up)
706 if(fTriggerPatchInfo)
707 FillTriggerPatchHistos();
709 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
712 //_______________________________________________________________________
713 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
715 // Called once at the end of the analysis.
717 //________________________________________________________________________
718 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet) const
720 // Get Z of constituent trk
721 return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
724 //________________________________________________________________________
725 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
728 // Get the z of a constituent inside of a jet
730 Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
732 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
734 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
739 //________________________________________________________________________
740 Int_t AliAnalysisTaskEmcalJetTriggerQA::GetLeadingCellId(const AliVCluster *clus) const
742 //Get energy of leading cell in cluster
748 Int_t iCellAbsIdMax = -1;
749 Int_t nCells = clus->GetNCells();
750 for(Int_t i = 0; i<nCells; i++) {
751 Int_t absId = clus->GetCellAbsId(i);
752 Double_t cellE = fCaloCells->GetCellAmplitude(absId);
755 iCellAbsIdMax = absId;
758 return iCellAbsIdMax;
761 //________________________________________________________________________
762 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetEnergyLeadingCell(const AliVCluster *clus) const
764 //Get energy of leading cell in cluster
768 Int_t absID = GetLeadingCellId(clus);
770 return fCaloCells->GetCellAmplitude(absID);
775 //________________________________________________________________________
776 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetECross(Int_t absID) const {
778 //Get Ecross = sum of energy of neighbouring cells (using uncalibrated energy)
783 Double_t ecross = -1.;
790 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
791 fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
792 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
794 if( iphi < AliEMCALGeoParams::fgkEMCALRows-1)
795 absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
797 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
799 if( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) {
800 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
801 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
803 else if( ieta == 0 && imod%2 ) {
804 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
805 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
808 if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
809 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
811 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
814 Double_t ecell1 = fCaloCells->GetCellAmplitude(absID1);
815 Double_t ecell2 = fCaloCells->GetCellAmplitude(absID2);
816 Double_t ecell3 = fCaloCells->GetCellAmplitude(absID3);
817 Double_t ecell4 = fCaloCells->GetCellAmplitude(absID4);
819 ecross = ecell1+ecell2+ecell3+ecell4;
824 //_________________________________________________________________________
825 Float_t AliAnalysisTaskEmcalJetTriggerQA::RelativeEP(Double_t objAng, Double_t EPAng) const
827 // function to calculate angle between object and EP in the 1st quadrant (0,Pi/2)
828 Double_t dphi = EPAng - objAng;
830 // ran into trouble with a few dEP<-Pi so trying this...
831 if( dphi<-1*TMath::Pi() )
832 dphi = dphi + 1*TMath::Pi();
833 if( dphi>1*TMath::Pi())
834 dphi = dphi - 1*TMath::Pi();
836 if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
837 // Do nothing! we are in quadrant 1
838 }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
839 dphi = 1*TMath::Pi() - dphi;
840 }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
842 }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
843 dphi = dphi + 1*TMath::Pi();
846 return dphi; // dphi in [0, Pi/2]