b6af76ac9405e6bc9fbef0674224d46c1f75ea59
[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 {
46   // default ctor
47
48   DefineOutput(1, TList::Class());
49 }
50
51 AliAnalysisTaskJetsTriggerTRD::~AliAnalysisTaskJetsTriggerTRD()
52 {
53   // dtor
54
55 }
56
57 void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects()
58 {
59   // create user output objects
60
61   // setup list
62   OpenFile(1);
63   fOutputList = new TList();
64   fOutputList->SetOwner();
65
66   // setup histograms
67   TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
68                                kStatLast-1, .5, kStatLast-.5);
69   histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
70   histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
71   histStat->GetXaxis()->SetBinLabel(ID(kStatMB));
72
73   AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts",
74                100, -.5, 99.5);
75
76   AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts",
77                100, 0., 25.);
78
79   AddHistogram(ID(kHistJetPt), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
80                fNoJetPtBins, 0., fJetPtBinMax);
81   AddHistogram(ID(kHistJetPtITS), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
82                fNoJetPtBins, 0., fJetPtBinMax);
83   AddHistogram(ID(kHistJetPt3x3), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
84                fNoJetPtBins, 0., fJetPtBinMax);
85
86   AddHistogram(ID(kHistJetPtEMC), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
87                fNoJetPtBins, 0., fJetPtBinMax);
88   AddHistogram(ID(kHistJetPtHJT), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
89                fNoJetPtBins, 0., fJetPtBinMax);
90
91   AddHistogram(ID(kHistJetPtNoTracks3), "number of tracks above 3 GeV;p_{T}^{jet};no. of tracks",
92                fNoJetPtBins, 0., fJetPtBinMax,
93                20, -.5, 19.5);
94
95   PostData(1, fOutputList);
96 }
97
98 Bool_t AliAnalysisTaskJetsTriggerTRD::Notify()
99 {
100   // actions to be taken upon notification about input file change
101
102   return AliAnalysisTaskSE::Notify();
103 }
104
105 void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
106 {
107   // actual work
108
109   // setup pointers to input data (null if unavailable)
110   // mcEvent:  MC input
111   // esdEvent: ESD input
112   // outEvent: AOD output
113   // aodEvent: AOD input if available, otherwise AOD output
114   // AliMCEvent  *mcEvent   = this->MCEvent();
115   AliESDEvent *esdEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
116   AliAODEvent* outEvent  = this->AODEvent();
117   AliAODEvent *aodEvent  = outEvent;
118   if (dynamic_cast<AliAODEvent*>(this->InputEvent()))
119     aodEvent = (AliAODEvent*) (this->InputEvent());
120
121   if ((fDebug > 0) && esdEvent)
122     printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile());
123
124   // Int_t nTracksMC  = 0; // no. of MC tracks
125   // Int_t nTracksESD = 0; // no. of global ESD tracks
126   Int_t nTracksGTU = esdEvent ? esdEvent->GetNumberOfTrdTracks() : 0; // no. of GTU tracks
127
128   // Int_t nTracks[6][90]; // tracks above lower pt threshold, counted stack-wise
129   // memset(nTracks, 0, sizeof(Int_t)*6*90);
130   // Int_t nMax[6] = { 0 };
131
132   TList partMC;
133   TList partEsd;
134   TList partGtu;
135
136   // Float_t leadingJetPtMC = 0.; // leading jet energy from MC information
137   Float_t leadingJetPtRec = 0.; // leading jet energy from AOD information
138
139   // record number of sampled events and detect trigger contributions
140   FillH1(kHistStat, kStatSeen);
141   if (!DetectTriggers()) {
142     AliError("Failed to detect the triggers");
143     return;
144   }
145
146   if (!IsTrigger(kTrgInt))
147     return;
148   FillH1(kHistStat, kStatUsed);
149
150   for (Int_t iTrack = 0; iTrack < nTracksGTU; ++iTrack) {
151     AliVTrdTrack *trk = esdEvent->GetTrdTrack(iTrack);
152     FillH1(kHistTrackGTU, trk->Pt());
153   }
154
155   AliAODJet *leadJet = 0x0;
156   AliAODJet *subleadJet = 0x0;
157
158   if (aodEvent) {
159     TClonesArray *jetArray = dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
160     if (jetArray) {
161       Int_t nJets = jetArray->GetEntriesFast();
162       FillH1(kHistNoJets, nJets);
163       for (Int_t iJet = 0; iJet < nJets; ++iJet) {
164         AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet];
165
166         // check contributing tracks
167         Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
168         Int_t nJetTracks3 = 0;
169         Int_t iLeadingTrack = -1;
170         Float_t ptLeadingTrack = 0.;
171         for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
172           AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
173           // count tracks above 3 GeV/c
174           if (track->Pt() > 3.)
175             ++nJetTracks3;
176           // find the leading track
177           if (track->Pt() > ptLeadingTrack) {
178             ptLeadingTrack = track->Pt();
179             iLeadingTrack = iTrack;
180           }
181         }
182
183         // retrieve the leading track
184         AliAODTrack *leadingTrack = (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
185
186         if (TMath::Abs(jet->Eta()) < 0.5) {
187
188           // find leading jet
189           if (TMath::Abs(jet->Pt()) > leadingJetPtRec) {
190             leadingJetPtRec = TMath::Abs(jet->Pt());
191             subleadJet = leadJet;
192             leadJet = jet;
193           }
194
195           FillH1(kHistJetPt, jet->Pt());
196
197           // discard the jet if the leading track has no ITS contribution with a hit in any SPD layer
198           if (leadingTrack &&
199               (leadingTrack->GetFlags() & AliVTrack::kITSrefit) &&
200               (leadingTrack->HasPointOnITSLayer(0) || leadingTrack->HasPointOnITSLayer(1)))
201             FillH1(kHistJetPtITS, jet->Pt());
202
203           // discard the jet if it does not have 3 tracks above 3 GeV/c
204           FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
205           if (nJetTracks3 > 2)
206             FillH1(kHistJetPt3x3, jet->Pt());
207
208           if (IsTrigger(kTrgEMC))
209             FillH1(kHistJetPtEMC, jet->Pt());
210           if (IsTrigger(kTrgHJT))
211             FillH1(kHistJetPtHJT, jet->Pt());
212         }
213       }
214     }
215     else {
216       printf("no jet array found\n");
217     }
218   }
219
220   PostData(1, fOutputList);
221 }
222
223 void AliAnalysisTaskJetsTriggerTRD::Terminate(const Option_t * /* option */)
224 {
225   // actions at task termination
226
227 }
228
229 Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
230 {
231   fTriggerMask = 0;
232
233   AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
234   TString trgClasses = InputEvent()->GetFiredTriggerClasses();
235
236   if (fDebug > 2)
237     printf("trg: %8s %8s %8s %8s %8s %8s %8s (%s)\n",
238            (physSel & AliVEvent::kAnyINT)  ? "kAnyINT"  : " ",
239            (physSel & AliVEvent::kCINT5)   ? "kCINT5" : " ",
240            (physSel & AliVEvent::kINT7)    ? "kINT7" : " ",
241            (physSel & AliVEvent::kINT8)    ? "kINT8" : " ",
242            (physSel & AliVEvent::kMB)      ? "kMB"   : " ",
243            (physSel & AliVEvent::kEMC7)    ? "kEMC7" : " ",
244            (physSel & AliVEvent::kTRD)     ? "kTRD" : " ",
245            trgClasses.Data()
246            );
247
248   // physics selection
249   if (physSel & (AliVEvent::kAnyINT))
250     MarkTrigger(kTrgInt);
251
252   // trigger classes
253   if (trgClasses.Contains("CINT7WUHJT") || trgClasses.Contains("CINT8WUHJT"))
254     MarkTrigger(kTrgHJT);
255   if (trgClasses.Contains("CEMC7WU-") || trgClasses.Contains("CEMC8WU-"))
256     MarkTrigger(kTrgEMC);
257
258   return kTRUE;
259 }
260
261 // ----- histogram management -----
262 TH1* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
263                                                  Int_t xbins, Float_t xmin, Float_t xmax,
264                                                  Int_t binType)
265 {
266   TString hName;
267   hName.Form("%s_%s", fShortTaskId, hid);
268   hName.ToLower();
269   TH1 *h = 0x0;
270   if (binType == 0)
271     h = new TH1I(hName.Data(), title,
272                  xbins, xmin, xmax);
273   else
274     h = new TH1F(hName.Data(), title,
275                  xbins, xmin, xmax);
276   GetHistogram(hist) = h;
277   fOutputList->Add(h);
278   return h;
279 }
280
281 TH2* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
282                                                  Int_t xbins, Float_t xmin, Float_t xmax,
283                                                  Int_t ybins, Float_t ymin, Float_t ymax,
284                                                  Int_t binType)
285 {
286   TString hName;
287   hName.Form("%s_%s", fShortTaskId, hid);
288   hName.ToLower();
289   TH1 *h = 0x0;
290   if (binType == 0)
291     h = GetHistogram(hist) = new TH2I(hName.Data(), title,
292                                      xbins, xmin, xmax,
293                                      ybins, ymin, ymax);
294   else
295     h = GetHistogram(hist) = new TH2F(hName.Data(), title,
296                                      xbins, xmin, xmax,
297                                      ybins, ymin, ymax);
298   fOutputList->Add(h);
299   return (TH2*) h;
300 }
301
302 TH3* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
303                                                  Int_t xbins, Float_t xmin, Float_t xmax,
304                                                  Int_t ybins, Float_t ymin, Float_t ymax,
305                                                  Int_t zbins, Float_t zmin, Float_t zmax,
306                                                  Int_t binType)
307 {
308   TString hName;
309   hName.Form("%s_%s", fShortTaskId, hid);
310   hName.ToLower();
311   TH1 *h = 0x0;
312   if (binType == 0)
313     h = GetHistogram(hist) = new TH3I(hName.Data(), title,
314                                      xbins, xmin, xmax,
315                                      ybins, ymin, ymax,
316                                      zbins, zmin, zmax);
317   else
318     h = GetHistogram(hist) = new TH3F(hName.Data(), title,
319                                      xbins, xmin, xmax,
320                                      ybins, ymin, ymax,
321                                      zbins, zmin, zmax);
322   fOutputList->Add(h);
323   return (TH3*) h;
324 }