Include L0 trigger in the trigger maker and the trigger patch object Add monitoring...
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetResponseV2.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //
17 // task compares jets in two or three branches,
18 // written for analysis of jet embedding in HI events
19 //
20
21
22 #include "TChain.h"
23 #include "TTree.h"
24 #include "TMath.h"
25 #include "TH1F.h"
26 #include "TH1I.h"
27 #include "TH2F.h"
28 #include "TH3F.h"
29 #include "THnSparse.h"
30 #include "TCanvas.h"
31
32 #include "AliLog.h"
33
34 #include "AliAnalysisTask.h"
35 #include "AliAnalysisManager.h"
36
37 #include "AliVEvent.h"
38 #include "AliESDEvent.h"
39 #include "AliESDInputHandler.h"
40 #include "AliCentrality.h"
41 #include "AliAnalysisHelperJetTasks.h"
42 #include "AliInputEventHandler.h"
43 #include "AliAODJetEventBackground.h"
44 #include "AliAnalysisTaskFastEmbedding.h"
45
46 #include "AliAODEvent.h"
47 #include "AliAODJet.h"
48 #include "AliAODHandler.h"
49
50 #include "AliAnalysisTaskJetResponseV2.h"
51
52 using std::cout;
53 using std::endl;
54
55 ClassImp(AliAnalysisTaskJetResponseV2)
56
57 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
58 AliAnalysisTaskSE(),
59   fESD(0x0),
60   fAOD(0x0),
61   fAODOut(0x0),
62   fAODExtension(0x0),
63   fNonStdFile(""),
64   fBackgroundBranch(""),
65   fIsPbPb(kTRUE),
66   fOfflineTrgMask(AliVEvent::kAny),
67   fMinContribVtx(1),
68   fVtxZMin(-8.),
69   fVtxZMax(8.),
70   fEvtClassMin(0),
71   fEvtClassMax(4),
72   fCentMin(0.),
73   fCentMax(100.),
74   fNInputTracksMin(0),
75   fNInputTracksMax(-1),
76   fJetEtaMin(-.5),
77   fJetEtaMax(.5),
78   fJetPtMin(20.),
79   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
80   fJetPtFractionMin(0.5),
81   fNMatchJets(4),
82   fMatchMaxDist(0.8),
83   fKeepJets(kFALSE),
84   fkNbranches(2),
85   fkEvtClasses(12),
86   fOutputList(0x0),
87   fbEvent(kTRUE),
88   fbJetsMismatch1(kTRUE),
89   fbJetsMismatch2(kTRUE),
90   fbJetsRp(kTRUE),
91   fbJetsDeltaPt(kTRUE),
92   fbJetsEta(kTRUE),
93   fbJetsPhi(kTRUE),
94   fbJetsArea(kTRUE),
95   fbJets3Branches(kFALSE),
96   fbJetsBeforeCut1(kTRUE),
97   fbJetsBeforeCut2(kTRUE),
98   fHistEvtSelection(0x0),
99   fHistJetSelection(0x0),
100   fh2JetSelection(0x0),
101   fhnEvent(0x0),
102   fhnJetsMismatch1(0x0),
103   fhnJetsMismatch2(0x0),
104   fhnJetsRp(0x0),
105   fhnJetsDeltaPt(0x0),
106   fhnJetsEta(0x0),
107   fhnJetsPhi(0x0),
108   fhnJetsArea(0x0),
109   fhnJets3Branches(0x0),
110   fhnJetsBeforeCut1(0x0),
111   fhnJetsBeforeCut2(0x0)
112 {
113   // default Constructor
114
115   fJetBranchName[0] = "";
116   fJetBranchName[1] = "";
117   fJetBranchName[2] = "";
118
119   fListJets[0] = new TList;
120   fListJets[1] = new TList;
121   fListJets[2] = new TList;
122 }
123
124 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
125   AliAnalysisTaskSE(name),
126   fESD(0x0),
127   fAOD(0x0),
128   fAODOut(0x0),
129   fAODExtension(0x0),
130   fNonStdFile(""),
131   fBackgroundBranch(""),
132   fIsPbPb(kTRUE),
133   fOfflineTrgMask(AliVEvent::kAny),
134   fMinContribVtx(1),
135   fVtxZMin(-8.),
136   fVtxZMax(8.),
137   fEvtClassMin(0),
138   fEvtClassMax(4),
139   fCentMin(0.),
140   fCentMax(100.),
141   fNInputTracksMin(0),
142   fNInputTracksMax(-1),
143   fJetEtaMin(-.5),
144   fJetEtaMax(.5),
145   fJetPtMin(20.),
146   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
147   fJetPtFractionMin(0.5),
148   fNMatchJets(4),
149   fMatchMaxDist(0.8),
150   fKeepJets(kFALSE),
151   fkNbranches(2),
152   fkEvtClasses(12),
153   fOutputList(0x0),
154   fbEvent(kTRUE),
155   fbJetsMismatch1(kTRUE),
156   fbJetsMismatch2(kTRUE),
157   fbJetsRp(kTRUE),
158   fbJetsDeltaPt(kTRUE),
159   fbJetsEta(kTRUE),
160   fbJetsPhi(kTRUE),
161   fbJetsArea(kTRUE),
162   fbJets3Branches(kFALSE),
163   fbJetsBeforeCut1(kTRUE),
164   fbJetsBeforeCut2(kTRUE),
165   fHistEvtSelection(0x0),
166   fHistJetSelection(0x0),
167   fh2JetSelection(0x0),
168   fhnEvent(0x0),
169   fhnJetsMismatch1(0x0),
170   fhnJetsMismatch2(0x0),
171   fhnJetsRp(0x0),
172   fhnJetsDeltaPt(0x0),
173   fhnJetsEta(0x0),
174   fhnJetsPhi(0x0),
175   fhnJetsArea(0x0),
176   fhnJets3Branches(0x0),
177   fhnJetsBeforeCut1(0x0),
178   fhnJetsBeforeCut2(0x0)
179 {
180   // Constructor
181
182   fJetBranchName[0] = "";
183   fJetBranchName[1] = "";
184   fJetBranchName[2] = "";
185
186   fListJets[0] = new TList;
187   fListJets[1] = new TList;
188   fListJets[2] = new TList;
189
190   DefineOutput(1, TList::Class());
191 }
192
193 AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
194 {
195   delete fListJets[0];
196   delete fListJets[1];
197   delete fListJets[2];
198 }
199
200 void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3)
201 {
202   fJetBranchName[0] = branch1;
203   fJetBranchName[1] = branch2;
204   fJetBranchName[2] = branch3;
205
206   if(strlen(fJetBranchName[2].Data()) ) {
207     fkNbranches = 3;
208     fbJets3Branches = kTRUE;
209   }
210 }
211
212 void AliAnalysisTaskJetResponseV2::Init()
213 {
214
215   // check for jet branches
216   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
217     AliError("Jet branch name not set.");
218   }
219
220 }
221
222 void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
223 {
224   // Create histograms
225   // Called once
226   OpenFile(1);
227   if(!fOutputList) fOutputList = new TList;
228   fOutputList->SetOwner(kTRUE);
229
230   Bool_t oldStatus = TH1::AddDirectoryStatus();
231   TH1::AddDirectory(kFALSE);
232
233
234   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
235   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
236   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
237   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
238   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
239   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
240   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
241
242   fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
243   fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
244   fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
245   fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
246   fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
247   fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
248   fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
249   fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
250   fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
251
252   fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
253   fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
254   fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
255   fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
256   fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
257   fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
258   fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
259   fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
260   fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
261
262
263   UInt_t entries = 0; // bit coded, see GetDimParams() below
264   UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
265
266   if(fbEvent){
267     entries = 1<<0 | 1<<1 | 1<<2 | 1<<26;  // cent : nInpTrks : rp psi : pT hard bin
268     opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution
269     fhnEvent = NewTHnSparseF("fhnEvent", entries, opt);
270   }
271    
272   if(fbJetsMismatch1){ // real mismatch (no related rec jet found)
273     // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin
274     entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26;
275     opt =  1<<6 | 1<<8 | 1<<10;
276     fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt);
277   }
278
279   if(fbJetsMismatch2){  // acceptance + fraction cut
280     // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin
281     entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26;
282     opt = 1<<6 | 1<<7 | 1<<8 | 1<<9;
283     fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt);
284
285   }
286    
287    
288   if(fbJetsRp){
289     // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho
290     /*
291       entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<12 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<24 | 1<<25;
292       opt = 1<<4 | 1<<5;*/
293     // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin
294     entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26;
295     opt = 1<<4 | 1<<5;
296     fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt);
297   }
298
299   // cent : nInpTrks : rp bins:  deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt)
300   if(fbJetsDeltaPt){
301     entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26;
302     opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
303     fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt);
304   }
305
306   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta)
307   if(fbJetsEta){
308     entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26;
309     opt = 1<<15 | 1<<8 | 1<<9;
310     fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt);
311   }
312
313   // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin
314   if(fbJetsPhi){
315     entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26;
316     opt = 1<<10 | 1<<11;
317     fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt);
318   }
319
320   // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) : pT hard bin (hr for area)
321   if(fbJetsArea){
322     entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7 | 1<<26;
323     opt = 1<<18 | 1<<12 | 1<<13;
324     fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
325   }
326
327   // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : deltaR(1x) : pT hard bin
328   if(fbJets3Branches){
329     entries = 1<<0 | 1<<1 | 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28 | 1<<12 | 1<<13 | 1<<29 | 1<<19 | 1<<30 | 1<<17 | 1<<26;
330     opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
331     fhnJets3Branches = NewTHnSparseF("fhnJets3Branches", entries, opt);
332   }
333
334
335
336
337   //before cut
338
339   // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
340   if(fbJetsBeforeCut1){
341     entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
342     opt = 0;
343     fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt);
344   }
345
346   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
347   if(fbJetsBeforeCut2){
348     entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
349     opt = 0;
350     fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt);
351   }
352
353   fOutputList->Add(fHistEvtSelection);
354   fOutputList->Add(fHistJetSelection);
355   fOutputList->Add(fh2JetSelection);
356   if(fbEvent)           fOutputList->Add(fhnEvent);
357   if(fbJetsMismatch1)   fOutputList->Add(fhnJetsMismatch1);
358   if(fbJetsMismatch2)   fOutputList->Add(fhnJetsMismatch2);
359   if(fbJetsRp)          fOutputList->Add(fhnJetsRp);
360   if(fbJetsDeltaPt)     fOutputList->Add(fhnJetsDeltaPt);
361   if(fbJetsEta)         fOutputList->Add(fhnJetsEta);
362   if(fbJetsPhi)         fOutputList->Add(fhnJetsPhi);
363   if(fbJetsArea)        fOutputList->Add(fhnJetsArea);
364   if(fbJets3Branches)   fOutputList->Add(fhnJets3Branches);
365   if(fbJetsBeforeCut1)  fOutputList->Add(fhnJetsBeforeCut1);
366   if(fbJetsBeforeCut2)  fOutputList->Add(fhnJetsBeforeCut2);
367
368   // =========== Switch on Sumw2 for all histos ===========
369   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
370     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
371     if (h1){
372       h1->Sumw2();
373       continue;
374     }
375     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
376     if (hn){
377       hn->Sumw2();
378     }     
379   }
380   TH1::AddDirectory(oldStatus);
381
382   PostData(1, fOutputList);
383 }
384
385 void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
386 {
387   // load events, apply event cuts, then compare leading jets
388
389   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
390     AliError("Jet branch name not set.");
391     return;
392   }
393
394   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
395   if (!fESD) {
396     AliError("ESD not available");
397
398     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
399     //  assume that the AOD is in the general output...
400     fAODOut  = AODEvent();
401   } else {
402     fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
403   }
404   if (!fAOD) {
405     AliError("AOD not available");
406     return;
407   }
408
409   if(fNonStdFile.Length()!=0){
410     // case that we have an AOD extension we need can fetch the jets from the extended output
411     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
412     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
413     if(!fAODExtension){
414       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
415     }
416   }
417
418   // -- event selection --
419   fHistEvtSelection->Fill(1); // number of events before event selection
420
421   // physics selection
422   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
423     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
424   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
425     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
426     fHistEvtSelection->Fill(2);
427     PostData(1, fOutputList);
428     return;
429   }
430
431   // vertex selection
432   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
433   Int_t nTracksPrim = primVtx->GetNContributors();
434   if ((nTracksPrim < fMinContribVtx) ||
435       (primVtx->GetZ() < fVtxZMin) ||
436       (primVtx->GetZ() > fVtxZMax) ){
437     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
438     fHistEvtSelection->Fill(3);
439     PostData(1, fOutputList);
440     return;
441   }
442
443   // event class selection (from jet helper task)
444   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
445   if(fDebug) Printf("Event class %d", eventClass);
446   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
447     fHistEvtSelection->Fill(4);
448     PostData(1, fOutputList);
449     return;
450   }
451
452   // centrality selection
453   AliCentrality *cent = 0x0;
454   Float_t centValue = 0.; 
455   if(fESD) cent = fESD->GetCentrality();
456   if(cent) centValue = cent->GetCentralityPercentile("V0M");
457   if(fDebug) printf("centrality: %f\n", centValue);
458   if (centValue < fCentMin || centValue > fCentMax){
459     fHistEvtSelection->Fill(4);
460     PostData(1, fOutputList);
461     return;
462   }
463
464
465   // multiplicity due to input tracks
466   Int_t nInputTracks = GetNInputTracks();
467
468   if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
469     fHistEvtSelection->Fill(5);
470     PostData(1, fOutputList);
471     return;
472   }
473
474    
475   fHistEvtSelection->Fill(0); // accepted events  
476   // -- end event selection --
477
478   // pt hard
479   Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
480   Int_t pthardbin = GetPtHardBin(pthard);
481
482   // reaction plane
483   Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
484   if(fbEvent){
485     Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
486     fhnEvent->Fill(eventEntries);
487   }
488
489
490   // get background
491   AliAODJetEventBackground* externalBackground = 0;
492   if(!externalBackground&&fBackgroundBranch.Length()){
493     externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
494     if(!externalBackground && fAODOut) externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
495     if(!externalBackground && fAODExtension)  externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
496     //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
497   }
498   Float_t rho = 0;
499   if(externalBackground) rho = externalBackground->GetBackground(0);
500
501
502   // fetch jets
503   TClonesArray *aodJets[3];
504   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
505   if(!aodJets[0] && fAODOut) aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
506   if(!aodJets[0] && fAODExtension) aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
507   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1
508   if(!aodJets[1] && fAODOut) aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet
509   if(!aodJets[1] && fAODExtension) aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
510   if( strlen(fJetBranchName[2].Data()) ) {
511     aodJets[2] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2
512     if(!aodJets[2] && fAODOut) aodJets[2] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet
513     if(!aodJets[2] && fAODExtension) aodJets[2] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[2].Data()));
514     if(aodJets[2]) fkNbranches=3;
515     if(fDebug>10) printf("3rd branch: %s\n",fJetBranchName[2].Data());
516   }
517
518   if(fDebug>10)  printf("fkNbranches %d\n",fkNbranches);
519   for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
520     fListJets[iJetType]->Clear();
521     if (!aodJets[iJetType]) {
522         if(fDebug) Printf("%s: no jet branch\n",fJetBranchName[iJetType].Data());
523       continue;
524     }
525     if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
526       
527     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
528       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
529       if (jet) fListJets[iJetType]->Add(jet);
530     }
531     fListJets[iJetType]->Sort();
532   }
533    
534    
535   // jet matching 
536   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
537   static TArrayF aPtFraction(fListJets[0]->GetEntries());
538   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
539   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
540
541   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
542   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
543                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
544                                             aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
545   static TArrayI aMatchIndexv2(fListJets[0]->GetEntries());
546   static TArrayF aPtFractionv2(fListJets[0]->GetEntries());
547   if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
548   if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
549   if( strlen(fJetBranchName[2].Data()) ) {
550     AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()), 
551                                               fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()), 
552                                               aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
553   }
554    
555   // loop over matched jets
556   Int_t ir = -1;  // index of matched reconstruced jet
557   Int_t ir2 = -1; // index of matched reconstruced jet
558   Float_t fraction = -1.;
559   Float_t fraction2 = -1.;
560   AliAODJet *jet[3]  = { 0x0, 0x0, 0x0 };
561   Float_t jetEta[3]  = { -990., -990., -990. };
562   Float_t jetPhi[3]  = { -990., -990., -990. };
563   Float_t jetPt[3]   = { -990., -990., -990. };
564   Float_t jetArea[3] = { -990., -990., -990. };
565   Float_t rpJet[3]   = { -990., -990., -990. };
566   Int_t rpBin = -990;
567    
568   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
569     ir = aMatchIndex[ig];
570     if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
571     else ir2 =0;
572
573     //fetch jets
574     jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
575     if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
576     else      jet[1] = 0x0;
577     if(ir2>=0) jet[2] = (AliAODJet*)(fListJets[2]->At(ir2));
578     else      jet[2] = 0x0;
579
580     for(Int_t i=0; i<fkNbranches; ++i){
581       if(!jet[i]){
582         jetEta[i]  = -990;
583         jetPhi[i]  = -990.;
584         jetPt[i]   = -990.;
585         jetArea[i] = -990.;
586         rpJet[i]   = -990.;
587         if(i==1) rpBin = -990;
588       } 
589       else {
590         jetEta[i]  = jet[i]->Eta();
591         jetPhi[i]  = jet[i]->Phi();
592         jetPt[i]   = GetPt(jet[i], i);
593         jetArea[i] = jet[i]->EffectiveAreaCharged();
594         rpJet[i]   = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
595         if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3);
596       }
597     }
598     fraction  = aPtFraction[ig];
599     if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
600     else fraction2 = 0.;
601
602     // jet statistics
603     fHistJetSelection->Fill(1); // all probe jets
604     if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
605
606     if(ir<0){
607       fHistJetSelection->Fill(2);
608       if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
609          
610       if(fbJetsMismatch1){
611         if(!jet[0]) continue;
612         Double_t jetEntriesMismatch1[7] = {
613           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
614           (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
615         };               
616         fhnJetsMismatch1->Fill(jetEntriesMismatch1);
617       }
618          
619       continue;
620     }
621       
622     if(!jet[0] || !jet[1]){
623       fHistJetSelection->Fill(3);
624       if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
625       continue;
626     }
627       
628     // look for distance to next rec jet
629     Float_t distNextJet = -0.01; // no neighbor
630     Float_t ptNextJet = -1.;
631     for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
632       if(r==ir) continue;
633       Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
634       if(distNextJet<0. || distNextJet>tmpDeltaR){
635         distNextJet = tmpDeltaR;
636         if(fKeepJets) ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0);
637         else          ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->Pt();
638       }
639     }
640       
641
642
643     //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
644     //Float_t relRho = rho>0. ? localRho / rho : 0.;
645
646       
647     // calculate parameters of associated jets
648     /* from Leticia, kT clusterizer
649        Float_t par0[4] = { 0.00409,       0.01229,      0.05,         0.26 };
650        Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
651        Float_t par2[4] = { 1.02865e-01,   1.49039e-01,  1.53910e-01,  1.51109e-01 };
652     */
653     // own, from embedded tracks
654     Float_t par0[4] = { 0.02841,       0.05039,      0.09092,      0.24089     };
655     Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
656     Float_t par2[4] = { 4.95415e-02,   9.79538e-02,  1.32814e-01,  1.71743e-01 };
657       
658     Float_t rpCorr = 0.;         
659       
660     if(eventClass>0&&eventClass<4){
661       rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
662       rpCorr *= rho * jetArea[1];
663     }
664
665     Float_t delta      = jetPt[2]-jetPt[1];
666     Float_t deltaPt    = jetPt[1]-jetPt[0];
667     Float_t deltaEta   = jetEta[1]-jetEta[0];
668     Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
669     Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
670     Float_t deltaArea  = jetArea[1]-jetArea[0];
671       
672       
673     // fill thnsparse before acceptance cut
674     if(fbJetsBeforeCut1){
675       Double_t jetBeforeCutEntries1[10] = { 
676         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
677         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], 
678         (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
679       };
680       fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
681     }
682       
683     if(fbJetsBeforeCut2){
684       Double_t jetBeforeCutEntries2[10] = {
685         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin, 
686         (Double_t)jetPt[0], (Double_t)jetPt[1],
687         (Double_t)jetEta[0], (Double_t)jetEta[1], 
688         (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
689       };
690       fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
691     }
692       
693     // jet selection
694     Bool_t jetAccepted = kTRUE;
695     // minimum fraction required
696     if(fraction<fJetPtFractionMin){
697       fHistJetSelection->Fill(4);
698       fh2JetSelection->Fill(4.,jetPt[0]);
699       jetAccepted = kFALSE;
700     }
701       
702     if(jetAccepted){
703       // jet acceptance + minimum pT check
704       if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
705          jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
706             
707         if(fDebug){
708           Printf("Jet not in eta acceptance.");
709           Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
710           Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
711         }
712         fHistJetSelection->Fill(5);
713         fh2JetSelection->Fill(5.,jetPt[0]);
714         jetAccepted = kFALSE;
715       }
716     }
717     if(jetAccepted){
718       if(jetPt[0] < fJetPtMin){
719         if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[0]);
720         fHistJetSelection->Fill(6);
721         fh2JetSelection->Fill(6.,jetPt[0]);
722         jetAccepted = kFALSE;
723       }
724     }
725
726     if(jetAccepted){
727       if(jet[1]->Trigger()&fJetTriggerExcludeMask){
728         fHistJetSelection->Fill(7);
729         fh2JetSelection->Fill(7.,jetPt[0]);
730         jetAccepted = kFALSE;
731       }
732          
733     }
734       
735     if(!jetAccepted){
736       if(fbJetsMismatch2){
737         Double_t jetEntriesMismatch2[10] = {
738           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
739           (Double_t)jetPt[0], (Double_t)jetPt[1],
740           (Double_t)jetEta[0], (Double_t)jetEta[1],
741           (Double_t)deltaEta, (Double_t)deltaR,
742           (Double_t)fraction
743         };
744         fhnJetsMismatch2->Fill(jetEntriesMismatch2);     
745       }
746       continue;
747     }
748       
749     // all accepted jets
750     fHistJetSelection->Fill(0);
751     fh2JetSelection->Fill(0.,jetPt[0]);
752       
753     // fill thnsparse
754     if(fbJetsRp){
755       Double_t jetEntriesRp[11] = {
756         (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
757         (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1], 
758         (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr,
759         (Double_t)pthardbin
760       };
761       fhnJetsRp->Fill(jetEntriesRp);
762     }
763       
764     if(fbJetsDeltaPt){
765       Double_t jetEntriesDeltaPt[8] = {
766         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
767         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea,
768         (Double_t)pthardbin
769       };                 
770       fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
771     }
772       
773     if(fbJetsEta){
774       Double_t jetEntriesEta[11] = {
775         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
776         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1], 
777         (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR,
778         (Double_t)pthardbin
779       };                                 
780       fhnJetsEta->Fill(jetEntriesEta);
781     }
782       
783     if(fbJetsPhi){
784       Double_t jetEntriesPhi[10] = {
785         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
786         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1],
787         (Double_t)deltaPt, (Double_t)deltaPhi,
788         (Double_t)pthardbin
789       };
790       fhnJetsPhi->Fill(jetEntriesPhi);
791     }
792       
793     if(fbJetsArea){
794       Double_t jetEntriesArea[14] = {
795         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
796         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1], 
797         (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea, 
798         (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
799         (Double_t)pthardbin
800       };                                 
801       fhnJetsArea->Fill(jetEntriesArea);
802     }
803
804     if(fbJets3Branches){
805       Double_t jetEntries3Branches[14] = {
806         (Double_t)centValue, (Double_t)nInputTracks, 
807         (Double_t)jetPt[0], (Double_t)jetPt[1],
808         (Double_t)jetArea[0], (Double_t)jetArea[1],
809         (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin,
810         (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2, (Double_t)deltaR
811       };                                 
812       fhnJets3Branches->Fill(jetEntries3Branches);
813     }
814       
815   }
816
817   PostData(1, fOutputList);
818 }
819
820 void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
821 {
822   // Draw result to the screen
823   // Called once at the end of the query
824
825   if (!GetOutputData(1))
826     return;
827 }
828
829 Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
830 {
831   // number of tracks in the event, obtained via jet finder
832
833   Int_t nInputTracks = 0;
834
835   TString jbname(fJetBranchName[1]);
836   //needs complete event, use jets without background subtraction
837   for(Int_t i=1; i<=3; ++i){
838     if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
839   }
840   // use only HI event
841   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
842   else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
843   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
844   else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD");
845
846   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
847   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
848   if(!tmpAODjets && fAODOut) tmpAODjets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(jbname.Data()));
849   if(!tmpAODjets && fAODExtension) tmpAODjets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(jbname.Data()));
850   if(!tmpAODjets){
851     Printf("Jet branch %s not found", jbname.Data());
852     Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
853     return -1;
854   }
855    
856   for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
857     AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
858     if(!jet) continue;
859     TRefArray *trackList = jet->GetRefTracks();
860     Int_t nTracks = trackList->GetEntriesFast();
861     nInputTracks += nTracks;
862     if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
863   }
864   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
865
866   return nInputTracks;  
867 }
868
869 THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
870 {
871   // generate new THnSparseF, axes are defined in GetDimParams()
872
873   Int_t count = 0;
874   UInt_t tmp = entries;
875   while(tmp!=0){
876     count++;
877     tmp = tmp &~ -tmp;  // clear lowest bit
878   }
879
880   TString hnTitle(name);
881   const Int_t dim = count;
882   Int_t nbins[dim];
883   Double_t xmin[dim];
884   Double_t xmax[dim];
885
886   Int_t i=0;
887   Int_t c=0;
888   while(c<dim && i<32){
889     if(entries&(1<<i)){
890       Bool_t highres = opt&(1<<i);
891       TString label("");
892       GetDimParams(i, highres, label, nbins[c], xmin[c], xmax[c]);
893       hnTitle += Form(";%s",label.Data());
894       c++;
895     }
896       
897     i++;
898   }
899   hnTitle += ";";
900
901   return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
902 }
903
904 void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
905 {
906   // stores label and binning of axis for THnSparse
907
908   const Double_t pi = TMath::Pi();
909    
910   switch(iEntry){
911       
912   case 0:
913     label = "V0 centrality (%)";
914     if(hr){
915       nbins = 100;
916       xmin = 0.;
917       xmax = 100.;
918     } else {
919       nbins = 8;
920       xmin = 0.;
921       xmax = 80.;
922     }
923     break;
924       
925       
926   case 1:
927     label = "nb. of input tracks";
928     if(fIsPbPb){
929       if(hr){
930         nbins = 400;
931         xmin = 0.;
932         xmax = 4000.;
933       } else {
934         nbins = 40;
935         xmin = 0.;
936         xmax = 4000.;
937       }
938     } else {
939       nbins = 40;
940       xmin = 0.;
941       xmax = 400.;
942     }
943     break;
944       
945       
946   case 2:
947     label = "event plane #psi";
948     if(hr){
949       nbins = 30;
950       xmin = 0.;
951       xmax = pi;
952     } else {
953       nbins = 30;
954       xmin = 0.;
955       xmax = pi;
956     }
957     break;
958       
959       
960   case 3:
961     label = "event plane bin";
962     nbins = 3;
963     xmin = -.5;
964     xmax = 2.5;
965     break;
966       
967       
968   case 4:
969   case 5:
970     if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
971     if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
972     nbins = 48;
973     xmin = -pi;
974     xmax =  pi;
975     break;
976       
977       
978   case 6:
979   case 7:
980   case 27:
981     if(iEntry==6)label = "probe p_{T} (GeV/c)";
982     if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
983     if(hr){
984       nbins = 300;
985       xmin = -50.;
986       xmax = 250.;
987     } else {
988       nbins = 50;
989       xmin = 0.;
990       xmax = 250.;
991     }
992     break;
993       
994       
995   case 8:
996   case 9:
997     if(iEntry==8)label = "probe #eta";
998     if(iEntry==9)label = "rec #eta";
999     if(hr){
1000       nbins = 56;
1001       xmin = -.7;
1002       xmax =  .7;
1003     } else {
1004       nbins = 28;
1005       xmin = -.7;
1006       xmax =  .7;
1007     }
1008     break;
1009       
1010       
1011   case 10:
1012   case 11:
1013     if(iEntry==10)label = "probe #phi";
1014     if(iEntry==11)label = "rec #phi";
1015     if(hr){
1016       nbins = 90;  // modulo 18 (sectors)
1017       xmin = 0.;
1018       xmax = 2*pi;
1019     } else {
1020       nbins = 25;
1021       xmin = 0.;
1022       xmax = 2*pi;
1023     }
1024     break;
1025       
1026       
1027   case 12:
1028   case 13:
1029   case 29:
1030     if(iEntry==12)label = "probe area";
1031     if(iEntry==13 || iEntry==29)label = "rec area";
1032     if(hr){
1033       nbins = 100;
1034       xmin = 0.;
1035       xmax = 1.;
1036     } else {
1037       nbins = 25;
1038       xmin = 0.;
1039       xmax = 1.;
1040     }
1041     break;
1042       
1043   case 14:
1044   case 28:
1045     if(iEntry==14) label = "#delta p_{T}";
1046     if(iEntry==28) label = "#delta";
1047     if(hr){
1048       nbins = 241;
1049       xmin = -120.5;
1050       xmax =  120.5;
1051     } else {
1052       nbins = 101;
1053       xmin = -101.;
1054       xmax =  101.;
1055     }
1056     break;
1057       
1058   case 15:
1059     label = "#Delta#eta";
1060     if(hr){
1061       nbins = 51;
1062       xmin = -1.02;
1063       xmax =  1.02;
1064     } else {
1065       nbins = 51;
1066       xmin = -1.02;
1067       xmax =  1.02;
1068     }
1069     break;
1070       
1071       
1072   case 16:
1073     label = "#Delta#phi";
1074     if(hr){
1075       nbins = 45;
1076       xmin = -pi;
1077       xmax =  pi;
1078     } else {
1079       nbins = 45;
1080       xmin = -pi;
1081       xmax =  pi;
1082     }
1083     break;
1084       
1085       
1086   case 17:
1087     label = "#DeltaR";
1088     if(hr){
1089       nbins = 50;
1090       xmin = 0.;
1091       xmax = 1.;
1092     } else {
1093       nbins = 50;
1094       xmin = 0.;
1095       xmax = 1.;
1096     }
1097     break;
1098       
1099       
1100   case 18:
1101     label = "#Deltaarea";
1102     if(hr){
1103       nbins = 81;
1104       xmin = -.81;
1105       xmax =  .81;
1106     } else {
1107       nbins = 33;
1108       xmin = -.825;
1109       xmax =  .825;
1110     }
1111     break;
1112       
1113       
1114   case 19:
1115   case 30:
1116     label = "fraction";
1117     if(hr){
1118       nbins = 52;
1119       xmin = 0.;
1120       xmax = 1.04;
1121     } else {
1122       nbins = 52;
1123       xmin = 0.;
1124       xmax = 1.04;
1125     }
1126     break;
1127       
1128       
1129   case 20:
1130     label = "distance to closest rec jet";
1131     if(hr){
1132       nbins = 51;
1133       xmin = -0.02;
1134       xmax =  1.;
1135     } else {
1136       nbins = 51;
1137       xmin = -0.02;
1138       xmax = 1.;
1139     }
1140     break;
1141       
1142       
1143   case 21:
1144     label = "p_{T} of closest rec jet";
1145     nbins = 100;
1146     xmin =  0.;
1147     xmax =  200.;
1148     break;
1149       
1150   case 22:
1151     label = "#rho";
1152     nbins = 125;
1153     xmin  = 0.;
1154     xmax  = 250.;
1155     break;
1156       
1157   case 23:
1158     label = "abs. correction of #rho for RP";
1159     nbins =  51;
1160     xmin  = -51.;
1161     xmax  =  51.;
1162     break;
1163       
1164   case 24:
1165     label = "local #rho";
1166     nbins = 125;
1167     xmin  = 0.;
1168     xmax  = 250.;
1169     break;
1170       
1171   case 25:
1172     label = "local #rho / #rho";
1173     nbins = 500;
1174     xmin  = 0.;
1175     xmax  = 5.;
1176     break;
1177       
1178   case 26:
1179     label = "p_{T,hard} bin";
1180     nbins = 10;
1181     xmin  = -.5;
1182     xmax  = 9.5;
1183     break;
1184       
1185   }
1186
1187 }
1188
1189 //____________________________________________________________________________
1190 Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1191
1192   // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1193
1194   const Int_t nBins = 10;
1195   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1196    
1197   Int_t bin = -1;
1198   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1199     bin++;
1200   }
1201    
1202   return bin;
1203 }
1204
1205 //____________________________________________________________________________
1206 Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1207
1208   // returns jet pt, also negative pt after background subtraction if available
1209
1210   Double_t pt = 0.;
1211
1212   if(fKeepJets && mode>0){  // background subtracted pt, can be negative
1213     pt = j->GetPtSubtracted(0);
1214   }
1215   else{
1216     pt = j->Pt();  // if negative, pt=0.01
1217   }
1218
1219   return pt;
1220 }