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),
45 fHistRhovsCentFull(0),
46 fHistRhovsCentCharged(0),
48 fh3PtEtaPhiTracksOnEmcal(0),
49 fh3PtEtaPhiTracksToProp(0),
50 fh3PtEtaPhiTracksProp(0),
51 fh3PtEtaPhiTracksNoProp(0),
53 fh2CentPtJetCharged(0),
54 fh3PtEtaPhiJetFull(0),
55 fh3PtEtaPhiJetCharged(0),
58 fh3PtEtaAreaJetFull(0),
59 fh3PtEtaAreaJetCharged(0),
60 fh2PtNConstituentsCharged(0),
61 fh2PtNConstituents(0),
62 fh2PtMeanPtConstituentsCharged(0),
63 fh2PtMeanPtConstituentsNeutral(0),
66 fh2NEFNConstituentsCharged(0),
67 fh2NEFNConstituentsNeutral(0),
70 fh2PtLeadJet1VsLeadJet2(0),
72 fh3PtLeadJet1VsPatchEnergy(0),
73 fh3PtLeadJet2VsPatchEnergy(0),
74 fh3PatchEnergyEtaPhiCenterJ1(0),
75 fh3PatchEnergyEtaPhiCenterJ2(0),
76 fh3PatchEnergyEtaPhiCenterJ1J2(0),
77 fh3PatchADCEnergyEtaPhiCenterJ1(0),
78 fh3PatchADCEnergyEtaPhiCenterJ2(0),
79 fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
80 fh3PatchADCEnergyEtaPhiCenterAll(0),
82 fh2CellEnergyVsTime(0),
83 fh3EClusELeadingCellVsTime(0),
86 // Default constructor.
88 SetMakeGeneralHistograms(kTRUE);
91 //________________________________________________________________________
92 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
93 AliAnalysisTaskEmcalJet(name, kTRUE),
102 fHistRhovsCentFull(0),
103 fHistRhovsCentCharged(0),
104 fh3PtEtaPhiTracks(0),
105 fh3PtEtaPhiTracksOnEmcal(0),
106 fh3PtEtaPhiTracksToProp(0),
107 fh3PtEtaPhiTracksProp(0),
108 fh3PtEtaPhiTracksNoProp(0),
110 fh2CentPtJetCharged(0),
111 fh3PtEtaPhiJetFull(0),
112 fh3PtEtaPhiJetCharged(0),
114 fh2NJetsPtCharged(0),
115 fh3PtEtaAreaJetFull(0),
116 fh3PtEtaAreaJetCharged(0),
117 fh2PtNConstituentsCharged(0),
118 fh2PtNConstituents(0),
119 fh2PtMeanPtConstituentsCharged(0),
120 fh2PtMeanPtConstituentsNeutral(0),
123 fh2NEFNConstituentsCharged(0),
124 fh2NEFNConstituentsNeutral(0),
127 fh2PtLeadJet1VsLeadJet2(0),
128 fh3EEtaPhiCluster(0),
129 fh3PtLeadJet1VsPatchEnergy(0),
130 fh3PtLeadJet2VsPatchEnergy(0),
131 fh3PatchEnergyEtaPhiCenterJ1(0),
132 fh3PatchEnergyEtaPhiCenterJ2(0),
133 fh3PatchEnergyEtaPhiCenterJ1J2(0),
134 fh3PatchADCEnergyEtaPhiCenterJ1(0),
135 fh3PatchADCEnergyEtaPhiCenterJ2(0),
136 fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
137 fh3PatchADCEnergyEtaPhiCenterAll(0),
139 fh2CellEnergyVsTime(0),
140 fh3EClusELeadingCellVsTime(0),
143 // Standard constructor.
145 SetMakeGeneralHistograms(kTRUE);
148 //________________________________________________________________________
149 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
153 delete fOutput; // delete output object list
159 //________________________________________________________________________
160 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
162 // Decide if event should be selected for analysis
165 fhNEvents->Fill(3.5);
167 if(!fTriggerClass.IsNull()) {
168 //Check if requested trigger was fired
169 TString trigType1 = "J1";
170 TString trigType2 = "J2";
171 if(fTriggerClass.Contains("G")) {
176 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
177 if(fTriggerClass.Contains(trigType1.Data()) && fTriggerClass.Contains(trigType2.Data())) { //if events with J1&&J2 are requested
178 if(!firedTrigClass.Contains(trigType1.Data()) || !firedTrigClass.Contains(trigType2.Data()) ) //check if both are fired
182 if(!firedTrigClass.Contains(fTriggerClass))
184 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
188 fhNEvents->Fill(1.5);
193 //________________________________________________________________________
194 void AliAnalysisTaskEmcalJetTriggerQA::FillTriggerPatchHistos() {
196 //Fill trigger patch histos for main trigger
198 AliEmcalTriggerPatchInfo *patch = GetMainTriggerPatch();
200 fMaxPatchEnergy = patch->GetPatchE();
201 Double_t patchADCGeV = patch->GetADCAmpGeVRough();
202 fh3PatchADCEnergyEtaPhiCenterAll->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
203 if(patch->IsJetLow() && !patch->IsJetHigh()) { //main patch only fired low threshold trigger
204 fh3PatchEnergyEtaPhiCenterJ2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
205 fh3PatchADCEnergyEtaPhiCenterJ2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
207 else if(patch->IsJetHigh() && !patch->IsJetLow()) { //main patch only fired high threshold trigger - should never happen
208 fh3PatchEnergyEtaPhiCenterJ1->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
209 fh3PatchADCEnergyEtaPhiCenterJ1->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
211 else if(patch->IsJetHigh() && patch->IsJetLow()) { //main patch fired both triggers
212 fh3PatchEnergyEtaPhiCenterJ1J2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
213 fh3PatchADCEnergyEtaPhiCenterJ1J2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
218 //________________________________________________________________________
219 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
221 // Create user output.
223 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
225 Bool_t oldStatus = TH1::AddDirectoryStatus();
226 TH1::AddDirectory(kFALSE);
228 fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
229 fOutput->Add(fhNEvents);
231 fHistRhovsCentFull = new TH2F("fHistRhovsCentFull", "fHistRhovsCentFull", 101, -1, 100, 300, 0., 300.);
232 fHistRhovsCentFull->GetXaxis()->SetTitle("Centrality (%)");
233 fHistRhovsCentFull->GetYaxis()->SetTitle("s#rho_{ch} (GeV/c * rad^{-1})");
234 fOutput->Add(fHistRhovsCentFull);
236 fHistRhovsCentCharged = new TH2F("fHistRhovsCentCharged", "fHistRhovsCentCharged", 101, -1, 100, 300, 0., 300.);
237 fHistRhovsCentCharged->GetXaxis()->SetTitle("Centrality (%)");
238 fHistRhovsCentCharged->GetYaxis()->SetTitle("#rho_{ch} (GeV/c * rad^{-1})");
239 fOutput->Add(fHistRhovsCentCharged);
241 Int_t fgkNCentBins = 21;
242 Float_t kMinCent = 0.;
243 Float_t kMaxCent = 105.;
244 Double_t *binsCent = new Double_t[fgkNCentBins+1];
245 for(Int_t i=0; i<=fgkNCentBins; i++) binsCent[i]=(Double_t)kMinCent + (kMaxCent-kMinCent)/fgkNCentBins*(Double_t)i ;
246 binsCent[fgkNCentBins-1] = 100.5;
247 binsCent[fgkNCentBins] = 101.5;
249 Int_t fgkNdEPBins = 18*8;
250 Float_t kMindEP = 0.;
251 Float_t kMaxdEP = 1.*TMath::Pi()/2.;
252 Double_t *binsdEP = new Double_t[fgkNdEPBins+1];
253 for(Int_t i=0; i<=fgkNdEPBins; i++) binsdEP[i]=(Double_t)kMindEP + (kMaxdEP-kMindEP)/fgkNdEPBins*(Double_t)i ;
255 Int_t fgkNPtBins = 200;
256 Float_t kMinPt = -50.;
257 Float_t kMaxPt = 150.;
258 Double_t *binsPt = new Double_t[fgkNPtBins+1];
259 for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
261 Int_t fgkNPhiBins = 18*8;
262 Float_t kMinPhi = 0.;
263 Float_t kMaxPhi = 2.*TMath::Pi();
264 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
265 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
267 Int_t fgkNEtaBins = 100;
268 Float_t fgkEtaMin = -1.;
269 Float_t fgkEtaMax = 1.;
270 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
271 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
273 Int_t fgkNAreaBins = 100;
274 Float_t kMinArea = 0.;
275 Float_t kMaxArea = 1.;
276 Double_t *binsArea = new Double_t[fgkNAreaBins+1];
277 for(Int_t i=0; i<=fgkNAreaBins; i++) binsArea[i]=(Double_t)kMinArea + (kMaxArea-kMinArea)/fgkNAreaBins*(Double_t)i ;
279 Int_t fgkNConstBins = 100;
280 Float_t kMinConst = 0.;
281 Float_t kMaxConst = 100.;
282 Double_t *binsConst = new Double_t[fgkNConstBins+1];
283 for(Int_t i=0; i<=fgkNConstBins; i++) binsConst[i]=(Double_t)kMinConst + (kMaxConst-kMinConst)/fgkNConstBins*(Double_t)i ;
285 Int_t fgkNMeanPtBins = 100;
286 Float_t kMinMeanPt = 0.;
287 Float_t kMaxMeanPt = 20.;
288 Double_t *binsMeanPt = new Double_t[fgkNMeanPtBins+1];
289 for(Int_t i=0; i<=fgkNMeanPtBins; i++) binsMeanPt[i]=(Double_t)kMinMeanPt + (kMaxMeanPt-kMinMeanPt)/fgkNMeanPtBins*(Double_t)i ;
291 Int_t fgkNNEFBins = 101;
292 Float_t kMinNEF = 0.;
293 Float_t kMaxNEF = 1.01;
294 Double_t *binsNEF = new Double_t[fgkNNEFBins+1];
295 for(Int_t i=0; i<=fgkNNEFBins; i++) binsNEF[i]=(Double_t)kMinNEF + (kMaxNEF-kMinNEF)/fgkNNEFBins*(Double_t)i ;
297 Int_t fgkNzBins = 101;
299 Float_t kMaxz = 1.01;
300 Double_t *binsz = new Double_t[fgkNzBins+1];
301 for(Int_t i=0; i<=fgkNzBins; i++) binsz[i]=(Double_t)kMinz + (kMaxz-kMinz)/fgkNzBins*(Double_t)i ;
303 Int_t fgkNJetTypeBins = 2;
304 Float_t kMinJetType = -0.5;
305 Float_t kMaxJetType = 1.5;
306 Double_t *binsJetType = new Double_t[fgkNJetTypeBins+1];
307 for(Int_t i=0; i<=fgkNJetTypeBins; i++) binsJetType[i]=(Double_t)kMinJetType + (kMaxJetType-kMinJetType)/fgkNJetTypeBins*(Double_t)i ;
309 Int_t fgkNTimeBins = 100;
310 Float_t kMinTime = -200.;
311 Float_t kMaxTime = 200;
312 Double_t *binsTime = new Double_t[fgkNTimeBins+1];
313 for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
315 Double_t enBinEdges[3][2];
316 enBinEdges[0][0] = 1.; //10 bins
317 enBinEdges[0][1] = 0.1;
318 enBinEdges[1][0] = 5.; //8 bins
319 enBinEdges[1][1] = 0.5;
320 enBinEdges[2][0] = 100.;//95 bins
321 enBinEdges[2][1] = 1.;
323 const Float_t enmin1 = 0;
324 const Float_t enmax1 = enBinEdges[0][0];
325 const Float_t enmin2 = enmax1 ;
326 const Float_t enmax2 = enBinEdges[1][0];
327 const Float_t enmin3 = enmax2 ;
328 const Float_t enmax3 = enBinEdges[2][0];//fgkEnMax;
329 const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
330 const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
331 const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
333 Int_t fgkNEnBins=nbin13;
334 Double_t *binsEn=new Double_t[fgkNEnBins+1];
335 for(Int_t i=0; i<=fgkNEnBins; i++) {
336 if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
337 if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
338 if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
341 fh3PtEtaPhiTracks = new TH3F("fh3PtEtaPhiTracks","fh3PtEtaPhiTracks;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
342 fOutput->Add(fh3PtEtaPhiTracks);
344 fh3PtEtaPhiTracksOnEmcal = new TH3F("fh3PtEtaPhiTracksOnEmcal","fh3PtEtaPhiTracksOnEmcal;#it{p}_{T}^{track}_{emc};#eta_{emc};#varphi_{emc}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
345 fOutput->Add(fh3PtEtaPhiTracksOnEmcal);
347 fh3PtEtaPhiTracksToProp = new TH3F("fh3PtEtaPhiTracksToProp","fh3PtEtaPhiTracksToProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
348 fOutput->Add(fh3PtEtaPhiTracksToProp);
350 fh3PtEtaPhiTracksProp = new TH3F("fh3PtEtaPhiTracksProp","fh3PtEtaPhiTracksProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
351 fOutput->Add(fh3PtEtaPhiTracksProp);
353 fh3PtEtaPhiTracksNoProp = new TH3F("fh3PtEtaPhiTracksNoProp","fh3PtEtaPhiTracksNoProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
354 fOutput->Add(fh3PtEtaPhiTracksNoProp);
356 fh2CentPtJetFull = new TH2F("fh2CentPtJetFull","fh2CentPtJetFull;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
357 fOutput->Add(fh2CentPtJetFull);
359 fh2CentPtJetCharged = new TH2F("fh2CentPtJetCharged","fh2CentPtJetCharged;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
360 fOutput->Add(fh2CentPtJetCharged);
362 fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
363 fOutput->Add(fh3PtEtaPhiJetFull);
365 fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
366 fOutput->Add(fh3PtEtaPhiJetCharged);
368 fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
369 fOutput->Add(fh2NJetsPtFull);
371 fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
372 fOutput->Add(fh2NJetsPtCharged);
374 fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
375 fOutput->Add(fh3PtEtaAreaJetFull);
377 fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
378 fOutput->Add(fh3PtEtaAreaJetCharged);
380 fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
381 fOutput->Add(fh2PtNConstituentsCharged);
383 fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
384 fOutput->Add(fh2PtNConstituents);
386 fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
387 fOutput->Add(fh2PtMeanPtConstituentsCharged);
389 fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
390 fOutput->Add(fh2PtMeanPtConstituentsNeutral);
392 fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",fgkNPtBins,binsPt,fgkNNEFBins,binsNEF);
393 fOutput->Add(fh2PtNEF);
395 fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",fgkNNEFBins,binsNEF,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
396 fOutput->Add(fh3NEFEtaPhi);
398 fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
399 fOutput->Add(fh2NEFNConstituentsCharged);
401 fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
402 fOutput->Add(fh2NEFNConstituentsNeutral);
404 fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
405 fOutput->Add(fh2Ptz);
407 fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
408 fOutput->Add(fh2PtzCharged);
410 fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",fgkNPtBins,binsPt,fgkNPtBins,binsPt);
411 fOutput->Add(fh2PtLeadJet1VsLeadJet2);
413 fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
414 fOutput->Add(fh3EEtaPhiCluster);
416 fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
417 fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
418 fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
419 fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
421 fh3PatchEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1","fh3PatchEnergyEtaPhiCenterJ1;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
422 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1);
424 fh3PatchEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ2","fh3PatchEnergyEtaPhiCenterJ2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
425 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ2);
427 fh3PatchEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1J2","fh3PatchEnergyEtaPhiCenterJ1J2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
428 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1J2);
430 fh3PatchADCEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1","fh3PatchADCEnergyEtaPhiCenterJ1;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
431 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1);
433 fh3PatchADCEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ2","fh3PatchADCEnergyEtaPhiCenterJ2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
434 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ2);
436 fh3PatchADCEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1J2","fh3PatchADCEnergyEtaPhiCenterJ1J2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
437 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1J2);
439 fh3PatchADCEnergyEtaPhiCenterAll = new TH3F("fh3PatchADCEnergyEtaPhiCenterAll","fh3PatchADCEnergyEtaPhiCenterAll;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
440 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterAll);
442 fh3EEtaPhiCell = new TH3F("fh3EEtaPhiCell","fh3EEtaPhiCell;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
443 fOutput->Add(fh3EEtaPhiCell);
445 fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
446 fOutput->Add(fh2CellEnergyVsTime);
448 fh3EClusELeadingCellVsTime = new TH3F("fh3EClusELeadingCellVsTime","fh3EClusELeadingCellVsTime;E_{cluster};E_{leading cell};time_{leading cell}",fgkNEnBins,binsEn,fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
449 fOutput->Add(fh3EClusELeadingCellVsTime);
451 fh3JetReacCent = new TH3F("fh3JetReacCent","fh3JetReacCent;E_{Jet};Centrality;dEP",fgkNEnBins,binsEn,fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
452 fOutput->Add(fh3JetReacCent);
454 // =========== Switch on Sumw2 for all histos ===========
455 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
456 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
461 TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
466 TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
471 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
475 TH1::AddDirectory(oldStatus);
477 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
479 if(binsCent) delete [] binsCent;
480 if(binsdEP) delete [] binsdEP;
481 if(binsEn) delete [] binsEn;
482 if(binsPt) delete [] binsPt;
483 if(binsPhi) delete [] binsPhi;
484 if(binsEta) delete [] binsEta;
485 if(binsArea) delete [] binsArea;
486 if(binsConst) delete [] binsConst;
487 if(binsMeanPt) delete [] binsMeanPt;
488 if(binsNEF) delete [] binsNEF;
489 if(binsz) delete [] binsz;
490 if(binsJetType) delete [] binsJetType;
491 if(binsTime) delete [] binsTime;
495 //________________________________________________________________________
496 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
501 AliParticleContainer *partCont = GetParticleContainer(0);
503 AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle(0));
505 Double_t trkphi = track->Phi()*TMath::RadToDeg();
506 fh3PtEtaPhiTracks->Fill(track->Pt(),track->Eta(),track->Phi());
507 //Select tracks which should be propagated
508 if(track->Pt()>=0.350) {
509 if (TMath::Abs(track->Eta())<=0.9 && trkphi > 10 && trkphi < 250) {
510 fh3PtEtaPhiTracksOnEmcal->Fill(track->GetTrackPtOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPhiOnEMCal());
511 fh3PtEtaPhiTracksToProp->Fill(track->Pt(),track->Eta(),track->Phi());
512 if(track->GetTrackPtOnEMCal()>=0)
513 fh3PtEtaPhiTracksProp->Fill(track->Pt(),track->Eta(),track->Phi());
515 fh3PtEtaPhiTracksNoProp->Fill(track->Pt(),track->Eta(),track->Phi());
518 track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle());
523 AliClusterContainer *clusCont = GetClusterContainer(0);
525 Int_t nclusters = clusCont->GetNClusters();
526 for (Int_t ic = 0; ic < nclusters; ic++) {
527 AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
529 AliDebug(2,Form("Could not receive cluster %d", ic));
532 if (!cluster->IsEMCAL()) {
533 AliDebug(2,Form("%s: Cluster is not emcal",GetName()));
538 cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
539 fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
541 Double_t leadCellE = GetEnergyLeadingCell(cluster);
542 Double_t leadCellT = cluster->GetTOF();
543 fh3EClusELeadingCellVsTime->Fill(lp.E(),leadCellE,leadCellT*1e9);
550 const Short_t nCells = fCaloCells->GetNumberOfCells();
552 for(Int_t iCell=0; iCell<nCells; ++iCell) {
553 Short_t cellId = fCaloCells->GetCellNumber(iCell);
554 Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
555 Double_t cellT = fCaloCells->GetCellTime(cellId);
557 fGeom->GetGlobal(cellId, pos);
558 TLorentzVector lv(pos,cellE);
559 Double_t cellEta = lv.Eta();
560 Double_t cellPhi = lv.Phi();
561 if(cellPhi<0.) cellPhi+=TMath::TwoPi();
562 if(cellPhi>TMath::TwoPi()) cellPhi-=TMath::TwoPi();
564 AliDebug(2,Form("cell energy = %f time = %f",cellE,cellT*1e9));
565 fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
566 fh3EEtaPhiCell->Fill(cellE,cellEta,cellPhi);
571 Double_t ptLeadJet1 = 0.;
572 Double_t ptLeadJet2 = 0.;
574 fHistRhovsCentFull->Fill(fCent,GetRhoVal(fContainerFull));
575 fHistRhovsCentCharged->Fill(fCent,GetRhoVal(fContainerCharged));
577 TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
579 nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
581 if (GetJetContainer(fContainerFull)) {
582 const Int_t njets = GetNJets(fContainerFull);
583 for (Int_t ij = 0; ij < njets; ij++) {
585 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
587 continue; //jet not selected
589 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerFull)*jet->Area();
590 if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
592 Double_t dEPJetFull = RelativeEP(jet->Phi() , fEPV0);
593 fh3JetReacCent->Fill(jet->E(),fCent,dEPJetFull);
595 fh2CentPtJetFull->Fill(fCent,jetPt);
596 fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
597 fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
599 //count jets above certain pT threshold
600 Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
601 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
602 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
604 fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
605 fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
606 fh2PtNEF->Fill(jetPt,jet->NEF());
607 fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
608 fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
609 fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
612 Double_t sumPtCh = 0.;
613 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
614 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
616 fh2Ptz->Fill(jetPt,GetZ(vp,jet));
620 if(jet->GetNumberOfTracks()>0)
621 fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
623 AliVCluster *vc = 0x0;
624 Double_t sumPtNe = 0.;
626 for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
627 vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
630 vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
634 if(jet->GetNumberOfClusters()>0)
635 fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
639 for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
640 Int_t nJetsInEvent = nJetsArr->At(i);
641 fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
645 //Reset array to zero to also count charged jets
648 //Loop over charged jets
649 if (GetJetContainer(fContainerCharged)) {
650 const Int_t njets = GetNJets(fContainerCharged);
651 for (Int_t ij = 0; ij < njets; ij++) {
653 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
655 continue; //jet not selected
657 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
658 if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
659 fh2CentPtJetCharged->Fill(fCent,jetPt);
660 fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
661 fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
664 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
665 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
667 fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
670 //count jets above certain pT threshold
671 Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
672 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
673 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
676 for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
677 Int_t nJetsInEvent = nJetsArr->At(i);
678 fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
682 if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged)) {
683 fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
686 fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
687 fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
689 if(nJetsArr) delete nJetsArr;
694 //________________________________________________________________________
695 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
697 // Run analysis code here, if needed. It will be executed before FillHistograms().
699 //Check if event is selected (vertex & pile-up)
703 if(fTriggerPatchInfo)
704 FillTriggerPatchHistos();
706 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
709 //_______________________________________________________________________
710 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
712 // Called once at the end of the analysis.
714 //________________________________________________________________________
715 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet) const
717 // Get Z of constituent trk
718 return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
721 //________________________________________________________________________
722 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
725 // Get the z of a constituent inside of a jet
727 Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
729 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
731 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
736 //________________________________________________________________________
737 Int_t AliAnalysisTaskEmcalJetTriggerQA::GetLeadingCellId(const AliVCluster *clus) const
739 //Get energy of leading cell in cluster
745 Int_t iCellAbsIdMax = -1;
746 Int_t nCells = clus->GetNCells();
747 for(Int_t i = 0; i<nCells; i++) {
748 Int_t absId = clus->GetCellAbsId(i);
749 Double_t cellE = fCaloCells->GetCellAmplitude(absId);
752 iCellAbsIdMax = absId;
756 return iCellAbsIdMax;
759 //________________________________________________________________________
760 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetEnergyLeadingCell(const AliVCluster *clus) const
762 //Get energy of leading cell in cluster
766 Int_t absID = GetLeadingCellId(clus);
768 return fCaloCells->GetCellAmplitude(absID);
773 //________________________________________________________________________
774 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetECross(Int_t absID) const {
776 //Get Ecross = sum of energy of neighbouring cells (using uncalibrated energy)
781 Double_t ecross = -1.;
788 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
789 fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
790 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
792 if( iphi < AliEMCALGeoParams::fgkEMCALRows-1)
793 absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
795 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
797 if( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) {
798 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
799 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
801 else if( ieta == 0 && imod%2 ) {
802 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
803 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
806 if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
807 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
809 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
812 Double_t ecell1 = fCaloCells->GetCellAmplitude(absID1);
813 Double_t ecell2 = fCaloCells->GetCellAmplitude(absID2);
814 Double_t ecell3 = fCaloCells->GetCellAmplitude(absID3);
815 Double_t ecell4 = fCaloCells->GetCellAmplitude(absID4);
817 ecross = ecell1+ecell2+ecell3+ecell4;
822 //_________________________________________________________________________
823 Float_t AliAnalysisTaskEmcalJetTriggerQA::RelativeEP(Double_t objAng, Double_t EPAng) const
825 // function to calculate angle between object and EP in the 1st quadrant (0,Pi/2)
826 Double_t dphi = EPAng - objAng;
828 // ran into trouble with a few dEP<-Pi so trying this...
829 if( dphi<-1*TMath::Pi() )
830 dphi = dphi + 1*TMath::Pi();
831 if( dphi>1*TMath::Pi())
832 dphi = dphi - 1*TMath::Pi();
834 if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
835 // Do nothing! we are in quadrant 1
836 }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
837 dphi = 1*TMath::Pi() - dphi;
838 }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
840 }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
841 dphi = dphi + 1*TMath::Pi();
844 return dphi; // dphi in [0, Pi/2]