]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/UserTasks/AliAnalysisTaskJetsTriggerTRD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetsTriggerTRD.cxx
CommitLineData
7032c043 1// ROOT
2#include "TFile.h"
3#include "TList.h"
4#include "TH1.h"
5#include "TH2.h"
6#include "TH3.h"
ebdf1dbb 7#include "TRandom.h"
7032c043 8
9// analysis framework
10#include "AliAnalysisManager.h"
11#include "AliInputEventHandler.h"
12#include "AliVEvent.h"
e7808ddf 13#include "AliVTrdTrack.h"
7032c043 14
15// MC stuff
16#include "AliMCEvent.h"
17#include "AliGenPythiaEventHeader.h"
18
19// ESD stuff
20#include "AliESDEvent.h"
21#include "AliESDInputHandler.h"
22#include "AliESDtrack.h"
23#include "AliESDTrdTrack.h"
24#include "AliESDTrdTracklet.h"
25#include "AliESDTrdTrigger.h"
26
27// AOD stuff
28#include "AliAODEvent.h"
29#include "AliAODJet.h"
30#include "AliAODTrack.h"
31
32// jet tasks
33#include "AliAnalysisTaskJetServices.h"
34#include "AliAnalysisHelperJetTasks.h"
35
36#include "AliAnalysisTaskJetsTriggerTRD.h"
37
7032c043 38AliAnalysisTaskJetsTriggerTRD::AliAnalysisTaskJetsTriggerTRD(const char *name) :
39 AliAnalysisTaskSE(name),
e7808ddf 40 fTriggerMask(0),
ebdf1dbb 41 fTrdTrg(),
7032c043 42 fOutputList(),
43 fHist(),
44 fShortTaskId("jets_trg_trd"),
ebdf1dbb 45 fNoTriggers(kTrgLast),
e7808ddf 46 fNoJetPtBins(80),
a26539a7 47 fJetPtBinMax(400),
ebdf1dbb 48 fAvgXsection(0.),
a26539a7 49 fAvgTrials(0.),
ebdf1dbb 50 fPtHard(0.),
51 fNTrials(0),
52 fGlobalEfficiencyGTU(.8),
53 fGtuLabel(-1) // -3 for hw, -1821 for re-simulation
7032c043 54{
55 // default ctor
56
57 DefineOutput(1, TList::Class());
58}
59
60AliAnalysisTaskJetsTriggerTRD::~AliAnalysisTaskJetsTriggerTRD()
61{
62 // dtor
63
64}
65
66void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects()
67{
68 // create user output objects
69
70 // setup list
71 OpenFile(1);
72 fOutputList = new TList();
73 fOutputList->SetOwner();
74
75 // setup histograms
e7808ddf 76 TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
77 kStatLast-1, .5, kStatLast-.5);
78 histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
a26539a7 79 histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
e7808ddf 80 histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
a26539a7 81 histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
82
ebdf1dbb 83 if (HasMC()) {
84 fNoTriggers = kTrgMCLast;
85 AddHistogram(ID(kHistXsection), "xsection stats",
86 1, 0., 1.);
87 AddHistogram(ID(kHistPtHard), "pt hard;#hat{p}_{T} (GeV/c);counts",
88 fNoJetPtBins, 0., fJetPtBinMax);
89 AddHistogram(ID(kHistJetPtMC), "leading jet spectrum (MC, |#eta| < 0.5);p_{T}^{jet} (GeV/c);counts",
90 fNoJetPtBins, 0., fJetPtBinMax);
91 }
92
93 AddHistogram(ID(kHistJetEtaAvg), "average eta",
94 100, -1., 1.);
e7808ddf 95
96 AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts",
a26539a7 97 400, -.5, 399.5);
e7808ddf 98
99 AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts",
100 100, 0., 25.);
101
ebdf1dbb 102 AddHistogram(ID(kHistTrackEffGTU), "found GTU tracks;p_{T} (GeV/c);#eta;#varphi",
103 100, 0., 25.,
104 100, -1., 1.,
105 100, 0., 2.*TMath::Pi());
106 AddHistogram(ID(kHistTrackEffMC), "wanted GTU tracks;p_{T} (GeV/c);#eta;#varphi",
107 100, 0., 25.,
108 100, -1., 1.,
109 100, 0., 2.*TMath::Pi());
110
a26539a7 111 AddHistogram(ID(kHistNPtMin),
112 "rejection;p_{T}^{min};N_{trk};trigger",
113 100, 0., 10.,
114 20, 0., 20.,
ebdf1dbb 115 fNoTriggers - 1, .5, fNoTriggers - .5);
a26539a7 116
ebdf1dbb 117 // leading jet
a26539a7 118 AddHistogram(ID(kHistLeadJetPt),
ebdf1dbb 119 "leading jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);counts",
120 fNoJetPtBins, 0., fJetPtBinMax,
121 fNoTriggers - 1, .5, fNoTriggers - .5);
122 AddHistogram(ID(kHistLeadJetPtEta),
123 "leading jet #eta (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#eta",
a26539a7 124 fNoJetPtBins, 0., fJetPtBinMax,
ebdf1dbb 125 fNoTriggers - 1, .5, fNoTriggers - .5,
126 40, -.5, .5);
127 AddHistogram(ID(kHistLeadJetPtPhi),
128 "leading jet #phi (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#varphi",
129 fNoJetPtBins, 0., fJetPtBinMax,
130 fNoTriggers - 1, .5, fNoTriggers - .5,
131 40, 0., 2. * TMath::Pi());
132 AddHistogram(ID(kHistLeadJetEtaPhi),
133 "leading jet #eta - #varphi (|#eta| < 0.5);#eta;trigger;#varphi",
134 40, -.5, .5,
135 fNoTriggers - 1, .5, fNoTriggers - .5,
136 40, 0., 2. * TMath::Pi());
137
138 AddHistogram(ID(kHistLeadJetPtTrackPt),
139 "leading jet pt vs track pt;p_{T} (GeV/c);trigger;p_{T}^{jet,ch} (GeV/c)",
140 100., 0., 20.,
141 fNoTriggers - 1, .5, fNoTriggers - .5,
142 fNoJetPtBins, 0., fJetPtBinMax);
143 AddHistogram(ID(kHistLeadJetPtZ),
144 "leading jet pt vs z;z;trigger;p_{T}^{jet,ch} (GeV/c)",
145 100., 0., 1.,
146 fNoTriggers - 1, .5, fNoTriggers - .5,
147 fNoJetPtBins, 0., fJetPtBinMax);
148 AddHistogram(ID(kHistLeadJetPtXi),
149 "leading jet pt vs #xi;#xi;trigger;p_{T}^{jet,ch} (GeV/c)",
150 100., 0., 10.,
151 fNoTriggers - 1, .5, fNoTriggers - .5,
152 fNoJetPtBins, 0., fJetPtBinMax);
153
154 // inclusive jets
a26539a7 155 AddHistogram(ID(kHistJetPt),
ebdf1dbb 156 "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger",
157 fNoJetPtBins, 0., fJetPtBinMax,
158 fNoTriggers - 1, .5, fNoTriggers - .5);
159 AddHistogram(ID(kHistJetPtEta),
160 "jet #eta (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#eta",
a26539a7 161 fNoJetPtBins, 0., fJetPtBinMax,
ebdf1dbb 162 fNoTriggers - 1, .5, fNoTriggers - .5,
163 40, -.5, .5);
164 AddHistogram(ID(kHistJetPtPhi),
165 "jet #phi (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#varphi",
166 fNoJetPtBins, 0., fJetPtBinMax,
167 fNoTriggers - 1, .5, fNoTriggers - .5,
168 40, 0., 2. * TMath::Pi());
169 AddHistogram(ID(kHistJetEtaPhi),
170 "jet #eta - #varphi (|#eta| < 0.5);#eta;trigger;#varphi",
171 40, -.5, .5,
172 fNoTriggers - 1, .5, fNoTriggers - .5,
173 40, 0., 2. * TMath::Pi());
174
a26539a7 175 AddHistogram(ID(kHistJetPtITS),
ebdf1dbb 176 "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger",
a26539a7 177 fNoJetPtBins, 0., fJetPtBinMax,
ebdf1dbb 178 fNoTriggers - 1, .5, fNoTriggers - .5);
a26539a7 179 AddHistogram(ID(kHistJetPt3x3),
ebdf1dbb 180 "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger",
a26539a7 181 fNoJetPtBins, 0., fJetPtBinMax,
ebdf1dbb 182 fNoTriggers - 1, .5, fNoTriggers - .5);
183
184 AddHistogram(ID(kHistJetPtTrackPt),
185 "jet pt vs track pt;p_{T} (GeV/c);trigger;p_{T}^{jet,ch} (GeV/c)",
186 100., 0., 20.,
187 fNoTriggers - 1, .5, fNoTriggers - .5,
188 fNoJetPtBins, 0., fJetPtBinMax);
189 AddHistogram(ID(kHistJetPtZ),
190 "jet pt vs z;z;trigger;p_{T}^{jet,ch} (GeV/c)",
191 100., 0., 1.,
192 fNoTriggers - 1, .5, fNoTriggers - .5,
193 fNoJetPtBins, 0., fJetPtBinMax);
194 AddHistogram(ID(kHistJetPtXi),
195 "jet pt vs #xi;#xi;trigger;p_{T}^{jet,ch} (GeV/c)",
196 100., 0., 10.,
197 fNoTriggers - 1, .5, fNoTriggers - .5,
198 fNoJetPtBins, 0., fJetPtBinMax);
199
200 for (Int_t iHist = kHistLeadJetPt; iHist <= kHistJetPtXi; ++iHist) {
a26539a7 201 TH1 *h = GetHistogram(Hist_t (iHist));
202 h->GetYaxis()->SetBinLabel(ID(kTrgMinBias));
203 h->GetYaxis()->SetBinLabel(ID(kTrgInt7));
204 h->GetYaxis()->SetBinLabel(ID(kTrgInt8));
205 h->GetYaxis()->SetBinLabel(ID(kTrgEMC7));
206 h->GetYaxis()->SetBinLabel(ID(kTrgEMC8));
ebdf1dbb 207 h->GetYaxis()->SetBinLabel(ID(kTrgPbPb));
208 h->GetYaxis()->SetBinLabel(ID(kTrgCentral));
209 h->GetYaxis()->SetBinLabel(ID(kTrgSemiCentral));
a26539a7 210 h->GetYaxis()->SetBinLabel(ID(kTrgInt7WUHJT));
211 h->GetYaxis()->SetBinLabel(ID(kTrgInt8WUHJT));
212 h->GetYaxis()->SetBinLabel(ID(kTrgEMC7WUHJT));
213 h->GetYaxis()->SetBinLabel(ID(kTrgEMC8WUHJT));
6883db3e 214 h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE));
215 h->GetYaxis()->SetBinLabel(ID(kTrgEMCEGA));
ebdf1dbb 216 h->GetYaxis()->SetBinLabel(ID(kTrgInt7_WU));
217 h->GetYaxis()->SetBinLabel(ID(kTrgInt7_WUHJT));
218 h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE_WU));
219 h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE_WUHJT));
220 if (HasMC()) {
221 h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3Vtx));
222 h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRD));
223 h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRDeff));
224 h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRDeffmap));
225 }
a26539a7 226 }
7032c043 227
a26539a7 228 AddHistogram(ID(kHistJetPtNoTracks3),
ebdf1dbb 229 "number of tracks above 3 GeV;p_{T}^{jet,ch};no. of tracks",
e7808ddf 230 fNoJetPtBins, 0., fJetPtBinMax,
a26539a7 231 40, -.5, 39.5);
e7808ddf 232
7032c043 233 PostData(1, fOutputList);
234}
235
236Bool_t AliAnalysisTaskJetsTriggerTRD::Notify()
237{
238 // actions to be taken upon notification about input file change
239
a26539a7 240 // ??? check ???
ebdf1dbb 241 // we should only count the cross section we see in the analysis,
242 // i.e. fXsection / nTrials
243
244 if (HasMC()) {
245 fAvgXsection = 0.;
246 Float_t xSection = 0.;
247 Float_t nTrials = 0.;
248 Float_t nEntries = 0.;
249
250 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
251
252 if (tree) {
253 TFile *curfile = tree->GetCurrentFile();
254 if (!curfile) {
255 AliError("No current file");
256 return kFALSE;
257 }
a26539a7 258
ebdf1dbb 259 if (!AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(), xSection, nTrials)) {
260 AliError("retrieval of cross section failed");
261 }
a26539a7 262
ebdf1dbb 263 nEntries = (Float_t) tree->GetTree()->GetEntries();
264 if (nEntries > 0.) {
265 fAvgTrials = nTrials / nEntries;
266 if (nTrials > 0.)
267 fAvgXsection = xSection / (nTrials / nEntries);
268 }
a26539a7 269 }
ebdf1dbb 270 printf("n_trials = %f\n", nTrials);
a26539a7 271 }
272
7032c043 273 return AliAnalysisTaskSE::Notify();
274}
275
276void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
277{
278 // actual work
279
280 // setup pointers to input data (null if unavailable)
281 // mcEvent: MC input
282 // esdEvent: ESD input
283 // outEvent: AOD output
284 // aodEvent: AOD input if available, otherwise AOD output
a26539a7 285 AliMCEvent *mcEvent = this->MCEvent();
7032c043 286 AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
287 AliAODEvent* outEvent = this->AODEvent();
288 AliAODEvent *aodEvent = outEvent;
289 if (dynamic_cast<AliAODEvent*>(this->InputEvent()))
290 aodEvent = (AliAODEvent*) (this->InputEvent());
291
292 if ((fDebug > 0) && esdEvent)
293 printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile());
294
ebdf1dbb 295 Int_t nTracksMC = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks
a26539a7 296 // Int_t nTracks = InputEvent()->GetNumberOfTracks(); // no. of global tracks
297 Int_t nTracksTRD = InputEvent()->GetNumberOfTrdTracks(); // no. of GTU tracks
7032c043 298
299 TList partMC;
300 TList partEsd;
301 TList partGtu;
302
a26539a7 303 Float_t leadingJetPtMC = 0.; // leading jet energy from MC information
7032c043 304 Float_t leadingJetPtRec = 0.; // leading jet energy from AOD information
305
e7808ddf 306 // record number of sampled events and detect trigger contributions
307 FillH1(kHistStat, kStatSeen);
308 if (!DetectTriggers()) {
309 AliError("Failed to detect the triggers");
310 return;
311 }
7032c043 312
a26539a7 313 // only continue for events from interesting triggers
314 if (fTriggerMask == 0)
7032c043 315 return;
a26539a7 316 FillH1(kHistStat, kStatTrg);
317
318 // no further technical requirements for the event at the moment
e7808ddf 319 FillH1(kHistStat, kStatUsed);
ebdf1dbb 320 if (HasMC()) {
321 if (fAvgXsection == 0.) {
322 AliError("zero cross section");
323 return;
324 }
325 FillH1(kHistXsection, .5, fAvgXsection);
326 }
7032c043 327
a26539a7 328 // apply event cuts
329 const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
330 if (!vtx ||
331 (vtx->GetNContributors() < 3.) ||
332 (vtx->GetZ() > 10.))
333 return;
334 FillH1(kHistStat, kStatEvCuts);
335
336 // extract MC information
337 if (mcEvent) {
338 // check for PYTHIA event header
339 AliGenPythiaEventHeader *pythiaHeader =
340 dynamic_cast<AliGenPythiaEventHeader*> (mcEvent->GenEventHeader());
341 if (!pythiaHeader) {
342 AliWarning("MC event without PYTHIA event header!\n");
343 }
344 else {
345 fPtHard = pythiaHeader->GetPtHard();
ebdf1dbb 346 Int_t nTrials = pythiaHeader->Trials();
347 fNTrials += nTrials;
348
349 FillH1(kHistPtHard, fPtHard);
a26539a7 350
351 // loop over jets from PYTHIA
ebdf1dbb 352 Float_t eta1 = -10.;
353 Float_t eta2 = -10.;
a26539a7 354 for (Int_t iJet = 0; iJet < pythiaHeader->NTriggerJets(); iJet++) {
355 Float_t p[4];
356 pythiaHeader->TriggerJet(iJet, p);
357 TLorentzVector pJet(p);
358 Float_t pt = pJet.Pt();
359 // only consider jets with |eta| < 0.5
360 Float_t eta = pJet.Eta();
361 if (TMath::Abs(eta) > 0.5)
362 continue;
363 if (pt > leadingJetPtMC)
364 leadingJetPtMC = pt;
ebdf1dbb 365
366 // for eta averge determination consider only jets above 60 GeV
367 if (pt < 60.)
368 continue;
369 if (eta1 < -1.)
370 eta1 = eta;
371 else if (eta2 < -1.)
372 eta2 = eta;
a26539a7 373 }
374 // fill histogram for leading jet pt spectrum
ebdf1dbb 375 FillH1(kHistJetPtMC, leadingJetPtMC);
376 // fill histogram for eta average
377 if ((eta1 > -1.) && (eta2 > -1.))
378 FillH1(kHistJetEtaAvg, (eta1 + eta2)/2.);
379 }
380
381 for (Int_t iTrack = 0; iTrack < nTracksMC; ++iTrack) {
382 AliVParticle *part = mcEvent->GetTrack(iTrack);
383 if (AcceptTrackMC(iTrack)) {
384 FillH3(kHistTrackEffMC, part->Pt(), part->Eta(), part->Phi());
385 }
a26539a7 386 }
e7808ddf 387 }
7032c043 388
a26539a7 389 // loop over GTU tracks
390 for (Int_t iTrack = 0; iTrack < nTracksTRD; ++iTrack) {
391 AliVTrdTrack *trk = InputEvent()->GetTrdTrack(iTrack);
ebdf1dbb 392 // printf("trk %p has pt %5.2f and label %i\n",
393 // trk, trk->Pt(), trk->GetLabel());
394 if ((fGtuLabel != -1) && (trk->GetLabel() != fGtuLabel))
395 continue;
a26539a7 396 FillH1(kHistTrackGTU, TMath::Abs(trk->Pt()));
397 partGtu.Add(trk);
398 }
399 partGtu.Sort(kSortAscending);
400
401 Int_t nTracksPerStack[90] = { 0 };
402 Int_t nTracksPerStackMax = 0;
403
404 TIter nextPartGtu(&partGtu);
405 while (AliVTrdTrack *trdTrack = (AliVTrdTrack*) nextPartGtu()) {
406 // count no. of tracks in stack,
407 // check whether this number was reached before,
408 // if not store pt^min(n),
409 // i.e. pt of current track because of sorting
410
411 Int_t sec = trdTrack->GetSector();
412 Int_t stack = trdTrack->GetStack();
413
414 if ((sec > -1) && (sec < 18) &&
415 (stack > -1) && (stack < 5)) {
416 ++nTracksPerStack[5*sec + stack];
417 if (nTracksPerStack[5*sec + stack] > nTracksPerStackMax) {
418 ++nTracksPerStackMax;
419
ebdf1dbb 420 for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
a26539a7 421 if (IsTrigger(Trigger_t (iTrigger)))
422 FillH3(kHistNPtMin,
423 TMath::Abs(trdTrack->Pt()), nTracksPerStackMax, iTrigger);
424 }
ebdf1dbb 425 if (HasMC()) {
426 Int_t label = trdTrack->GetLabel();
427 if (label > -1) {
428 if (AcceptTrackMC(label)) {
429 AliVParticle *part = MCEvent()->GetTrack(label);
430 FillH3(kHistTrackEffGTU, part->Pt(), part->Eta(), part->Phi());
431 }
432 }
433 else {
434 AliWarning(Form("GTU track at %p with no label", trdTrack));
435 const Int_t nLayers = 6;
436 for (Int_t iLayer = 0; iLayer < nLayers; ++iLayer) {
437 AliVTrdTracklet *trkl = trdTrack->GetTracklet(iLayer);
438 if (trkl)
439 AliWarning(Form("tracklet in layer %i has label %i\n",
440 iLayer, trkl->GetLabel()));
441 }
442 }
443 }
a26539a7 444 }
445 else {
446 AliError(Form("Invalid sector or stack: %i %i",
447 sec, stack));
448 }
449
450 }
7032c043 451
a26539a7 452 // loop over jets from AOD event
7032c043 453 if (aodEvent) {
a26539a7 454 TClonesArray *jetArray =
455 dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
7032c043 456 if (jetArray) {
457 Int_t nJets = jetArray->GetEntriesFast();
e7808ddf 458 FillH1(kHistNoJets, nJets);
7032c043 459
a26539a7 460 AliAODJet *leadJet = 0x0;
461 // AliAODJet *subleadJet = 0x0;
7032c043 462
a26539a7 463 for (Int_t iJet = 0; iJet < nJets; ++iJet) {
464 AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet];
7032c043 465 if (TMath::Abs(jet->Eta()) < 0.5) {
a26539a7 466 // check contributing tracks
467 Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
468 Int_t nJetTracks3 = 0;
469 AliAODTrack *leadingTrack = 0x0;
470 for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
471 AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
472
473 // count constituents above 3 GeV/c
474 if (track->Pt() > 3.)
475 ++nJetTracks3;
476
477 // find the leading track
478 if (!leadingTrack ||
479 (track->Pt() > leadingTrack->Pt()))
480 leadingTrack = track;
481 }
e7808ddf 482
483 // find leading jet
7032c043 484 if (TMath::Abs(jet->Pt()) > leadingJetPtRec) {
485 leadingJetPtRec = TMath::Abs(jet->Pt());
a26539a7 486 // subleadJet = leadJet;
7032c043 487 leadJet = jet;
488 }
489
ebdf1dbb 490 // jet pt spectrum and
491 // fragmentation function
492 for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
493 if (IsTrigger(Trigger_t (iTrigger))) {
494 Float_t jetPt = jet->Pt();
495
496 FillH1(kHistJetPt, jetPt, iTrigger);
497
498 FillH3(kHistJetPtEta, jet->Pt(), iTrigger, jet->Eta());
499 FillH3(kHistJetPtPhi, jet->Pt(), iTrigger, jet->Phi());
500 if (jet->Pt() > 50.)
501 FillH3(kHistJetEtaPhi, jet->Eta(), iTrigger, jet->Phi());
502
503 Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
504 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
505 AliVParticle *track = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(iTrack));
506 if (track) {
507 Float_t trackPt = track->Pt();
508 Float_t z = trackPt / jetPt;
509 Float_t xi = - TMath::Log(z);
510 FillH3(kHistJetPtTrackPt, trackPt, iTrigger, jet->Pt());
511 FillH3(kHistJetPtZ, z, iTrigger, jet->Pt());
512 FillH3(kHistJetPtXi, xi, iTrigger, jet->Pt());
513 }
514 }
515 }
a26539a7 516
517 // integrated over all triggers
518 FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
e7808ddf 519
a26539a7 520 // limit to jets with leading track having an ITS contribution
521 // with a hit in any SPD layer
e7808ddf 522 if (leadingTrack &&
523 (leadingTrack->GetFlags() & AliVTrack::kITSrefit) &&
524 (leadingTrack->HasPointOnITSLayer(0) || leadingTrack->HasPointOnITSLayer(1)))
ebdf1dbb 525 for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
a26539a7 526 if (IsTrigger(Trigger_t (iTrigger)))
527 FillH1(kHistJetPtITS, jet->Pt(), iTrigger);
e7808ddf 528
a26539a7 529 // limit to jets having 3 tracks above 3 GeV/c
e7808ddf 530 if (nJetTracks3 > 2)
ebdf1dbb 531 for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
a26539a7 532 if (IsTrigger(Trigger_t (iTrigger)))
533 FillH1(kHistJetPt3x3, jet->Pt(), iTrigger);
7032c043 534 }
535 }
ebdf1dbb 536
a26539a7 537 // fill leading jet information
ebdf1dbb 538 for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
539 if (IsTrigger(Trigger_t (iTrigger))) {
540 FillH2(kHistLeadJetPt, leadJet ? leadJet->Pt() : 0., iTrigger);
541 if (leadJet) {
542 Float_t leadJetPt = leadJet->Pt();
543
544 FillH3(kHistLeadJetPtEta, leadJet->Pt(), iTrigger, leadJet->Eta());
545 FillH3(kHistLeadJetPtPhi, leadJet->Pt(), iTrigger, leadJet->Phi());
546 if (leadJet->Pt() > 50.)
547 FillH3(kHistLeadJetEtaPhi, leadJet->Eta(), iTrigger, leadJet->Phi());
548
549 Int_t nTracks = leadJet->GetRefTracks()->GetEntriesFast();
550 for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
551 AliVParticle *track = dynamic_cast<AliVParticle*>(leadJet->GetRefTracks()->At(iTrack));
552 if (track) {
553 Float_t trackPt = track->Pt();
554 Float_t z = trackPt / leadJetPt;
555 Float_t xi = - TMath::Log(z);
556 FillH3(kHistLeadJetPtTrackPt, trackPt, iTrigger, leadJet->Pt());
557 FillH3(kHistLeadJetPtZ, z, iTrigger, leadJet->Pt());
558 FillH3(kHistLeadJetPtXi, xi, iTrigger, leadJet->Pt());
559 }
560 }
561 }
562 }
7032c043 563 }
564 else {
a26539a7 565 printf("no jet array found as branch %s\n", fJetBranchName);
566 aodEvent->Print();
7032c043 567 }
a26539a7 568 } else {
569 printf("no AOD event found\n");
7032c043 570 }
571
572 PostData(1, fOutputList);
573}
574
575void AliAnalysisTaskJetsTriggerTRD::Terminate(const Option_t * /* option */)
576{
577 // actions at task termination
578
ebdf1dbb 579 printf("total trials: %d\n", fNTrials);
7032c043 580}
581
e7808ddf 582Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
583{
584 fTriggerMask = 0;
585
a26539a7 586 AliInputEventHandler *inputHandler =
587 (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
588
589 AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) inputHandler->IsEventSelected();
e7808ddf 590 TString trgClasses = InputEvent()->GetFiredTriggerClasses();
591
ebdf1dbb 592 fTrdTrg.CalcTriggers(InputEvent());
593
594 UInt_t inputMaskL1 = 0;
595 if (AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*> (InputEvent()))
596 inputMaskL1 = esdEvent->GetHeader()->GetL1TriggerInputs();
597 else if (AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*> (InputEvent()))
598 inputMaskL1 = aodEvent->GetHeader()->GetL1TriggerInputs();
599
e7808ddf 600 if (fDebug > 2)
601 printf("trg: %8s %8s %8s %8s %8s %8s %8s (%s)\n",
602 (physSel & AliVEvent::kAnyINT) ? "kAnyINT" : " ",
603 (physSel & AliVEvent::kCINT5) ? "kCINT5" : " ",
604 (physSel & AliVEvent::kINT7) ? "kINT7" : " ",
605 (physSel & AliVEvent::kINT8) ? "kINT8" : " ",
606 (physSel & AliVEvent::kMB) ? "kMB" : " ",
607 (physSel & AliVEvent::kEMC7) ? "kEMC7" : " ",
608 (physSel & AliVEvent::kTRD) ? "kTRD" : " ",
609 trgClasses.Data()
610 );
611
ebdf1dbb 612 if (fDebug > 2)
613 printf("trg: %8s %8s %8s %8s (%s)\n",
614 fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHJT) ? "HJT" : " ",
615 fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHSE) ? "HSE" : " ",
616 fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHQU) ? "HQU" : " ",
617 fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHEE) ? "HEE" : " ",
618 trgClasses.Data()
619 );
620
e7808ddf 621 // physics selection
a26539a7 622 if ((physSel & (AliVEvent::kMB)))
623 MarkTrigger(kTrgMinBias);
624
ebdf1dbb 625 if ((physSel & (AliVEvent::kAnyINT | AliVEvent::kCentral | AliVEvent::kSemiCentral)))
626 MarkTrigger(kTrgPbPb);
627
628 if ((physSel & (AliVEvent::kINT7))) {
a26539a7 629 MarkTrigger(kTrgInt7);
ebdf1dbb 630 if (trgClasses.Contains("WU")) {
631 MarkTrigger(kTrgInt7_WU);
632 if (inputMaskL1 & (1 << 9))
633 MarkTrigger(kTrgInt7_WUHJT);
634 }
635 }
a26539a7 636
637 if ((physSel & (AliVEvent::kINT8)))
638 MarkTrigger(kTrgInt8);
639
6883db3e 640 if ((physSel & (AliVEvent::kEMC7)) &&
641 trgClasses.Contains("CEMC7"))
a26539a7 642 MarkTrigger(kTrgEMC7);
643
6883db3e 644 if ((physSel & (AliVEvent::kEMC8)) &&
645 trgClasses.Contains("CEMC8"))
a26539a7 646 MarkTrigger(kTrgEMC8);
647
ebdf1dbb 648 if ((physSel & (AliVEvent::kEMCEJE))) {
6883db3e 649 MarkTrigger(kTrgEMCEJE);
ebdf1dbb 650 if (trgClasses.Contains("WU")) {
651 MarkTrigger(kTrgEMCEJE_WU);
652 if (inputMaskL1 & (1 << 9))
653 MarkTrigger(kTrgEMCEJE_WUHJT);
654 }
655 }
6883db3e 656
657 if ((physSel & (AliVEvent::kEMCEGA)))
658 MarkTrigger(kTrgEMCEGA);
659
660 // for the TRD-triggered events we use the classes
a26539a7 661 if (trgClasses.Contains("CINT7WUHJT-"))
662 MarkTrigger(kTrgInt7WUHJT);
663
664 if (trgClasses.Contains("CINT8WUHJT-"))
665 MarkTrigger(kTrgInt8WUHJT);
666
667 if (trgClasses.Contains("CEMC7WUHJT-"))
668 MarkTrigger(kTrgEMC7WUHJT);
669
670 if (trgClasses.Contains("CEMC8WUHJT-"))
671 MarkTrigger(kTrgEMC8WUHJT);
e7808ddf 672
ebdf1dbb 673 if (HasMC())
674 DetectMCTriggers();
675
676 return kTRUE;
677}
678
679Bool_t AliAnalysisTaskJetsTriggerTRD::DetectMCTriggers()
680{
681 AliMCEvent *mcEvent = MCEvent();
682
683 Int_t nTracksMC = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks
684 TList partMC; // list of particles
685
686 // fill tracks passing cuts to list
687 for (Int_t iTrackMC = 0; iTrackMC < nTracksMC; ++iTrackMC) {
688 AliVParticle *part = mcEvent->GetTrack(iTrackMC);
689
690 // only consider primaries
691 if (!mcEvent->IsPhysicalPrimary(iTrackMC))
692 continue;
693
694 // only consider charged particles
695 if (part->Charge() == 0)
696 continue;
697
698 // only look at particles in eta-acceptance of the central barrel
699 if (TMath::Abs(part->Eta()) >= 0.9)
700 continue;
701
702 partMC.Add(part);
703 }
704
705 // sort, starting from highest pt
706 partMC.Sort(kSortAscending);
707
708 Int_t nTracksInStack[4][90] = { { 0 } };
709
710 // iterate over accepted tracks
711 TIter nextMcPart(&partMC);
712 while (AliVParticle *part = (AliMCParticle*) (nextMcPart())) {
713 Float_t pt = part->Pt();
714 Float_t eta = part->Eta();
715 Float_t phi = part->Phi();
716
717 Int_t phiBin = (Int_t) (phi / (2 / 18. *TMath::Pi()));
718 // ??? use actual geometry
719 Int_t etaBin = (Int_t) ((eta + .9) / (1.8 / 5.));
720
721 Int_t pdgCode = ((AliMCParticle*) part)->Particle()->GetPdgCode();
722 if (fDebug > 2)
723 printf("pt = %f, eta = %f, phi = %f, phiBin = %d, etaBin = %d, pdgCode = %d\n",
724 pt, eta, phi, phiBin, etaBin, pdgCode);
725
726 // increment number of tracks in stack
727 Int_t bin = phiBin * 5 + etaBin;
728 ++nTracksInStack[0][bin];
729
730 // if minimum number reached and
731 // the track has pt above threshold
732 // trigger fired
733 if ((nTracksInStack[0][bin] >= 3) &&
734 (pt >= 3.))
735 MarkTrigger(kTrgMC3x3Vtx);
736
737 // only propagate particles that do not decay before ???
738
739 // propagate to TRD
740 // why not use track references?
741 Float_t rTrack = pt / .15;
742 Float_t rTrd = 3.;
743 Float_t phiTrd = phi + TMath::ASin(rTrd / 2 / rTrack);
744 if (phiTrd < 0)
745 phiTrd += 2*TMath::Pi();
746 else if (phiTrd > 2*TMath::Pi())
747 phiTrd -= 2*TMath::Pi();
748 Float_t etaTrd = eta;
749
750 Int_t secTrd = (Int_t) (phiTrd / (2.*TMath::Pi()/18.));
751 Int_t stackTrd = (Int_t) ((etaTrd+0.9) / (1.8/5.));
752 Int_t binTrd = secTrd * 5 + stackTrd;
753
754 // increment number of tracks in stack
755 ++nTracksInStack[1][binTrd];
756 if (gRandom->Uniform() < fGlobalEfficiencyGTU)
757 ++nTracksInStack[2][binTrd];
758 if (gRandom->Uniform() < GetEfficiencyTRD(pt, eta, phi))
759 ++nTracksInStack[3][binTrd];
760
761 // if minimum number reached and
762 // the track has pt above threshold
763 // trigger fired
764 if (pt >= 3.) {
765 if (nTracksInStack[1][binTrd] >= 3)
766 MarkTrigger(kTrgMC3x3TRD);
767 if (nTracksInStack[2][binTrd] >= 3)
768 MarkTrigger(kTrgMC3x3TRDeff);
769 if (nTracksInStack[3][binTrd] >= 3)
770 MarkTrigger(kTrgMC3x3TRDeffmap);
771 }
772 }
773
774 return kTRUE;
775}
776
777Bool_t AliAnalysisTaskJetsTriggerTRD::AcceptTrackMC(Int_t track) const
778{
779 // only consider primaries
780 if (!MCEvent()->IsPhysicalPrimary(track))
781 return kFALSE;
782
783 const AliVParticle *part = MCEvent()->GetTrack(track);
784
785 // only consider charged particles
786 if (part->Charge() == 0)
787 return kFALSE;
788
789 // only look at particles in eta-acceptance of the central barrel
790 if (TMath::Abs(part->Eta()) >= 0.9)
791 return kFALSE;
792
e7808ddf 793 return kTRUE;
794}
795
ebdf1dbb 796
7032c043 797// ----- histogram management -----
798TH1* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
799 Int_t xbins, Float_t xmin, Float_t xmax,
800 Int_t binType)
801{
802 TString hName;
803 hName.Form("%s_%s", fShortTaskId, hid);
804 hName.ToLower();
805 TH1 *h = 0x0;
806 if (binType == 0)
807 h = new TH1I(hName.Data(), title,
808 xbins, xmin, xmax);
809 else
810 h = new TH1F(hName.Data(), title,
811 xbins, xmin, xmax);
812 GetHistogram(hist) = h;
813 fOutputList->Add(h);
814 return h;
815}
816
817TH2* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
818 Int_t xbins, Float_t xmin, Float_t xmax,
819 Int_t ybins, Float_t ymin, Float_t ymax,
820 Int_t binType)
821{
822 TString hName;
823 hName.Form("%s_%s", fShortTaskId, hid);
824 hName.ToLower();
825 TH1 *h = 0x0;
826 if (binType == 0)
827 h = GetHistogram(hist) = new TH2I(hName.Data(), title,
828 xbins, xmin, xmax,
829 ybins, ymin, ymax);
830 else
831 h = GetHistogram(hist) = new TH2F(hName.Data(), title,
832 xbins, xmin, xmax,
833 ybins, ymin, ymax);
834 fOutputList->Add(h);
835 return (TH2*) h;
836}
837
838TH3* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
839 Int_t xbins, Float_t xmin, Float_t xmax,
840 Int_t ybins, Float_t ymin, Float_t ymax,
841 Int_t zbins, Float_t zmin, Float_t zmax,
842 Int_t binType)
843{
844 TString hName;
845 hName.Form("%s_%s", fShortTaskId, hid);
846 hName.ToLower();
847 TH1 *h = 0x0;
848 if (binType == 0)
849 h = GetHistogram(hist) = new TH3I(hName.Data(), title,
850 xbins, xmin, xmax,
851 ybins, ymin, ymax,
852 zbins, zmin, zmax);
853 else
854 h = GetHistogram(hist) = new TH3F(hName.Data(), title,
855 xbins, xmin, xmax,
856 ybins, ymin, ymax,
857 zbins, zmin, zmax);
858 fOutputList->Add(h);
859 return (TH3*) h;
860}