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