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