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