909a4ac15e37f19df72e028a4da9384580c5bef0
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetResponse.cxx
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
25 ClassImp(AliAnalysisTaskJetResponse)
26
27 AliAnalysisTaskJetResponse::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]),
49   fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
50   fHistDeltaEtaDeltaPhiLeadingJet(0x0),
51   fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
52   fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
53   fHistPtPtExtra(0x0),
54   fHistPtResponse(new TH2F*[fkEvtClasses]),
55   fHistPtSmearing(new TH2F*[fkEvtClasses]),
56   fHistdR(new TH2F*[fkEvtClasses]),
57   fHistArea(new TH2F*[fkEvtClasses])
58 {
59   // default Constructor
60
61   fJetBranchName[0] = "";
62   fJetBranchName[1] = "";
63
64   fListJets[0] = new TList;
65   fListJets[1] = new TList;
66 }
67
68 AliAnalysisTaskJetResponse::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]),
90   fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
91   fHistDeltaEtaDeltaPhiLeadingJet(0x0),
92   fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
93   fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
94   fHistPtPtExtra(0x0),
95   fHistPtResponse(new TH2F*[fkEvtClasses]),
96   fHistPtSmearing(new TH2F*[fkEvtClasses]),
97   fHistdR(new TH2F*[fkEvtClasses]),
98   fHistArea(new TH2F*[fkEvtClasses])
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
111 AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
112 {
113   delete fListJets[0];
114   delete fListJets[1];
115 }
116
117 void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
118 {
119   fJetBranchName[0] = branch1;
120   fJetBranchName[1] = branch2;
121 }
122
123 void 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
133 void 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.);
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);
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]);
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]);
199   }
200 }
201
202 void 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];
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
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]) {
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     }
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
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 {
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);
320
321           fHistPtPtExtra->Fill(pt0, pt1);
322
323           if (eventClass > -1 && eventClass < fkEvtClasses){
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);
333           }
334        }
335     }
336   }
337
338   PostData(1, fOutputList);
339 }
340
341 void 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 }