]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTriggerQA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Possible set of pt weight flag as macro parameter (R.Bala)
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetTriggerQA.cxx
CommitLineData
0531c107 1//
2// Jet trigger QA analysis task.
3//
4// Author: M.Verweij
5
6#include <TClonesArray.h>
7#include <TH1F.h>
8#include <TH2F.h>
9#include <TH3F.h>
10#include <THnSparse.h>
11#include <TList.h>
12#include <TLorentzVector.h>
13
14#include "AliVCluster.h"
15#include "AliVTrack.h"
16#include "AliEmcalJet.h"
17#include "AliRhoParameter.h"
18#include "AliLog.h"
19#include "AliAnalysisUtils.h"
20#include "AliEmcalParticle.h"
21#include "AliAODCaloTrigger.h"
22#include "AliEMCALGeometry.h"
23#include "AliVCaloCells.h"
24
25
26#include "AliAnalysisTaskEmcalJetTriggerQA.h"
27
28ClassImp(AliAnalysisTaskEmcalJetTriggerQA)
29
30//________________________________________________________________________
31AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA() :
32 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTriggerQA", kTRUE),
33 fDebug(kFALSE),
34 fUseAnaUtils(kTRUE),
35 fAnalysisUtils(0),
609c7a0d 36 fJetsName2("Jet_AKTChargedR040_PicoTracks_pT0150_caloClustersCorr_ET0300"),
37 fJets2(0),
38 fEtaMinJet2(0.5),
39 fEtaMaxJet2(0.5),
40 fPhiMinJet2(-10.),
41 fPhiMaxJet2(10.),
42 fMaxTrackPtJet2(100.),
0531c107 43 fRhoChName(""),
44 fRhoCh(0),
45 fRhoChVal(0),
46 fTriggerClass(""),
47 fBitJ1((1<<8)),
48 fBitJ2((1<<11)),
49 fMaxPatchEnergy(0),
50 fTriggerType(-1),
51 fNFastOR(16),
52 fhNEvents(0),
53 fh3PtEtaPhiJetFull(0),
54 fh3PtEtaPhiJetCharged(0),
55 fh2NJetsPtFull(0),
56 fh2NJetsPtCharged(0),
57 fh3PtEtaAreaJetFull(0),
58 fh3PtEtaAreaJetCharged(0),
59 fh2PtNConstituentsCharged(0),
60 fh2PtNConstituents(0),
61 fh2PtMeanPtConstituentsCharged(0),
62 fh2PtMeanPtConstituentsNeutral(0),
63 fh2PtNEF(0),
64 fh2Ptz(0),
65 fh2PtLeadJet1VsLeadJet2(0),
66 fh3EEtaPhiCluster(0),
67 fh3PtLeadJet1VsPatchEnergy(0),
68 fh3PtLeadJet2VsPatchEnergy(0),
69 fh3PatchEnergyEtaPhiCenter(0)
70{
71 // Default constructor.
72
73
74 SetMakeGeneralHistograms(kTRUE);
75}
76
77//________________________________________________________________________
78AliAnalysisTaskEmcalJetTriggerQA::AliAnalysisTaskEmcalJetTriggerQA(const char *name) :
79 AliAnalysisTaskEmcalJet(name, kTRUE),
80 fDebug(kFALSE),
81 fUseAnaUtils(kTRUE),
82 fAnalysisUtils(0),
609c7a0d 83 fJetsName2("Jet_AKTChargedR040_PicoTracks_pT0150_caloClustersCorr_ET0300"),
84 fJets2(0),
85 fEtaMinJet2(0.5),
86 fEtaMaxJet2(0.5),
87 fPhiMinJet2(-10.),
88 fPhiMaxJet2(10.),
89 fMaxTrackPtJet2(100.),
0531c107 90 fRhoChName(""),
91 fRhoCh(0),
92 fRhoChVal(0),
93 fTriggerClass(""),
94 fBitJ1((1<<8)),
95 fBitJ2((1<<11)),
96 fMaxPatchEnergy(0),
97 fTriggerType(-1),
98 fNFastOR(16),
99 fhNEvents(0),
100 fh3PtEtaPhiJetFull(0),
101 fh3PtEtaPhiJetCharged(0),
102 fh2NJetsPtFull(0),
103 fh2NJetsPtCharged(0),
104 fh3PtEtaAreaJetFull(0),
105 fh3PtEtaAreaJetCharged(0),
106 fh2PtNConstituentsCharged(0),
107 fh2PtNConstituents(0),
108 fh2PtMeanPtConstituentsCharged(0),
109 fh2PtMeanPtConstituentsNeutral(0),
110 fh2PtNEF(0),
111 fh2Ptz(0),
112 fh2PtLeadJet1VsLeadJet2(0),
113 fh3EEtaPhiCluster(0),
114 fh3PtLeadJet1VsPatchEnergy(0),
115 fh3PtLeadJet2VsPatchEnergy(0),
116 fh3PatchEnergyEtaPhiCenter(0)
117{
118 // Standard constructor.
119
120 SetMakeGeneralHistograms(kTRUE);
121}
122
123//________________________________________________________________________
124AliAnalysisTaskEmcalJetTriggerQA::~AliAnalysisTaskEmcalJetTriggerQA()
125{
126 // Destructor.
127}
128
129//----------------------------------------------------------------------------------------------
130void AliAnalysisTaskEmcalJetTriggerQA::InitOnce() {
131 //
132 // only initialize once
133 //
134
135 // Initialize analysis util class for vertex selection
136 if(fUseAnaUtils) {
137 fAnalysisUtils = new AliAnalysisUtils();
138 fAnalysisUtils->SetMinVtxContr(2);
139 fAnalysisUtils->SetMaxVtxZ(10.);
140 }
141}
142
143//________________________________________________________________________
144Bool_t AliAnalysisTaskEmcalJetTriggerQA::SelectEvent() {
145 //
146 // Decide if event should be selected for analysis
147 //
148
149 fhNEvents->Fill(0.5);
150
151 if(fAnalysisUtils) {
152 if(!fAnalysisUtils->IsVertexSelected2013pA(InputEvent()))
153 return kFALSE;
154
155 fhNEvents->Fill(2.5);
156
157 if(fAnalysisUtils->IsPileUpEvent(InputEvent()))
158 return kFALSE;
159 }
160 else{
161 if(fUseAnaUtils)
162 AliError("fAnalysisUtils not initialized. Call AliAnalysisTaskEmcalJetTriggerQA::InitOnce()");
163 }
164
165 fhNEvents->Fill(3.5);
166
167 if(!fTriggerClass.IsNull()) {
168 //Check if requested trigger was fired
169 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
170 if(!firedTrigClass.Contains(fTriggerClass))
171 return kFALSE;
172 }
173
174 fhNEvents->Fill(1.5);
175
176 return kTRUE;
177
178}
179
180//________________________________________________________________________
181void AliAnalysisTaskEmcalJetTriggerQA::FindTriggerPatch() {
182
183 //Code to get position of trigger
184
185 if(!fGeom)
186 fGeom = AliEMCALGeometry::GetInstance();
187
188
189 TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
190
191 AliAODCaloTrigger *trg = dynamic_cast<AliAODCaloTrigger*>(InputEvent()->GetCaloTrigger("EMCAL"));
192 trg->Reset();
193
194 int col, row; //FASTOR position
195
196 Int_t nPatchNotEmpty = 0; //counter number of patches which are not empty
197 fMaxPatchEnergy = 0.;
198 while (trg->Next())
199 {
200
201 trg->GetPosition(col, row); //col (0 to 63), row (0 to 47)
202
203 if (col > -1 && row > -1)
204 {
205 Int_t id = -1; //FASTOR index
206 Int_t cellIndex[4] = {-1};
207 fGeom->GetAbsFastORIndexFromPositionInEMCAL(col,row,id); //phi is column, eta is row
208 if(id<0) continue;
209
210 fGeom->GetCellIndexFromFastORIndex(id,cellIndex);
211
212 Int_t ts = 0;
213 trg->GetL1TimeSum(ts);
214
215 Float_t ampTrg = 0.;
216 trg->GetAmplitude(ampTrg);
217
218 //L1 analysis
219 Int_t bit = 0;
220 trg->GetTriggerBits(bit);
221
222 Bool_t bTrigJ = TestFilterBit(bit,fBitJ1|fBitJ2);
223 if(!bTrigJ)
224 continue;
225
226 Bool_t bTrigJ1 = TestFilterBit(bit,fBitJ1);
227 Bool_t bTrigJ2 = TestFilterBit(bit,fBitJ2);
228
229 if(bTrigJ1) fTriggerType = 0;
230 if(bTrigJ2) fTriggerType = 1;
231 if(!bTrigJ1 && !bTrigJ2) fTriggerType = -1;
232
233 Double_t minPhiPatch = 10.;
234 Double_t maxPhiPatch = -10.;
235 Double_t minEtaPatch = 10.;
236 Double_t maxEtaPatch = -10.;
237
238 //Get energy in trigger patch 8x8 FASTOR
239 Double_t patchEnergy = 0.;
240 Double_t sumAmp = 0.;
241 // const Int_t nFastOR = 8;//16;
242 for(Int_t fastrow = 0; fastrow<fNFastOR; fastrow++) {
243 for(Int_t fastcol = 0; fastcol<fNFastOR; fastcol++) {
244 Int_t nrow = row+fastrow;
245 Int_t ncol = col+fastcol;
246 fGeom->GetAbsFastORIndexFromPositionInEMCAL(ncol,nrow,id);
247
248 if(id<0) {
249 AliWarning(Form("%s: id smaller than 0 %d",GetName(),id));
250 continue;
251 }
252
253 fGeom->GetCellIndexFromFastORIndex(id,cellIndex);
254 for(Int_t icell=0; icell<4; icell++) {
255
256 if(!fGeom->CheckAbsCellId(cellIndex[icell])) continue;
257
258 Double_t amp =0., time = 0., efrac = 0;
259 Int_t mclabel = -1;
260 Short_t absId = -1;
261 Int_t nSupMod = -1, nModule = -1, nIphi = -1, nIeta = -1;
262
263 fGeom->GetCellIndex(cellIndex[icell], nSupMod, nModule, nIphi, nIeta );
264
265 fCaloCells->GetCell(cellIndex[icell], absId, amp, time,mclabel,efrac);
266
267 Double_t eta, phi;
268 fGeom->EtaPhiFromIndex(cellIndex[icell], eta, phi);
269
270 if(phi<minPhiPatch) minPhiPatch = phi;
271 if(phi>maxPhiPatch) maxPhiPatch = phi;
272 if(eta<minEtaPatch) minEtaPatch = eta;
273 if(eta>maxEtaPatch) maxEtaPatch = eta;
274
275 sumAmp+=fCaloCells->GetAmplitude(cellIndex[icell]);
276 patchEnergy+=amp;
277
278 }//cells in fastor loop
279 }//fastor col loop
280 }//fastor row loop
281
282 Double_t etaCent = minEtaPatch + 0.5*(maxEtaPatch-minEtaPatch);
283 Double_t phiCent = minPhiPatch + 0.5*(maxPhiPatch-minPhiPatch);
284 fh3PatchEnergyEtaPhiCenter->Fill(patchEnergy,etaCent,phiCent);
285
286 if(patchEnergy>0)
287 if(fDebug>2) AliInfo(Form("%s: patch edges eta: %f - %f phi: %f - %f ampTrg: %f patchEnergy: %f sumAmp: %f",GetName(),minEtaPatch,maxEtaPatch,minPhiPatch,maxPhiPatch,ampTrg,patchEnergy,sumAmp));
288
289 if(patchEnergy>0) nPatchNotEmpty++;
290
291 if(patchEnergy>fMaxPatchEnergy) fMaxPatchEnergy=patchEnergy;
292
293 }
294 }
295
296}
297
298//________________________________________________________________________
299void AliAnalysisTaskEmcalJetTriggerQA::LoadExtraBranches() {
300 //
301 // get charged jets
302 //
303
609c7a0d 304 if (!fJetsName2.IsNull() && !fJets2) {
305 fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName2.Data()));
306 if (!fJets2) {
307 AliError(Form("%s: Could not retrieve charged jets %s!", GetName(), fJetsName2.Data()));
0531c107 308 return;
309 }
609c7a0d 310 else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
311 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetsName2.Data()));
312 fJets2 = 0;
0531c107 313 return;
314 }
315 }
316
317 //Get Charged Rho
318 if (!fRhoChName.IsNull() && !fRhoCh) {
319 fRhoCh = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoChName.Data()));
320 if (!fRhoCh) {
321 AliError(Form("%s: Could not retrieve charged rho %s!", GetName(), fRhoChName.Data()));
322 return;
323 }
324 else
325 fRhoChVal = 0.;
326 }
327 if (fRhoChName.IsNull() ) {
328 if(fDebug>10) AliInfo(Form("%s: fRhoChName empty %s",GetName(),fRhoChName.Data()));
329 fRhoChVal = 0.;
330 }
331 else if(!fRhoChName.IsNull() && fRhoCh) //do not use rho if rhoType==0
332 fRhoChVal = fRhoCh->GetVal();
333
334}
335
336//________________________________________________________________________
609c7a0d 337Bool_t AliAnalysisTaskEmcalJetTriggerQA::AcceptJet2(const AliEmcalJet *jet) const {
0531c107 338
609c7a0d 339 // Accept jet in 2nd branch
0531c107 340
341 Bool_t accept = kFALSE;
342
343 if(jet->Pt()<0.15) //reject ghost jets
344 return accept;
345
346 //do area cut
347 if(jet->Area()<fJetAreaCut)
348 return accept;
349
609c7a0d 350 //Check if jet is within fiducial acceptance
0531c107 351 Double_t eta = jet->Eta();
352 Double_t phi = jet->Phi();
609c7a0d 353
354 if(eta<fEtaMinJet2 || eta>fEtaMaxJet2)
0531c107 355 return accept;
356
609c7a0d 357 if(phi<fPhiMinJet2 || phi>fPhiMaxJet2)
0531c107 358 return accept;
359
609c7a0d 360 if (jet->MaxTrackPt() > fMaxTrackPtJet2)
0531c107 361 return kFALSE;
362
363 return kTRUE;
364
365}
366
367//________________________________________________________________________
368void AliAnalysisTaskEmcalJetTriggerQA::UserCreateOutputObjects()
369{
370 // Create user output.
371
372 InitOnce();
373
374 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
375
376 Bool_t oldStatus = TH1::AddDirectoryStatus();
377 TH1::AddDirectory(kFALSE);
378
379 fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
380 fOutput->Add(fhNEvents);
381
382 Int_t nBinsPt = 120;
383 Double_t minPt = -20.;
384 Double_t maxPt = 100.;
385 Int_t nBinsEta = 40;
386 Double_t minEta = -1.;
387 Double_t maxEta = 1.;
388 Int_t nBinsPhi = 18*6;
389 Double_t minPhi = 0.;
390 Double_t maxPhi = TMath::TwoPi();
391 Int_t nBinsArea = 100;
392 Double_t minArea = 0.;
393 Double_t maxArea = 1.;
394
395 fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
396 fOutput->Add(fh3PtEtaPhiJetFull);
397
398 fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
399 fOutput->Add(fh3PtEtaPhiJetCharged);
400
401 fh2NJetsPtFull = new TH2F("fh2NJetsPtFull","fh2NJetsPtFull;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,nBinsPt,minPt,maxPt);
402 fOutput->Add(fh2NJetsPtFull);
403
404 fh2NJetsPtCharged = new TH2F("fh2NJetsPtCharged","fh2NJetsPtCharged;N_{jets};#it{p}_{T}^{jet}",20,-0.5,19.5,nBinsPt,minPt,maxPt);
405 fOutput->Add(fh2NJetsPtCharged);
406
407 fh3PtEtaAreaJetFull = new TH3F("fh3PtEtaAreaJetFull","fh3PtEtaAreaJetFull;#it{p}_{T}^{jet};#eta;A",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsArea,minArea,maxArea);
408 fOutput->Add(fh3PtEtaAreaJetFull);
409
410 fh3PtEtaAreaJetCharged = new TH3F("fh3PtEtaAreaJetCharged","fh3PtEtaAreaJetCharged;#it{p}_{T}^{jet};#eta;A",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsArea,minArea,maxArea);
411 fOutput->Add(fh3PtEtaAreaJetCharged);
412
413 Int_t nBinsConst =100;
414 Double_t minConst = 0.;
415 Double_t maxConst = 100.;
416
417 Int_t nBinsMeanPt = 200;
418 Double_t minMeanPt = 0.;
419 Double_t maxMeanPt = 20.;
420
421 Int_t nBinsNEF = 101;
422 Double_t minNEF = 0.;
423 Double_t maxNEF = 1.01;
424
425 Int_t nBinsz = 101;
426 Double_t minz = 0.;
427 Double_t maxz = 1.01;
428
429 Int_t nBinsECluster =100;
430 Double_t minECluster = 0.;
431 Double_t maxECluster = 100.;
432
433
434 fh2PtNConstituentsCharged = new TH2F("fh2PtNConstituentsCharged","fh2PtNConstituentsCharged;#it{p}_{T}^{jet};N_{charged constituents}",nBinsPt,minPt,maxPt,nBinsConst,minConst,maxConst);
435 fOutput->Add(fh2PtNConstituentsCharged);
436
437 fh2PtNConstituents = new TH2F("fh2PtNConstituents","fh2PtNConstituents;#it{p}_{T}^{jet};N_{constituents}",nBinsPt,minPt,maxPt,nBinsConst,minConst,maxConst);
438 fOutput->Add(fh2PtNConstituents);
439
440 fh2PtMeanPtConstituentsCharged = new TH2F("fh2PtMeanPtConstituentsCharged","fh2PtMeanPtConstituentsCharged;#it{p}_{T}^{jet};charged #langle #it{p}_{T} #rangle",nBinsPt,minPt,maxPt,nBinsMeanPt,minMeanPt,maxMeanPt);
441 fOutput->Add(fh2PtMeanPtConstituentsCharged);
442
443 fh2PtMeanPtConstituentsNeutral = new TH2F("fh2PtMeanPtConstituentsNeutral","fh2PtMeanPtConstituentsNeutral;#it{p}_{T}^{jet};neutral langle #it{p}_{T} #rangle",nBinsPt,minPt,maxPt,nBinsMeanPt,minMeanPt,maxMeanPt);
444 fOutput->Add(fh2PtMeanPtConstituentsNeutral);
445
446 fh2PtNEF = new TH2F("fh2PtNEF","fh2PtNEF;#it{p}_{T}^{jet};NEF",nBinsPt,minPt,maxPt,nBinsNEF,minNEF,maxNEF);
447 fOutput->Add(fh2PtNEF);
448
449 fh2Ptz = new TH2F("fh2Ptz","fh2Ptz;#it{p}_{T}^{jet};z=p_{t,trk}^{proj}/p_{jet}",nBinsPt,minPt,maxPt,nBinsz,minz,maxz);
450 fOutput->Add(fh2Ptz);
451
452 fh2PtLeadJet1VsLeadJet2 = new TH2F("fh2PtLeadJet1VsLeadJet2","fh2PtLeadJet1VsLeadJet2;#it{p}_{T}^{jet 1};#it{p}_{T}^{jet 2}",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
453 fOutput->Add(fh2PtLeadJet1VsLeadJet2);
454
455 fh3EEtaPhiCluster = new TH3F("fh3EEtaPhiCluster","fh3EEtaPhiCluster;E_{clus};#eta,#phi",nBinsECluster,minECluster,maxECluster,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
456 fOutput->Add(fh3EEtaPhiCluster);
457
458 fh3PtLeadJet1VsPatchEnergy = new TH3F("fh3PtLeadJet1VsPatchEnergy","fh3PtLeadJet1VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,2,-0.5,1.5);
459 fOutput->Add(fh3PtLeadJet1VsPatchEnergy);
460 fh3PtLeadJet2VsPatchEnergy = new TH3F("fh3PtLeadJet2VsPatchEnergy","fh3PtLeadJet2VsPatchEnergy;#it{p}_{T}^{jet 1};Amplitude_{patch};trig type",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,2,-0.5,1.5);
461 fOutput->Add(fh3PtLeadJet2VsPatchEnergy);
462
463 fh3PatchEnergyEtaPhiCenter = new TH3F("fh3PatchEnergyEtaPhiCenter","fh3PatchEnergyEtaPhiCenter;E_{patch};#eta;#phi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
464 fOutput->Add(fh3PatchEnergyEtaPhiCenter);
465
466
467 // =========== Switch on Sumw2 for all histos ===========
468 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
469 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
470 if (h1){
471 h1->Sumw2();
472 continue;
473 }
474 TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
475 if (h2){
476 h2->Sumw2();
477 continue;
478 }
479 TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
480 if (h3){
481 h3->Sumw2();
482 continue;
483 }
484 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
485 if(hn)hn->Sumw2();
486 }
487
488 TH1::AddDirectory(oldStatus);
489
490 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
491}
492
493//________________________________________________________________________
494Bool_t AliAnalysisTaskEmcalJetTriggerQA::FillHistograms()
495{
496 // Fill histograms.
497
498
499 if (fCaloClusters) {
500 const Int_t nclusters = fCaloClusters->GetEntriesFast();
501
502 for (Int_t ic = 0; ic < nclusters; ic++) {
503 AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
504 if (!cluster) {
505 AliError(Form("Could not receive cluster %d", ic));
506 continue;
507 }
508 if (!cluster->IsEMCAL())
509 continue;
510
511 TLorentzVector lp;
512 cluster->GetMomentum(lp, const_cast<Double_t*>(fVertex));
513
514 //Fill eta,phi,E of clusters here
515 fh3EEtaPhiCluster->Fill(lp.E(),lp.Eta(),lp.Phi());
516 }
517 }
518
519 Double_t ptLeadJet1 = 0.;
520 Double_t ptLeadJet2 = 0.;
521
522 TArrayI *nJetsArr = new TArrayI(fh2NJetsPtFull->GetNbinsY()+1);
523 nJetsArr->Reset(0);
524 nJetsArr->Set(fh2NJetsPtFull->GetNbinsY()+1);
525
526 if (fJets) {
527 const Int_t njets = fJets->GetEntriesFast();
528 for (Int_t ij = 0; ij < njets; ij++) {
529
530 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
531 if (!jet) {
532 AliError(Form("Could not receive jet %d", ij));
533 continue;
534 }
535
536 if (!AcceptJet(jet))
537 continue;
538
539 Double_t jetPt = jet->Pt();
540 if(jetPt>ptLeadJet1) ptLeadJet1=jetPt;
541 fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
542 fh3PtEtaAreaJetFull->Fill(jetPt,jet->Eta(),jet->Area());
543
544 //count jets above certain pT threshold
545 Int_t ptbin = fh2NJetsPtFull->GetYaxis()->FindBin(jetPt);
546 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtFull->GetNbinsY(); iptbin++)
547 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
548
549 fh2PtNConstituentsCharged->Fill(jetPt,jet->GetNumberOfTracks());
550 fh2PtNConstituents->Fill(jetPt,jet->GetNumberOfConstituents());
551
552 fh2PtNEF->Fill(jetPt,jet->NEF());
553
554 AliVParticle *vp;
555 Double_t sumPtCh = 0.;
556 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
557 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
558 if(!vp) continue;
559 fh2Ptz->Fill(jetPt,GetZ(vp,jet));
560
561 sumPtCh+=vp->Pt();
562
563 }
564
565 if(jet->GetNumberOfTracks()>0)
566 fh2PtMeanPtConstituentsCharged->Fill(jetPt,sumPtCh/(double)(jet->GetNumberOfTracks()) );
567
568
569 AliVCluster *vc = 0x0;
570 Double_t sumPtNe = 0.;
571 for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
572 vc = static_cast<AliVCluster*>(jet->ClusterAt(icc, fCaloClusters));
573 if(!vc) continue;
574
575 TLorentzVector lp;
576 vc->GetMomentum(lp, const_cast<Double_t*>(fVertex));
577
578 sumPtNe+=lp.Pt();
579
580 }
581
582 if(jet->GetNumberOfClusters()>0)
583 fh2PtMeanPtConstituentsNeutral->Fill(jetPt,sumPtNe/(double)(jet->GetNumberOfClusters()) );
584
585 } //full jet loop
586
587 for(Int_t i=1; i<=fh2NJetsPtFull->GetNbinsY(); i++) {
588 Int_t nJetsInEvent = nJetsArr->At(i);
589 fh2NJetsPtFull->Fill(nJetsInEvent,fh2NJetsPtFull->GetYaxis()->GetBinCenter(i));
590 }
591
592 }
593
594 //Reset array to zero to also count charged jets
595 nJetsArr->Reset(0);
596
597 //Loop over charged jets
609c7a0d 598 if (fJets2) {
599 const Int_t njetsCh = fJets2->GetEntriesFast();
0531c107 600 for (Int_t ij = 0; ij < njetsCh; ij++) {
601
609c7a0d 602 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets2->At(ij));
0531c107 603 if (!jet) {
604 AliError(Form("Could not receive charged jet %d", ij));
605 continue;
606 }
607
609c7a0d 608 if(!AcceptJet2(jet))
0531c107 609 continue;
610
611 Double_t jetPt = jet->Pt();
612 if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
613 fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
614 fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
615
616 //count jets above certain pT threshold
617 Int_t ptbin = fh2NJetsPtCharged->GetYaxis()->FindBin(jetPt);
618 for(Int_t iptbin = ptbin; iptbin<=fh2NJetsPtCharged->GetNbinsY(); iptbin++)
619 nJetsArr->AddAt(nJetsArr->At(iptbin)+1,iptbin);
620
621 }//ch jet loop
622 for(Int_t i=1; i<=fh2NJetsPtCharged->GetNbinsY(); i++) {
623 Int_t nJetsInEvent = nJetsArr->At(i);
624 fh2NJetsPtCharged->Fill(nJetsInEvent,fh2NJetsPtCharged->GetYaxis()->GetBinCenter(i));
625 }
626 }
627
609c7a0d 628 if(fJets && fJets2) {
0531c107 629 fh2PtLeadJet1VsLeadJet2->Fill(ptLeadJet1,ptLeadJet2);
630 }
631
632 fh3PtLeadJet1VsPatchEnergy->Fill(ptLeadJet1,fMaxPatchEnergy,fTriggerType);
633 fh3PtLeadJet2VsPatchEnergy->Fill(ptLeadJet2,fMaxPatchEnergy,fTriggerType);
634
635 if(nJetsArr) delete nJetsArr;
636
637 return kTRUE;
638}
639
640//________________________________________________________________________
641Bool_t AliAnalysisTaskEmcalJetTriggerQA::Run()
642{
643 // Run analysis code here, if needed. It will be executed before FillHistograms().
644
645 //Check if event is selected (vertex & pile-up)
646 if(!SelectEvent())
647 return kFALSE;
648
649 LoadExtraBranches();
650
651 if(fRhoChName.IsNull())
652 fRhoChVal = 0.;
653 else
654 fRhoChVal = fRhoCh->GetVal();
655
609c7a0d 656 if(!fTriggerClass.IsNull())
657 FindTriggerPatch();
0531c107 658
659 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
660}
661
662//_______________________________________________________________________
663void AliAnalysisTaskEmcalJetTriggerQA::Terminate(Option_t *)
664{
665 // Called once at the end of the analysis.
666}
667//________________________________________________________________________
668Double_t AliAnalysisTaskEmcalJetTriggerQA::GetZ(const AliVParticle *trk, const AliEmcalJet *jet) const
669{
670 // Get Z of constituent trk
671
672 return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
673}
674
675//________________________________________________________________________
676Double_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
677{
678 //
679 // Get the z of a constituent inside of a jet
680 //
681
682 Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
683
684 if(pJetSq>0.)
685 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
686 else {
687 AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
688 return 0;
689 }
690}