From Jochen:
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetsTriggerTRD.cxx
1 // ROOT
2 #include "TFile.h"
3 #include "TList.h"
4 #include "TH1.h"
5 #include "TH2.h"
6 #include "TH3.h"
7
8 // analysis framework
9 #include "AliAnalysisManager.h"
10 #include "AliInputEventHandler.h"
11 #include "AliVEvent.h"
12 #include "AliVTrdTrack.h"
13
14 // MC stuff
15 #include "AliMCEvent.h"
16 #include "AliGenPythiaEventHeader.h"
17
18 // ESD stuff
19 #include "AliESDEvent.h"
20 #include "AliESDInputHandler.h"
21 #include "AliESDtrack.h"
22 #include "AliESDTrdTrack.h"
23 #include "AliESDTrdTracklet.h"
24 #include "AliESDTrdTrigger.h"
25
26 // AOD stuff
27 #include "AliAODEvent.h"
28 #include "AliAODJet.h"
29 #include "AliAODTrack.h"
30
31 // jet tasks
32 #include "AliAnalysisTaskJetServices.h"
33 #include "AliAnalysisHelperJetTasks.h"
34
35 #include "AliAnalysisTaskJetsTriggerTRD.h"
36
37 AliAnalysisTaskJetsTriggerTRD::AliAnalysisTaskJetsTriggerTRD(const char *name) :
38   AliAnalysisTaskSE(name),
39   fTriggerMask(0),
40   fOutputList(),
41   fHist(),
42   fShortTaskId("jets_trg_trd"),
43   fNoJetPtBins(80),
44   fJetPtBinMax(400),
45   fXsection(0.),
46   fAvgTrials(0.),
47   fPtHard(0.)
48 {
49   // default ctor
50
51   DefineOutput(1, TList::Class());
52 }
53
54 AliAnalysisTaskJetsTriggerTRD::~AliAnalysisTaskJetsTriggerTRD()
55 {
56   // dtor
57
58 }
59
60 void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects()
61 {
62   // create user output objects
63
64   // setup list
65   OpenFile(1);
66   fOutputList = new TList();
67   fOutputList->SetOwner();
68
69   // setup histograms
70   TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
71                                kStatLast-1, .5, kStatLast-.5);
72   histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
73   histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
74   histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
75   histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
76
77   AddHistogram(ID(kHistJetPtMC), "leading jet spectrum (MC, |#eta| < 0.5);p_{T} (GeV/c);counts",
78                fNoJetPtBins, 0., fJetPtBinMax);
79
80   AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts",
81                400, -.5, 399.5);
82
83   AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts",
84                100, 0., 25.);
85
86   AddHistogram(ID(kHistNPtMin),
87                "rejection;p_{T}^{min};N_{trk};trigger",
88                100, 0., 10.,
89                20, 0., 20.,
90                kTrgLast - 1, .5, kTrgLast - .5);
91
92   AddHistogram(ID(kHistLeadJetPt),
93                "leading jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
94                fNoJetPtBins, 0., fJetPtBinMax,
95                kTrgLast - 1, .5, kTrgLast - .5);
96   AddHistogram(ID(kHistJetPt),
97                "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
98                fNoJetPtBins, 0., fJetPtBinMax,
99                kTrgLast - 1, .5, kTrgLast - .5);
100   AddHistogram(ID(kHistJetPtITS),
101                "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
102                fNoJetPtBins, 0., fJetPtBinMax,
103                kTrgLast - 1, .5, kTrgLast - .5);
104   AddHistogram(ID(kHistJetPt3x3),
105                "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
106                fNoJetPtBins, 0., fJetPtBinMax,
107                kTrgLast - 1, .5, kTrgLast - .5);
108
109   for (Int_t iHist = kHistLeadJetPt; iHist <= kHistJetPt3x3; ++iHist) {
110     TH1 *h = GetHistogram(Hist_t (iHist));
111     h->GetYaxis()->SetBinLabel(ID(kTrgMinBias));
112     h->GetYaxis()->SetBinLabel(ID(kTrgInt7));
113     h->GetYaxis()->SetBinLabel(ID(kTrgInt8));
114     h->GetYaxis()->SetBinLabel(ID(kTrgEMC7));
115     h->GetYaxis()->SetBinLabel(ID(kTrgEMC8));
116     h->GetYaxis()->SetBinLabel(ID(kTrgInt7WUHJT));
117     h->GetYaxis()->SetBinLabel(ID(kTrgInt8WUHJT));
118     h->GetYaxis()->SetBinLabel(ID(kTrgEMC7WUHJT));
119     h->GetYaxis()->SetBinLabel(ID(kTrgEMC8WUHJT));
120   }
121
122   AddHistogram(ID(kHistJetPtNoTracks3),
123                "number of tracks above 3 GeV;p_{T}^{jet};no. of tracks",
124                fNoJetPtBins, 0., fJetPtBinMax,
125                40, -.5, 39.5);
126
127   PostData(1, fOutputList);
128 }
129
130 Bool_t AliAnalysisTaskJetsTriggerTRD::Notify()
131 {
132   // actions to be taken upon notification about input file change
133
134   // ??? check ???
135
136   fXsection = 0.;
137   fAvgTrials = 1.;
138
139   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
140
141   if (tree) {
142     TFile *curfile = tree->GetCurrentFile();
143     if (!curfile) {
144       Error("Notify","No current file");
145       return kFALSE;
146     }
147
148     Float_t nEntries = (Float_t) tree->GetTree()->GetEntries();
149     if (!AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(), fXsection, fAvgTrials)) {
150       AliError("retrieval of cross section failed");
151       // ??? what to set as cross section?
152     }
153
154     if (nEntries > 0.)
155       fAvgTrials /= nEntries;
156   }
157
158   return AliAnalysisTaskSE::Notify();
159 }
160
161 void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
162 {
163   // actual work
164
165   // setup pointers to input data (null if unavailable)
166   // mcEvent:  MC input
167   // esdEvent: ESD input
168   // outEvent: AOD output
169   // aodEvent: AOD input if available, otherwise AOD output
170   AliMCEvent  *mcEvent   = this->MCEvent();
171   AliESDEvent *esdEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
172   AliAODEvent* outEvent  = this->AODEvent();
173   AliAODEvent *aodEvent  = outEvent;
174   if (dynamic_cast<AliAODEvent*>(this->InputEvent()))
175     aodEvent = (AliAODEvent*) (this->InputEvent());
176
177   if ((fDebug > 0) && esdEvent)
178     printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile());
179
180   // Int_t nTracksMC  = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks
181   // Int_t nTracks    = InputEvent()->GetNumberOfTracks(); // no. of global tracks
182   Int_t nTracksTRD = InputEvent()->GetNumberOfTrdTracks(); // no. of GTU tracks
183
184   TList partMC;
185   TList partEsd;
186   TList partGtu;
187
188   Float_t leadingJetPtMC = 0.; // leading jet energy from MC information
189   Float_t leadingJetPtRec = 0.; // leading jet energy from AOD information
190
191   // record number of sampled events and detect trigger contributions
192   FillH1(kHistStat, kStatSeen);
193   if (!DetectTriggers()) {
194     AliError("Failed to detect the triggers");
195     return;
196   }
197
198   // only continue for events from interesting triggers
199   if (fTriggerMask == 0)
200     return;
201   FillH1(kHistStat, kStatTrg);
202
203   // no further technical requirements for the event at the moment
204   FillH1(kHistStat, kStatUsed);
205
206   // apply event cuts
207   const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
208   if (!vtx ||
209       (vtx->GetNContributors() < 3.) ||
210       (vtx->GetZ() > 10.))
211     return;
212   FillH1(kHistStat, kStatEvCuts);
213
214   // extract MC information
215   if (mcEvent) {
216     // check for PYTHIA event header
217     AliGenPythiaEventHeader *pythiaHeader =
218       dynamic_cast<AliGenPythiaEventHeader*> (mcEvent->GenEventHeader());
219     if (!pythiaHeader) {
220       AliWarning("MC event without PYTHIA event header!\n");
221     }
222     else {
223       fPtHard  = pythiaHeader->GetPtHard();
224       // Int_t nTrials = pythiaHeader->Trials();
225
226       // loop over jets from PYTHIA
227       for (Int_t iJet = 0; iJet < pythiaHeader->NTriggerJets(); iJet++) {
228         Float_t p[4];
229         pythiaHeader->TriggerJet(iJet, p);
230         TLorentzVector pJet(p);
231         Float_t pt  = pJet.Pt();
232         // only consider jets with |eta| < 0.5
233         Float_t eta = pJet.Eta();
234         if (TMath::Abs(eta) > 0.5)
235           continue;
236         if (pt > leadingJetPtMC)
237           leadingJetPtMC = pt;
238       }
239       // fill histogram for leading jet pt spectrum
240       FillH1(kHistJetPtMC, leadingJetPtMC, fXsection);
241     }
242   }
243
244   // loop over GTU tracks
245   for (Int_t iTrack = 0; iTrack < nTracksTRD; ++iTrack) {
246     AliVTrdTrack *trk = InputEvent()->GetTrdTrack(iTrack);
247     FillH1(kHistTrackGTU, TMath::Abs(trk->Pt()));
248     partGtu.Add(trk);
249   }
250   partGtu.Sort(kSortAscending);
251
252   Int_t nTracksPerStack[90] = { 0 };
253   Int_t nTracksPerStackMax = 0;
254
255   TIter nextPartGtu(&partGtu);
256   while (AliVTrdTrack *trdTrack = (AliVTrdTrack*) nextPartGtu()) {
257     // count no. of tracks in stack,
258     // check whether this number was reached before,
259     // if not store pt^min(n),
260     // i.e. pt of current track because of sorting
261
262     Int_t sec    = trdTrack->GetSector();
263     Int_t stack  = trdTrack->GetStack();
264
265     if ((sec > -1) && (sec < 18) &&
266         (stack > -1) && (stack < 5)) {
267       ++nTracksPerStack[5*sec + stack];
268       if (nTracksPerStack[5*sec + stack] > nTracksPerStackMax) {
269         ++nTracksPerStackMax;
270
271         for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
272           if (IsTrigger(Trigger_t (iTrigger)))
273               FillH3(kHistNPtMin,
274                      TMath::Abs(trdTrack->Pt()), nTracksPerStackMax, iTrigger);
275       }
276     }
277     else {
278       AliError(Form("Invalid sector or stack: %i %i",
279                     sec, stack));
280     }
281     
282   }
283
284   // loop over jets from AOD event
285   if (aodEvent) {
286     TClonesArray *jetArray =
287       dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
288     if (jetArray) {
289       Int_t nJets = jetArray->GetEntriesFast();
290       FillH1(kHistNoJets, nJets);
291
292       AliAODJet *leadJet = 0x0;
293       // AliAODJet *subleadJet = 0x0;
294
295       for (Int_t iJet = 0; iJet < nJets; ++iJet) {
296         AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet];
297         if (TMath::Abs(jet->Eta()) < 0.5) {
298           // check contributing tracks
299           Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
300           Int_t nJetTracks3 = 0;
301           AliAODTrack *leadingTrack = 0x0;
302           for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
303             AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
304
305             // count constituents above 3 GeV/c
306             if (track->Pt() > 3.)
307               ++nJetTracks3;
308
309             // find the leading track
310             if (!leadingTrack ||
311                 (track->Pt() > leadingTrack->Pt()))
312               leadingTrack = track;
313           }
314
315           // find leading jet
316           if (TMath::Abs(jet->Pt()) > leadingJetPtRec) {
317             leadingJetPtRec = TMath::Abs(jet->Pt());
318             // subleadJet = leadJet;
319             leadJet = jet;
320           }
321
322           // jet pt spectrum
323           for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
324             if (IsTrigger(Trigger_t (iTrigger)))
325               FillH1(kHistJetPt, jet->Pt(), iTrigger);
326
327           // integrated over all triggers
328           FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
329
330           // limit to jets with leading track having an ITS contribution
331           // with a hit in any SPD layer
332           if (leadingTrack &&
333               (leadingTrack->GetFlags() & AliVTrack::kITSrefit) &&
334               (leadingTrack->HasPointOnITSLayer(0) || leadingTrack->HasPointOnITSLayer(1)))
335             for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
336               if (IsTrigger(Trigger_t (iTrigger)))
337                 FillH1(kHistJetPtITS, jet->Pt(), iTrigger);
338
339           // limit to jets having 3 tracks above 3 GeV/c
340           if (nJetTracks3 > 2)
341             for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
342               if (IsTrigger(Trigger_t (iTrigger)))
343                 FillH1(kHistJetPt3x3, jet->Pt(), iTrigger);
344         }
345       }
346       // fill leading jet information
347       for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
348         if (IsTrigger(Trigger_t (iTrigger)))
349           FillH1(kHistLeadJetPt, leadJet ? leadJet->Pt() : 0., iTrigger);
350     }
351     else {
352       printf("no jet array found as branch %s\n", fJetBranchName);
353       aodEvent->Print();
354     }
355   } else {
356     printf("no AOD event found\n");
357   }
358
359   PostData(1, fOutputList);
360 }
361
362 void AliAnalysisTaskJetsTriggerTRD::Terminate(const Option_t * /* option */)
363 {
364   // actions at task termination
365
366 }
367
368 Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
369 {
370   fTriggerMask = 0;
371
372   AliInputEventHandler *inputHandler =
373     (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
374
375   AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) inputHandler->IsEventSelected();
376   TString trgClasses = InputEvent()->GetFiredTriggerClasses();
377
378   if (fDebug > 2)
379     printf("trg: %8s %8s %8s %8s %8s %8s %8s (%s)\n",
380            (physSel & AliVEvent::kAnyINT)  ? "kAnyINT"  : " ",
381            (physSel & AliVEvent::kCINT5)   ? "kCINT5" : " ",
382            (physSel & AliVEvent::kINT7)    ? "kINT7" : " ",
383            (physSel & AliVEvent::kINT8)    ? "kINT8" : " ",
384            (physSel & AliVEvent::kMB)      ? "kMB"   : " ",
385            (physSel & AliVEvent::kEMC7)    ? "kEMC7" : " ",
386            (physSel & AliVEvent::kTRD)     ? "kTRD" : " ",
387            trgClasses.Data()
388            );
389
390   // physics selection
391   if ((physSel & (AliVEvent::kMB)))
392     MarkTrigger(kTrgMinBias);
393
394   if ((physSel & (AliVEvent::kINT7)))
395     MarkTrigger(kTrgInt7);
396
397   if ((physSel & (AliVEvent::kINT8)))
398     MarkTrigger(kTrgInt8);
399
400   if ((physSel & (AliVEvent::kEMC7)))
401     MarkTrigger(kTrgEMC7);
402
403   if ((physSel & (AliVEvent::kEMC8)))
404     MarkTrigger(kTrgEMC8);
405
406   // for the triggered events we use the classes
407   if (trgClasses.Contains("CINT7WUHJT-"))
408     MarkTrigger(kTrgInt7WUHJT);
409
410   if (trgClasses.Contains("CINT8WUHJT-"))
411     MarkTrigger(kTrgInt8WUHJT);
412
413   if (trgClasses.Contains("CEMC7WUHJT-"))
414     MarkTrigger(kTrgEMC7WUHJT);
415
416   if (trgClasses.Contains("CEMC8WUHJT-"))
417     MarkTrigger(kTrgEMC8WUHJT);
418
419   return kTRUE;
420 }
421
422 // ----- histogram management -----
423 TH1* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
424                                                  Int_t xbins, Float_t xmin, Float_t xmax,
425                                                  Int_t binType)
426 {
427   TString hName;
428   hName.Form("%s_%s", fShortTaskId, hid);
429   hName.ToLower();
430   TH1 *h = 0x0;
431   if (binType == 0)
432     h = new TH1I(hName.Data(), title,
433                  xbins, xmin, xmax);
434   else
435     h = new TH1F(hName.Data(), title,
436                  xbins, xmin, xmax);
437   GetHistogram(hist) = h;
438   fOutputList->Add(h);
439   return h;
440 }
441
442 TH2* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
443                                                  Int_t xbins, Float_t xmin, Float_t xmax,
444                                                  Int_t ybins, Float_t ymin, Float_t ymax,
445                                                  Int_t binType)
446 {
447   TString hName;
448   hName.Form("%s_%s", fShortTaskId, hid);
449   hName.ToLower();
450   TH1 *h = 0x0;
451   if (binType == 0)
452     h = GetHistogram(hist) = new TH2I(hName.Data(), title,
453                                      xbins, xmin, xmax,
454                                      ybins, ymin, ymax);
455   else
456     h = GetHistogram(hist) = new TH2F(hName.Data(), title,
457                                      xbins, xmin, xmax,
458                                      ybins, ymin, ymax);
459   fOutputList->Add(h);
460   return (TH2*) h;
461 }
462
463 TH3* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
464                                                  Int_t xbins, Float_t xmin, Float_t xmax,
465                                                  Int_t ybins, Float_t ymin, Float_t ymax,
466                                                  Int_t zbins, Float_t zmin, Float_t zmax,
467                                                  Int_t binType)
468 {
469   TString hName;
470   hName.Form("%s_%s", fShortTaskId, hid);
471   hName.ToLower();
472   TH1 *h = 0x0;
473   if (binType == 0)
474     h = GetHistogram(hist) = new TH3I(hName.Data(), title,
475                                      xbins, xmin, xmax,
476                                      ybins, ymin, ymax,
477                                      zbins, zmin, zmax);
478   else
479     h = GetHistogram(hist) = new TH3F(hName.Data(), title,
480                                      xbins, xmin, xmax,
481                                      ybins, ymin, ymax,
482                                      zbins, zmin, zmax);
483   fOutputList->Add(h);
484   return (TH3*) h;
485 }