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