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