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