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),
85 fh2CellEnergyVsTime(0),
86 fh3EClusELeadingCellVsTime(0),
89 // Default constructor.
91 SetMakeGeneralHistograms(kTRUE);
94 //________________________________________________________________________
95 AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
96 AliAnalysisTaskEmcalJet(name, kTRUE),
100 fContainerCharged(1),
106 fHistRhovsCentFull(0),
107 fHistRhovsCentCharged(0),
108 fh3PtEtaPhiTracks(0),
109 fh3PtEtaPhiTracksOnEmcal(0),
110 fh3PtEtaPhiTracksToProp(0),
111 fh3PtEtaPhiTracksProp(0),
112 fh3PtEtaPhiTracksNoProp(0),
114 fh2CentPtJetCharged(0),
115 fh3PtEtaPhiJetFull(0),
116 fh3PtEtaPhiJetCharged(0),
118 fh2NJetsPtCharged(0),
119 fh3PtEtaAreaJetFull(0),
120 fh3PtEtaAreaJetCharged(0),
121 fh2PtNConstituentsCharged(0),
122 fh2PtNConstituents(0),
123 fh2PtMeanPtConstituentsCharged(0),
124 fh2PtMeanPtConstituentsNeutral(0),
127 fh2NEFNConstituentsCharged(0),
128 fh2NEFNConstituentsNeutral(0),
131 fh2PtLeadJet1VsLeadJet2(0),
132 fh3EEtaPhiCluster(0),
133 fh3PtLeadJet1VsPatchEnergy(0),
134 fh3PtLeadJet2VsPatchEnergy(0),
135 fh3PatchEnergyEtaPhiCenterJ1(0),
136 fh3PatchEnergyEtaPhiCenterJ2(0),
137 fh3PatchEnergyEtaPhiCenterJ1J2(0),
138 fh3PatchADCEnergyEtaPhiCenterJ1(0),
139 fh3PatchADCEnergyEtaPhiCenterJ2(0),
140 fh3PatchADCEnergyEtaPhiCenterJ1J2(0),
141 fh3PatchADCEnergyEtaPhiCenterAll(0),
144 fh2CellEnergyVsTime(0),
145 fh3EClusELeadingCellVsTime(0),
148 // Standard constructor.
150 SetMakeGeneralHistograms(kTRUE);
153 //________________________________________________________________________
154 AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
160 //________________________________________________________________________
161 Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
163 // Decide if event should be selected for analysis
166 fhNEvents->Fill(3.5);
168 if(!fTriggerClass.IsNull()) {
169 //Check if requested trigger was fired
170 TString trigType1 = "J1";
171 TString trigType2 = "J2";
172 if(fTriggerClass.Contains("G")) {
177 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
178 if(fTriggerClass.Contains(trigType1.Data()) && fTriggerClass.Contains(trigType2.Data())) { //if events with J1&&J2 are requested
179 if(!firedTrigClass.Contains(trigType1.Data()) || !firedTrigClass.Contains(trigType2.Data()) ) //check if both are fired
183 if(!firedTrigClass.Contains(fTriggerClass))
185 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
189 fhNEvents->Fill(1.5);
194 //________________________________________________________________________
195 void AliAnalysisTaskEmcalJetTriggerQA::FillTriggerPatchHistos() {
197 //Fill trigger patch histos for main trigger
199 AliEmcalTriggerPatchInfo *patch = GetMainTriggerPatch();
201 fMaxPatchEnergy = patch->GetPatchE();
202 Double_t patchADCGeV = patch->GetADCAmpGeVRough();
203 fh3PatchADCEnergyEtaPhiCenterAll->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
204 if(patch->IsJetLow() && !patch->IsJetHigh()) { //main patch only fired low threshold trigger
205 fh3PatchEnergyEtaPhiCenterJ2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
206 fh3PatchADCEnergyEtaPhiCenterJ2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
208 else if(patch->IsJetHigh() && !patch->IsJetLow()) { //main patch only fired high threshold trigger - should never happen
209 fh3PatchEnergyEtaPhiCenterJ1->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
210 fh3PatchADCEnergyEtaPhiCenterJ1->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
212 else if(patch->IsJetHigh() && patch->IsJetLow()) { //main patch fired both triggers
213 fh3PatchEnergyEtaPhiCenterJ1J2->Fill(patch->GetPatchE(),patch->GetEtaGeo(),patch->GetPhiGeo());
214 fh3PatchADCEnergyEtaPhiCenterJ1J2->Fill(patchADCGeV,patch->GetEtaGeo(),patch->GetPhiGeo());
219 //________________________________________________________________________
220 void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
222 // Create user output.
224 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
226 Bool_t oldStatus = TH1::AddDirectoryStatus();
227 TH1::AddDirectory(kFALSE);
229 fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
230 fOutput->Add(fhNEvents);
232 fhTriggerbit = new TProfile("fhTriggerbit","fhTriggerbit;;TriggerBit",1,0,1);
233 fOutput->Add(fhTriggerbit);
235 fHistRhovsCentFull = new TH2F("fHistRhovsCentFull", "fHistRhovsCentFull", 101, -1, 100, 300, 0., 300.);
236 fHistRhovsCentFull->GetXaxis()->SetTitle("Centrality (%)");
237 fHistRhovsCentFull->GetYaxis()->SetTitle("s#rho_{ch} (GeV/c * rad^{-1})");
238 fOutput->Add(fHistRhovsCentFull);
240 fHistRhovsCentCharged = new TH2F("fHistRhovsCentCharged", "fHistRhovsCentCharged", 101, -1, 100, 300, 0., 300.);
241 fHistRhovsCentCharged->GetXaxis()->SetTitle("Centrality (%)");
242 fHistRhovsCentCharged->GetYaxis()->SetTitle("#rho_{ch} (GeV/c * rad^{-1})");
243 fOutput->Add(fHistRhovsCentCharged);
245 Int_t fgkNCentBins = 21;
246 Float_t kMinCent = 0.;
247 Float_t kMaxCent = 105.;
248 Double_t *binsCent = new Double_t[fgkNCentBins+1];
249 for(Int_t i=0; i<=fgkNCentBins; i++) binsCent[i]=(Double_t)kMinCent + (kMaxCent-kMinCent)/fgkNCentBins*(Double_t)i ;
250 binsCent[fgkNCentBins-1] = 100.5;
251 binsCent[fgkNCentBins] = 101.5;
253 Int_t fgkNdEPBins = 18*8;
254 Float_t kMindEP = 0.;
255 Float_t kMaxdEP = 1.*TMath::Pi()/2.;
256 Double_t *binsdEP = new Double_t[fgkNdEPBins+1];
257 for(Int_t i=0; i<=fgkNdEPBins; i++) binsdEP[i]=(Double_t)kMindEP + (kMaxdEP-kMindEP)/fgkNdEPBins*(Double_t)i ;
259 Int_t fgkNPtBins = 200;
260 Float_t kMinPt = -50.;
261 Float_t kMaxPt = 150.;
262 Double_t *binsPt = new Double_t[fgkNPtBins+1];
263 for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
265 Int_t fgkNPhiBins = 18*8;
266 Float_t kMinPhi = 0.;
267 Float_t kMaxPhi = 2.*TMath::Pi();
268 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
269 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
271 Int_t fgkNEtaBins = 100;
272 Float_t fgkEtaMin = -1.;
273 Float_t fgkEtaMax = 1.;
274 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
275 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
277 Int_t fgkNAreaBins = 100;
278 Float_t kMinArea = 0.;
279 Float_t kMaxArea = 1.;
280 Double_t *binsArea = new Double_t[fgkNAreaBins+1];
281 for(Int_t i=0; i<=fgkNAreaBins; i++) binsArea[i]=(Double_t)kMinArea + (kMaxArea-kMinArea)/fgkNAreaBins*(Double_t)i ;
283 Int_t fgkNConstBins = 100;
284 Float_t kMinConst = 0.;
285 Float_t kMaxConst = 100.;
286 Double_t *binsConst = new Double_t[fgkNConstBins+1];
287 for(Int_t i=0; i<=fgkNConstBins; i++) binsConst[i]=(Double_t)kMinConst + (kMaxConst-kMinConst)/fgkNConstBins*(Double_t)i ;
289 Int_t fgkNMeanPtBins = 100;
290 Float_t kMinMeanPt = 0.;
291 Float_t kMaxMeanPt = 20.;
292 Double_t *binsMeanPt = new Double_t[fgkNMeanPtBins+1];
293 for(Int_t i=0; i<=fgkNMeanPtBins; i++) binsMeanPt[i]=(Double_t)kMinMeanPt + (kMaxMeanPt-kMinMeanPt)/fgkNMeanPtBins*(Double_t)i ;
295 Int_t fgkNNEFBins = 101;
296 Float_t kMinNEF = 0.;
297 Float_t kMaxNEF = 1.01;
298 Double_t *binsNEF = new Double_t[fgkNNEFBins+1];
299 for(Int_t i=0; i<=fgkNNEFBins; i++) binsNEF[i]=(Double_t)kMinNEF + (kMaxNEF-kMinNEF)/fgkNNEFBins*(Double_t)i ;
301 Int_t fgkNzBins = 101;
303 Float_t kMaxz = 1.01;
304 Double_t *binsz = new Double_t[fgkNzBins+1];
305 for(Int_t i=0; i<=fgkNzBins; i++) binsz[i]=(Double_t)kMinz + (kMaxz-kMinz)/fgkNzBins*(Double_t)i ;
307 Int_t fgkNJetTypeBins = 2;
308 Float_t kMinJetType = -0.5;
309 Float_t kMaxJetType = 1.5;
310 Double_t *binsJetType = new Double_t[fgkNJetTypeBins+1];
311 for(Int_t i=0; i<=fgkNJetTypeBins; i++) binsJetType[i]=(Double_t)kMinJetType + (kMaxJetType-kMinJetType)/fgkNJetTypeBins*(Double_t)i ;
313 Int_t fgkNTimeBins = 100;
314 Float_t kMinTime = -200.;
315 Float_t kMaxTime = 200;
316 Double_t *binsTime = new Double_t[fgkNTimeBins+1];
317 for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
319 Double_t enBinEdges[3][2];
320 enBinEdges[0][0] = 1.; //10 bins
321 enBinEdges[0][1] = 0.1;
322 enBinEdges[1][0] = 5.; //8 bins
323 enBinEdges[1][1] = 0.5;
324 enBinEdges[2][0] = 100.;//95 bins
325 enBinEdges[2][1] = 1.;
327 const Float_t enmin1 = 0;
328 const Float_t enmax1 = enBinEdges[0][0];
329 const Float_t enmin2 = enmax1 ;
330 const Float_t enmax2 = enBinEdges[1][0];
331 const Float_t enmin3 = enmax2 ;
332 const Float_t enmax3 = enBinEdges[2][0];//fgkEnMax;
333 const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
334 const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
335 const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
337 Int_t fgkNEnBins=nbin13;
338 Double_t *binsEn=new Double_t[fgkNEnBins+1];
339 for(Int_t i=0; i<=fgkNEnBins; i++) {
340 if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
341 if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
342 if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
345 fh3PtEtaPhiTracks = new TH3F("fh3PtEtaPhiTracks","fh3PtEtaPhiTracks;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
346 fOutput->Add(fh3PtEtaPhiTracks);
348 fh3PtEtaPhiTracksOnEmcal = new TH3F("fh3PtEtaPhiTracksOnEmcal","fh3PtEtaPhiTracksOnEmcal;#it{p}_{T}^{track}_{emc};#eta_{emc};#varphi_{emc}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
349 fOutput->Add(fh3PtEtaPhiTracksOnEmcal);
351 fh3PtEtaPhiTracksToProp = new TH3F("fh3PtEtaPhiTracksToProp","fh3PtEtaPhiTracksToProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
352 fOutput->Add(fh3PtEtaPhiTracksToProp);
354 fh3PtEtaPhiTracksProp = new TH3F("fh3PtEtaPhiTracksProp","fh3PtEtaPhiTracksProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
355 fOutput->Add(fh3PtEtaPhiTracksProp);
357 fh3PtEtaPhiTracksNoProp = new TH3F("fh3PtEtaPhiTracksNoProp","fh3PtEtaPhiTracksNoProp;#it{p}_{T}^{track}_{vtx};#eta_{vtx};#varphi_{vtx}",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
358 fOutput->Add(fh3PtEtaPhiTracksNoProp);
360 fh2CentPtJetFull = new TH2F("fh2CentPtJetFull","fh2CentPtJetFull;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
361 fOutput->Add(fh2CentPtJetFull);
363 fh2CentPtJetCharged = new TH2F("fh2CentPtJetCharged","fh2CentPtJetCharged;cent;#it{p}_{T}^{jet}",fgkNCentBins,binsCent,fgkNPtBins,binsPt);
364 fOutput->Add(fh2CentPtJetCharged);
366 fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
367 fOutput->Add(fh3PtEtaPhiJetFull);
369 fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
370 fOutput->Add(fh3PtEtaPhiJetCharged);
372 fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
373 fOutput->Add(fh2NJetsPtFull);
375 fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,fgkNPtBins,binsPt);
376 fOutput->Add(fh2NJetsPtCharged);
378 fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
379 fOutput->Add(fh3PtEtaAreaJetFull);
381 fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",fgkNPtBins,binsPt,fgkNEtaBins,binsEta,fgkNAreaBins,binsArea);
382 fOutput->Add(fh3PtEtaAreaJetCharged);
384 fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
385 fOutput->Add(fh2PtNConstituentsCharged);
387 fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",fgkNPtBins,binsPt,fgkNConstBins,binsConst);
388 fOutput->Add(fh2PtNConstituents);
390 fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
391 fOutput->Add(fh2PtMeanPtConstituentsCharged);
393 fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",fgkNPtBins,binsPt,fgkNMeanPtBins,binsMeanPt);
394 fOutput->Add(fh2PtMeanPtConstituentsNeutral);
396 fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",fgkNPtBins,binsPt,fgkNNEFBins,binsNEF);
397 fOutput->Add(fh2PtNEF);
399 fh3NEFEtaPhi = new TH3F("fh3NEFEtaPhi","fh3NEFEtaPhi;NEF;#eta;#varphi",fgkNNEFBins,binsNEF,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
400 fOutput->Add(fh3NEFEtaPhi);
402 fh2NEFNConstituentsCharged = new TH2F("fh2NEFNConstituentsCharged","fh2NEFNConstituentsCharged;NEF;N_{charged constituents}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
403 fOutput->Add(fh2NEFNConstituentsCharged);
405 fh2NEFNConstituentsNeutral = new TH2F("fh2NEFNConstituentsNeutral","fh2NEFNConstituentsNeutral;NEF;N_{clusters}",fgkNNEFBins,binsNEF,fgkNConstBins,binsConst);
406 fOutput->Add(fh2NEFNConstituentsNeutral);
408 fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
409 fOutput->Add(fh2Ptz);
411 fh2PtzCharged = new TH2F("fh2PtzCharged","fh2Ptz;#it{p}_{T}^{ch jet};z=p_{t,trk}^{proj}/p_{ch jet}",fgkNPtBins,binsPt,fgkNzBins,binsz);
412 fOutput->Add(fh2PtzCharged);
414 fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",fgkNPtBins,binsPt,fgkNPtBins,binsPt);
415 fOutput->Add(fh2PtLeadJet1VsLeadJet2);
417 fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
418 fOutput->Add(fh3EEtaPhiCluster);
420 fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
421 fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
422 fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",fgkNPtBins,binsPt,fgkNPtBins,binsPt,fgkNJetTypeBins,binsJetType);
423 fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
425 fh3PatchEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1","fh3PatchEnergyEtaPhiCenterJ1;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
426 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1);
428 fh3PatchEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ2","fh3PatchEnergyEtaPhiCenterJ2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
429 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ2);
431 fh3PatchEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchEnergyEtaPhiCenterJ1J2","fh3PatchEnergyEtaPhiCenterJ1J2;E_{patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
432 fOutput->Add(fh3PatchEnergyEtaPhiCenterJ1J2);
434 fh3PatchADCEnergyEtaPhiCenterJ1 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1","fh3PatchADCEnergyEtaPhiCenterJ1;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
435 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1);
437 fh3PatchADCEnergyEtaPhiCenterJ2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ2","fh3PatchADCEnergyEtaPhiCenterJ2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
438 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ2);
440 fh3PatchADCEnergyEtaPhiCenterJ1J2 = new TH3F("fh3PatchADCEnergyEtaPhiCenterJ1J2","fh3PatchADCEnergyEtaPhiCenterJ1J2;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
441 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterJ1J2);
443 fh3PatchADCEnergyEtaPhiCenterAll = new TH3F("fh3PatchADCEnergyEtaPhiCenterAll","fh3PatchADCEnergyEtaPhiCenterAll;E_{ADC,patch};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
444 fOutput->Add(fh3PatchADCEnergyEtaPhiCenterAll);
446 fh3EEtaPhiCell = new TH3F("fh3EEtaPhiCell","fh3EEtaPhiCell;E_{cell};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
447 fOutput->Add(fh3EEtaPhiCell);
449 fh2ECellVsCent = new TH2F("fh2ECellVsCent","fh2ECellVsCent;centrality;E_{cell}",101,-1,100,500,0.,5.);
450 fOutput->Add(fh2ECellVsCent);
452 fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
453 fOutput->Add(fh2CellEnergyVsTime);
455 fh3EClusELeadingCellVsTime = new TH3F("fh3EClusELeadingCellVsTime","fh3EClusELeadingCellVsTime;E_{cluster};E_{leading cell};time_{leading cell}",fgkNEnBins,binsEn,fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
456 fOutput->Add(fh3EClusELeadingCellVsTime);
458 fh3JetReacCent = new TH3F("fh3JetReacCent","fh3JetReacCent;E_{Jet};Centrality;dEP",fgkNEnBins,binsEn,fgkNCentBins,binsCent,fgkNdEPBins,binsdEP);
459 fOutput->Add(fh3JetReacCent);
461 // =========== Switch on Sumw2 for all histos ===========
462 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
463 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
468 TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
473 TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
478 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
482 TH1::AddDirectory(oldStatus);
484 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
487 if(binsdEP) delete [] binsdEP;
488 if(binsEn) delete [] binsEn;
489 if(binsPt) delete [] binsPt;
490 if(binsPhi) delete [] binsPhi;
491 if(binsEta) delete [] binsEta;
492 if(binsArea) delete [] binsArea;
493 if(binsConst) delete [] binsConst;
494 if(binsMeanPt) delete [] binsMeanPt;
495 if(binsNEF) delete [] binsNEF;
496 if(binsz) delete [] binsz;
497 if(binsJetType) delete [] binsJetType;
498 if(binsTime) delete [] binsTime;
502 //________________________________________________________________________
503 Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
508 AliParticleContainer *partCont = GetParticleContainer(0);
510 AliPicoTrack *track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle(0));
512 Double_t trkphi = track->Phi()*TMath::RadToDeg();
513 fh3PtEtaPhiTracks->Fill(track->Pt(),track->Eta(),track->Phi());
514 //Select tracks which should be propagated
515 if(track->Pt()>=0.350) {
516 if (TMath::Abs(track->Eta())<=0.9 && trkphi > 10 && trkphi < 250) {
517 fh3PtEtaPhiTracksOnEmcal->Fill(track->GetTrackPtOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPhiOnEMCal());
518 fh3PtEtaPhiTracksToProp->Fill(track->Pt(),track->Eta(),track->Phi());
519 if(track->GetTrackPtOnEMCal()>=0)
520 fh3PtEtaPhiTracksProp->Fill(track->Pt(),track->Eta(),track->Phi());
522 fh3PtEtaPhiTracksNoProp->Fill(track->Pt(),track->Eta(),track->Phi());
525 track = dynamic_cast<AliPicoTrack*>(partCont->GetNextAcceptParticle());
530 AliClusterContainer *clusCont = GetClusterContainer(0);
532 Int_t nclusters = clusCont->GetNClusters();
533 for (Int_t ic = 0; ic < nclusters; ic++) {
534 AliVCluster *cluster = static_cast<AliVCluster*>(clusCont->GetCluster(ic));
536 AliDebug(2,Form("Could not receive cluster %d", ic));
539 if (!cluster->IsEMCAL()) {
540 AliDebug(2,Form("%s: Cluster is not emcal",GetName()));
545 cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
546 fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
548 Double_t leadCellE = GetEnergyLeadingCell(cluster);
549 Double_t leadCellT = cluster->GetTOF();
550 fh3EClusELeadingCellVsTime->Fill(lp.E(),leadCellE,leadCellT*1e9);
557 const Short_t nCells = fCaloCells->GetNumberOfCells();
559 for(Int_t iCell=0; iCell<nCells; ++iCell) {
560 Short_t cellId = fCaloCells->GetCellNumber(iCell);
561 Double_t cellE = fCaloCells->GetCellAmplitude(cellId);
562 Double_t cellT = fCaloCells->GetCellTime(cellId);
564 fGeom->GetGlobal(cellId, pos);
565 TLorentzVector lv(pos,cellE);
566 Double_t cellEta = lv.Eta();
567 Double_t cellPhi = lv.Phi();
568 if(cellPhi<0.) cellPhi+=TMath::TwoPi();
569 if(cellPhi>TMath::TwoPi()) cellPhi-=TMath::TwoPi();
571 AliDebug(2,Form("cell energy = %f time = %f",cellE,cellT*1e9));
572 fh2CellEnergyVsTime->Fill(cellE,cellT*1e9);
573 fh3EEtaPhiCell->Fill(cellE,cellEta,cellPhi);
574 fh2ECellVsCent->Fill(fCent,cellE);
579 Double_t ptLeadJet1 = 0.;
580 Double_t ptLeadJet2 = 0.;
582 fHistRhovsCentFull->Fill(fCent,GetRhoVal(fContainerFull));
583 fHistRhovsCentCharged->Fill(fCent,GetRhoVal(fContainerCharged));
585 TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
587 nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
589 if (GetJetContainer(fContainerFull)) {
590 const Int_t njets = GetNJets(fContainerFull);
591 for (Int_t ij = 0; ij < njets; ij++) {
593 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerFull);
595 continue; //jet not selected
597 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerFull)*jet->Area();
598 if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
600 Double_t dEPJetFull = RelativeEP(jet->Phi() , fEPV0);
601 fh3JetReacCent->Fill(jet->E(),fCent,dEPJetFull);
603 fh2CentPtJetFull->Fill(fCent,jetPt);
604 fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
605 fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
607 //count jets above certain pT threshold
608 Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
609 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
610 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
612 fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
613 fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
614 fh2PtNEF->Fill(jetPt,jet->NEF());
615 fh3NEFEtaPhi->Fill(jet->NEF(),jet->Eta(),jet->Phi());
616 fh2NEFNConstituentsCharged->Fill(jet->NEF(),jet->GetNumberOfTracks());
617 fh2NEFNConstituentsNeutral->Fill(jet->NEF(),jet->GetNumberOfClusters());
620 Double_t sumPtCh = 0.;
621 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
622 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
624 fh2Ptz->Fill(jetPt,GetZ(vp,jet));
628 if(jet->GetNumberOfTracks()>0)
629 fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
631 AliVCluster *vc = 0x0;
632 Double_t sumPtNe = 0.;
634 for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
635 vc = static_cast<AliVCluster*>(clusCont->GetCluster(icc));
638 vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
642 if(jet->GetNumberOfClusters()>0)
643 fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
647 for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
648 Int_t nJetsInEvent = nJetsArr->At(i);
649 fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
653 //Reset array to zero to also count charged jets
656 //Loop over charged jets
657 if (GetJetContainer(fContainerCharged)) {
658 const Int_t njets = GetNJets(fContainerCharged);
659 for (Int_t ij = 0; ij < njets; ij++) {
661 AliEmcalJet* jet = GetAcceptJetFromArray(ij,fContainerCharged);
663 continue; //jet not selected
665 Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
666 if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
667 fh2CentPtJetCharged->Fill(fCent,jetPt);
668 fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
669 fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
672 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
673 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
675 fh2PtzCharged->Fill(jetPt,GetZ(vp,jet));
678 //count jets above certain pT threshold
679 Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
680 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
681 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
684 for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
685 Int_t nJetsInEvent = nJetsArr->At(i);
686 fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
690 if(GetJetContainer(fContainerFull) && GetJetContainer(fContainerCharged))
691 fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
693 fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
694 fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
701 //________________________________________________________________________
702 Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
704 // Run analysis code here, if needed. It will be executed before FillHistograms().
706 fhTriggerbit->Fill(0.5,GetCollisionCandidates());
708 //Check if event is selected (vertex & pile-up)
712 if(fTriggerPatchInfo)
713 FillTriggerPatchHistos();
715 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
718 //_______________________________________________________________________
719 void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
721 // Called once at the end of the analysis.
723 //________________________________________________________________________
724 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet) const
726 // Get Z of constituent trk
727 return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
730 //________________________________________________________________________
731 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(Double_t trkPx, Double_t trkPy, Double_t trkPz, Double_t jetPx, Double_t jetPy, Double_t jetPz) const
734 // Get the z of a constituent inside of a jet
736 Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
738 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
740 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
745 //________________________________________________________________________
746 Int_t AliAnalysisTaskEmcalJetTriggerQA::GetLeadingCellId(const AliVCluster *clus) const
748 //Get energy of leading cell in cluster
754 Int_t iCellAbsIdMax = -1;
755 Int_t nCells = clus->GetNCells();
756 for(Int_t i = 0; i<nCells; i++) {
757 Int_t absId = clus->GetCellAbsId(i);
758 Double_t cellE = fCaloCells->GetCellAmplitude(absId);
761 iCellAbsIdMax = absId;
764 return iCellAbsIdMax;
767 //________________________________________________________________________
768 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetEnergyLeadingCell(const AliVCluster *clus) const
770 //Get energy of leading cell in cluster
774 Int_t absID = GetLeadingCellId(clus);
776 return fCaloCells->GetCellAmplitude(absID);
781 //________________________________________________________________________
782 Double_t AliAnalysisTaskEmcalJetTriggerQA::GetECross(Int_t absID) const {
784 //Get Ecross = sum of energy of neighbouring cells (using uncalibrated energy)
789 Double_t ecross = -1.;
796 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
797 fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta);
798 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
800 if( iphi < AliEMCALGeoParams::fgkEMCALRows-1)
801 absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
803 absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
805 if( ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2) ) {
806 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
807 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
809 else if( ieta == 0 && imod%2 ) {
810 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
811 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
814 if( ieta < AliEMCALGeoParams::fgkEMCALCols-1 )
815 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
817 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
820 Double_t ecell1 = fCaloCells->GetCellAmplitude(absID1);
821 Double_t ecell2 = fCaloCells->GetCellAmplitude(absID2);
822 Double_t ecell3 = fCaloCells->GetCellAmplitude(absID3);
823 Double_t ecell4 = fCaloCells->GetCellAmplitude(absID4);
825 ecross = ecell1+ecell2+ecell3+ecell4;
830 //_________________________________________________________________________
831 Float_t AliAnalysisTaskEmcalJetTriggerQA::RelativeEP(Double_t objAng, Double_t EPAng) const
833 // function to calculate angle between object and EP in the 1st quadrant (0,Pi/2)
834 Double_t dphi = EPAng - objAng;
836 // ran into trouble with a few dEP<-Pi so trying this...
837 if( dphi<-1*TMath::Pi() )
838 dphi = dphi + 1*TMath::Pi();
839 if( dphi>1*TMath::Pi())
840 dphi = dphi - 1*TMath::Pi();
842 if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
843 // Do nothing! we are in quadrant 1
844 }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
845 dphi = 1*TMath::Pi() - dphi;
846 }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
848 }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
849 dphi = dphi + 1*TMath::Pi();
852 return dphi; // dphi in [0, Pi/2]