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