Adding ability to use multiple test functions
[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)aod = fAODIn;
960      else aod = fAODOut;   
961     
962     for(int it = 0;it < aod->GetNumberOfTracks();++it){
963       AliAODTrack *tr = aod->GetTrack(it);
964       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
965       if(TMath::Abs(tr->Eta())>0.9)continue;
966       if(tr->Pt()<0.15)continue;
967       list->Add(tr);
968       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
969       iCount++;
970     }
971   
972    
973   return iCount;
974  
975 }
976
977    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
978  
979     AliAODEvent *aod = 0;
980     if(!fESD)aod = fAODIn;
981     else aod = fAODOut;     
982     Int_t index=-1;
983     Double_t ptmax=-10;
984     Double_t dphi=0;
985     Double_t dif=0;
986     Int_t iCount=0;
987     for(int it = 0;it < aod->GetNumberOfTracks();++it){
988       AliAODTrack *tr = aod->GetTrack(it);
989       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
990       if(TMath::Abs(tr->Eta())>0.9)continue;
991       if(tr->Pt()<0.15)continue;
992       iCount=iCount+1;
993       dphi=RelativePhi(tr->Phi(),jetbig->Phi());  
994       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
995       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
996       index=iCount-1;
997       dif=dphi;  }}
998   
999       return index;
1000
1001    }
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1012
1013     Int_t iCount = 0;
1014       AliAODEvent *aod = 0;
1015      if(!fESD)aod = fAODIn;
1016      else aod = fAODOut;   
1017   
1018     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1019       AliAODTrack *tr = aod->GetTrack(it);
1020       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1021       if(TMath::Abs(tr->Eta())>0.9)continue;
1022       if(tr->Pt()<0.15)continue;
1023       Double_t disR=jetbig->DeltaR(tr);
1024       if(disR>0.8)  continue;
1025       list->Add(tr);
1026       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1027       iCount++;
1028     }
1029   
1030    list->Sort();
1031    return iCount;
1032
1033 }
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1046 {
1047
1048    Int_t nInputTracks = 0;
1049      AliAODEvent *aod = 0;
1050      if(!fESD)aod = fAODIn;
1051      else aod = fAODOut;   
1052    TString jbname(fJetBranchName[1]);
1053    //needs complete event, use jets without background subtraction
1054    for(Int_t i=1; i<=3; ++i){
1055       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1056    }
1057    // use only HI event
1058    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1059    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1060
1061    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1062    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1063    if(!tmpAODjets){
1064       Printf("Jet branch %s not found", jbname.Data());
1065       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1066       return -1;
1067    }
1068    
1069    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1070       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1071       if(!jet) continue;
1072       TRefArray *trackList = jet->GetRefTracks();
1073       Int_t nTracks = trackList->GetEntriesFast();
1074       nInputTracks += nTracks;
1075       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1076    }
1077    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1078
1079    return nInputTracks;  
1080 }
1081
1082
1083
1084 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1085
1086   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1087   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1088   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1089   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1090   double dphi = mphi-vphi;
1091   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1092   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1093   return dphi;//dphi in [-Pi, Pi]
1094 }
1095
1096
1097
1098 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1099 {
1100    // generate new THnSparseF, axes are defined in GetDimParams()
1101
1102    Int_t count = 0;
1103    UInt_t tmp = entries;
1104    while(tmp!=0){
1105       count++;
1106       tmp = tmp &~ -tmp;  // clear lowest bit
1107    }
1108
1109    TString hnTitle(name);
1110    const Int_t dim = count;
1111    Int_t nbins[dim];
1112    Double_t xmin[dim];
1113    Double_t xmax[dim];
1114
1115    Int_t i=0;
1116    Int_t c=0;
1117    while(c<dim && i<32){
1118       if(entries&(1<<i)){
1119       
1120          TString label("");
1121          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1122          hnTitle += Form(";%s",label.Data());
1123          c++;
1124       }
1125       
1126       i++;
1127    }
1128    hnTitle += ";";
1129
1130    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1131 }
1132
1133 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1134 {
1135    // stores label and binning of axis for THnSparse
1136
1137    const Double_t pi = TMath::Pi();
1138    
1139    switch(iEntry){
1140       
1141    case 0:
1142       label = "V0 centrality (%)";
1143      
1144          nbins = 10;
1145          xmin = 0.;
1146          xmax = 100.;
1147          break;
1148       
1149       
1150    case 1:
1151       label = "corrected jet pt";
1152          nbins = 20;
1153          xmin = 0.;
1154          xmax = 200.;
1155           break;
1156       
1157       
1158    case 2:
1159       label = "track pT";
1160      
1161          nbins = 10;
1162          xmin = 0.;
1163          xmax = 150;
1164          break;
1165       
1166       
1167     case 3:
1168       label = "deltaR";
1169       nbins = 15;
1170       xmin = 0.;
1171       xmax = 1.5;
1172       break;
1173
1174
1175
1176    case 4:
1177       label = "deltaEta";
1178       nbins = 15;
1179       xmin = -1.5;
1180       xmax = 1.5;
1181       break;
1182
1183
1184   case 5:
1185       label = "deltaPhi";
1186       nbins = 90;
1187       xmin = -0.5*pi;
1188       xmax = 1.5*pi;
1189       break;   
1190    
1191       
1192         
1193     case 6:
1194       label = "leading track";
1195       nbins = 14;
1196       xmin = 0;
1197       xmax = 50;
1198       break;
1199            
1200      case 7:
1201     
1202       label = "trigger track";
1203       nbins =5;
1204       xmin = 0;
1205       xmax = 50;
1206       break;
1207
1208
1209    
1210   
1211
1212
1213
1214
1215    }
1216
1217 }
1218