e16e40db3bef801c763a8997e784aefe32e0dab5
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
1
2 // ******************************************
3 // This task computes several jet observables like 
4 // the fraction of energy in inner and outer coronnas,
5 // jet-track correlations,triggered jet shapes and 
6 // correlation strength distribution of particles inside jets.    
7 // Author: lcunquei@cern.ch
8 // *******************************************
9
10
11 /**************************************************************************
12  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
13  *                                                                        *
14  * Author: The ALICE Off-line Project.                                    *
15  * Contributors are mentioned in the code where appropriate.              *
16  *                                                                        *
17  * Permission to use, copy, modify and distribute this software and its   *
18  * documentation strictly for non-commercial purposes is hereby granted   *
19  * without fee, provided that the above copyright notice appears in all   *
20  * copies and that both the copyright notice and this permission notice   *
21  * appear in the supporting documentation. The authors make no claims     *
22  * about the suitability of this software for any purpose. It is          *
23  * provided "as is" without express or implied warranty.                  *
24  **************************************************************************/
25
26
27 #include "TChain.h"
28 #include "TTree.h"
29 #include "TMath.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TH3F.h"
33 #include "THnSparse.h"
34 #include "TCanvas.h"
35
36 #include "AliLog.h"
37
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
40
41 #include "AliVEvent.h"
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliCentrality.h"
45 #include "AliAnalysisHelperJetTasks.h"
46 #include "AliInputEventHandler.h"
47 #include "AliAODJetEventBackground.h"
48 #include "AliAnalysisTaskFastEmbedding.h"
49 #include "AliAODEvent.h"
50 #include "AliAODHandler.h"
51 #include "AliAODJet.h"
52
53 #include "AliAnalysisTaskJetCore.h"
54
55 ClassImp(AliAnalysisTaskJetCore)
56
57 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
58 AliAnalysisTaskSE(),
59 fESD(0x0),
60 fAOD(0x0),
61 fAODExtension(0x0),
62 fBackgroundBranch(""),
63 fNonStdFile(""),
64 fIsPbPb(kTRUE),
65 fOfflineTrgMask(AliVEvent::kAny),
66 fMinContribVtx(1),
67 fVtxZMin(-10.),
68 fVtxZMax(10.),
69 fEvtClassMin(0),
70 fEvtClassMax(4),
71 fFilterMask(0),
72 fRadioFrac(0.2),
73 fMinDist(0.1),
74 fCentMin(0.),
75 fCentMax(100.),
76 fNInputTracksMin(0),
77 fNInputTracksMax(-1),
78 fAngStructCloseTracks(0),
79 fCheckMethods(0),
80 fDoEventMixing(0), 
81 fJetEtaMin(-.5),
82 fJetEtaMax(.5),
83 fNevents(0x0),
84 fTindex(0x0),
85 fTrigBufferIndex(0x0),
86 fJetPtMin(20.),
87 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
88 fJetPtFractionMin(0.5),
89 fNMatchJets(4),
90 fMatchMaxDist(0.8),
91 fKeepJets(kFALSE),
92 fkNbranches(2),
93 fkEvtClasses(12),
94 fOutputList(0x0),
95 fbEvent(kTRUE),
96 fHistEvtSelection(0x0),
97 fhnDeltaR(0x0),
98 fhnMixedEvents(0x0),
99 fh2JetCoreMethod1C10(0x0),
100 fh2JetCoreMethod2C10(0x0),
101 fh2JetCoreMethod1C20(0x0),
102 fh2JetCoreMethod2C20(0x0),
103 fh2JetCoreMethod1C30(0x0),
104 fh2JetCoreMethod2C30(0x0),
105 fh2JetCoreMethod1C60(0x0),
106 fh2JetCoreMethod2C60(0x0),
107 fh2AngStructpt1C10(0x0),
108 fh2AngStructpt2C10(0x0),
109 fh2AngStructpt3C10(0x0),
110 fh2AngStructpt4C10(0x0),
111 fh2AngStructpt1C20(0x0),
112 fh2AngStructpt2C20(0x0),
113 fh2AngStructpt3C20(0x0),
114 fh2AngStructpt4C20(0x0),    
115 fh2AngStructpt1C30(0x0),
116 fh2AngStructpt2C30(0x0),
117 fh2AngStructpt3C30(0x0),
118 fh2AngStructpt4C30(0x0),   
119 fh2AngStructpt1C60(0x0),
120 fh2AngStructpt2C60(0x0),
121 fh2AngStructpt3C60(0x0),
122 fh2AngStructpt4C60(0x0),
123 fh2JetsumHT1R2(0x0),
124 fh2JetsumHT4R2(0x0),
125 fh2JetsumHT8R2(0x0),
126 fh2JetsumHT1R4(0x0),
127 fh2JetsumHT4R4(0x0),
128 fh2JetsumHT8R4(0x0),
129 fh2JetsumHT1R6(0x0),
130 fh2JetsumHT4R6(0x0),
131 fh2JetsumHT8R6(0x0),
132 fh2JetsumHT1R8(0x0),
133 fh2JetsumHT4R8(0x0),
134 fh2JetsumHT8R8(0x0),
135 fh2JetsumHT1R10(0x0),
136 fh2JetsumHT4R10(0x0),
137 fh2JetsumHT8R10(0x0),
138 fh3spectriggered(0x0),
139 fh3specbiased(0x0),
140 fh3spectot(0x0),
141 fh3spectotb(0x0)
142  
143 {
144    // default Constructor
145
146
147  // Trigger buffer.
148    for(Int_t i=0; i<10; i++) {
149                 for(Int_t j=0; j<7; j++) {
150                         fTrigBuffer[i][j]=0;
151                 }                               
152         }       
153
154
155
156
157
158    fJetBranchName[0] = "";
159    fJetBranchName[1] = "";
160
161    fListJets[0] = new TList;
162    fListJets[1] = new TList;
163 }
164
165 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
166 AliAnalysisTaskSE(name),
167 fESD(0x0),
168 fAOD(0x0),
169 fAODExtension(0x0),
170 fBackgroundBranch(""),
171 fNonStdFile(""),
172 fIsPbPb(kTRUE),
173 fOfflineTrgMask(AliVEvent::kAny),
174 fMinContribVtx(1),
175 fVtxZMin(-10.),
176 fVtxZMax(10.),
177 fEvtClassMin(0),
178 fEvtClassMax(4),
179 fFilterMask(0),
180 fRadioFrac(0.2),
181 fMinDist(0.1),
182 fCentMin(0.),
183 fCentMax(100.),
184 fNInputTracksMin(0),
185 fNInputTracksMax(-1),
186 fAngStructCloseTracks(0),
187 fCheckMethods(0),
188 fDoEventMixing(0),
189 fJetEtaMin(-.5),
190 fJetEtaMax(.5),
191 fNevents(0x0),
192 fTindex(0x0),
193 fTrigBufferIndex(0x0),
194 fJetPtMin(20.),
195 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
196 fJetPtFractionMin(0.5),
197 fNMatchJets(4),
198 fMatchMaxDist(0.8),
199 fKeepJets(kFALSE),
200 fkNbranches(2),
201 fkEvtClasses(12),
202 fOutputList(0x0),
203 fbEvent(kTRUE),
204 fHistEvtSelection(0x0),
205 fhnDeltaR(0x0),
206 fhnMixedEvents(0x0),
207 fh2JetCoreMethod1C10(0x0),
208 fh2JetCoreMethod2C10(0x0),
209 fh2JetCoreMethod1C20(0x0),
210 fh2JetCoreMethod2C20(0x0),
211 fh2JetCoreMethod1C30(0x0),
212 fh2JetCoreMethod2C30(0x0),
213 fh2JetCoreMethod1C60(0x0),
214 fh2JetCoreMethod2C60(0x0),
215 fh2AngStructpt1C10(0x0),
216 fh2AngStructpt2C10(0x0),
217 fh2AngStructpt3C10(0x0),
218 fh2AngStructpt4C10(0x0),
219 fh2AngStructpt1C20(0x0),
220 fh2AngStructpt2C20(0x0),
221 fh2AngStructpt3C20(0x0),
222 fh2AngStructpt4C20(0x0),    
223 fh2AngStructpt1C30(0x0),
224 fh2AngStructpt2C30(0x0),
225 fh2AngStructpt3C30(0x0),
226 fh2AngStructpt4C30(0x0),   
227 fh2AngStructpt1C60(0x0),
228 fh2AngStructpt2C60(0x0),
229 fh2AngStructpt3C60(0x0),
230 fh2AngStructpt4C60(0x0),    
231 fh2JetsumHT1R2(0x0),
232 fh2JetsumHT4R2(0x0),
233 fh2JetsumHT8R2(0x0),
234 fh2JetsumHT1R4(0x0),
235 fh2JetsumHT4R4(0x0),
236 fh2JetsumHT8R4(0x0),
237 fh2JetsumHT1R6(0x0),
238 fh2JetsumHT4R6(0x0),
239 fh2JetsumHT8R6(0x0),
240 fh2JetsumHT1R8(0x0),
241 fh2JetsumHT4R8(0x0),
242 fh2JetsumHT8R8(0x0),
243 fh2JetsumHT1R10(0x0),
244 fh2JetsumHT4R10(0x0),
245 fh2JetsumHT8R10(0x0),
246 fh3spectriggered(0x0),
247 fh3specbiased(0x0),
248 fh3spectot(0x0),
249 fh3spectotb(0x0)
250  {
251    // Constructor
252
253
254     for(Int_t i=0; i<10; i++) {
255        for(Int_t j=0; j<7; j++) {
256             fTrigBuffer[i][j]=0;
257                 }                               
258     }   
259
260
261
262    fJetBranchName[0] = "";
263    fJetBranchName[1] = "";
264
265    fListJets[0] = new TList;
266    fListJets[1] = new TList;
267
268    DefineOutput(1, TList::Class());
269 }
270
271 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
272 {
273    delete fListJets[0];
274    delete fListJets[1];
275 }
276
277 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
278 {
279    fJetBranchName[0] = branch1;
280    fJetBranchName[1] = branch2;
281 }
282
283 void AliAnalysisTaskJetCore::Init()
284 {
285
286    // check for jet branches
287    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
288       AliError("Jet branch name not set.");
289    }
290
291 }
292
293 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
294 {
295    // Create histograms
296    // Called once
297    OpenFile(1);
298    if(!fOutputList) fOutputList = new TList;
299    fOutputList->SetOwner(kTRUE);
300
301    Bool_t oldStatus = TH1::AddDirectoryStatus();
302    TH1::AddDirectory(kFALSE);
303
304
305    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
306    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
307    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
308    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
309    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
310    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
311    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
312
313      UInt_t entries = 0; // bit coded, see GetDimParams() below 
314      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7; 
315      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
316
317      //change binning in pTtrack
318      Double_t *xPt3 = new Double_t[10];
319      xPt3[0] = 0.;
320      for(int i = 1; i<=9;i++){
321       if(xPt3[i-1]<1)xPt3[i] = xPt3[i-1] + 0.2; // 1 - 5
322       else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
323       else xPt3[i] = xPt3[i-1] + 150.; // 18 
324      }
325     fhnDeltaR->SetBinEdges(2,xPt3);
326     delete [] xPt3;
327
328     //change binning in HTI
329      Double_t *xPt4 = new Double_t[14];
330      xPt4[0] = 0.;
331      for(int i = 1; i<=13;i++){
332       if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 1; // 1 - 10
333       else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 5; // 10 - 12
334       else xPt4[i] = xPt4[i-1] + 30.; // 13
335      }
336     fhnDeltaR->SetBinEdges(6,xPt4);
337     delete [] xPt4;
338
339     
340
341
342
343
344    
345      if(fDoEventMixing){    
346      UInt_t cifras = 0; // bit coded, see GetDimParams() below 
347      cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7; 
348      fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
349
350     if(fCheckMethods){
351
352
353     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
354     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
355     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
356     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
357     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
358     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
359     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
360     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
361
362    
363     if(fAngStructCloseTracks>0){
364     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
365     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
366     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
367     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
368     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
369     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
370     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
371     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
372     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
373     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
374     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
375     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
376     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
377     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
378     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
379     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
380
381     fh2JetsumHT1R2 = new TH2F("Pt sum R02 HT1","",20,0.,200.,100,0.,10.);
382     fh2JetsumHT4R2 = new TH2F("Pt sum R02 HT4","",20,0.,200.,100,0.,10.);
383     fh2JetsumHT8R2 = new TH2F("Pt sum R02 HT8","",20,0.,200.,100,0.,10.);
384     fh2JetsumHT1R4 = new TH2F("Pt sum R04 HT1","",20,0.,200.,100,0.,10.);
385     fh2JetsumHT4R4 = new TH2F("Pt sum R04 HT4","",20,0.,200.,100,0.,10.);
386     fh2JetsumHT8R4 = new TH2F("Pt sum R04 HT8","",20,0.,200.,100,0.,10.);
387     fh2JetsumHT1R6 = new TH2F("Pt sum R06 HT1","",20,0.,200.,100,0.,10.);
388     fh2JetsumHT4R6 = new TH2F("Pt sum R06 HT4","",20,0.,200.,100,0.,10.);
389     fh2JetsumHT8R6 = new TH2F("Pt sum R06 HT8","",20,0.,200.,100,0.,10.); 
390     fh2JetsumHT1R8 = new TH2F("Pt sum R08 HT1","",20,0.,200.,100,0.,10.);
391     fh2JetsumHT4R8 = new TH2F("Pt sum R08 HT4","",20,0.,200.,100,0.,10.);
392     fh2JetsumHT8R8 = new TH2F("Pt sum R08 HT8","",20,0.,200.,100,0.,10.); 
393     fh2JetsumHT1R10 = new TH2F("Pt sum R10 HT1","",20,0.,200.,100,0.,10.);
394     fh2JetsumHT4R10 = new TH2F("Pt sum R10 HT4","",20,0.,200.,100,0.,10.);
395     fh2JetsumHT8R10 = new TH2F("Pt sum R10 HT8","",20,0.,200.,100,0.,10.);  
396     fh3spectriggered = new TH3F("Triggered spectrum","",10,0,100,50,0.,200,50,0.,50.);
397     fh3specbiased = new TH3F("Biased spectrum","",10,0,100,50,0.,200.,50,0.,50.);
398     fh3spectot = new TH3F("Total spectrum 0-10","",50,0.,200.,50,0.,50.,50,0.,50.);
399     fh3spectotb = new TH3F("Total spectrum 0-20","",50,0.,200.,50,0.,50.,50,0.,50.);    
400    fOutputList->Add(fHistEvtSelection);
401
402    fOutputList->Add(fhnDeltaR);
403    
404    fOutputList->Add(fhnMixedEvents);
405          
406      
407   
408       if(fCheckMethods){
409
410       fOutputList->Add(fh2JetCoreMethod1C10);
411       fOutputList->Add(fh2JetCoreMethod2C10);
412       fOutputList->Add(fh2JetCoreMethod1C20);
413       fOutputList->Add(fh2JetCoreMethod2C20);
414       fOutputList->Add(fh2JetCoreMethod1C30);
415       fOutputList->Add(fh2JetCoreMethod2C30);
416       fOutputList->Add(fh2JetCoreMethod1C60);
417       fOutputList->Add(fh2JetCoreMethod2C60);}
418       
419       
420      
421
422
423         if(fAngStructCloseTracks>0){
424        fOutputList->Add(fh2AngStructpt1C10);
425        fOutputList->Add(fh2AngStructpt2C10);
426        fOutputList->Add(fh2AngStructpt3C10);
427        fOutputList->Add(fh2AngStructpt4C10); 
428        fOutputList->Add(fh2AngStructpt1C20);
429        fOutputList->Add(fh2AngStructpt2C20);
430        fOutputList->Add(fh2AngStructpt3C20);
431        fOutputList->Add(fh2AngStructpt4C20); 
432        fOutputList->Add(fh2AngStructpt1C30);
433        fOutputList->Add(fh2AngStructpt2C30);
434        fOutputList->Add(fh2AngStructpt3C30);
435        fOutputList->Add(fh2AngStructpt4C30);
436        fOutputList->Add(fh2AngStructpt1C60);
437        fOutputList->Add(fh2AngStructpt2C60);
438        fOutputList->Add(fh2AngStructpt3C60);
439        fOutputList->Add(fh2AngStructpt4C60);}  
440
441
442         fOutputList->Add(fh2JetsumHT1R2);
443         fOutputList->Add(fh2JetsumHT4R2);
444         fOutputList->Add(fh2JetsumHT8R2);
445
446         fOutputList->Add(fh2JetsumHT1R4);
447         fOutputList->Add(fh2JetsumHT4R4);
448         fOutputList->Add(fh2JetsumHT8R4);
449
450         fOutputList->Add(fh2JetsumHT1R6);
451         fOutputList->Add(fh2JetsumHT4R6);
452         fOutputList->Add(fh2JetsumHT8R6);
453
454         fOutputList->Add(fh2JetsumHT1R8);
455         fOutputList->Add(fh2JetsumHT4R8);
456         fOutputList->Add(fh2JetsumHT8R8);
457
458         fOutputList->Add(fh2JetsumHT1R10);
459         fOutputList->Add(fh2JetsumHT4R10);
460         fOutputList->Add(fh2JetsumHT8R10); 
461
462        fOutputList->Add(fh3spectriggered);
463        fOutputList->Add(fh3specbiased);
464        fOutputList->Add(fh3spectot);
465        fOutputList->Add(fh3spectotb); 
466    // =========== Switch on Sumw2 for all histos ===========
467    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
468       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
469       if (h1){
470          h1->Sumw2();
471          continue;
472       }
473     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
474       if (hn){
475          hn->Sumw2();
476       }   
477    }
478    TH1::AddDirectory(oldStatus);
479
480    PostData(1, fOutputList);
481 }
482
483 void AliAnalysisTaskJetCore::UserExec(Option_t *)
484 {
485    
486
487    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
488       AliError("Jet branch name not set.");
489       return;
490    }
491
492    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
493    if (!fESD) {
494       AliError("ESD not available");
495       fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
496    } else {
497       fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
498    }
499  
500     if(fNonStdFile.Length()!=0){
501     // case that we have an AOD extension we need can fetch the jets from the extended output
502     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
503     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
504     if(!fAODExtension){
505       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
506     }}
507     
508
509
510
511
512    // -- event selection --
513    fHistEvtSelection->Fill(1); // number of events before event selection
514
515    // physics selection
516    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
517    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
518    cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
519    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
520       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
521       fHistEvtSelection->Fill(2);
522       PostData(1, fOutputList);
523       return;
524    }
525
526    // vertex selection
527    if(!fAOD){
528      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
529      fHistEvtSelection->Fill(3);
530       PostData(1, fOutputList);
531    }
532    AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
533
534    if(!primVtx){
535      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
536      fHistEvtSelection->Fill(3);
537      PostData(1, fOutputList);
538      return;
539    }
540
541    Int_t nTracksPrim = primVtx->GetNContributors();
542    if ((nTracksPrim < fMinContribVtx) ||
543          (primVtx->GetZ() < fVtxZMin) ||
544          (primVtx->GetZ() > fVtxZMax) ){
545       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
546       fHistEvtSelection->Fill(3);
547       PostData(1, fOutputList);
548       return;
549    }
550
551    // event class selection (from jet helper task)
552    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
553    if(fDebug) Printf("Event class %d", eventClass);
554    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
555       fHistEvtSelection->Fill(4);
556       PostData(1, fOutputList);
557       return;
558    }
559
560    // centrality selection
561    AliCentrality *cent = 0x0;
562    Double_t centValue = 0.; 
563    if(fESD) {cent = fESD->GetCentrality();
564      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
565    else     centValue=fAOD->GetHeader()->GetCentrality();
566    
567    if(fDebug) printf("centrality: %f\n", centValue);
568    if (centValue < fCentMin || centValue > fCentMax){
569       fHistEvtSelection->Fill(4);
570       PostData(1, fOutputList);
571       return;
572    }
573
574
575    fHistEvtSelection->Fill(0); 
576    // accepted events  
577    // -- end event selection --
578   
579    // get background
580    AliAODJetEventBackground* externalBackground = 0;
581    if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
582       externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
583       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
584    }
585    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
586      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
587       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
588    }
589    
590    Float_t rho = 0;
591    if(externalBackground)rho = externalBackground->GetBackground(0);
592
593
594    // fetch jets
595    TClonesArray *aodJets[2];
596    aodJets[0]=0;
597    if(fAOD&&!aodJets[0]){
598    aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); 
599    aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));  }
600    if(fAODExtension && !aodJets[0]){ 
601    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); 
602    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
603
604    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
605    //Int_t inord[aodJets[0]->GetEntriesFast()];
606    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
607    //  ptsub[n]=0;
608    //  inord[n]=0;}   
609
610    TList ParticleList;
611    Int_t nT = GetListOfTracks(&ParticleList);
612      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
613       fListJets[iJetType]->Clear();
614       if (!aodJets[iJetType]) continue;
615
616       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
617       
618    
619       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
620          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
621          if (jet) fListJets[iJetType]->Add(jet);
622          // if(iJetType==0){
623          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
624       }}
625    
626    Double_t etabig=0;
627    Double_t ptbig=0;
628    Double_t areabig=0;
629    Double_t phibig=0.;
630    Double_t etasmall=0;
631    Double_t ptsmall=0;
632    Double_t areasmall=0;
633    //Double_t distr=0.;
634    Double_t phismall=0.;
635          
636   
637    // Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
638    // Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
639    // Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
640    // Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
641    // Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
642    // Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
643    // Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
644    // Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
645    Int_t iCount=0; 
646    Int_t trigJet=-1;
647    Int_t trigBBTrack=-1;
648    Int_t trigInTrack=-1;
649      
650    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
651            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
652            etabig  = jetbig->Eta();
653            phibig  = jetbig->Phi();
654            ptbig   = jetbig->Pt();
655            if(ptbig==0) continue; 
656            areabig = jetbig->EffectiveAreaCharged();
657            Double_t ptcorr=ptbig-rho*areabig;
658            if(ptcorr<=0) continue;
659            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
660                    Double_t dismin=100.;
661                    Double_t ptmax=-10.; 
662                    Int_t index1=-1;
663                    Int_t index2=-1;
664                   
665                Int_t point=GetHardestTrackBackToJet(jetbig);    
666                AliVParticle *partback = (AliVParticle*)ParticleList.At(point);                            
667                if(!partback) continue; 
668                    fh3spectriggered->Fill(centValue,ptcorr,partback->Pt());
669                    //if(partback->Pt()<6.) continue;
670                        AliAODTrack* leadtrack; 
671                        Int_t ippt=0;
672                        Double_t ppt=-10;   
673                        TRefArray *genTrackList = jetbig->GetRefTracks();
674                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
675                        AliAODTrack* genTrack;
676                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
677                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
678                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
679                        ippt=ir;}}
680                         leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
681                         if(!leadtrack) continue;
682                         fh3specbiased->Fill(centValue,ptcorr,leadtrack->Pt());
683                         if(centValue<10)fh3spectot->Fill(ptcorr,leadtrack->Pt(),partback->Pt());
684                         if(centValue<20)fh3spectotb->Fill(ptcorr,leadtrack->Pt(),partback->Pt()); 
685                         //store one trigger info                   
686                         if((partback->Pt()>10.)&&(iCount==0)){                        
687                           trigJet=i;
688                           trigBBTrack=point;
689                           trigInTrack=ippt;
690                           iCount=iCount+1;} 
691
692
693                   if(fCheckMethods){
694                   for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
695                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
696                   etasmall  = jetsmall->Eta();
697                   phismall = jetsmall->Phi();
698                   ptsmall   = jetsmall->Pt();
699                   areasmall = jetsmall->EffectiveAreaCharged();
700                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
701                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
702                      //Fraction in the jet core  
703                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
704                     index2=j;}  
705                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
706                       index1=j;}} //en of loop over R=0.2 jets
707                 //method1:most concentric jet=core 
708                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
709                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
710                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
711                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
712                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
713                 //method2:hardest contained jet=core   
714                 if(index2!=-1){ 
715                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
716                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
717                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
718                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
719                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
720                   Double_t sumpt2a=0.;
721                   Double_t sumpt2b=0.;
722                   Double_t sumpt2d=0.;
723                   Double_t sumpt4a=0.;
724                   Double_t sumpt4b=0.;
725                   Double_t sumpt4d=0.;
726                   Double_t sumpt6a=0.;
727                   Double_t sumpt6b=0.;
728                   Double_t sumpt6d=0.;
729                   Double_t sumpt8a=0.;
730                   Double_t sumpt8b=0.;
731                   Double_t sumpt8d=0.;  
732                   Double_t sumpt10a=0.;
733                   Double_t sumpt10b=0.;
734                   Double_t sumpt10d=0.;  
735
736           for(int it = 0;it<nT;++it){
737           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
738           Double_t deltaR = jetbig->DeltaR(part);
739           if(partback->Pt()>10.){
740             if(centValue<10){
741
742                              if(deltaR<0.2){if(leadtrack->Pt()>1.) sumpt2a=sumpt2a+part->Pt();
743                                             if(leadtrack->Pt()>4.) sumpt2b=sumpt2b+part->Pt(); 
744                                             if(leadtrack->Pt()>8.) sumpt2d=sumpt2d+part->Pt();}    
745                if(deltaR>=0.2 && deltaR<0.4){if(leadtrack->Pt()>1.) sumpt4a=sumpt4a+part->Pt();
746                                             if(leadtrack->Pt()>4.) sumpt4b=sumpt4b+part->Pt(); 
747                                             if(leadtrack->Pt()>8.) sumpt4d=sumpt4d+part->Pt();}                                                    if(deltaR>=0.4 && deltaR<0.6){if(leadtrack->Pt()>1.) sumpt6a=sumpt6a+part->Pt();
748                                             if(leadtrack->Pt()>4.) sumpt6b=sumpt6b+part->Pt(); 
749                                             if(leadtrack->Pt()>8.) sumpt6d=sumpt6d+part->Pt();}      
750                if(deltaR>=0.6 && deltaR<0.8){if(leadtrack->Pt()>1.) sumpt8a=sumpt8a+part->Pt();
751                                             if(leadtrack->Pt()>4.) sumpt8b=sumpt8b+part->Pt(); 
752                                             if(leadtrack->Pt()>8.) sumpt8d=sumpt8d+part->Pt();}      
753  
754                if(deltaR>=0.8 && deltaR<1.2){if(leadtrack->Pt()>1.)sumpt10a=sumpt10a+part->Pt();
755                                             if(leadtrack->Pt()>4.) sumpt10b=sumpt10b+part->Pt(); 
756                                             if(leadtrack->Pt()>8.) sumpt10d=sumpt10d+part->Pt();}   }}
757  
758
759           Double_t deltaEta = etabig-part->Eta();
760           Double_t deltaPhi=phibig-part->Phi();
761           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
762           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
763           Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,leadtrack->Pt(),partback->Pt()};                     fhnDeltaR->Fill(jetEntries);
764           }
765           //end of track loop
766             Double_t rhoin2=rho*TMath::Pi()*0.2*0.2; 
767             Double_t rhoin4=rho*TMath::Pi()*(0.4*0.4-0.2*0.2);
768             Double_t rhoin6=rho*TMath::Pi()*(0.6*0.6-0.4*0.4);
769             Double_t rhoin8=rho*TMath::Pi()*(0.8*0.8-0.6*0.6);
770             Double_t rhoin10=rho*TMath::Pi()*(1.2*1.2-0.8*0.8); 
771            
772               
773           if(rho!=0){
774          
775           fh2JetsumHT1R2->Fill(ptcorr,sumpt2a/rhoin2);
776           fh2JetsumHT4R2->Fill(ptcorr,sumpt2b/rhoin2);
777           fh2JetsumHT8R2->Fill(ptcorr,sumpt2d/rhoin2);
778          
779           fh2JetsumHT1R4->Fill(ptcorr,sumpt4a/rhoin4);
780           fh2JetsumHT4R4->Fill(ptcorr,sumpt4b/rhoin4);
781           fh2JetsumHT8R4->Fill(ptcorr,sumpt4d/rhoin4);
782          
783           fh2JetsumHT1R6->Fill(ptcorr,sumpt6a/rhoin6);
784           fh2JetsumHT4R6->Fill(ptcorr,sumpt6b/rhoin6);
785           fh2JetsumHT8R6->Fill(ptcorr,sumpt6d/rhoin6);
786          
787           fh2JetsumHT1R8->Fill(ptcorr,sumpt8a/rhoin8);
788           fh2JetsumHT4R8->Fill(ptcorr,sumpt8b/rhoin8);
789           fh2JetsumHT8R8->Fill(ptcorr,sumpt8d/rhoin8);
790          
791           fh2JetsumHT1R10->Fill(ptcorr,sumpt10a/rhoin10);
792           fh2JetsumHT4R10->Fill(ptcorr,sumpt10b/rhoin10);
793           fh2JetsumHT8R10->Fill(ptcorr,sumpt10d/rhoin10);}
794
795    }
796           //end of jet loop
797
798
799
800
801           if(fDoEventMixing){
802             //check before if the trigger exists
803             // fTrigBuffer[i][0] = zvtx
804             // fTrigBuffer[i][1] = phi
805             // fTrigBuffer[i][2] = eta
806             // fTrigBuffer[i][3] = pt_jet
807             // fTrigBuffer[i][4] = pt_trig
808             // fTrigBuffer[i][5]= pt_track_in
809             // fTrigBuffer[i][6]= centrality
810             if(fTindex==11) fTindex=0;
811             if(fTrigBuffer[fTindex][3]>0){
812             if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
813             if (TMath::Abs(fTrigBuffer[fTindex][6]-centValue<10)){  
814               
815                         for(int it = 0;it<nT;++it){
816                         AliVParticle *part = (AliVParticle*)ParticleList.At(it);         
817                         Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
818                         Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
819                         Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
820                         if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
821                         if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
822                         Double_t triggerEntries[8] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4],fTrigBuffer[fTindex][5]};                      
823                         fhnMixedEvents->Fill(triggerEntries);
824                         }
825                         fNevents=fNevents+1;  
826                         if(fNevents==9) {fTindex=fTindex+1;
827                         fNevents=0;} 
828             }}}
829         
830
831                // Copy the triggers from the current event into the buffer.
832                //again, only if the trigger exists:
833                 if(trigJet>-1){
834                 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));                 
835                 AliVParticle *partL = (AliVParticle*)ParticleList.At(trigInTrack);
836                 AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);         
837                 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
838                 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
839                 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
840                 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
841                 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
842                 fTrigBuffer[fTrigBufferIndex][5] = partL->Pt();
843                 fTrigBuffer[fTrigBufferIndex][6] = centValue;
844                 fTrigBufferIndex++;
845                 if(fTrigBufferIndex==9) fTrigBufferIndex=0;
846                 }
847           }
848           
849
850
851
852
853      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
854
855      //tracks up to R=0.8 distant from the jet axis
856    //   if(fAngStructCloseTracks==1){
857    //    TList CloseTrackList;
858    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
859    //    Double_t difR=0.04;
860    //    for(Int_t l=0;l<15;l++){
861    //    Double_t rr=l*0.1+0.1;
862    //     for(int it = 0;it<nn;++it){
863    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
864    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){      
865    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);  
866    //         Double_t ptm=part1->Pt();
867    //         Double_t ptn=part2->Pt(); 
868    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
869    //                    Rnm=TMath::Sqrt(Rnm);
870    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));      
871    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
872    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
873    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
874    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
875    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
876    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
877    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
878    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
879    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
880    //   }
881     
882    //   //only jet constituents
883    //    if(fAngStructCloseTracks==2){
884
885    //    Double_t difR=0.04;
886    //    for(Int_t l=0;l<15;l++){
887    //    Double_t rr=l*0.1+0.1;
888
889     
890    //    AliAODTrack* part1;
891    //    AliAODTrack* part2;
892
893    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
894    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
895           
896              
897
898    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
899    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
900    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
901    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
902    //         Double_t ptm=part1->Pt();
903    //         Double_t ptn=part2->Pt(); 
904    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
905    //                    Rnm=TMath::Sqrt(Rnm);
906    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
907    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR))); 
908    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
909    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
910    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
911    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
912    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
913    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
914    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
915    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
916    // }
917      // //end loop over R=0.4 jets      
918      // if(fAngStructCloseTracks>0){
919      // for(Int_t l=0;l<15;l++){
920      // Double_t rr=l*0.1+0.1;
921      //    if(down1[l]!=0){  
922      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
923      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
924      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
925      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
926      //    if(down2[l]!=0){  
927      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
928      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
929      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
930      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
931      //    if(down3[l]!=0){  
932      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
933      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
934      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
935      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
936      //    if(down4[l]!=0){  
937      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
938      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
939      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
940      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
941
942     
943
944
945
946
947
948    PostData(1, fOutputList);
949 }
950
951 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
952 {
953    // Draw result to the screen
954    // Called once at the end of the query
955
956    if (!GetOutputData(1))
957    return;
958 }
959
960
961
962
963
964
965
966
967
968
969
970 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
971
972     Int_t iCount = 0;
973  
974     
975     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
976       AliAODTrack *tr = fAOD->GetTrack(it);
977       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
978       if(TMath::Abs(tr->Eta())>0.9)continue;
979       if(tr->Pt()<0.15)continue;
980       list->Add(tr);
981       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
982       iCount++;
983     }
984   
985    
986   return iCount;
987  
988 }
989
990    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
991
992    
993     Int_t index=-1;
994     Double_t ptmax=-10;
995     Double_t dphi=0;
996     Double_t dif=0;
997     Int_t iCount=0;
998     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
999       AliAODTrack *tr = fAOD->GetTrack(it);
1000       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1001       if(TMath::Abs(tr->Eta())>0.9)continue;
1002       if(tr->Pt()<0.15)continue;
1003       iCount=iCount+1;
1004       dphi=RelativePhi(tr->Phi(),jetbig->Phi());  
1005       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
1006       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1007       index=iCount-1;
1008       dif=dphi;  }}
1009   
1010       return index;
1011
1012    }
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1023
1024     Int_t iCount = 0;
1025  
1026   
1027     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1028       AliAODTrack *tr = fAOD->GetTrack(it);
1029       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1030       if(TMath::Abs(tr->Eta())>0.9)continue;
1031       if(tr->Pt()<0.15)continue;
1032       Double_t disR=jetbig->DeltaR(tr);
1033       if(disR>0.8)  continue;
1034       list->Add(tr);
1035       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1036       iCount++;
1037     }
1038   
1039    list->Sort();
1040    return iCount;
1041
1042 }
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1055 {
1056
1057    Int_t nInputTracks = 0;
1058
1059    TString jbname(fJetBranchName[1]);
1060    //needs complete event, use jets without background subtraction
1061    for(Int_t i=1; i<=3; ++i){
1062       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1063    }
1064    // use only HI event
1065    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1066    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1067
1068    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1069    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
1070    if(!tmpAODjets){
1071       Printf("Jet branch %s not found", jbname.Data());
1072       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1073       return -1;
1074    }
1075    
1076    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1077       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1078       if(!jet) continue;
1079       TRefArray *trackList = jet->GetRefTracks();
1080       Int_t nTracks = trackList->GetEntriesFast();
1081       nInputTracks += nTracks;
1082       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1083    }
1084    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1085
1086    return nInputTracks;  
1087 }
1088
1089
1090
1091 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1092
1093   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1094   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1095   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1096   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1097   double dphi = mphi-vphi;
1098   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1099   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1100   return dphi;//dphi in [-Pi, Pi]
1101 }
1102
1103
1104
1105 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1106 {
1107    // generate new THnSparseF, axes are defined in GetDimParams()
1108
1109    Int_t count = 0;
1110    UInt_t tmp = entries;
1111    while(tmp!=0){
1112       count++;
1113       tmp = tmp &~ -tmp;  // clear lowest bit
1114    }
1115
1116    TString hnTitle(name);
1117    const Int_t dim = count;
1118    Int_t nbins[dim];
1119    Double_t xmin[dim];
1120    Double_t xmax[dim];
1121
1122    Int_t i=0;
1123    Int_t c=0;
1124    while(c<dim && i<32){
1125       if(entries&(1<<i)){
1126       
1127          TString label("");
1128          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1129          hnTitle += Form(";%s",label.Data());
1130          c++;
1131       }
1132       
1133       i++;
1134    }
1135    hnTitle += ";";
1136
1137    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1138 }
1139
1140 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1141 {
1142    // stores label and binning of axis for THnSparse
1143
1144    const Double_t pi = TMath::Pi();
1145    
1146    switch(iEntry){
1147       
1148    case 0:
1149       label = "V0 centrality (%)";
1150      
1151          nbins = 10;
1152          xmin = 0.;
1153          xmax = 100.;
1154          break;
1155       
1156       
1157    case 1:
1158       label = "corrected jet pt";
1159          nbins = 20;
1160          xmin = 0.;
1161          xmax = 200.;
1162           break;
1163       
1164       
1165    case 2:
1166       label = "track pT";
1167      
1168          nbins = 9;
1169          xmin = 0.;
1170          xmax = 150.;
1171          break;
1172       
1173       
1174     case 3:
1175       label = "deltaR";
1176       nbins = 15;
1177       xmin = 0.;
1178       xmax = 1.5;
1179       break;
1180
1181
1182
1183    case 4:
1184       label = "deltaEta";
1185       nbins = 8;
1186       xmin = -1.6;
1187       xmax = 1.6;
1188       break;
1189
1190
1191     case 5:
1192       label = "deltaPhi";
1193       nbins = 90;
1194       xmin = -0.5*pi;
1195       xmax = 1.5*pi;
1196       break;   
1197    
1198       
1199         
1200     case 6:
1201       label = "leading track";
1202       nbins = 13;
1203       xmin = 0;
1204       xmax = 50;
1205       break;
1206            
1207      case 7:
1208     
1209       label = "trigger track";
1210       nbins =5;
1211       xmin = 0;
1212       xmax = 50;
1213       break;
1214
1215
1216    
1217   
1218
1219
1220
1221
1222    }
1223
1224 }
1225