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