f1f8fce64df12c9b7f4eba1b252360110cbf3868
[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   fHistDeltaEtaDeltaPhiLeadingJet(0x0),
50   fHistPtPtExtra(0x0),
51   fHistPtResponse(new TH2F*[fkEvtClasses]),
52   fHistPtSmearing(new TH2F*[fkEvtClasses])
53 {
54   // default Constructor
55
56   fJetBranchName[0] = "";
57   fJetBranchName[1] = "";
58
59   fListJets[0] = new TList;
60   fListJets[1] = new TList;
61 }
62
63 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
64   AliAnalysisTaskSE(name),
65   fESD(0x0),
66   fAOD(0x0),
67   fOfflineTrgMask(AliVEvent::kAny),
68   fMinContribVtx(1),
69   fVtxZMin(-8.),
70   fVtxZMax(8.),
71   fEvtClassMin(0),
72   fEvtClassMax(10),
73   fCentMin(0.),
74   fCentMax(100.),
75   fJetEtaMin(-.5),
76   fJetEtaMax(.5),
77   fJetDeltaEta(.4),
78   fJetDeltaPhi(.4),
79   fkNbranches(2),
80   fkEvtClasses(10),
81   fOutputList(0x0),
82   fHistEvtSelection(0x0),
83   fHistPtLeadingJet(new TH1F*[fkNbranches]),
84   fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
85   fHistDeltaEtaDeltaPhiLeadingJet(0x0),
86   fHistPtPtExtra(0x0),
87   fHistPtResponse(new TH2F*[fkEvtClasses]),
88   fHistPtSmearing(new TH2F*[fkEvtClasses])
89 {
90   // Constructor
91
92   fJetBranchName[0] = "";
93   fJetBranchName[1] = "";
94
95   fListJets[0] = new TList;
96   fListJets[1] = new TList;
97
98   DefineOutput(1, TList::Class());
99 }
100
101 AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
102 {
103   delete fListJets[0];
104   delete fListJets[1];
105 }
106
107 void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
108 {
109   fJetBranchName[0] = branch1;
110   fJetBranchName[1] = branch2;
111 }
112
113 void AliAnalysisTaskJetResponse::Init()
114 {
115
116    // check for jet branches
117    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
118       AliError("Jet branch name not set.");
119    }
120
121 }
122
123 void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
124 {
125   // Create histograms
126   // Called once
127
128   fHistEvtSelection = new TH1F("fHistEvtSelection", "event selection", 5, -0.5, 5.5);
129   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
130   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
131   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
132   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
133   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
134
135   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
136     fHistPtLeadingJet[iBranch] = new TH1F(Form("fHistPtLeadingJet%i", iBranch),
137                                           Form("p_{T} of leading jet from branch %i;p_{T} (GeV/c);counts", iBranch),
138                                           300, 0., 300.);
139     fHistPtLeadingJet[iBranch]->SetMarkerStyle(kFullCircle);
140
141     fHistEtaPhiLeadingJet[iBranch] = new TH2F(Form("fHistEtaPhiLeadingJet%i", iBranch),
142                                               Form("#eta - #phi of leading jet from branch %i;#eta;#phi", iBranch),
143                                               100, -2., 2., 100, 0., 2*TMath::Pi());
144
145   }
146   fHistDeltaEtaDeltaPhiLeadingJet = new TH2F("fHistDeltaEtaDeltaPhiLeadingJet",
147                                              "#Delta#eta - #Delta#phi of leading jet;#Delta#eta;#Delta#phi",
148                                              100, -4., 4., 100, -2.*TMath::Pi(), 2*TMath::Pi());
149
150   fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
151                             300, 0., 300., 300, 0., 300.);
152
153   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
154     fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
155                                           300, 0., 300., 300, 0., 300.);
156     fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
157                                           200, -50., 150., 300, 0., 300.);
158   }
159
160   OpenFile(1);
161   fOutputList = new TList;
162   fOutputList->Add(fHistEvtSelection);
163   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
164     fOutputList->Add(fHistPtLeadingJet[iBranch]);
165     fOutputList->Add(fHistEtaPhiLeadingJet[iBranch]);
166   }
167   fOutputList->Add(fHistDeltaEtaDeltaPhiLeadingJet);
168   fOutputList->Add(fHistPtPtExtra);
169   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
170     fOutputList->Add(fHistPtResponse[iEvtClass]);
171     fOutputList->Add(fHistPtSmearing[iEvtClass]);
172   }
173 }
174
175 void AliAnalysisTaskJetResponse::UserExec(Option_t *)
176 {
177   // load events, apply event cuts, then compare leading jets
178
179   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
180       AliError("Jet branch name not set.");
181       return;
182    }
183
184   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
185   if (!fESD) {
186     AliError("ESD not available");
187     return;
188   }
189   fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
190   if (!fAOD) {
191     AliError("AOD not available");
192     return;
193   }
194
195   fHistEvtSelection->Fill(1); // number of events before event selection
196
197   // physics selection
198   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
199     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
200   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
201     Printf(" Trigger Selection: event REJECTED ... ");
202     fHistEvtSelection->Fill(2);
203     PostData(1, fOutputList);
204     return;
205   }
206
207   // vertex selection
208   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
209   Int_t nTracksPrim = primVtx->GetNContributors();
210   if ((nTracksPrim < fMinContribVtx) ||
211       (primVtx->GetZ() < fVtxZMin) ||
212       (primVtx->GetZ() > fVtxZMax) ){
213     Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
214     fHistEvtSelection->Fill(3);
215     PostData(1, fOutputList);
216     return;
217   }
218
219   // event class selection (from jet helper task)
220   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
221   Printf("Event class %d", eventClass);
222   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
223     fHistEvtSelection->Fill(4);
224     PostData(1, fOutputList);
225     return;
226   }
227
228   // centrality selection
229   AliCentrality *cent = fESD->GetCentrality();
230   Float_t centValue = cent->GetCentralityPercentile("TRK");
231   printf("centrality: %f\n", centValue);
232   if (centValue < fCentMin || centValue > fCentMax){
233     fHistEvtSelection->Fill(5);
234     PostData(1, fOutputList);
235     return;
236   }
237
238   fHistEvtSelection->Fill(0); // accepted events
239
240
241   TClonesArray *aodJets[2];
242   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data()));
243   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));
244   AliAODJet *leadingJet[2] = { 0x0, 0x0 };
245
246   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
247     fListJets[iJetType]->Clear();
248     if (!aodJets[iJetType])
249       continue;
250     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
251       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
252       if (jet)
253         fListJets[iJetType]->Add(jet);
254     }
255     fListJets[iJetType]->Sort();
256     leadingJet[iJetType] = (AliAODJet*) fListJets[iJetType]->First();
257
258     if (leadingJet[iJetType])
259       fHistEtaPhiLeadingJet[iJetType]->Fill(leadingJet[iJetType]->Eta(), leadingJet[iJetType]->Phi());
260   }
261
262   // check if leading jets are close in eta-phi and compare their Pt
263   if (leadingJet[0] && leadingJet[1]) {
264
265     // leading jets in "jet" acceptance
266     if(leadingJet[0]->Eta()>fJetEtaMax || leadingJet[0]->Eta()<fJetEtaMin ||
267        leadingJet[1]->Eta()>fJetEtaMax || leadingJet[1]->Eta()<fJetEtaMin){
268        Printf("Jet not in eta acceptance.");
269     }
270     else{
271        // check association of jets
272        Float_t deltaEta = leadingJet[0]->Eta() - leadingJet[1]->Eta();
273        Float_t deltaPhi = leadingJet[0]->Phi() - leadingJet[1]->Phi();
274        fHistDeltaEtaDeltaPhiLeadingJet->Fill(deltaEta, deltaPhi);
275
276        if (TMath::Abs(deltaEta) > fJetDeltaEta || (TMath::Pi() - TMath::Abs(TMath::Abs(deltaPhi) - TMath::Pi())) > fJetDeltaPhi)
277           printf("leading jets two far apart\n");
278        else {
279           fHistPtLeadingJet[0]->Fill(leadingJet[0]->Pt());
280           fHistPtLeadingJet[1]->Fill(leadingJet[1]->Pt());
281
282           fHistPtPtExtra->Fill(leadingJet[0]->Pt(), leadingJet[1]->Pt());
283
284           if (eventClass > -1 && eventClass < fkEvtClasses){
285               fHistPtResponse[eventClass]->Fill(leadingJet[1]->Pt(), leadingJet[0]->Pt());
286               fHistPtSmearing[eventClass]->Fill(leadingJet[1]->Pt()-leadingJet[0]->Pt(), leadingJet[0]->Pt());
287           }
288        }
289     }
290   }
291
292   PostData(1, fOutputList);
293 }
294
295 void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
296 {
297   // Draw result to the screen
298   // Called once at the end of the query
299
300   if (!GetOutputData(1))
301     return;
302 }