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