ae84c0f8dbb7a45f392595b5701a207426846e63
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCorePP.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 "TList.h"
30 #include "TMath.h"
31 #include "TH1F.h"
32 #include "TH1D.h"
33 #include "TH2F.h"
34 #include "TH3F.h"
35 #include "THnSparse.h"
36 #include "TCanvas.h"
37 #include "TArrayI.h" 
38 #include "TProfile.h"
39 #include "TFile.h"
40 #include "TKey.h"
41 #include "TRandom3.h"
42
43 #include "AliLog.h"
44
45 #include "AliAnalysisTask.h"
46 #include "AliAnalysisManager.h"
47
48 #include "AliVEvent.h"
49 #include "AliESDEvent.h"
50 #include "AliMCEvent.h"
51 #include "AliESDInputHandler.h"
52 #include "AliCentrality.h"
53 #include "AliAnalysisHelperJetTasks.h"
54 #include "AliInputEventHandler.h"
55 #include "AliAODJetEventBackground.h"
56 #include "AliGenPythiaEventHeader.h"
57 #include "AliAODMCParticle.h"
58 #include "AliMCParticle.h"
59 #include "AliAODEvent.h"
60 #include "AliAODHandler.h"
61 #include "AliAODJet.h"
62
63 #include "AliAnalysisTaskJetCorePP.h"
64
65 using std::cout;
66 using std::endl;
67
68 ClassImp(AliAnalysisTaskJetCorePP)
69
70 //Filip Krizek 1st March 2013
71
72 //---------------------------------------------------------------------
73 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
74 AliAnalysisTaskSE(),
75 fESD(0x0),
76 fAODIn(0x0),
77 fAODOut(0x0),
78 fAODExtension(0x0),
79 fJetBranchName(""),
80 fJetBranchNameChargMC(""),
81 fJetBranchNameFullMC(""),
82 fJetBranchNameBg(""),
83 fJetBranchNameBgChargMC(""),
84 fListJets(0x0),
85 fListJetsGen(0x0),
86 fListJetsGenFull(0x0),
87 fListJetsBg(0x0),
88 fListJetsBgGen(0x0),
89 fNonStdFile(""),
90 fSystem(0), //pp=0  pPb=1
91 fJetParamR(0.4),
92 fBgJetParamR(0.3),
93 fBgMaxJetPt(8.0),
94 fBgConeR(0.4),
95 fOfflineTrgMask(AliVEvent::kAny),
96 fMinContribVtx(1),
97 fVtxZMin(-10.0),
98 fVtxZMax(10.0),
99 fFilterMask(0),
100 fCentMin(0.0),
101 fCentMax(100.0),
102 fJetEtaMin(-0.5),
103 fJetEtaMax(0.5),
104 fTriggerEtaCut(0.9),
105 fTrackEtaCut(0.9),
106 fTrackLowPtCut(0.15),
107 fOutputList(0x0),
108 fHistEvtSelection(0x0),
109 fh2Ntriggers(0x0),
110 fHJetSpec(0x0),
111 fHJetSpecSubUeMedian(0x0),
112 fHJetSpecSubUeCone(0x0),
113 fHJetUeMedian(0x0),
114 fHJetUeCone(0x0),
115 fHRhoUeMedianVsCone(0x0),
116 //fHJetDensity(0x0),
117 //fHJetDensityA4(0x0),
118 fhJetPhi(0x0),
119 fhTriggerPhi(0x0),
120 fhJetEta(0x0),
121 fhTriggerEta(0x0),
122 fhVertexZ(0x0),
123 fhVertexZAccept(0x0),
124 fhContribVtx(0x0),
125 fhContribVtxAccept(0x0),
126 fhDphiTriggerJet(0x0),
127 fhDphiTriggerJetAccept(0x0),
128 fhCentrality(0x0),
129 fhCentralityAccept(0x0),
130 //fHJetPtRaw(0x0),
131 //fHLeadingJetPtRaw(0x0), 
132 //fHDphiVsJetPtAll(0x0), 
133 fhJetPtGenVsJetPtRec(0x0),
134 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
135 fhJetPtGenVsJetPtRecSubUeCone(0x0),
136 fhJetPtGen(0x0), 
137 fhJetPtSubUeMedianGen(0x0), 
138 fhJetPtSubUeConeGen(0x0), 
139 fhJetPtGenChargVsJetPtGenFull(0x0),
140 fhJetPtGenFull(0x0),
141 fh2NtriggersGen(0x0),
142 fHJetSpecGen(0x0),
143 fHJetSpecSubUeMedianGen(0x0),
144 fHJetSpecSubUeConeGen(0x0),
145 fHJetUeMedianGen(0x0),
146 fHJetUeConeGen(0x0),
147 fhPtTrkTruePrimRec(0x0),
148 fhPtTrkTruePrimGen(0x0),
149 fhPtTrkSecOrFakeRec(0x0),
150 fHRhoUeMedianVsConeGen(0x0),
151 fIsChargedMC(0),
152 fIsFullMC(0),
153 faGenIndex(0),
154 faRecIndex(0),
155 fkAcceptance(2.0*TMath::Pi()*1.8),
156 fkDeltaPhiCut(TMath::Pi()-0.8),
157 fh1Xsec(0x0),
158 fh1Trials(0x0),
159 fh1AvgTrials(0x0),
160 fh1PtHard(0x0),
161 fh1PtHardNoW(0x0),  
162 fh1PtHardTrials(0x0),
163 fAvgTrials(1),
164 fHardest(0),
165 fEventNumberRangeLow(0),
166 fEventNumberRangeHigh(99),
167 fTriggerPtRangeLow(0.0),
168 fTriggerPtRangeHigh(10000.0),
169 fFillRespMx(0),
170 fRandom(0x0),
171 fnTrials(1000),
172 fnPhi(9),
173 fnEta(2),
174 fEtaSize(0.7),
175 fPhiSize(2*TMath::Pi()/fnPhi),
176 fCellArea(fPhiSize*fEtaSize),
177 fSafetyMargin(1.1)
178 {
179    // default Constructor
180 }
181
182 //---------------------------------------------------------------------
183
184 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
185 AliAnalysisTaskSE(name),
186 fESD(0x0),
187 fAODIn(0x0),
188 fAODOut(0x0),
189 fAODExtension(0x0),
190 fJetBranchName(""),
191 fJetBranchNameChargMC(""),
192 fJetBranchNameFullMC(""),
193 fJetBranchNameBg(""),
194 fJetBranchNameBgChargMC(""),
195 fListJets(0x0),
196 fListJetsGen(0x0),
197 fListJetsGenFull(0x0),
198 fListJetsBg(0x0),
199 fListJetsBgGen(0x0),
200 fNonStdFile(""),
201 fSystem(0),  //pp=0   pPb=1
202 fJetParamR(0.4),
203 fBgJetParamR(0.3),
204 fBgMaxJetPt(8.0),
205 fBgConeR(0.4),
206 fOfflineTrgMask(AliVEvent::kAny),
207 fMinContribVtx(1),
208 fVtxZMin(-10.0),
209 fVtxZMax(10.0),
210 fFilterMask(0),
211 fCentMin(0.0),
212 fCentMax(100.0),
213 fJetEtaMin(-0.5),
214 fJetEtaMax(0.5),
215 fTriggerEtaCut(0.9),
216 fTrackEtaCut(0.9),
217 fTrackLowPtCut(0.15),
218 fOutputList(0x0),
219 fHistEvtSelection(0x0),
220 fh2Ntriggers(0x0),
221 fHJetSpec(0x0),
222 fHJetSpecSubUeMedian(0x0),
223 fHJetSpecSubUeCone(0x0),
224 fHJetUeMedian(0x0),
225 fHJetUeCone(0x0),
226 fHRhoUeMedianVsCone(0x0),
227 //fHJetDensity(0x0),
228 //fHJetDensityA4(0x0),
229 fhJetPhi(0x0),
230 fhTriggerPhi(0x0),
231 fhJetEta(0x0),
232 fhTriggerEta(0x0),
233 fhVertexZ(0x0),
234 fhVertexZAccept(0x0),
235 fhContribVtx(0x0),
236 fhContribVtxAccept(0x0),
237 fhDphiTriggerJet(0x0),
238 fhDphiTriggerJetAccept(0x0),
239 fhCentrality(0x0),
240 fhCentralityAccept(0x0),
241 //fHJetPtRaw(0x0),
242 //fHLeadingJetPtRaw(0x0), 
243 //fHDphiVsJetPtAll(0x0), 
244 fhJetPtGenVsJetPtRec(0x0),
245 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
246 fhJetPtGenVsJetPtRecSubUeCone(0x0),
247 fhJetPtGen(0x0),
248 fhJetPtSubUeMedianGen(0x0), 
249 fhJetPtSubUeConeGen(0x0), 
250 fhJetPtGenChargVsJetPtGenFull(0x0),
251 fhJetPtGenFull(0x0),
252 fh2NtriggersGen(0x0),
253 fHJetSpecGen(0x0),
254 fHJetSpecSubUeMedianGen(0x0),
255 fHJetSpecSubUeConeGen(0x0),
256 fHJetUeMedianGen(0x0),
257 fHJetUeConeGen(0x0),
258 fhPtTrkTruePrimRec(0x0),
259 fhPtTrkTruePrimGen(0x0),
260 fhPtTrkSecOrFakeRec(0x0),
261 fHRhoUeMedianVsConeGen(0x0),
262 fIsChargedMC(0),
263 fIsFullMC(0),
264 faGenIndex(0),
265 faRecIndex(0),
266 fkAcceptance(2.0*TMath::Pi()*1.8),
267 fkDeltaPhiCut(TMath::Pi()-0.8),
268 fh1Xsec(0x0),
269 fh1Trials(0x0), 
270 fh1AvgTrials(0x0),
271 fh1PtHard(0x0),
272 fh1PtHardNoW(0x0),  
273 fh1PtHardTrials(0x0),
274 fAvgTrials(1),
275 fHardest(0),
276 fEventNumberRangeLow(0),
277 fEventNumberRangeHigh(99),
278 fTriggerPtRangeLow(0.0),
279 fTriggerPtRangeHigh(10000.0),
280 fFillRespMx(0),
281 fRandom(0x0),
282 fnTrials(1000),
283 fnPhi(9),
284 fnEta(2),
285 fEtaSize(0.7),
286 fPhiSize(2*TMath::Pi()/fnPhi),
287 fCellArea(fPhiSize*fEtaSize),
288 fSafetyMargin(1.1)
289 {
290 // Constructor
291
292    DefineOutput(1, TList::Class());
293 }
294
295 //--------------------------------------------------------------
296 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
297 AliAnalysisTaskSE(a.GetName()),
298 fESD(a.fESD),
299 fAODIn(a.fAODIn),
300 fAODOut(a.fAODOut),
301 fAODExtension(a.fAODExtension),
302 fJetBranchName(a.fJetBranchName),
303 fJetBranchNameChargMC(a.fJetBranchNameChargMC),
304 fJetBranchNameFullMC(a.fJetBranchNameFullMC),
305 fJetBranchNameBg(a.fJetBranchNameBg),
306 fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
307 fListJets(a.fListJets),
308 fListJetsGen(a.fListJetsGen),
309 fListJetsGenFull(a.fListJetsGenFull),
310 fListJetsBg(a.fListJetsBg),
311 fListJetsBgGen(a.fListJetsBgGen),
312 fNonStdFile(a.fNonStdFile),
313 fSystem(a.fSystem),  
314 fJetParamR(a.fJetParamR),
315 fBgJetParamR(a.fBgJetParamR),
316 fBgMaxJetPt(a.fBgMaxJetPt),
317 fBgConeR(a.fBgConeR),
318 fOfflineTrgMask(a.fOfflineTrgMask),
319 fMinContribVtx(a.fMinContribVtx),
320 fVtxZMin(a.fVtxZMin),
321 fVtxZMax(a.fVtxZMax),
322 fFilterMask(a.fFilterMask),
323 fCentMin(a.fCentMin),
324 fCentMax(a.fCentMax),
325 fJetEtaMin(a.fJetEtaMin),
326 fJetEtaMax(a.fJetEtaMax),
327 fTriggerEtaCut(a.fTriggerEtaCut),
328 fTrackEtaCut(a.fTrackEtaCut),
329 fTrackLowPtCut(a.fTrackLowPtCut),
330 fOutputList(a.fOutputList),
331 fHistEvtSelection(a.fHistEvtSelection),
332 fh2Ntriggers(a.fh2Ntriggers),
333 fHJetSpec(a.fHJetSpec),
334 fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
335 fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
336 fHJetUeMedian(a.fHJetUeMedian),
337 fHJetUeCone(a.fHJetUeCone),
338 fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone), 
339 //fHJetDensity(a.fHJetDensity),
340 //fHJetDensityA4(a.fHJetDensityA4),
341 fhJetPhi(a.fhJetPhi),
342 fhTriggerPhi(a.fhTriggerPhi),
343 fhJetEta(a.fhJetEta),
344 fhTriggerEta(a.fhTriggerEta),
345 fhVertexZ(a.fhVertexZ),
346 fhVertexZAccept(a.fhVertexZAccept),
347 fhContribVtx(a.fhContribVtx),
348 fhContribVtxAccept(a.fhContribVtxAccept),
349 fhDphiTriggerJet(a.fhDphiTriggerJet),
350 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
351 fhCentrality(a.fhCentrality),
352 fhCentralityAccept(a.fhCentralityAccept),
353 //fHJetPtRaw(a.fHJetPtRaw),
354 //fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
355 //fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
356 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
357 fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
358 fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
359 fhJetPtGen(a.fhJetPtGen),
360 fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen), 
361 fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen), 
362 fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
363 fhJetPtGenFull(a.fhJetPtGenFull),
364 fh2NtriggersGen(a.fh2NtriggersGen),
365 fHJetSpecGen(a.fHJetSpecGen),
366 fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
367 fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
368 fHJetUeMedianGen(a.fHJetUeMedianGen),
369 fHJetUeConeGen(a.fHJetUeConeGen),
370 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
371 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
372 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
373 fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
374 fIsChargedMC(a.fIsChargedMC),
375 fIsFullMC(a.fIsFullMC),
376 faGenIndex(a.faGenIndex),
377 faRecIndex(a.faRecIndex),
378 fkAcceptance(a.fkAcceptance),
379 fkDeltaPhiCut(a.fkDeltaPhiCut),
380 fh1Xsec(a.fh1Xsec),
381 fh1Trials(a.fh1Trials),
382 fh1AvgTrials(a.fh1AvgTrials),
383 fh1PtHard(a.fh1PtHard),
384 fh1PtHardNoW(a.fh1PtHardNoW),  
385 fh1PtHardTrials(a.fh1PtHardTrials),
386 fAvgTrials(a.fAvgTrials),
387 fHardest(a.fHardest),
388 fEventNumberRangeLow(a.fEventNumberRangeLow),
389 fEventNumberRangeHigh(a.fEventNumberRangeHigh),
390 fTriggerPtRangeLow(a.fTriggerPtRangeLow),
391 fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
392 fFillRespMx(a.fFillRespMx),
393 fRandom(a.fRandom),
394 fnTrials(a.fnTrials),
395 fnPhi(a.fnPhi),
396 fnEta(a.fnEta),
397 fEtaSize(a.fEtaSize),
398 fPhiSize(a.fPhiSize),
399 fCellArea(a.fCellArea),
400 fSafetyMargin(a.fSafetyMargin)
401 {
402    //Copy Constructor
403    fRandom->SetSeed(0);
404 }
405 //--------------------------------------------------------------
406
407 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
408   // assignment operator
409   this->~AliAnalysisTaskJetCorePP();
410   new(this) AliAnalysisTaskJetCorePP(a);
411   return *this;
412 }
413 //--------------------------------------------------------------
414
415 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
416 {
417    //Destructor 
418    delete fListJets;
419    delete fListJetsGen;
420    delete fListJetsGenFull;
421    delete fListJetsBg;
422    delete fListJetsBgGen;
423    delete fOutputList; // ????
424    delete fRandom;
425 }
426
427 //--------------------------------------------------------------
428
429
430 Bool_t AliAnalysisTaskJetCorePP::Notify()
431 {
432    //Implemented Notify() to read the cross sections
433    //and number of trials from pyxsec.root
434    //inspired by AliAnalysisTaskJetSpectrum2::Notify()
435    if(!fIsChargedMC) return kFALSE; 
436
437    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
438    if(!fESD){
439       if(fDebug>1) AliError("ESD not available");
440       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
441    } 
442  
443    fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
444
445
446    if(fNonStdFile.Length()!=0){
447       // case that we have an AOD extension we can fetch the jets from the extended output
448       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
449       fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
450       if(!fAODExtension){
451          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
452       } 
453    }
454  
455    TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
456    Float_t xsection = 0;
457    Float_t ftrials  = 1;
458
459    fAvgTrials = 1;
460    if(tree){
461       TFile *curfile = tree->GetCurrentFile();
462       if(!curfile) {
463          Error("Notify","No current file");
464          return kFALSE;
465       }
466       if(!fh1Xsec || !fh1Trials){
467          Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
468          return kFALSE;
469       }
470       AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
471       fh1Xsec->Fill("<#sigma>",xsection);
472       // construct a poor man average trials
473       Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
474       if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
475       fh1Trials->Fill("#sum{ntrials}",ftrials);
476    }  
477
478    if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
479
480    return kTRUE;
481 }
482 //--------------------------------------------------------------
483
484 void AliAnalysisTaskJetCorePP::Init()
485 {
486    // check for jet branches
487    if(!strlen(fJetBranchName.Data())){
488       AliError("Jet branch name not set.");
489    }
490 }
491
492 //--------------------------------------------------------------
493
494 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
495 {
496    // Create histograms and initilize variables
497    fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
498   
499    // Called once
500    fListJets   = new TList();  //reconstructed level antikt jets
501    fListJetsBg = new TList();  //reconstructed jets to be removed from bg
502
503    fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
504    fIsFullMC    = (fJetBranchNameFullMC.Length()>0)  ? kTRUE : kFALSE;
505
506    fRandom = new TRandom3(0);
507
508    if(fIsChargedMC){
509       fListJetsGen   = new TList(); //generator level charged antikt jets
510       fListJetsBgGen = new TList(); //generator level jets to be removed from bg 
511
512       if(fIsFullMC)
513          fListJetsGenFull = new TList(); //generator level full jets
514    }
515    OpenFile(1);
516    if(!fOutputList) fOutputList = new TList();
517    fOutputList->SetOwner(kTRUE);
518
519    Bool_t oldStatus = TH1::AddDirectoryStatus();
520    TH1::AddDirectory(kFALSE);
521
522    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
523    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
524    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
525    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
526    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
527    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
528    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
529    
530    fOutputList->Add(fHistEvtSelection);
531
532    Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
533    //trigger pt spectrum (reconstructed) 
534    fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
535                              nBinsCentrality,0.0,100.0,50,0.0,50.0);
536    fOutputList->Add(fh2Ntriggers);
537
538    //Centrality, A, pTjet, pTtrigg, dphi
539    const Int_t dimSpec   = 5;
540    const Int_t nBinsSpec[dimSpec]     = {nBinsCentrality, 100,   220,  50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
541    const Double_t lowBinSpec[dimSpec] = {0.0,             0.0,   -20, 0.0, fkDeltaPhiCut};
542    const Double_t hiBinSpec[dimSpec]  = {100.0,           1.0, 200.0,50.0, TMath::Pi()};
543    fHJetSpec = new THnSparseF("fHJetSpec",
544                    "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
545                    dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
546    fOutputList->Add(fHJetSpec);  
547
548    //background estimated as  median of kT jets 
549    fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
550    fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
551    fOutputList->Add(fHJetSpecSubUeMedian); 
552    //background estimated as weighted  median of kT jets  ala Cone
553    fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
554    fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
555    fOutputList->Add(fHJetSpecSubUeCone); 
556
557
558
559    //------------------- HISTOS FOR DIAGNOSTIC ----------------------
560    //A, pTjet, pTjet-pTUe, pTUe, rhoUe     bg estimate from kT median
561    const Int_t    dimSpecMed   = 5;
562    const Int_t    nBinsSpecMed[dimSpecMed]  = {50,     50,  120,     50,   50};
563    const Double_t lowBinSpecMed[dimSpecMed] = {0.0,   0.0, -20.0,   0.0,  0.0};
564    const Double_t hiBinSpecMed[dimSpecMed]  = {1.0, 100.0, 100.0,  20.0, 20.0};
565    fHJetUeMedian = new THnSparseF("fHJetUeMedian",
566                    "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
567                    dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
568    fOutputList->Add(fHJetUeMedian);  
569    
570    //A, pTjet, pTjet-pTUe, pTUe, rhoUe     bg estimate from kT median Cone
571    fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
572    fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
573    fOutputList->Add(fHJetUeCone); 
574
575    //rho bacground reconstructed data
576    const Int_t    dimRho   = 2;
577    const Int_t    nBinsRho[dimRho]  = {200  , 200};
578    const Double_t lowBinRho[dimRho] = {0.0  , 0.0};
579    const Double_t hiBinRho[dimRho]  = {20.0 , 20.0};
580
581    fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",  
582                                       dimRho, nBinsRho, lowBinRho, hiBinRho);
583    fOutputList->Add(fHRhoUeMedianVsCone);
584
585    //Jet number density histos [Trk Mult, jet density, pT trigger]
586    /*const Int_t    dimJetDens   = 3;
587    const Int_t    nBinsJetDens[dimJetDens]  = {100,   100,   10};
588    const Double_t lowBinJetDens[dimJetDens] = {0.0,   0.0,  0.0};
589    const Double_t hiBinJetDens[dimJetDens]  = {500.0, 5.0, 50.0};
590
591    fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
592                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
593
594    fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
595                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
596
597    fOutputList->Add(fHJetDensity);
598    fOutputList->Add(fHJetDensityA4);
599    */      
600
601    //inclusive azimuthal and pseudorapidity histograms
602    fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
603                         50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
604    fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
605                         50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
606    fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
607                         50,0, 100, 40,-0.9,0.9);
608    fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
609                         50, 0, 50, 40,-0.9,0.9);
610
611    fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);  
612    fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);  
613    fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);   
614    fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);   
615    fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi()); 
616    fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi()); 
617    fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
618    fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
619
620    fOutputList->Add(fhJetPhi);
621    fOutputList->Add(fhTriggerPhi);
622    fOutputList->Add(fhJetEta);
623    fOutputList->Add(fhTriggerEta);
624    fOutputList->Add(fhVertexZ);    
625    fOutputList->Add(fhVertexZAccept);    
626    fOutputList->Add(fhContribVtx); 
627    fOutputList->Add(fhContribVtxAccept); 
628    fOutputList->Add(fhDphiTriggerJet);
629    fOutputList->Add(fhDphiTriggerJetAccept);
630    fOutputList->Add(fhCentrality); 
631    fOutputList->Add(fhCentralityAccept);
632
633    // raw spectra of INCLUSIVE jets  
634    //Centrality, pTjet, A
635    /*const Int_t dimRaw   = 3;
636    const Int_t nBinsRaw[dimRaw]     = {nBinsCentrality,  50,   100};
637    const Double_t lowBinRaw[dimRaw] = {0.0,             0.0,   0.0};
638    const Double_t hiBinRaw[dimRaw]  = {100.0,           100,   1.0};
639    fHJetPtRaw = new THnSparseF("fHJetPtRaw",
640                                 "Incl. jet spectrum [cent,pTjet,A]",
641                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
642    fOutputList->Add(fHJetPtRaw);  
643
644    // raw spectra of LEADING jets  
645    //Centrality, pTjet, A
646    fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
647                                 "Leading jet spectrum [cent,pTjet,A]",
648                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
649    fOutputList->Add(fHLeadingJetPtRaw);  
650
651    // Dphi versus pT jet 
652    //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg 
653    const Int_t dimDp   = 4;
654    const Int_t nBinsDp[dimDp]     = {nBinsCentrality,  50,     50,    50};
655    const Double_t lowBinDp[dimDp] = {0.0,       -TMath::Pi(),   0.0,   0.0};
656    const Double_t hiBinDp[dimDp]  = {100.0,      TMath::Pi(), 100.0, 100.0};
657    fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
658                                 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
659                                 dimDp,nBinsDp,lowBinDp,hiBinDp);
660    fOutputList->Add(fHDphiVsJetPtAll);  
661    */
662
663    //analyze MC generator level 
664    if(fIsChargedMC){   
665       if(fFillRespMx){
666          //Fill response matrix only once 
667          fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 200,0,200, 200,0,200); 
668          fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
669          //....
670          fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", 220,-20,200, 220,-20,200); 
671          fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr 
672          //....
673          fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
674          fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
675          fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr 
676          //....
677          fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",200,0,200); //MC generator charged jet pt spectrum
678          fOutputList->Add(fhJetPtGen);  
679          //....
680          fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",220,-20,200); 
681          fOutputList->Add(fhJetPtSubUeMedianGen);  // with kT median bg subtr
682          //....
683          fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
684          fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
685       
686          //....
687          if(fIsFullMC){
688             fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 200,0,200, 200,0,200);
689             fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
690          //....
691             fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",200,0,200); //MC generator full jet pt spectrum
692             fOutputList->Add(fhJetPtGenFull); 
693          }
694       }
695       //....
696       fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
697       fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
698       fOutputList->Add(fh2NtriggersGen);
699
700       //Centrality, A, pT, pTtrigg
701       fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
702       fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
703       fOutputList->Add(fHJetSpecGen); 
704
705       fHJetSpecSubUeMedianGen = (THnSparseF*)  fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
706       fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
707       fOutputList->Add(fHJetSpecSubUeMedianGen);  
708
709       fHJetSpecSubUeConeGen =  (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen"); 
710       fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
711       fOutputList->Add(fHJetSpecSubUeConeGen);
712       //---
713       fHJetUeMedianGen  =  (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");  
714       fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
715       fOutputList->Add(fHJetUeMedianGen);
716
717       fHJetUeConeGen =  (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen"); 
718       fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle())); 
719       fOutputList->Add(fHJetUeConeGen);
720
721       if(fFillRespMx){
722          //track efficiency/contamination histograms  eta versus pT
723          fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
724          fOutputList->Add(fhPtTrkTruePrimRec); 
725       
726          fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
727          fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");    
728          fOutputList->Add(fhPtTrkTruePrimGen);
729       
730          fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");    
731          fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");    
732          fOutputList->Add(fhPtTrkSecOrFakeRec);
733       }
734
735       fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
736       fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle())); 
737       fOutputList->Add(fHRhoUeMedianVsConeGen);
738    }
739    //-------------------------------------
740    //     pythia histograms
741    const Int_t nBinPt = 150;
742    Double_t binLimitsPt[nBinPt+1];
743    for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
744       if(iPt == 0){
745          binLimitsPt[iPt] = -50.0;
746       }else{// 1.0
747          binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
748       }
749    }
750    
751    fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
752    fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
753    fOutputList->Add(fh1Xsec);
754    fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
755    fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
756    fOutputList->Add(fh1Trials);
757    fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
758    fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
759    fOutputList->Add(fh1AvgTrials);
760    fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
761    fOutputList->Add(fh1PtHard);
762    fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
763    fOutputList->Add(fh1PtHardNoW);
764    fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
765    fOutputList->Add(fh1PtHardTrials);      
766    
767
768    // =========== Switch on Sumw2 for all histos ===========
769    for(Int_t i=0; i<fOutputList->GetEntries(); i++){
770       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
771       if(h1){
772          h1->Sumw2();
773          continue;
774       }
775       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
776       if(hn){
777          hn->Sumw2();
778       }   
779    }
780    TH1::AddDirectory(oldStatus);
781
782    PostData(1, fOutputList);
783 }
784
785 //--------------------------------------------------------------------
786
787 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
788 {
789    //Event loop
790    Double_t eventW  = 1.0;
791    Double_t ptHard  = 0.0;
792    Double_t nTrials = 1.0; // Trials for MC trigger
793    if(fIsChargedMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); 
794
795    if(TMath::Abs((Float_t) fJetParamR)<0.00001){
796       AliError("ANTIKT Cone radius is set to zero.");  
797       return;
798    }
799    if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
800       AliError("KT Cone radius is set to zero.");  
801       return;
802    }
803    if(!strlen(fJetBranchName.Data())){
804       AliError("Jet branch name not set.");
805       return;
806    }
807    if(!strlen(fJetBranchNameBg.Data())){
808       AliError("Jet bg branch name not set.");
809       return;
810    }
811
812
813    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
814    if(!fESD){
815       if(fDebug>1) AliError("ESD not available");
816       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
817    } 
818  
819    fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
820    AliAODEvent* aod = NULL;
821    // take all other information from the aod we take the tracks from
822    if(!aod){
823       if(!fESD) aod = fAODIn;
824       else      aod = fAODOut;
825    }
826
827    if(fNonStdFile.Length()!=0){
828       // case that we have an AOD extension we can fetch the jets from the extended output
829       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
830       fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
831       if(!fAODExtension){
832          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
833       } 
834    }
835     
836    // ----------------- EVENT SELECTION  --------------------------
837    fHistEvtSelection->Fill(1); // number of events before event selection
838
839    // physics selection
840    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
841            ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
842
843    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
844       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
845       fHistEvtSelection->Fill(2);
846       PostData(1, fOutputList);
847       return;
848    }
849   
850    //check AOD pointer
851    if(!aod){
852       if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
853       fHistEvtSelection->Fill(3);
854       PostData(1, fOutputList);
855       return;
856    }
857    
858    // vertex selection
859    AliAODVertex* primVtx = aod->GetPrimaryVertex();
860
861    if(!primVtx){
862       if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
863       fHistEvtSelection->Fill(3);
864       PostData(1, fOutputList);
865       return;
866    }
867
868    Int_t nTracksPrim = primVtx->GetNContributors();
869    Float_t vtxz = primVtx->GetZ();
870    //Input events
871    fhContribVtx->Fill(nTracksPrim);
872    if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
873
874    if((nTracksPrim < fMinContribVtx) ||
875       (vtxz < fVtxZMin) ||
876       (vtxz > fVtxZMax)){
877       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
878                          (char*)__FILE__,__LINE__,vtxz);
879       fHistEvtSelection->Fill(3);
880       PostData(1, fOutputList);
881       return;
882    }else{
883       //Accepted events
884       fhContribVtxAccept->Fill(nTracksPrim);
885       fhVertexZAccept->Fill(vtxz);
886    }
887    //FK// No event class selection imposed
888    // event class selection (from jet helper task)
889    //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
890    //if(fDebug) Printf("Event class %d", eventClass);
891    //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
892    //   fHistEvtSelection->Fill(4);
893    //   PostData(1, fOutputList);
894    //   return;
895    //}
896
897    //------------------ CENTRALITY SELECTION ---------------
898    AliCentrality *cent = 0x0;
899    Double_t centValue  = 0.0; 
900    if(fSystem){  //fSystem=0 for pp,   fSystem=1 for pPb
901       if(fESD){
902          cent = fESD->GetCentrality();
903          if(cent) centValue = cent->GetCentralityPercentile("V0M");
904       }else{
905          centValue = aod->GetHeader()->GetCentrality();
906       }   
907       if(fDebug) printf("centrality: %f\n", centValue);
908       //Input events
909       fhCentrality->Fill(centValue); 
910
911       if(centValue < fCentMin || centValue > fCentMax){
912          fHistEvtSelection->Fill(4);
913          PostData(1, fOutputList);
914          return;
915       }else{
916          //Accepted events
917          fhCentralityAccept->Fill( centValue );
918       }
919    }
920  
921    //-----------------select disjunct event subsamples ----------------
922    Int_t eventnum  = aod->GetHeader()->GetEventNumberESDFile();
923    Int_t lastdigit = eventnum % 10;
924    if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
925       fHistEvtSelection->Fill(5);
926       PostData(1, fOutputList);
927       return;
928    } 
929
930    if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
931    fHistEvtSelection->Fill(0); // accepted events 
932    // ==================== end event selection ============================ 
933  
934    Double_t tmpArrayFive[5];
935  
936    //=============== Reconstructed antikt jets =============== 
937    ReadTClonesArray(fJetBranchName.Data()  , fListJets); 
938    ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg); 
939
940    //============ Estimate background in reconstructed events ===========
941    Double_t rhoFromCellMedian=0.0,    rhoCone=0.0;
942
943    //Find Hadron trigger and Estimate rho from cone
944    TList particleList; //list of tracks
945    Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
946    EstimateBgCone(fListJets, &particleList, rhoCone);
947
948    //Estimate rho from cell median minus jets
949    EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian);
950    //fListJetsBg->Clear(); //this list is further not needed
951
952    //==============  analyze generator level MC  ================ 
953    TList particleListGen; //list of tracks in MC
954
955    if(fIsChargedMC){       
956       //================= generated charged antikt jets ================
957       ReadTClonesArray(fJetBranchNameChargMC.Data(),   fListJetsGen); 
958       ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen); 
959
960       if(fIsFullMC){ //generated full jets
961          ReadTClonesArray(fJetBranchNameFullMC.Data(),  fListJetsGenFull); 
962       }
963
964       //========================================================
965       //serarch for charged trigger at the MC generator level
966
967       Int_t    indexTriggGen = -1;
968       Double_t ptTriggGen    = -1;
969       Int_t    iCounterGen   =  0; //number of entries in particleListGen array
970       Int_t    triggersMC[200];//list of trigger candidates
971       Int_t    ntriggersMC   = 0; //index in triggers array
972
973       if(fESD){//ESD input
974
975          AliMCEvent* mcEvent = MCEvent();
976          if(!mcEvent){
977             PostData(1, fOutputList);
978             return;
979          }
980          
981          AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
982          if(pythiaGenHeader){
983             nTrials = pythiaGenHeader->Trials();
984             ptHard  = pythiaGenHeader->GetPtHard();
985             fh1PtHard->Fill(ptHard,eventW);
986             fh1PtHardNoW->Fill(ptHard,1);
987             fh1PtHardTrials->Fill(ptHard,nTrials);
988          }
989
990          for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
991             if(!mcEvent->IsPhysicalPrimary(it)) continue;
992             AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
993             if(!part) continue;  
994             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){ 
995
996                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
997                   if(indexTriggGen > -1){//trigger candidate was found
998                      triggersMC[ntriggersMC] = indexTriggGen;
999                      ntriggersMC++; 
1000                   }
1001                }
1002
1003                iCounterGen++;//index in particleListGen array 
1004             }
1005          }
1006       }else{  //AOD input
1007          if(!fAODIn){
1008             PostData(1, fOutputList);
1009             return;
1010          }
1011          TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1012          if(!tca){ 
1013             PostData(1, fOutputList);
1014             return;
1015          }
1016          for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1017             AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1018             if(!part) continue;  
1019             if(!part->IsPhysicalPrimary()) continue;
1020             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1021
1022                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1023                   if(indexTriggGen > -1){ //trigger candidater was found
1024                      triggersMC[ntriggersMC] = indexTriggGen;
1025                      ntriggersMC++; 
1026                   }
1027                }
1028
1029                iCounterGen++;//number of entries in particleListGen array
1030             }
1031          }
1032       }
1033
1034       //==============  Estimate bg in generated events ==============
1035       Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1036
1037       //Estimate rho from cone
1038       EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1039
1040       //Estimate rho from cell median minus jets
1041       EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen);
1042       //fListJetsBgGen->Clear();
1043
1044       //============  Generator trigger+jet ==================
1045       if(fHardest==0){ //single inclusive trigger
1046          if(ntriggersMC>0){ //there is at least one trigger 
1047             Int_t rnd     = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1048             indexTriggGen = triggersMC[rnd];
1049          }else{
1050             indexTriggGen = -1; //trigger not found
1051          }
1052       }
1053
1054       //-----------
1055       Int_t ilowGen  = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1056       Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
1057       Bool_t fillOnceGen = kTRUE;
1058       //-----------
1059       
1060       for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1061          indexTriggGen = igen; //trigger hadron 
1062
1063          if(indexTriggGen == -1) continue;  
1064          AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);  
1065          if(!triggerGen) continue;
1066  
1067          if(fHardest >= 2){ 
1068             if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6  
1069          }
1070          if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1071
1072          ptTriggGen = triggerGen->Pt(); //count triggers
1073          fh2NtriggersGen->Fill(centValue, ptTriggGen);
1074
1075          //Count jets and trigger-jet pairs at MC  generator level
1076          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1077             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1078             if(!jet) continue;
1079             Double_t etaJetGen = jet->Eta();
1080             Double_t areaJetGen = jet->EffectiveAreaCharged();
1081  
1082             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ 
1083
1084                Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi()); 
1085                if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1086
1087                //Centrality, A, pT, pTtrigg
1088                tmpArrayFive[0] = centValue;
1089                tmpArrayFive[1] = areaJetGen;
1090                tmpArrayFive[2] = jet->Pt();
1091                tmpArrayFive[3] = ptTriggGen;
1092                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1093                fHJetSpecGen->Fill(tmpArrayFive);
1094
1095
1096                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1097                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1098
1099                //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1100                tmpArrayFive[0] = centValue;
1101                tmpArrayFive[1] = areaJetGen;
1102                tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1103                tmpArrayFive[3] = ptTriggGen;
1104                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1105                fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1106
1107                //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1108                tmpArrayFive[0] = centValue;
1109                tmpArrayFive[1] = areaJetGen;
1110                tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1111                tmpArrayFive[3] = ptTriggGen;
1112                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1113                fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
1114
1115                //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1116                tmpArrayFive[0] = areaJetGen;
1117                tmpArrayFive[1] = jet->Pt();
1118                tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1119                tmpArrayFive[3] = ptUeFromCellMedianGen;
1120                tmpArrayFive[4] = rhoFromCellMedianGen;
1121                fHJetUeMedianGen->Fill(tmpArrayFive); 
1122
1123                //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
1124                tmpArrayFive[0] = areaJetGen;
1125                tmpArrayFive[1] = jet->Pt();
1126                tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1127                tmpArrayFive[3] = ptUeConeGen;
1128                tmpArrayFive[4] = rhoConeGen;
1129                fHJetUeConeGen->Fill(tmpArrayFive);
1130
1131                if(fillOnceGen){
1132                   Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1133                   fHRhoUeMedianVsConeGen->Fill(fillRhoGen); 
1134                   fillOnceGen = kFALSE;
1135                }
1136             }//back to back jet-trigger pair
1137          }//jet loop
1138       }//trigger loop
1139
1140       if(fFillRespMx){
1141          //================ RESPONSE MATRIX ===============
1142          //Count jets and trigger-jet pairs at MC  generator level
1143          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1144             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1145             if(!jet) continue;
1146             Double_t etaJetGen = jet->Eta();
1147             Double_t ptJetGen  = jet->Pt();
1148        
1149             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ 
1150                fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1151
1152                Double_t areaJetGen = jet->EffectiveAreaCharged();
1153                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1154                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1155
1156                fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen); 
1157                fhJetPtSubUeConeGen->Fill(ptJetGen   - ptUeConeGen);    
1158             }
1159          }
1160          if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1161          
1162             Int_t ng = (Int_t) fListJetsGen->GetEntries();
1163             Int_t nr = (Int_t) fListJets->GetEntries();
1164
1165             //Find closest MC generator - reconstructed jets
1166             if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1167             if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
1168
1169             if(fDebug){
1170                Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1171                Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
1172             }
1173             //matching of MC genrator level and reconstructed jets
1174             AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug); 
1175
1176             // Fill response matrix
1177             for(Int_t ir = 0; ir < nr; ir++){
1178                AliAODJet *recJet   = (AliAODJet*) fListJets->At(ir);
1179                Double_t etaJetRec  = recJet->Eta();
1180                Double_t ptJetRec   = recJet->Pt();
1181                Double_t areaJetRec = recJet->EffectiveAreaCharged();
1182                //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1183
1184                if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){ 
1185                   Int_t ig = faGenIndex[ir]; //associated generator level jet
1186                   if(ig >= 0 && ig < ng){
1187                      if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1188                      AliAODJet *genJet  = (AliAODJet*) fListJetsGen->At(ig);
1189                      Double_t ptJetGen  = genJet->Pt();
1190                      Double_t etaJetGen = genJet->Eta();
1191
1192                      //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1193                      if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1194                         fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1195
1196                         Double_t areaJetGen  = genJet->EffectiveAreaCharged();
1197                         Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1198                         Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1199                         Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1200                         Double_t ptUeConeRec         = rhoCone*areaJetRec;
1201                         fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec, 
1202                                                            ptJetGen-ptUeFromCellMedianGen);
1203                         fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
1204                      }
1205                   }//ig>=0
1206                }//rec jet in eta acceptance
1207             }//loop over reconstructed jets
1208          }// # of  rec jets >0
1209       
1210          //=========================== Full jet vs charged jet matrix  ==========================
1211          if(fIsFullMC){
1212             //Count full jets and charged-jet pairs at MC  generator level
1213             for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1214                AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1215                if(!jetFull) continue;
1216                Double_t etaJetFull = jetFull->Eta();
1217                Double_t ptJetFull  = jetFull->Pt();
1218
1219                if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1220                   fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets 
1221                }
1222             }
1223             if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1224                Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1225                Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1226  
1227                //Find closest MC generator full - charged jet
1228                if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1229                if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1230  
1231                if(fDebug){
1232                   Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1233                   Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1234                }
1235                //matching of MC genrator level and reconstructed jets
1236                AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1237
1238               // Fill response matrix
1239                for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1240                   AliAODJet *chJet  = (AliAODJet*) fListJetsGen->At(ichr);
1241                   Double_t etaJetCh = chJet->Eta();
1242                   Double_t ptJetCh  = chJet->Pt();
1243                   //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1244     
1245                   if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1246                      Int_t iful = faGenIndex[ichr]; //associated generator level jet
1247                      if(iful >= 0 && iful < nful){
1248                         if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1249                         AliAODJet *genJetFull  = (AliAODJet*) fListJetsGenFull->At(iful);
1250                         Double_t ptJetFull  = genJetFull->Pt();
1251                         Double_t etaJetFull = genJetFull->Eta();
1252
1253                         //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1254                         if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1255                            fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1256                         }
1257                      }//iful>=0
1258                   }//rec jet in eta acceptance
1259                }//loop over reconstructed jets
1260             }// # of  rec jets >0
1261          }//pointer MC generator jets
1262       } //fill resp mx only for bin 
1263    }//analyze generator level MC
1264     
1265     
1266     
1267    //=============  RECONSTRUCTED INCLUSIVE JETS ===============
1268
1269    Double_t etaJet  = 0.0;
1270    Double_t pTJet   = 0.0;
1271    Double_t areaJet = 0.0;
1272    Double_t phiJet  = 0.0;
1273    //Int_t indexLeadingJet     = -1;
1274    //Double_t pTLeadingJet     = -10.0; 
1275    //Double_t areaLeadingJet   = -10.0;
1276   
1277    for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1278       AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1279       if(!jet) continue;
1280       etaJet  = jet->Eta();
1281       phiJet  = jet->Phi();
1282       pTJet   = jet->Pt();
1283       if(pTJet==0) continue; 
1284      
1285       if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1286       /*areaJet = jet->EffectiveAreaCharged();*/
1287
1288       //Jet Diagnostics---------------------------------
1289       fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1290       fhJetEta->Fill(pTJet, etaJet);
1291       //search for leading jet
1292       /*if(pTJet > pTLeadingJet){
1293          indexLeadingJet  = ij; 
1294          pTLeadingJet     = pTJet; 
1295          areaLeadingJet   = areaJet; 
1296       } 
1297  
1298       // raw spectra of INCLUSIVE jets  
1299       //Centrality, pTjet, A
1300       Double_t fillraw[] = { centValue,
1301                              pTJet,
1302                              areaJet
1303                            };
1304       fHJetPtRaw->Fill(fillraw);*/
1305    }
1306    /*
1307    if(indexLeadingJet > -1){ 
1308       // raw spectra of LEADING jets  
1309       //Centrality, pTjet,  A
1310       Double_t fillleading[] = { centValue,
1311                                  pTLeadingJet,
1312                                  areaLeadingJet
1313                                };
1314       fHLeadingJetPtRaw->Fill(fillleading);
1315    } 
1316    */
1317  
1318    // ===============  RECONSTRUCTED TRIGGER-JET PAIRS ================
1319    if(fIsChargedMC && fFillRespMx){
1320      FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1321    }
1322    Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
1323    //set ranges of the trigger loop
1324    Int_t ilow  = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1325    Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1326
1327    for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1328       indexTrigg = itrk; //trigger hadron with pT >10 GeV 
1329  
1330       if(indexTrigg < 0) continue;
1331
1332       AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);     
1333       if(!triggerHadron){  
1334          PostData(1, fOutputList);
1335          return;
1336       }
1337       if(fHardest >= 2){ 
1338          if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6  
1339       }
1340       if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1341  
1342       //Fill trigger histograms
1343       fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1344       fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1345       fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1346
1347   
1348       //---------- make trigger-jet pairs ---------
1349       //Int_t injet4     = 0;
1350       //Int_t injet      = 0; 
1351
1352       for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1353          AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1354          if(!jet) continue;
1355          etaJet  = jet->Eta();
1356          phiJet  = jet->Phi();
1357          pTJet   = jet->Pt();
1358          if(pTJet==0) continue; 
1359      
1360          if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1361          areaJet = jet->EffectiveAreaCharged();
1362          //if(areaJet >= 0.07) injet++; 
1363          //if(areaJet >= 0.4)  injet4++;
1364
1365          Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet); 
1366          fhDphiTriggerJet->Fill(dphi); //Input
1367
1368          //Dphi versus jet pT   
1369          //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg 
1370          /*Double_t filldp[] = { centValue,
1371                                dphi,
1372                                pTJet,
1373                                triggerHadron->Pt()
1374                              };
1375          fHDphiVsJetPtAll->Fill(filldp);
1376          */ 
1377          // Select back to back trigger - jet pairs
1378          if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1379          fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1380           
1381          //NO bg subtraction
1382          //Centrality, A, pTjet, pTtrigg, dphi
1383          tmpArrayFive[0] = centValue;
1384          tmpArrayFive[1] = areaJet;
1385          tmpArrayFive[2] = pTJet;
1386          tmpArrayFive[3] = triggerHadron->Pt();
1387          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1388          fHJetSpec->Fill(tmpArrayFive);
1389
1390          //subtract bg using estinates base on median of kT jets
1391          Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1392          Double_t ptUeCone         = rhoCone*areaJet;
1393
1394          //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1395          tmpArrayFive[0] = centValue;
1396          tmpArrayFive[1] = areaJet;
1397          tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
1398          tmpArrayFive[3] = triggerHadron->Pt();
1399          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1400          fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1401
1402          //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1403          tmpArrayFive[0] = centValue;
1404          tmpArrayFive[1] = areaJet;
1405          tmpArrayFive[2] = pTJet - ptUeCone;
1406          tmpArrayFive[3] = triggerHadron->Pt();
1407          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1408          fHJetSpecSubUeCone->Fill(tmpArrayFive);
1409
1410          //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1411          tmpArrayFive[0] = areaJet;
1412          tmpArrayFive[1] = pTJet;
1413          tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1414          tmpArrayFive[3] = ptUeFromCellMedian;
1415          tmpArrayFive[4] = rhoFromCellMedian;
1416          fHJetUeMedian->Fill(tmpArrayFive);
1417  
1418          //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
1419          tmpArrayFive[0] = areaJet;
1420          tmpArrayFive[1] = pTJet;
1421          tmpArrayFive[2] = pTJet - ptUeCone;
1422          tmpArrayFive[3] = ptUeCone;
1423          tmpArrayFive[4] = rhoCone;
1424          fHJetUeCone->Fill(tmpArrayFive);
1425
1426          if(filledOnce){ //fill for each event only once
1427             Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1428             fHRhoUeMedianVsCone->Fill(fillRho);
1429             filledOnce = kFALSE;
1430          }
1431       }//jet loop
1432  
1433       //Fill Jet Density In the Event A>0.07
1434       /*if(injet>0){
1435          Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1436                             injet/fkAcceptance,
1437                             triggerHadron->Pt()
1438                           };
1439          fHJetDensity->Fill(filldens);
1440       } 
1441
1442       //Fill Jet Density In the Event A>0.4
1443       if(injet4>0){ 
1444          Double_t filldens4[]={ (Double_t) particleList.GetEntries(), 
1445                                 injet4/fkAcceptance,
1446                                 triggerHadron->Pt()
1447                               };
1448          fHJetDensityA4->Fill(filldens4);
1449       }*/
1450    }
1451
1452
1453    PostData(1, fOutputList);
1454 }
1455
1456 //----------------------------------------------------------------------------
1457 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1458 {
1459    // Draw result to the screen
1460    // Called once at the end of the query
1461
1462    if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1463    if(!GetOutputData(1)) return;
1464 }
1465
1466
1467 //----------------------------------------------------------------------------
1468 Int_t  AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1469    //Fill the list of accepted tracks (passed track cut)
1470    //return consecutive index of the hardest ch hadron in the list
1471    Int_t iCount        = 0;
1472    AliAODEvent *aodevt = NULL;
1473
1474    if(!fESD) aodevt = fAODIn;
1475    else      aodevt = fAODOut;   
1476
1477    if(!aodevt) return -1;
1478
1479    Int_t    index = -1; //index of the highest particle in the list
1480    Double_t ptmax = -10;
1481    Int_t    triggers[200];
1482    Int_t    ntriggers = 0; //index in triggers array
1483
1484    for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1485       AliAODTrack *tr = aodevt->GetTrack(it);
1486       
1487       //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1488       if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1489       if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1490       if(tr->Pt() < fTrackLowPtCut) continue;
1491       list->Add(tr);
1492
1493       //Search trigger:
1494       if(fHardest>0){ //leading particle 
1495          if(tr->Pt()>ptmax){ 
1496             ptmax = tr->Pt();   
1497             index = iCount;
1498          }
1499       }
1500     
1501       if(fHardest==0 && ntriggers<200){ //single inclusive 
1502          if(fTriggerPtRangeLow <= tr->Pt() && 
1503             tr->Pt() < fTriggerPtRangeHigh){ 
1504             triggers[ntriggers] = iCount;
1505             ntriggers++;
1506          }
1507       }
1508
1509       iCount++;
1510    }
1511
1512    if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1513       Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1514       index     = triggers[rnd];
1515    }
1516
1517    return index;
1518 }
1519
1520 //----------------------------------------------------------------------------
1521 /*
1522 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1523    //calculate sum of track pT in the cone perpendicular in phi to the jet 
1524    //jetR = cone radius
1525    // jetPhi, jetEta = direction of the jet 
1526    Int_t numberOfTrks = trkList->GetEntries();
1527    Double_t pTsum = 0.0;
1528    Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1529    for(Int_t it=0; it<numberOfTrks; it++){
1530       AliVParticle *trk = (AliVParticle*) trkList->At(it); 
1531       Double_t dphi = RelativePhi(perpConePhi,trk->Phi());     
1532       Double_t deta = trk->Eta()-jetEta;     
1533
1534       if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1535          pTsum += trk->Pt();
1536       } 
1537    }
1538
1539    return pTsum;
1540 }
1541 */
1542 //----------------------------------------------------------------------------
1543
1544 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1545    //Get relative azimuthal angle of two particles -pi to pi
1546    if      (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1547    else if (vphi > TMath::Pi())  vphi -= TMath::TwoPi();
1548
1549    if      (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1550    else if (mphi > TMath::Pi())  mphi -= TMath::TwoPi();
1551
1552    Double_t dphi = mphi - vphi;
1553    if      (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1554    else if (dphi > TMath::Pi())  dphi -= TMath::TwoPi();
1555
1556    return dphi;//dphi in [-Pi, Pi]
1557 }
1558
1559
1560 //----------------------------------------------------------------------------
1561 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1562    //fill track efficiency denominator
1563    if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1564    if(trk->Charge()==0)        return kFALSE;
1565    if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1566    trkList->Add(trk);
1567    if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
1568
1569    //Search trigger:
1570    if(fHardest>0){ //leading particle 
1571       if(ptLeading < trk->Pt()){
1572          index      = counter;
1573          ptLeading  = trk->Pt();
1574       }
1575    }
1576
1577    if(fHardest==0){ //single inclusive 
1578       index = -1;  
1579       if(fTriggerPtRangeLow <= trk->Pt() && 
1580             trk->Pt() < fTriggerPtRangeHigh){
1581          index      = counter;
1582       }
1583    }
1584
1585    return kTRUE;
1586
1587
1588 //----------------------------------------------------------------------------
1589 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1590    //fill track efficiency numerator 
1591    if(!recList) return;
1592    if(!genList) return;
1593    Int_t nRec = recList->GetEntries();
1594    Int_t nGen = genList->GetEntries();
1595    for(Int_t ir=0; ir<nRec; ir++){ 
1596       AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1597       if(!trkRec) continue; 
1598       Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1599       Bool_t matched = kFALSE;
1600  
1601       for(Int_t ig=0; ig<nGen; ig++){
1602         AliVParticle *trkGen = (AliVParticle*) genList->At(ig);  
1603         if(!trkGen) continue;
1604         Int_t genLabel = TMath::Abs(trkGen->GetLabel()); 
1605         if(recLabel==genLabel){
1606           fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1607           matched = kTRUE;
1608           break;
1609         }
1610       }
1611       if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta()); 
1612    }
1613
1614    return;
1615 }
1616 //________________________________________________________________
1617 void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart,  Double_t &rhoMedian){
1618    //Estimate background rho by means of integrating track pT outside identified jet cones
1619
1620    rhoMedian = 0;
1621
1622    //identify leading jet in the full track acceptance
1623    Int_t idxLeadingJet    = -1;
1624    Double_t pTleading     = 0.0; 
1625    AliAODJet* jet = NULL;
1626  
1627    for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
1628       jet = (AliAODJet*)(listJet->At(ij));
1629       if(!jet) continue;
1630       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
1631       if(jet->Pt() > pTleading){
1632          idxLeadingJet = ij; 
1633          pTleading     = jet->Pt();
1634       }
1635    } 
1636
1637    Int_t  njets = 0;
1638    static Double_t jphi[100];
1639    static Double_t jeta[100];
1640    static Double_t jRsquared[100];
1641
1642    for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
1643
1644       jet = (AliAODJet*)(listJet->At(ij));
1645       if(!jet) continue;
1646       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; 
1647
1648       if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
1649          jphi[njets]      = RelativePhi(jet->Phi(),0.0); //-pi,pi
1650          jeta[njets]      = jet->Eta();
1651          jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
1652          njets++;
1653       }
1654    }
1655
1656    
1657    static Double_t nOutCone[10][4];
1658    static Double_t sumPtOutOfCone[10][4];
1659
1660    for(Int_t ie=0; ie<fnEta; ie++){
1661       for(Int_t ip=0; ip<fnPhi; ip++){
1662          nOutCone[ip][ie]       = 0.0;     //initialize counter
1663          sumPtOutOfCone[ip][ie] = 0.0;
1664       }
1665    }
1666
1667    Double_t rndphi, rndeta;
1668    Double_t rndphishift, rndetashift;
1669    Double_t dphi, deta; 
1670    Bool_t   bIsInCone;
1671
1672    for(Int_t it=0; it<fnTrials; it++){
1673
1674       rndphi = fRandom->Uniform(0, fPhiSize);
1675       rndeta = fRandom->Uniform(0, fEtaSize);
1676                               
1677       for(Int_t ip=0; ip<fnPhi; ip++){  //move radom position to each cell
1678          rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
1679          for(Int_t ie=0; ie<fnEta; ie++){
1680             rndetashift = rndeta + ie*fEtaSize - fEtaSize;
1681
1682             bIsInCone = 0; //tag if trial is in the jet cone
1683             for(Int_t ij=0; ij<njets; ij++){
1684                deta = jeta[ij] - rndetashift;
1685                dphi = RelativePhi(rndphishift,jphi[ij]);
1686                if((dphi*dphi + deta*deta) < jRsquared[ij]){
1687                   bIsInCone = 1;
1688                   break;
1689                }
1690             }
1691             if(!bIsInCone) nOutCone[ip][ie]++;
1692          }
1693       }
1694    }
1695   
1696
1697    //Sum particle pT outside jets in cells
1698    Int_t npart = listPart->GetEntries();
1699    Int_t phicell,etacell; 
1700    AliVParticle *part = NULL; 
1701    for(Int_t ip=0; ip < npart; ip++){ 
1702          
1703       part = (AliVParticle*) listPart->At(ip);
1704       if(!part) continue;
1705       if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
1706  
1707       bIsInCone = 0; //init
1708       for(Int_t ij=0; ij<njets; ij++){
1709          dphi = RelativePhi(jphi[ij], part->Phi());
1710          deta = jeta[ij] - part->Eta();
1711          if((dphi*dphi + deta*deta) < jRsquared[ij]){
1712             bIsInCone = 1;
1713             break;
1714          }
1715       } 
1716       if(!bIsInCone){
1717          phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
1718          etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
1719          sumPtOutOfCone[phicell][etacell]+= part->Pt();
1720       }
1721    }
1722
1723    static Double_t rhoInCells[20];
1724    Double_t  relativeArea;
1725    Int_t  nCells=0;
1726
1727    for(Int_t ip=0; ip<fnPhi; ip++){
1728       for(Int_t ie=0; ie<fnEta; ie++){
1729          relativeArea = nOutCone[ip][ie]/fnTrials;
1730          
1731          if(relativeArea < 0.5) continue; 
1732          rhoInCells[nCells] =  sumPtOutOfCone[ip][ie]/(relativeArea*fCellArea);
1733          nCells++;
1734       }
1735    }
1736
1737    if(nCells>0)
1738      rhoMedian = TMath::Median(nCells, rhoInCells); 
1739
1740 }
1741
1742 //__________________________________________________________________
1743 void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart,  Double_t &rhoPerpCone){
1744
1745    rhoPerpCone = 0.0; //init
1746
1747    if(!listJet) return;
1748    Int_t njet  = listJet->GetEntries();
1749    Int_t npart = listPart->GetEntries();
1750    Double_t pTleading  = 0.0;
1751    Double_t phiLeading = 1000.;
1752    Double_t etaLeading = 1000.;
1753    
1754    for(Int_t jsig=0; jsig < njet; jsig++){ 
1755       AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
1756       if(!signaljet) continue;
1757       if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
1758       if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
1759          pTleading  = signaljet->Pt();
1760          phiLeading = signaljet->Phi();
1761          etaLeading = signaljet->Eta();
1762       }
1763    }   
1764   
1765    Double_t phileftcone  = phiLeading + TMath::Pi()/2; 
1766    Double_t phirightcone = phiLeading - TMath::Pi()/2; 
1767    Double_t dp, de;
1768
1769    for(Int_t ip=0; ip < npart; ip++){ 
1770          
1771       AliVParticle *part = (AliVParticle*) listPart->At(ip);
1772       if(!part){
1773          continue;
1774       }
1775       
1776       dp = RelativePhi(phileftcone, part->Phi());
1777       de = etaLeading - part->Eta();
1778       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
1779
1780       dp = RelativePhi(phirightcone, part->Phi());
1781       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
1782  
1783    }
1784
1785    //normalize total pT by two times cone are 
1786    rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
1787
1788 }
1789 //__________________________________________________________________
1790
1791 void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
1792    //Convert TClones array of jets to TList 
1793
1794    if(!strlen(bname.Data())){
1795       AliError(Form("Jet branch %s not set.", bname.Data()));
1796       return;
1797    }
1798
1799    TClonesArray *array=0x0;
1800
1801    if(fAODOut&&!array){
1802       array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
1803    }
1804    if(fAODExtension&&!array){
1805       array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
1806    }
1807    if(fAODIn&&!array){
1808       array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
1809    }
1810
1811    list->Clear(); //clear list beforehand
1812
1813    if(array){
1814       if(fDebug){
1815          Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
1816       }
1817       for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
1818          AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
1819          if (jet) list->Add(jet);
1820       }
1821    }
1822
1823    return;
1824 }