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