Updated to Jet embedding to add real jets (J. Klein)
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetResponse.cxx
CommitLineData
2a436fa1 1#include "TChain.h"
2#include "TTree.h"
3#include "TMath.h"
4#include "TH1F.h"
5#include "TH2F.h"
6#include "TCanvas.h"
7
8#include "AliLog.h"
9
10#include "AliAnalysisTask.h"
11#include "AliAnalysisManager.h"
12
13#include "AliVEvent.h"
14#include "AliESDEvent.h"
15#include "AliESDInputHandler.h"
16#include "AliCentrality.h"
17#include "AliAnalysisHelperJetTasks.h"
18#include "AliInputEventHandler.h"
19
20#include "AliAODEvent.h"
21#include "AliAODJet.h"
22
23#include "AliAnalysisTaskJetResponse.h"
24
25ClassImp(AliAnalysisTaskJetResponse)
26
27AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
28 AliAnalysisTaskSE(),
29 fESD(0x0),
30 fAOD(0x0),
31 fOfflineTrgMask(AliVEvent::kAny),
32 fMinContribVtx(1),
33 fVtxZMin(-8.),
34 fVtxZMax(8.),
35 fEvtClassMin(0),
36 fEvtClassMax(10),
37 fCentMin(0.),
38 fCentMax(100.),
39 fJetEtaMin(-.5),
40 fJetEtaMax(.5),
41 fJetDeltaEta(.4),
42 fJetDeltaPhi(.4),
43 fkNbranches(2),
44 fkEvtClasses(10),
45 fOutputList(0x0),
46 fHistEvtSelection(0x0),
47 fHistPtLeadingJet(new TH1F*[fkNbranches]),
48 fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
a11972e9 49 fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
2a436fa1 50 fHistDeltaEtaDeltaPhiLeadingJet(0x0),
a11972e9 51 fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
52 fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
2a436fa1 53 fHistPtPtExtra(0x0),
54 fHistPtResponse(new TH2F*[fkEvtClasses]),
a11972e9 55 fHistPtSmearing(new TH2F*[fkEvtClasses]),
56 fHistdR(new TH2F*[fkEvtClasses]),
57 fHistArea(new TH2F*[fkEvtClasses])
2a436fa1 58{
59 // default Constructor
60
61 fJetBranchName[0] = "";
62 fJetBranchName[1] = "";
63
64 fListJets[0] = new TList;
65 fListJets[1] = new TList;
66}
67
68AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
69 AliAnalysisTaskSE(name),
70 fESD(0x0),
71 fAOD(0x0),
72 fOfflineTrgMask(AliVEvent::kAny),
73 fMinContribVtx(1),
74 fVtxZMin(-8.),
75 fVtxZMax(8.),
76 fEvtClassMin(0),
77 fEvtClassMax(10),
78 fCentMin(0.),
79 fCentMax(100.),
80 fJetEtaMin(-.5),
81 fJetEtaMax(.5),
82 fJetDeltaEta(.4),
83 fJetDeltaPhi(.4),
84 fkNbranches(2),
85 fkEvtClasses(10),
86 fOutputList(0x0),
87 fHistEvtSelection(0x0),
88 fHistPtLeadingJet(new TH1F*[fkNbranches]),
89 fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
a11972e9 90 fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
2a436fa1 91 fHistDeltaEtaDeltaPhiLeadingJet(0x0),
a11972e9 92 fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
93 fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
2a436fa1 94 fHistPtPtExtra(0x0),
95 fHistPtResponse(new TH2F*[fkEvtClasses]),
a11972e9 96 fHistPtSmearing(new TH2F*[fkEvtClasses]),
97 fHistdR(new TH2F*[fkEvtClasses]),
98 fHistArea(new TH2F*[fkEvtClasses])
2a436fa1 99{
100 // Constructor
101
102 fJetBranchName[0] = "";
103 fJetBranchName[1] = "";
104
105 fListJets[0] = new TList;
106 fListJets[1] = new TList;
107
108 DefineOutput(1, TList::Class());
109}
110
111AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
112{
113 delete fListJets[0];
114 delete fListJets[1];
115}
116
117void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
118{
119 fJetBranchName[0] = branch1;
120 fJetBranchName[1] = branch2;
121}
122
123void AliAnalysisTaskJetResponse::Init()
124{
125
126 // check for jet branches
127 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
128 AliError("Jet branch name not set.");
129 }
130
131}
132
133void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
134{
135 // Create histograms
136 // Called once
137
138 fHistEvtSelection = new TH1F("fHistEvtSelection", "event selection", 5, -0.5, 5.5);
139 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
140 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
141 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
142 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
143 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
144
145 for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
146 fHistPtLeadingJet[iBranch] = new TH1F(Form("fHistPtLeadingJet%i", iBranch),
147 Form("p_{T} of leading jet from branch %i;p_{T} (GeV/c);counts", iBranch),
148 300, 0., 300.);
149 fHistPtLeadingJet[iBranch]->SetMarkerStyle(kFullCircle);
150
151 fHistEtaPhiLeadingJet[iBranch] = new TH2F(Form("fHistEtaPhiLeadingJet%i", iBranch),
152 Form("#eta - #phi of leading jet from branch %i;#eta;#phi", iBranch),
153 100, -2., 2., 100, 0., 2*TMath::Pi());
154
155 }
156 fHistDeltaEtaDeltaPhiLeadingJet = new TH2F("fHistDeltaEtaDeltaPhiLeadingJet",
157 "#Delta#eta - #Delta#phi of leading jet;#Delta#eta;#Delta#phi",
158 100, -4., 4., 100, -2.*TMath::Pi(), 2*TMath::Pi());
159
160 fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
161 300, 0., 300., 300, 0., 300.);
162
163 for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
164 fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
165 300, 0., 300., 300, 0., 300.);
166 fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
167 200, -50., 150., 300, 0., 300.);
a11972e9 168 fHistdR[iEvtClass] = new TH2F(Form("hist_dR%i",iEvtClass), "", 200, -50.,150., 200, -1.,1.);
169 fHistArea[iEvtClass] = new TH2F(Form("hist_Area%i",iEvtClass), "", 200, -50.,150., 100, 0.,1.);
170
171 fHistEtaPhiLeadingJetCut[iEvtClass] = new TH2F(Form("fHistEtaPhiLeadingJetCut%i", iEvtClass),
172 Form("#eta - #phi of leading jet from event class %i;#eta;#phi", iEvtClass),
173 100, -2., 2., 100, 0., 2*TMath::Pi());
174 fHistDeltaEtaEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaLeadingJet%i", iEvtClass),
175 Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#eta", iEvtClass),
176 100, -1., 1., 100, -.5, .5);
177 fHistDeltaPtEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaLeadingJet%i", iEvtClass),
178 Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#p_{T}", iEvtClass),
179 100, -1., 1., 200, -50., 150);
2a436fa1 180 }
181
182 OpenFile(1);
183 fOutputList = new TList;
184 fOutputList->Add(fHistEvtSelection);
185 for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
186 fOutputList->Add(fHistPtLeadingJet[iBranch]);
187 fOutputList->Add(fHistEtaPhiLeadingJet[iBranch]);
188 }
189 fOutputList->Add(fHistDeltaEtaDeltaPhiLeadingJet);
190 fOutputList->Add(fHistPtPtExtra);
191 for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
192 fOutputList->Add(fHistPtResponse[iEvtClass]);
193 fOutputList->Add(fHistPtSmearing[iEvtClass]);
a11972e9 194 fOutputList->Add(fHistdR[iEvtClass]);
195 fOutputList->Add(fHistArea[iEvtClass]);
196 fOutputList->Add(fHistEtaPhiLeadingJetCut[iEvtClass]);
197 fOutputList->Add(fHistDeltaEtaEtaLeadingJet[iEvtClass]);
198 fOutputList->Add(fHistDeltaPtEtaLeadingJet[iEvtClass]);
2a436fa1 199 }
200}
201
202void AliAnalysisTaskJetResponse::UserExec(Option_t *)
203{
204 // load events, apply event cuts, then compare leading jets
205
206 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
207 AliError("Jet branch name not set.");
208 return;
209 }
210
211 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
212 if (!fESD) {
213 AliError("ESD not available");
214 return;
215 }
216 fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
217 if (!fAOD) {
218 AliError("AOD not available");
219 return;
220 }
221
222 fHistEvtSelection->Fill(1); // number of events before event selection
223
224 // physics selection
225 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
226 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
227 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
228 Printf(" Trigger Selection: event REJECTED ... ");
229 fHistEvtSelection->Fill(2);
230 PostData(1, fOutputList);
231 return;
232 }
233
234 // vertex selection
235 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
236 Int_t nTracksPrim = primVtx->GetNContributors();
237 if ((nTracksPrim < fMinContribVtx) ||
238 (primVtx->GetZ() < fVtxZMin) ||
239 (primVtx->GetZ() > fVtxZMax) ){
240 Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
241 fHistEvtSelection->Fill(3);
242 PostData(1, fOutputList);
243 return;
244 }
245
246 // event class selection (from jet helper task)
247 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
248 Printf("Event class %d", eventClass);
249 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
250 fHistEvtSelection->Fill(4);
251 PostData(1, fOutputList);
252 return;
253 }
254
255 // centrality selection
256 AliCentrality *cent = fESD->GetCentrality();
257 Float_t centValue = cent->GetCentralityPercentile("TRK");
258 printf("centrality: %f\n", centValue);
259 if (centValue < fCentMin || centValue > fCentMax){
260 fHistEvtSelection->Fill(5);
261 PostData(1, fOutputList);
262 return;
263 }
264
265 fHistEvtSelection->Fill(0); // accepted events
266
267
268 TClonesArray *aodJets[2];
a11972e9 269 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
270 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
2a436fa1 271 AliAODJet *leadingJet[2] = { 0x0, 0x0 };
272
273 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
274 fListJets[iJetType]->Clear();
275 if (!aodJets[iJetType])
276 continue;
277 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
278 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
279 if (jet)
280 fListJets[iJetType]->Add(jet);
281 }
282 fListJets[iJetType]->Sort();
283 leadingJet[iJetType] = (AliAODJet*) fListJets[iJetType]->First();
284
285 if (leadingJet[iJetType])
286 fHistEtaPhiLeadingJet[iJetType]->Fill(leadingJet[iJetType]->Eta(), leadingJet[iJetType]->Phi());
287 }
288
289 // check if leading jets are close in eta-phi and compare their Pt
290 if (leadingJet[0] && leadingJet[1]) {
a11972e9 291 Float_t deltaEta = leadingJet[0]->Eta() - leadingJet[1]->Eta();
292 Float_t deltaPhi = TVector2::Phi_mpi_pi(leadingJet[0]->Phi() - leadingJet[1]->Phi());
293
294 if (eventClass > -1 && eventClass < fkEvtClasses){
295
296 if(leadingJet[1]->Eta()<fJetEtaMax && leadingJet[1]->Eta()>fJetEtaMin){
297 fHistEtaPhiLeadingJetCut[eventClass]->Fill(leadingJet[1]->Eta(), leadingJet[1]->Phi());
298 }
299 }
2a436fa1 300
301 // leading jets in "jet" acceptance
302 if(leadingJet[0]->Eta()>fJetEtaMax || leadingJet[0]->Eta()<fJetEtaMin ||
303 leadingJet[1]->Eta()>fJetEtaMax || leadingJet[1]->Eta()<fJetEtaMin){
304 Printf("Jet not in eta acceptance.");
305 }
306 else{
307 // check association of jets
2a436fa1 308 fHistDeltaEtaDeltaPhiLeadingJet->Fill(deltaEta, deltaPhi);
309
310 if (TMath::Abs(deltaEta) > fJetDeltaEta || (TMath::Pi() - TMath::Abs(TMath::Abs(deltaPhi) - TMath::Pi())) > fJetDeltaPhi)
311 printf("leading jets two far apart\n");
312 else {
a11972e9 313 Float_t pt0 = leadingJet[0]->Pt();
314 Float_t pt1 = leadingJet[1]->Pt();
315 Float_t dPt = pt1-pt0;
316 Float_t jetarea = leadingJet[1]->EffectiveAreaCharged();
317
318 fHistPtLeadingJet[0]->Fill(pt0);
319 fHistPtLeadingJet[1]->Fill(pt1);
2a436fa1 320
a11972e9 321 fHistPtPtExtra->Fill(pt0, pt1);
2a436fa1 322
323 if (eventClass > -1 && eventClass < fkEvtClasses){
a11972e9 324 fHistPtResponse[eventClass]->Fill(pt1, pt0);
325 fHistPtSmearing[eventClass]->Fill(dPt, pt0);
326
327 Float_t dR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
328 fHistdR[eventClass]->Fill(dPt, dR);
329 fHistArea[eventClass]->Fill(dPt, jetarea);
330
331 fHistDeltaEtaEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), deltaEta);
332 fHistDeltaPtEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), dPt);
2a436fa1 333 }
334 }
335 }
336 }
337
338 PostData(1, fOutputList);
339}
340
341void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
342{
343 // Draw result to the screen
344 // Called once at the end of the query
345
346 if (!GetOutputData(1))
347 return;
348}