AliPhysics  1148d43 (1148d43)
AliAnalysisTaskHJetSpectra.cxx
Go to the documentation of this file.
1 #ifndef ALIANALYSISTASKSE_H
2 
3 #include <Riostream.h>
4 #include <TROOT.h>
5 #include <TFile.h>
6 #include <TChain.h>
7 #include <TTree.h>
8 #include <TKey.h>
9 #include <TProfile.h>
10 #include <TProfile2D.h>
11 #include <TH1.h>
12 #include <TH1F.h>
13 #include <TF1.h>
14 #include <TH2F.h>
15 #include <TH1D.h>
16 #include <TH1I.h>
17 #include <TArrayF.h>
18 #include <TArrayD.h>
19 #include <THnSparse.h>
20 #include <TCanvas.h>
21 #include <TList.h>
22 #include <TClonesArray.h>
23 #include <TObject.h>
24 #include <TMath.h>
25 #include <TSystem.h>
26 #include <TInterpreter.h>
27 #include "AliAnalysisTask.h"
28 #include "AliCentrality.h"
29 #include "AliStack.h"
30 #include "AliESDEvent.h"
31 #include "AliESDInputHandler.h"
32 #include "AliAODEvent.h"
33 #include "AliAODHandler.h"
34 #include "AliAnalysisManager.h"
35 #include "AliAnalysisTaskSE.h"
36 #include "AliAnalysisHelperJetTasks.h"
37 #include "AliParticleContainer.h"
38 #include "AliInputEventHandler.h"
39 #endif
40 
41 #include <time.h>
42 #include <TRandom3.h>
43 #include "AliGenEventHeader.h"
44 #include "AliGenPythiaEventHeader.h"
45 #include "AliGenHijingEventHeader.h"
46 #include "AliAODMCHeader.h"
47 #include "AliMCEvent.h"
48 #include "AliLog.h"
49 #include <AliEmcalJet.h>
50 #include <AliPicoTrack.h>
51 #include "AliVEventHandler.h"
52 #include "AliVParticle.h"
53 #include "AliAODMCParticle.h"
54 #include "AliAnalysisUtils.h"
55 #include "AliRhoParameter.h"
56 #include "TVector3.h"
57 #include "AliVVertex.h"
58 #include "AliExternalTrackParam.h"
59 
60 #include <stdio.h>
61 #include <stdlib.h>
62 
63 #include "AliGenDPMjetEventHeader.h"
64 #include "AliJetContainer.h"
65 #include "AliAnalysisTaskEmcal.h"
66 #include "AliAnalysisTaskEmcalJet.h"
67 #include "AliAnalysisTaskHJetSpectra.h"
68 #include "AliHeader.h"
69 #include "AliRunLoader.h"
70 #include "AliVVZERO.h"
71 #include "AliAODZDC.h"
72 #include "AliVZDC.h"
73 
74 using namespace std;
75 
76 // ANALYSIS OF HIGH PT HADRON TRIGGER ASSOCIATED SPECTRUM OF RECOIL JETS IN P+PB
77 // Author Filip Krizek (7.Oct. 2015)
78 
79 ClassImp(AliAnalysisTaskHJetSpectra)
80 //________________________________________________________________________________________
81 
82 AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra():
83 AliAnalysisTaskEmcalJet("AliAnalysisTaskHJetSpectra", kTRUE),
84  fCollisionSystem(0), fTypeOfData(0), fTypeOfAnal(0),
85  fUseDefaultVertexCut(1), fUsePileUpCut(1),
86  fRhoTaskName(), fRhoTaskNameMC(),
87  fSignalJetRadius(0.4), fSignalJetRadiusSquared(fSignalJetRadius*fSignalJetRadius),
88 fSignalJetEtaWindow(0.9 - fSignalJetRadius), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.0),
89 fCentralityType("V0A"), fMinFractionShared(0.5),
90 fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitializedLocal(0),
91 fTTType(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0),
92 fHistEvtSelection(0x0),fhKTAreaPt(0x0),
93  fhVertexZ(0x0), fhVertexXAccept(0x0), fhVertexYAccept(0x0), fhVertexZAccept(0x0),
94  fhVertexXAcceptTT(0x0), fhVertexYAcceptTT(0x0), fhVertexZAcceptTT(0x0),
95 fhVertexZMC(0x0), fhVertexZAcceptMC(0x0),
96  fhDphiTriggerJetAccept(0x0), fhJetPhiIncl(0x0), fhJetEtaIncl(0x0),
97 fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0),
98 fhDiffPtVsPtTrackTrue(0x0),
99 fhTrackPhiCG(0x0), fhTrackPhiTPCG(0x0),
100 fhDCAinXVsPt(0x0),
101 fhDCAinYVsPt(0x0),
102 fhDCAinXVsPtStrange(0x0),
103 fhDCAinYVsPtStrange(0x0),
104 fhDCAinXVsPtNonStrange(0x0),
105 fhDCAinYVsPtNonStrange(0x0),
106 fCentralityBins(kCAll),
107 fNofRandomCones(1),
108 fZVertexCut(10.0),fCutPhi(0.6),
109 fpyVtx(3)
110 {
111  //default constructor
112  for(Int_t it=0; it<kTT;it++){
113  fRhoRec[it].Set(kRho);
114  fRhoMC[it].Set(kRho);
115  fhDphiTTTT[it]=NULL;
116 
117  fhCentralityTT[it] = NULL;
118  }
119 
120  for(Int_t iw=0; iw<21; iw++){
121  fhPtJetPrimVsPtJetRec[iw] = NULL;
122  }
123 
124  for(Int_t ic =0; ic<kCAll; ic++){
125  for(Int_t it=0; it<kTT;it++){
126  fh1Ntriggers[ic][it]=NULL;
127  fh1TriggerMult[ic][it]=NULL;
128  fh1NtriggersGen[ic][it]=NULL;
129  fh1TriggerMultGen[ic][it]=NULL;
130 
131 
132  fhImpactParameterTT[ic][it]=NULL;
133 
134  for(Int_t ir=0; ir<kRho; ir++){
135  fhDphiTriggerJet[ic][it][ir]=NULL;
136  fhDphiTriggerJetGen[ic][it][ir]=NULL;
137  fHJetSpec[ic][it][ir]=NULL;
138  fHJetSpecGen[ic][it][ir]=NULL;
139 
140  fhDeltaPt[ic][it][ir]=NULL;
141  fhDeltaPtEmb[ic][it][ir]=NULL;
142  fhDeltaPtEmb2D[ic][it][ir]=NULL;
143  fhDeltaPtEmbPerp[ic][it][ir]=NULL;
144  fhDeltaPtEmbPerp2D[ic][it][ir]=NULL;
145  fhDeltaPtEmbBc2Bc[ic][it][ir]=NULL;
146  fhDeltaPtEmbBc2Bc2D[ic][it][ir]=NULL;
147 
148  fhRhoTT[ic][it][ir]=NULL;
149  fARhoTT[ic][it][ir]=NULL;
150 
151  }
152  fhVzeroATotMultTT[ic][it]=NULL;
153  fhZNAEnergyTT[ic][it]=NULL;
154  fhTrackMultiplicityTT[ic][it]=NULL;
155  fhZNAVzeroATrackTT[ic][it]=NULL;
156  }
157  for(Int_t it=0; it<kTT; it++){
158  fhJetPhi[ic][it]=NULL;
159  fhJetPhiGen[ic][it]=NULL;
160  }
161  fhTrackPhi[ic]=NULL;
162 
163  for(Int_t it=0; it<kTT; it++){
164  fhJetEta[ic][it]=NULL;
165  fhJetEtaGen[ic][it]=NULL;
166  fhJetEtaRecoil[ic][it]=NULL;
167  fhJetEtaRecoilGen[ic][it]=NULL;
168  fhJetPhiRecoil[ic][it]=NULL;
169  fhJetPhiRecoilGen[ic][it]=NULL;
170  }
171  fhTrackEta[ic]=NULL;
172  fhTrackPt[ic]=NULL;
173  fhTrackPtGen[ic]=NULL;
174  fhCentrality[ic]=NULL;
175  fhVzeroATotMult[ic]=NULL;
176  fhZNAEnergy[ic]=NULL;
177  fhTrackMultiplicity[ic]=NULL;
178  fhZNAVzeroATrack[ic]=NULL;
179  fhImpactParameter[ic]=NULL;
180 
181  for(Int_t ir=0; ir<kRho; ir++){
182  fhJetPtGen[ic][ir]=NULL;
183  fhJetPtGenVsJetPtRec[ic][ir]=NULL;
184  fhJetPtResolutionVsPtGen[ic][ir]=NULL;
185 
186  }
187 
188  for(Int_t ir=0; ir<kRho-1; ir++){
189  fhRhoIncl[ic][ir]=NULL;
190  fhDeltaPtIncl[ic][ir]=NULL;
191  }
192 
193  fhPtTrkTruePrimRec[ic]=NULL;
194  fhPtTrkTruePrimGen[ic]=NULL;
195  fhPtTrkSecOrFakeRec[ic]=NULL;
196  }
197  //Centrality bin borders
198  fCentralityBins[0]=0.;
199  fCentralityBins[1]=20.;
200  fCentralityBins[2]=50.;
201  fCentralityBins[3]=100.;
202  fCentralityBins[4]=1e6; //centrality overflow bin
203 
204  ficb[0]=-1;
205  ficb[1]=-1;
206  ftmpArray[0]=0.;
207  ftmpArray[1]=0.;
208  ftmpArrayX[0]=0.;
209  ftmpArrayX[1]=0.;
210  ftmpArrayX[2]=0.;
211  fVtxArray[0]=0.;
212  fVtxArray[1]=0.;
213  fVtxArray[2]=0.;
214 
215 
216  for(Int_t i=0; i<999; i++){
217  frhovec[i] = 0.;
218 
219  for(Int_t it=0; it<kTT; it++){
220  fTrigTracksGen[it][i]=0x0;
221  fTrigTracks[it][i]=0x0;
222  }
223  }
224 
225  fTTlow[kRef] = 6.0;
226  fTThigh[kRef]= 7.0;
227  fTTlow[kSig] = 12.0;
228  fTThigh[kSig]= 50.0;
229 
230 
231  for(Int_t i=0; i<2; i++){
232  fhInvPtQVsPhi[i] = NULL;
233  fhInvPtQVsEta[i] = NULL;
234  fhInvPtQVsPhiASide[i] = NULL;
235  fhInvPtQVsPhiCSide[i] = NULL;
236  fhSigmaPtOverPtVsPt[i] = NULL;
237  }
238 }
239 
240 //________________________________________________________________________
241 AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra(const char *name) :
242 AliAnalysisTaskEmcalJet(name,kTRUE),
243 fCollisionSystem(0), fTypeOfData(0), fTypeOfAnal(0),
244  fUseDefaultVertexCut(1), fUsePileUpCut(1),
245  fRhoTaskName(), fRhoTaskNameMC(),
246 fSignalJetRadius(0.4), fSignalJetRadiusSquared(fSignalJetRadius*fSignalJetRadius),
247 fSignalJetEtaWindow(0.9 - fSignalJetRadius), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.0),
248 fCentralityType("V0A"), fMinFractionShared(0.5),
249 fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitializedLocal(0),
250 fTTType(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0),
251 fHistEvtSelection(0x0), fhKTAreaPt(0x0),
252  fhVertexZ(0x0), fhVertexXAccept(0x0), fhVertexYAccept(0x0), fhVertexZAccept(0x0),
253  fhVertexXAcceptTT(0x0), fhVertexYAcceptTT(0x0), fhVertexZAcceptTT(0x0),
254 fhVertexZMC(0x0), fhVertexZAcceptMC(0x0),
255 fhDphiTriggerJetAccept(0x0), fhJetPhiIncl(0x0), fhJetEtaIncl(0x0),
256 fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0),
257 /*fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0),*/
258 fhDiffPtVsPtTrackTrue(0x0),
259 fhTrackPhiCG(0x0), fhTrackPhiTPCG(0x0),
260 fhDCAinXVsPt(0x0),
261 fhDCAinYVsPt(0x0),
262 fhDCAinXVsPtStrange(0x0),
263 fhDCAinYVsPtStrange(0x0),
264 fhDCAinXVsPtNonStrange(0x0),
265 fhDCAinYVsPtNonStrange(0x0),
266 fCentralityBins(kCAll),
267 fNofRandomCones(1),
268 fZVertexCut(10.0),fCutPhi(0.6),
269 fpyVtx(3)
270 {
271 //Constructor
272  for(Int_t it=0; it<kTT;it++){
273  fRhoRec[it].Set(kRho);
274  fRhoMC[it].Set(kRho);
275  fhDphiTTTT[it] = NULL;
276  fhCentralityTT[it] = NULL;
277  }
278 
279  for(Int_t iw=0; iw<21; iw++){
280  fhPtJetPrimVsPtJetRec[iw] = NULL;
281  }
282 
283  for(Int_t ic =0; ic<kCAll; ic++){
284  for(Int_t it=0; it<kTT;it++){
285  fh1Ntriggers[ic][it]=NULL;
286  fh1TriggerMult[ic][it]=NULL;
287  fh1NtriggersGen[ic][it]=NULL;
288  fh1TriggerMultGen[ic][it]=NULL;
289 
290  fhImpactParameterTT[ic][it]=NULL;
291 
292 
293  for(Int_t ir=0; ir<kRho; ir++){
294  fhDphiTriggerJet[ic][it][ir]=NULL;
295  fhDphiTriggerJetGen[ic][it][ir]=NULL;
296  fHJetSpec[ic][it][ir]=NULL;
297  fHJetSpecGen[ic][it][ir]=NULL;
298 
299  fhDeltaPt[ic][it][ir]=NULL;
300  fhDeltaPtEmb[ic][it][ir]=NULL;
301  fhDeltaPtEmb2D[ic][it][ir]=NULL;
302  fhDeltaPtEmbPerp[ic][it][ir]=NULL;
303  fhDeltaPtEmbPerp2D[ic][it][ir]=NULL;
304  fhDeltaPtEmbBc2Bc[ic][it][ir]=NULL;
305  fhDeltaPtEmbBc2Bc2D[ic][it][ir]=NULL;
306 
307  fhRhoTT[ic][it][ir]=NULL;
308  fARhoTT[ic][it][ir]=NULL;
309  }
310  fhVzeroATotMultTT[ic][it]=NULL;
311  fhZNAEnergyTT[ic][it]=NULL;
312  fhTrackMultiplicityTT[ic][it]=NULL;
313  fhZNAVzeroATrackTT[ic][it]=NULL;
314  }
315  for(Int_t it=0; it<kTT; it++){
316  fhJetPhi[ic][it]=NULL;
317  fhJetPhiGen[ic][it]=NULL;
318  }
319  fhTrackPhi[ic]=NULL;
320  for(Int_t it=0; it<kTT; it++){
321  fhJetEta[ic][it]=NULL;
322  fhJetEtaGen[ic][it]=NULL;
323  fhJetEtaRecoil[ic][it]=NULL;
324  fhJetEtaRecoilGen[ic][it]=NULL;
325  fhJetPhiRecoil[ic][it]=NULL;
326  fhJetPhiRecoilGen[ic][it]=NULL;
327  }
328  fhTrackEta[ic]=NULL;
329  fhTrackPt[ic]=NULL;
330  fhTrackPtGen[ic]=NULL;
331  fhCentrality[ic]=NULL;
332  fhVzeroATotMult[ic]=NULL;
333  fhZNAEnergy[ic]=NULL;
334  fhTrackMultiplicity[ic]=NULL;
335  fhZNAVzeroATrack[ic]=NULL;
336  fhImpactParameter[ic]=NULL;
337 
338  for(Int_t ir=0; ir<kRho; ir++){
339  fhJetPtGen[ic][ir]=NULL;
340  fhJetPtGenVsJetPtRec[ic][ir]=NULL;
341  fhJetPtResolutionVsPtGen[ic][ir]=NULL;
342 
343  }
344 
345  for(Int_t ir=0; ir<kRho-1; ir++){
346  fhRhoIncl[ic][ir]=NULL;
347  fhDeltaPtIncl[ic][ir]=NULL;
348  }
349 
350  fhPtTrkTruePrimRec[ic]=NULL;
351  fhPtTrkTruePrimGen[ic]=NULL;
352  fhPtTrkSecOrFakeRec[ic]=NULL;
353 
354  }
355 
356  //Centrality bin borders
357  fCentralityBins[0]=0.;
358  fCentralityBins[1]=20.;
359  fCentralityBins[2]=50.;
360  fCentralityBins[3]=100.;
361  fCentralityBins[4]=1e6; //centrality overflow bin
362 
363  ficb[0]=-1;
364  ficb[1]=-1;
365  ftmpArray[0]=0.;
366  ftmpArray[1]=0.;
367  ftmpArrayX[0]=0.;
368  ftmpArrayX[1]=0.;
369  ftmpArrayX[2]=0.;
370  fVtxArray[0]=0.;
371  fVtxArray[1]=0.;
372  fVtxArray[2]=0.;
373 
374 
375  for(Int_t i=0; i<999; i++){
376  frhovec[i] = 0.;
377  for(Int_t it=0; it<kTT; it++){
378  fTrigTracksGen[it][i]=0x0;
379  fTrigTracks[it][i]=0x0;
380  }
381  }
382 
383  fTTlow[kRef] = 6.0;
384  fTThigh[kRef]= 7.0;
385  fTTlow[kSig] = 12.0;
386  fTThigh[kSig]= 50.0;
387 
388  for(Int_t i=0; i<2; i++){
389  fhInvPtQVsPhi[i] = NULL;
390  fhInvPtQVsEta[i] = NULL;
391  fhInvPtQVsPhiASide[i] = NULL;
392  fhInvPtQVsPhiCSide[i] = NULL;
393  fhSigmaPtOverPtVsPt[i] = NULL;
394  }
395 
396 
397  //inclusive pT spectrum times the boost function
398 
399  DefineOutput(1, TList::Class());
400 }
401 //________________________________________________________________________
402 void AliAnalysisTaskHJetSpectra::SetTT(Double_t tlr, Double_t thr,Double_t tls, Double_t ths){
403  //set trigger track pT bins
404  fTTlow[kRef] = tlr;
405  fTThigh[kRef] = thr;
406  fTTlow[kSig] = tls;
407  fTThigh[kSig] = ths;
408 }
409 //________________________________________________________________________
410 Double_t AliAnalysisTaskHJetSpectra::GetConePt(Double_t eta, Double_t phi, Double_t radius, AliParticleContainer *recTrkCont, Bool_t isGen){
411  //sum up pt inside a cone
412  Double_t tmpConePt = 0.0;
413 
414  if(!recTrkCont) return 0.0;
415 
416  AliVParticle* tmpTrack=NULL;
417 
418  if(recTrkCont){
419  recTrkCont->ResetCurrentID();
420  while((tmpTrack = recTrkCont->GetNextAcceptParticle())){
421  if(!tmpTrack) continue;
422  if(IsTrackInAcceptance(tmpTrack, isGen)){
423  if(GetDeltaR(tmpTrack->Phi(), phi, tmpTrack->Eta(), eta) < radius){
424  tmpConePt = tmpConePt + tmpTrack->Pt();
425  }
426  }
427  }
428  }
429  return tmpConePt;
430 }
431 
432 //________________________________________________________________________
433 Double_t AliAnalysisTaskHJetSpectra::GetImpactParameter(){
434  //get impact parameter from hijing
435  AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader());
436  if(MCEvent()){
437  if(!hijingHeader){
438  // Check if AOD
439  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
440 
441  if(aodMCH){
442  for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++){
443  hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i));
444  if(hijingHeader) break;
445  }
446  }
447  }
448  }
449  if(hijingHeader){
450  return (Double_t) (hijingHeader->ImpactParameter());
451  }
452  AliWarning(Form("In task %s: GetImpactParameter() failed!", GetName()));
453  return -1.0;
454 }
455 //________________________________________________________________________
456 Double_t AliAnalysisTaskHJetSpectra::GetSimPrimaryVertex(){
457  //get generator level primary vertex
458 
459  AliGenEventHeader* mcHeader = NULL;
460  AliAODMCHeader* aodMCH = NULL;
461 
462  if(fTypeOfAnal == kKine){ //KINE
463  AliRunLoader *rl = AliRunLoader::Instance();
464  if(rl) mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
465 
466  } else {
467 
468  if(MCEvent()){
469  if(fTypeOfData == kPythia){
470  mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
471  if(!mcHeader){
472  // Check if AOD
473  aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
474 
475  if(aodMCH){
476  for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){
477  mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
478  if(mcHeader) break;
479  }
480  }
481  }
482  }else if(fTypeOfData == kHijing){
483  mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader());
484  if(!mcHeader){
485  // Check if AOD
486  aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
487 
488  if(aodMCH){
489  for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){
490  mcHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i));
491  if(mcHeader) break;
492  }
493  }
494  }
495  }else if(fTypeOfData == kDmpjet){
496  mcHeader = dynamic_cast<AliGenDPMjetEventHeader*>(MCEvent()->GenEventHeader());
497  if(!mcHeader){
498  // Check if AOD
499  aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
500 
501  if(aodMCH){
502  for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){
503  mcHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aodMCH->GetCocktailHeader(i));
504  if(mcHeader) break;
505  }
506  }
507  }
508  }
509  }
510  }
511 
512 
513  if(mcHeader){
514 
515  mcHeader->PrimaryVertex(fpyVtx);
516  return (Double_t) (fpyVtx[2]);
517  }
518  AliWarning(Form("In task %s: Pythia Vertex failed!", GetName()));
519  return 9999.0;
520 }
521 
522 //________________________________________________________________________
523 Bool_t AliAnalysisTaskHJetSpectra::IsMCEventInAcceptance(AliVEvent* event){
524  //EVENT SELECTION PURE MC
525 
526  if(!event) return kFALSE;
527 
528  //___________________________________________________
529 
530  if(fTypeOfAnal!= kKine && !MCEvent()) return kFALSE;
531 
532 
533  Double_t vtxMC = GetSimPrimaryVertex();
534  fhVertexZMC->Fill(vtxMC); //Fill BEFORE vertex cut
535 
536  if(TMath::Abs(vtxMC) > fZVertexCut){
537  fHistEvtSelection->Fill(3); //count events rejected by vertex cut
538  return kFALSE;
539  }
540  fhVertexZAcceptMC->Fill(vtxMC);//Fill AFTER vertex cut
541 
542  return kTRUE;
543 
544 }
545 //________________________________________________________________________
546 Bool_t AliAnalysisTaskHJetSpectra::IsEventInAcceptance(AliVEvent* event){
547  //EVENT SELECTION RECONSTRUCTED DATA
548 
549 
550  if(!event) return kFALSE;
551 
552 
553  //___________________________________________________
554  //TEST PILE UP
555  if(fUsePileUpCut){
556  if(!fHelperClass || fHelperClass->IsPileUpEvent(event)){
557  fHistEvtSelection->Fill(2); //count events rejected by pileup
558  return kFALSE;
559  }
560  }
561  //___________________________________________________
562  //BEFORE VERTEX CUT
563  fhVertexZ->Fill(event->GetPrimaryVertex()->GetZ());
564 
565  if(fUseDefaultVertexCut){
566  if(!fHelperClass || !fHelperClass->IsVertexSelected2013pA(event)){
567  fHistEvtSelection->Fill(3); //count events rejected by vertex cut
568  return kFALSE;
569  }
570  }else{
571  if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > fZVertexCut){
572  fHistEvtSelection->Fill(3); //count events rejected by vertex cut
573  return kFALSE;
574  }
575  if(event->GetPrimaryVertex()->GetNContributors()<1){
576  fHistEvtSelection->Fill(3); //count events rejected by vertex cut
577  return kFALSE;
578  }
579  }
580  //___________________________________________________
581  //AFTER VERTEX CUT
582  fhVertexXAccept->Fill(event->GetPrimaryVertex()->GetX());
583  fhVertexYAccept->Fill(event->GetPrimaryVertex()->GetY());
584  fhVertexZAccept->Fill(event->GetPrimaryVertex()->GetZ());
585  fVtxArray[0]=event->GetPrimaryVertex()->GetX();
586  fVtxArray[1]=event->GetPrimaryVertex()->GetY();
587  fVtxArray[2]=event->GetPrimaryVertex()->GetZ();
588 
589 
590  return kTRUE;
591 }
592 
593 //________________________________________________________________________
594 Bool_t AliAnalysisTaskHJetSpectra::IsTrackInAcceptance(AliVParticle* track, Bool_t isGen){
595  // Check if the track pt and eta range
596  if(!track) return kFALSE;
597 
598  if(isGen){ //pure MC select charged primary tracks
599  //Apply only for kine level or MC containers
600  if(!track->Charge()) return kFALSE;
601  if(fTypeOfAnal == kEff){
602  if(!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) return kFALSE;
603  }
604  }
605  if(TMath::Abs(track->Eta()) <= fTrackEtaWindow){ //APPLY TRACK ETA CUT
606  if(track->Pt() >= fMinTrackPt){ //APPLY TRACK CUT
607  return kTRUE;
608  }
609  }
610  return kFALSE;
611 }
612 
613 
614 //________________________________________________________________________
615 Bool_t AliAnalysisTaskHJetSpectra::IsSignalJetInAcceptance(AliEmcalJet *jet, Bool_t suppressGhost){
616  //select jets in acceptance
617  if(!jet) return kFALSE;
618  if(TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow){
619  if(jet->Area() >= fMinJetArea){
620  if(suppressGhost){
621  if(jet->Pt() >= fMinTrackPt) return kTRUE;
622  }else{
623  return kTRUE;
624  }
625  }
626  }
627  return kFALSE;
628 }
629 
630 
631 //________________________________________________________________________
632 void AliAnalysisTaskHJetSpectra::ExecOnceLocal(){
633  // Initialization of jet containers done in AliAnalysisTaskEmcalJet::ExecOnce()
634  //Read arrays of jets and tracks
635  fInitializedLocal = kTRUE;
636 
637  // Initialize helper class (for vertex selection & pile up correction)
638  fHelperClass = new AliAnalysisUtils();
639  fHelperClass->SetCutOnZVertexSPD(kFALSE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm
640 
641  return;
642 }
643 
644 
645 //________________________________________________________________________
646 void AliAnalysisTaskHJetSpectra::GetDeltaPt(Int_t nrho, TArrayD &rho, Double_t *dpt,
647  Double_t ttPhi, Double_t ttEta, AliParticleContainer *recTrkCont, Bool_t isGen){
648  //Double_t leadingJetExclusionProbability)
649 
650  //delta pt = pT in random cone - rho * Area of random cone
651  // processes real reconstructed data. Exclude region around jet with TT
652 
653  for(Int_t ir=0;ir<nrho;ir++){
654  dpt[ir] = -10000.0; // Set an invalid delta pt
655  }
656 
657  // Define random cone Eta+Phi
658  Bool_t coneValid = kTRUE;
659  Double_t tmpRandConeEta = -fSignalJetEtaWindow + fRandom->Rndm()*2*fSignalJetEtaWindow;
660  Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
661 
662  if(ttEta > -2.0){ //make sure that one generates RC far away from the trigger track jet
663  while(GetDeltaR( tmpRandConePhi, ttPhi, tmpRandConeEta, ttEta)<fSignalJetRadius+0.1){
664  tmpRandConeEta = -fSignalJetEtaWindow + fRandom->Rndm()*2*fSignalJetEtaWindow;
665  tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
666  }
667  }
668 
669  //if(fDebug>20) Printf("RND CONE RCphi=%f RCeta=%f TTphi=%f TTeta=%f",
670  // tmpRandConePhi,tmpRandConeEta,ttPhi,ttEta);
671 
672  //Based on Ncoll probability that Reject a jet based on the probability check for overlap if demanded
673  /*
674 
675  //WHAT TO DO HERE ??? SHOULD RANDOM CONES BE CORRELAtED WITH JETs and rejected based on 1./Ncoll ?
676 
677  if(leadingJetExclusionProbability>0){ //skips pp
678 
679  AliEmcalJet* tmpLeading = NULL;
680  Double_t lpt = -1.0;
681  // Get leading jet (regardless of pT)
682 
683  AliEmcalJet* tmpJet = NULL;
684  AliJetContainer *jetContRec = GetJetContainer(kContainerOne);
685  if(jetContRec){
686  jetContRec->ResetCurrentID();
687  while((tmpJet = jetContRec->GetNextAcceptJet())) {
688  if(!tmpJet) continue;
689  if((TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) && (tmpJet->Area() >= fMinJetArea)){
690  if(tmpJet->Pt() > lpt){
691  tmpLeading = tmpJet;
692  lpt = tmpJet->Pt();
693  }
694  }
695  }
696 
697  if(tmpLeading){
698  Double_t tmpDeltaPhi = RelativePhi(tmpRandConePhi, tmpLeading->Phi());
699  Double_t tmpDeltaEta = tmpLeading->Eta()-tmpRandConeEta;
700 
701  // Check, if cone has overlap with jet
702  if(sqrt(tmpDeltaPhi*tmpDeltaPhi + tmpDeltaEta*tmpDeltaEta) <= fSignalJetRadius){
703  // probability to exclude the RC
704  // Only exclude cone with a given probability
705  if(fRandom->Rndm()<=leadingJetExclusionProbability) coneValid = kFALSE;
706  }
707  }
708  }
709  if(fRandom->Rndm()<=leadingJetExclusionProbability) coneValid = kFALSE;
710  }*/
711 
712  // Get the cones' pt and calculate delta pt
713  if(coneValid){
714  Double_t conePt = GetConePt(tmpRandConeEta,tmpRandConePhi, fSignalJetRadius, recTrkCont, isGen);
715  for(Int_t ir=0; ir < nrho; ir++){
716  dpt[ir] = conePt - (rho[ir]*fSignalJetRadiusSquared*TMath::Pi());
717  }
718  }
719 
720 }
721 
722 //________________________________________________________________________
723 Bool_t AliAnalysisTaskHJetSpectra::FillHistograms(){
724  // executed in each event
725  //called in AliAnalysisTaskEmcal::UserExec(Option_t *)
726  // Analyze the event and Fill histograms
727 
728  if(!InputEvent()){
729  AliError("??? Event pointer == 0 ???");
730  return kFALSE;
731  }
732 
733  //Execute only once: Get tracks, jets from arrays if not already given
734  if(!fInitializedLocal) ExecOnceLocal();
735 
736  //_________________________________________________________________
737  // fill MC information
738 
739  //if(fTypeOfAnal == kKine && fTypeOfData == kPythia){
740  // fh1PtHard->Fill(GetPtHard()); //Fills cross section
741  //}
742 
743  //_________________________________________________________________
744  // FILL EVENT STATISTICS
745  fHistEvtSelection->Fill(1); //Count input event
746 
747  if(fTypeOfData != kReal){ //Check MC event vertex
748  if(!IsMCEventInAcceptance(InputEvent())) return kFALSE; //post data is in UserExec
749  }
750 
751  //Check Reconstructed event vertex
752  if(fTypeOfAnal != kKine){ //Check MC event vertex
753  if(!IsEventInAcceptance(InputEvent())) return kFALSE; //post data is in UserExec
754  }
755 
756  //___________________
757  // Get centrality
758  Double_t centralityPercentile = -1.0; //KINE
759  Double_t centralityPercentileV0A = -1.0;
760  Double_t centralityPercentileV0C = -1.0;
761  Double_t centralityPercentileV0M = -1.0;
762  Double_t centralityPercentileZNA = -1.0;
763  Double_t multVzero = -1.0;
764  Double_t energyZdcNA = -1.0;
765  ficb[0] = 0; //0=min bias
766  ficb[1] = -1; //-1 assigned centrality bin
767 
768  if(InputEvent()->GetVZEROData()){
769  multVzero = InputEvent()->GetVZEROData()->GetMTotV0A();
770  }
771 
772  //if(InputEvent()->GetZDCData()){
773  // AliESDZDC.h : C --> 1 A--->2
774  //Double_t GetZNAEnergy() const {return (Double_t) fZDCN2Energy;}
775  //Double_t GetZNCEnergy() const {return (Double_t) fZDCN1Energy;}
776  // energyZdcNA = InputEvent()->GetZDCData()->GetZNAEnergy();
777  //energyZdcNA = InputEvent()->GetZDCN2Energy(); //should be ZNA
778  //}
779  AliVZDC *aodZDC = dynamic_cast<AliVZDC*> (InputEvent()->GetZDCData());
780  if(aodZDC){
781  const Double_t *ZNAtower = aodZDC->GetZNATowerEnergy();
782  energyZdcNA = ZNAtower[0]*4.*82./208./12.96; //ZNA energy in TeV
783  }
784 
785  if(fCollisionSystem != kpp){ //KINE Check MC event vertex
786  AliCentrality* tmpCentrality = InputEvent()->GetCentrality();
787  if(!tmpCentrality){
788  fHistEvtSelection->Fill(4);
789  return kFALSE; //post data is in UserExec
790  }else{
791  centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
792  centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
793  centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
794  centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
795  centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA");
796  }
797 
798  if((centralityPercentile < fCentralityBins[0]) || (centralityPercentile > fCentralityBins[4])){ //cut on centrality
799  AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile));
800  fHistEvtSelection->Fill(4);
801  return kFALSE;
802  }
803 
804  if(!(fTypeOfData == kPythia || fTypeOfData == kDmpjet || fTypeOfAnal == kKine)){
805  for(Int_t ic=0; ic<fCentralityBins.GetSize()-1;ic++){
806  if(fCentralityBins[ic] <= centralityPercentile &&
807  centralityPercentile < fCentralityBins[ic+1]){
808  ficb[1] = ic+1; //MB is ficb[0]=0; other centralities ficb[1]>0
809  }
810  }
811  }
812 
813  for(Int_t ic=0; ic<2; ic++){
814  if(ficb[ic]==-1) continue; //MB + CENT BIASED
815  fhCentrality[ficb[ic]]->Fill(centralityPercentile);
816  }
817 
818  //MB
819  fhCentralityV0M->Fill(centralityPercentileV0M);
820  fhCentralityV0A->Fill(centralityPercentileV0A);
821  fhCentralityV0C->Fill(centralityPercentileV0C);
822  fhCentralityZNA->Fill(centralityPercentileZNA);
823  }
824 
825  for(Int_t ic=0; ic<2; ic++){
826  if(ficb[ic]==-1) continue; //MB +CENT bias
827  fhVzeroATotMult[ficb[ic]]->Fill(multVzero);
828  fhZNAEnergy[ficb[ic]]->Fill(energyZdcNA);
829 
830  if(fTypeOfData == kHijing){
831  fImpParam = GetImpactParameter();
832  fhImpactParameter[ficb[ic]]->Fill(fImpParam);
833  }
834  }
835  //cout<<"energyZdcNA = "<<energyZdcNA<<endl;
836 
837  fHistEvtSelection->Fill(0); //Count Accepted input event
838 
839  // END EVENT SELECTION
840  //_________________________________________________________________
841  TClonesArray* arrayMC = 0; // array particles in the MC event
842  AliAODEvent* aodIn = 0;
843  Int_t iNTracksMC = 0; //number of particles in arrayMC
844  if(fTypeOfAnal == kEff){
845  aodIn = dynamic_cast<AliAODEvent*>(InputEvent()); // input AOD
846  if(!aodIn){
847  AliError("No input AOD found!");
848  return kFALSE;
849  }
850  arrayMC = (TClonesArray*) aodIn->FindListObject(AliAODMCParticle::StdBranchName());
851  if(!arrayMC){
852  AliError("No MC array found!");
853  return kFALSE;
854  }
855 
856  iNTracksMC = arrayMC->GetEntriesFast();
857  }
858 
859 
860  //_________________________________________________________________
861  // JET+TRACK CONTAINERS
862  AliJetContainer *jetContRec = NULL; //AKTjet container from reconstruced tracks
863  AliJetContainer *jetContGen = NULL; //AKT jet container from MC particles
864  AliJetContainer *jetContRecKT = NULL; //KT jet container from reconstruced tracks
865  AliJetContainer *jetContGenKT = NULL; //KT jet container from MC particles
866 
867  AliEmcalJet *jetGen = NULL; //jet pointer MC jet
868  AliEmcalJet *jetRec = NULL; //jet pointer real jet
869 
870  AliParticleContainer *trkContRec = NULL; //track contaier real reconstructed tracks
871  AliParticleContainer *parContGen = NULL; //particle container of MC particles
872 
873  AliVParticle *constTrackRec = NULL; //jet constituent
874  AliVParticle *constTrackGen = NULL; //jet constituent
875  //_________________________________________________________
876  //READ JET TRACK CONTAINERS
877  if(fTypeOfAnal == kRec){
878  //REAL DATA
879  jetContRec = GetJetContainer(kContainerOne); //AKT jet
880  trkContRec = GetParticleContainer(kContainerOne); //reconstructed particle container
881  jetContRecKT = GetJetContainer(kContainerTwo);//KT jet
882  }
883 
884  if(fTypeOfAnal == kEff || fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
885  //MC + EMB DATA
886  jetContRec = GetJetContainer(kContainerOne); //AKT
887  trkContRec = GetParticleContainer(kContainerOne); //reconstructed particle container
888  jetContGen = GetJetContainer(kContainerTwo); //AKT generator level jets
889  parContGen = GetParticleContainer(kContainerTwo); //true MC particle container
890 
891  jetContRecKT = GetJetContainer(kContainerThree); //KT reconstructed jets
892  jetContGenKT = GetJetContainer(kContainerFour); //KT generator level jets
893  }
894 
895  if(fTypeOfAnal == kKine){ // Kine written to the 0th container !!!!!!!!!!!!!!!!
896  //KINE
897  jetContGen = GetJetContainer(kContainerOne); //AKT jets
898  parContGen = GetParticleContainer(kContainerOne); //kine particle container
899 
900  jetContGenKT = GetJetContainer(kContainerTwo); //KT jets
901  }
902 
903  //if(fDebug>20) Printf("POINTER TO CONTAINERS JETrec=%p TRKrec=%p JETgen=%p TRKgen=%p",
904  // jetContRec,trkContRec,jetContGen, parContGen);
905 
906 
907  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
908  //LOOP OVER TRACKS SEARCH FOR TRIGGER CANDIDATES IN GEN TRACKS
909  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
910 
911  Int_t indexSingleRndTrigGen[kTT] = {-1,-1}; //index of single random trigger
912  AliVParticle* trackTTGen[kTT] = {NULL,NULL}; //pointer to the trigger hadron
913  Int_t ntriggersGen[kTT] = {0,0};
914 
915  Double_t sumFakeTrackPtInJetNotStrange = 0.;//sum of track pT from non strange secondaries
916  Double_t sumFakeTrackPtInJetStrange = 0.;//sum of track pT from strange secondaries
917  Double_t sumStrangeTrackPtInJet = 0.;//sum of track pT from strange secondaries
918  AliAODMCParticle* particleMC = NULL; //
919  AliAODMCParticle* particleMCMother = NULL;//
920  Int_t iIndexMother = 0; //
921  Int_t iPdgCode = 0;//
922  Int_t iPdgCodeMother = 0;//
923  Bool_t bStrange = 0;
924  Int_t lab = 0;
925  Double_t newJetPt=0., newSumFakePt=0., enhw =0.; //
926  Double_t pTofTT, dphiTTTT;
927 
928  if(fTypeOfAnal == kEff || fTypeOfAnal == kKine){
929 
930 
931  if(parContGen){
932  parContGen->ResetCurrentID();
933  while((constTrackGen = (AliVParticle*) parContGen->GetNextAcceptParticle())){
934  if(!constTrackGen) continue;
935 
936  if(IsTrackInAcceptance(constTrackGen, kTRUE)){
937  for(Int_t ic=0; ic<2; ic++){
938  if(ficb[ic]==-1) continue;
939  fhTrackPtGen[ficb[ic]]->Fill(constTrackGen->Pt()); //inclusive pT spectrum of tracks
940  }
941 
942  for(Int_t it=0; it<kTT; it++){
943  if((fTTlow[it] <= constTrackGen->Pt()) && (constTrackGen->Pt() < fTThigh[it])){
944  if(ntriggersGen[it]<999){
945  fTrigTracksGen[it][ntriggersGen[it]] = constTrackGen; //GEN trigger candidates
946  ntriggersGen[it]++;
947  }
948  }
949  }
950  }
951  }
952  }
953 
954 
955  for(Int_t it=0; it<kTT; it++){
956  for(Int_t ic=0; ic<2; ic++){
957  if(ficb[ic]==-1) continue;
958  fh1TriggerMultGen[ficb[ic]][it]->Fill(ntriggersGen[it]);
959  }
960 
961  if(ntriggersGen[it]==1){ //only one trigger candidate
962  indexSingleRndTrigGen[it] = 0;
963  trackTTGen[it] = (AliVParticle*) (fTrigTracksGen[it][indexSingleRndTrigGen[it]]);
964 
965  }else if(ntriggersGen[it]>1){ //more trigger candiates
966 
967  if(fTTType){ //use TT with the highest pT
968  pTofTT = -1.0;
969 
970  for(Int_t ik=0; ik<ntriggersGen[it]; ik++){
971  if(pTofTT < ((AliVParticle*) (fTrigTracksGen[it][ik]))->Pt()){
972  indexSingleRndTrigGen[it] = ik;
973  pTofTT = ((AliVParticle*) (fTrigTracksGen[it][ik]))->Pt();
974  }
975  }
976  trackTTGen[it] = (AliVParticle*) (fTrigTracksGen[it][indexSingleRndTrigGen[it]]);
977  }else{ //TT selection without weighting
978  indexSingleRndTrigGen[it] = fRandom->Integer(ntriggersGen[it]); //Integer 0 ... ntriggers-1
979  trackTTGen[it] = (AliVParticle*) (fTrigTracksGen[it][indexSingleRndTrigGen[it]]);
980  }
981  }
982  }
983  }
984  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
985  //LOOP OVER TRACKS SEARCH FOR TRIGGER CANDIDATES IN REC TRACKS
986  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
987 
988  Int_t indexSingleRndTrig[kTT] = {-1,-1}; //index of single random trigger
989  AliVParticle* trackTT[kTT] = {NULL,NULL};
990  Int_t ntriggers[kTT] = {0, 0};
991 
992  if(fTypeOfAnal != kKine){
993 
994  if(trkContRec){
995  trkContRec->ResetCurrentID();
996  while((constTrackRec = (AliVParticle*) (trkContRec->GetNextAcceptParticle()))){
997  if(!constTrackRec) continue;
998 
999  if(fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
1000  //if(fDebug>99) Printf("TRACK LABEL %d", track->GetLabel());
1001  //in embed evts search for TT only among real tracks; do not consider embedded tracks as trigger
1002  if(TMath::Abs(constTrackRec->GetLabel()) == 99999) continue;//9999 set in AddTaskHJetSpectra.C
1003  }
1004 
1005 
1006  if(IsTrackInAcceptance(constTrackRec, kFALSE)){ //rec
1007  //Fill some inclusive spectra (Same analysis for gen level?)
1008  for(Int_t ic=0; ic<2; ic++){
1009  if(ficb[ic]==-1) continue;
1010  fhTrackPhi[ficb[ic]]->Fill(constTrackRec->Pt(), RelativePhi(constTrackRec->Phi(),0.0)); // phi = -pi,pi
1011  fhTrackEta[ficb[ic]]->Fill(constTrackRec->Pt(), constTrackRec->Eta());
1012  fhTrackPt[ficb[ic]]->Fill(constTrackRec->Pt());
1013  }
1014 
1015  for(Int_t it=0; it<kTT; it++){
1016  if(fTTlow[it] <= constTrackRec->Pt() && constTrackRec->Pt() < fTThigh[it]){
1017  if(ntriggers[it]<999){
1018  fTrigTracks[it][ntriggers[it]] = constTrackRec; //trigger candidates
1019  ntriggers[it]++;
1020  }
1021  }
1022  }
1023  }
1024  }
1025  }
1026 
1027  // SELECT SINGLE INCLUSIVE REC LEVEL TRIGGER
1028  for(Int_t it=0; it<kTT; it++){
1029  for(Int_t ic=0; ic<2; ic++){
1030  if(ficb[ic]==-1) continue;
1031  fh1TriggerMult[ficb[ic]][it]->Fill(ntriggers[it]);
1032  }
1033 
1034  // if(ntriggers[it]>0){
1035  // indexSingleRndTrig[it] = fRandom->Integer(ntriggers[it]); //Integer 0 ... ntriggers-1
1036  // trackTT[it] = (AliVParticle*) (fTrigTracks[it][indexSingleRndTrig[it]]);
1037  //if(fDebug>20) Printf("TT index = %d size =%d", indexSingleRndTrig, (int)fTrigTracks.size());
1038  // }
1039  if(ntriggers[it]==1){ //only one trigger candidate
1040  indexSingleRndTrig[it] = 0;
1041  trackTT[it] = (AliVParticle*) (fTrigTracks[it][indexSingleRndTrig[it]]);
1042 
1043  }else if(ntriggers[it]>1){ //multiple trigger candiates
1044 
1045  if(fTTType){
1046 
1047  for(Int_t ik=0; ik<ntriggers[it]; ik++){
1048 
1049  if(ik>0){ //azimuthal angle between triggers
1050  dphiTTTT = RelativePhi(((AliVParticle*) (fTrigTracks[it][0]))->Phi(),
1051  ((AliVParticle*) (fTrigTracks[it][ik]))->Phi());
1052  if(dphiTTTT<-0.5*TMath::Pi()) dphiTTTT += TMath::TwoPi();
1053  if(dphiTTTT> 1.5*TMath::Pi()) dphiTTTT -= TMath::TwoPi();
1054 
1055  fhDphiTTTT[it]->Fill(dphiTTTT);
1056  }
1057  }
1058 
1059  pTofTT = -1.0; //select TT with the highest pT
1060  for(Int_t ik=0; ik<ntriggers[it]; ik++){
1061  if(pTofTT < ((AliVParticle*) (fTrigTracks[it][ik]))->Pt()){
1062  indexSingleRndTrig[it] = ik;
1063  pTofTT = ((AliVParticle*) (fTrigTracks[it][ik]))->Pt();
1064  }
1065  }
1066 
1067  trackTT[it] = (AliVParticle*) (fTrigTracks[it][indexSingleRndTrig[it]]);
1068  }else{ //TT selection without weighting
1069  indexSingleRndTrig[it] = fRandom->Integer(ntriggers[it]);
1070  trackTT[it] = (AliVParticle*) (fTrigTracks[it][indexSingleRndTrig[it]]);
1071  }
1072  }
1073  }
1074 
1075  if(ntriggers[kRef]>0 || ntriggers[kSig]>0){
1076  fhVertexXAcceptTT->Fill(fVtxArray[0]);
1077  fhVertexYAcceptTT->Fill(fVtxArray[1]);
1078  fhVertexZAcceptTT->Fill(fVtxArray[2]);
1079  }
1080  }
1081 
1082 
1083  //___________________________________________________________
1084  // CALCULATE RHO
1085  for(Int_t it=0;it<kTT;it++){
1086  fRhoRec[it].Reset(0.);
1087  fRhoMC[it].Reset(0.);
1088  }
1089 
1090  for(Int_t it=0;it<kTT;it++){
1091  if(fTypeOfAnal == kRec || fTypeOfAnal == kEff || fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
1092  fRhoRec[it][kConeRho] = EstimateBgCone( jetContRec, trkContRec, trackTT[it], kFALSE); //container ID=0 reconstructed tracks
1093  fRhoRec[it][kCMSRho] = EstimateBgKTcms(jetContRecKT, trkContRec, trackTT[it]);
1094  fRhoRec[it][kKtRho] = EstimateBgKT( jetContRecKT, trkContRec, trackTT[it]);
1095  fRhoRec[it][kZeroRho] = 0.0;
1096  }
1097 
1098  if(fTypeOfAnal == kEff || fTypeOfAnal == kEmb){ //rho in MC events
1099  fRhoMC[it][kConeRho] = EstimateBgCone( jetContGen, parContGen, trackTTGen[it], kTRUE); //container ID=1 mc particles
1100  fRhoMC[it][kCMSRho] = EstimateBgKTcms(jetContGenKT, parContGen, trackTTGen[it]);
1101  fRhoMC[it][kKtRho] = EstimateBgKT( jetContGenKT, parContGen, trackTTGen[it]);
1102  fRhoMC[it][kZeroRho] = 0.0;
1103  }
1104 
1105  if(fTypeOfAnal == kEmbSingl){ //embedding single track
1106  fRhoMC[it][kConeRho] = 0.0;
1107  fRhoMC[it][kCMSRho] = 0.0;
1108  fRhoMC[it][kKtRho] = 0.0;
1109  fRhoMC[it][kZeroRho] = 0.0;
1110  }
1111 
1112  if(fTypeOfAnal == kKine){ //rho in KINE MC events
1113  fRhoMC[it][kConeRho] = EstimateBgCone( jetContGen, parContGen, trackTTGen[it], kTRUE); //container ID=1 mc particles
1114  fRhoMC[it][kCMSRho] = EstimateBgKTcms(jetContGenKT, parContGen, trackTTGen[it]);
1115  fRhoMC[it][kKtRho] = EstimateBgKT( jetContGenKT, parContGen, trackTTGen[it]);
1116  fRhoMC[it][kZeroRho] = 0.0;
1117  }
1118  }
1119 
1120 
1121  //_________________________________________________________
1122  //EVALUATE SINGLE PARTICLE EFFICIENCY + FILL RESPONSE MATRIX
1123  Bool_t bRecPrim = kFALSE; //tags the reconstructed primary particles
1124 
1125  if(fTypeOfAnal == kEff){
1126 
1127  //1) FILL HISTOS FOR SINGLE PARTICLE EFFICIENCY
1128  if(parContGen){
1129  parContGen->ResetCurrentID();
1130  while((constTrackGen = (AliVParticle*)(parContGen->GetNextAcceptParticle()))){
1131  if(!constTrackGen) continue;
1132  if(IsTrackInAcceptance(constTrackGen, kTRUE)){
1133  //pT spectrum of generator level particles
1134  for(Int_t ic=0; ic<2; ic++){
1135  if(ficb[ic]==-1) continue;
1136  fhPtTrkTruePrimGen[ficb[ic]]->Fill(constTrackGen->Pt(),constTrackGen->Eta());
1137  }
1138  }
1139  }
1140 
1141  //single particle efficiency and contamination
1142 
1143  if(trkContRec && parContGen){
1144  trkContRec->ResetCurrentID();
1145  while((constTrackRec =(AliVParticle*) (trkContRec->GetNextAcceptParticle()))){
1146  if(!constTrackRec) continue;
1147  if(!IsTrackInAcceptance(constTrackRec, kFALSE)) continue; //reconstructed level tracks
1148  bRecPrim = kFALSE; //not yet matched to generator level physical primary
1149 
1150  parContGen->ResetCurrentID();
1151  while((constTrackGen = (AliVParticle*)(parContGen->GetNextAcceptParticle()))){
1152  if(!constTrackGen) continue;
1153  if(!IsTrackInAcceptance(constTrackGen, kTRUE)) continue; //gen level physical primary
1154  if(TMath::Abs(constTrackRec->GetLabel()) == TMath::Abs(constTrackGen->GetLabel())){
1155  //has the same label as reconstr track
1156 
1157  bRecPrim = kTRUE;
1158  for(Int_t ic=0; ic<2; ic++){
1159  if(ficb[ic]==-1) continue;
1160  fhPtTrkTruePrimRec[ficb[ic]]->Fill(constTrackGen->Pt(),constTrackGen->Eta()); //this is well recontr phys primary
1161  }
1162  fhDiffPtVsPtTrackTrue->Fill(constTrackGen->Pt(), constTrackRec->Pt()-constTrackGen->Pt());
1163 
1164  break;
1165  }//same label with rec particle
1166  }//loop over gen tracks
1167  if(!bRecPrim){
1168  for(Int_t ic=0; ic<2; ic++){
1169  if(ficb[ic]==-1) continue;
1170  fhPtTrkSecOrFakeRec[ficb[ic]]->Fill(constTrackRec->Pt(),constTrackRec->Eta()); //matchnig to phys primary not found, this is fake or second.
1171  }
1172  }
1173  }//loop over rec tracks
1174  }//rec track array exists
1175  }//gen particle array exists
1176 
1177  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1178  //2) FILL JET RESPONSE MATRIX
1179  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1180  Double_t ptGenCorr; //GEN jet pt corrected for rho
1181  Double_t ptRecCorr; //REC jet pt corrected for rho
1182 
1183  //Response matrix normalization - spectrum of all generator level jets in acceptance
1184  if(jetContGen){
1185  jetContGen->ResetCurrentID();
1186  while((jetGen = jetContGen->GetNextAcceptJet())){
1187  if(!jetGen) continue;
1188  if(!IsSignalJetInAcceptance(jetGen,kTRUE)) continue; //cuts on eta, pT ,area
1189 
1190  //if(fDebug>20) Printf("GEN JET phi=%f eta=%f pt=%f", jetGen->Phi(), jetGen->Eta(), jetGen->Pt());
1191 
1192  for(Int_t ir=0; ir< kRho; ir++){
1193  ptGenCorr = jetGen->Pt() - jetGen->Area()*fRhoMC[kRef][ir]; // correct for rho
1194 
1195  for(Int_t ic=0; ic<2; ic++){
1196  if(ficb[ic]==-1) continue;
1197  fhJetPtGen[ficb[ic]][ir]->Fill(ptGenCorr);
1198  }
1199  }
1200  }
1201  }
1202 
1203  //Find closest gen level+rec level jets
1204  //Get momentum shift due to fake tracks
1205  if(jetContRec){
1206  jetContRec->ResetCurrentID();
1207  while((jetRec = jetContRec->GetNextAcceptJet())) {
1208  if(!jetRec) continue;
1209  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue; //cuts on eta, pT ,area
1210 
1211 
1212  //Get momentum shift due to fake tracks
1213  sumFakeTrackPtInJetNotStrange = 0.;
1214  sumFakeTrackPtInJetStrange = 0.;
1215  sumStrangeTrackPtInJet = 0.;
1216  particleMC = NULL;
1217  particleMCMother = NULL;
1218  iIndexMother = 0;
1219  iPdgCode = 0;
1220  iPdgCodeMother = 0;
1221 
1222 
1223  for(Int_t iq=0; iq < jetRec->GetNumberOfTracks(); iq++) {
1224  constTrackRec = static_cast<AliVParticle*> (jetRec->TrackAt(iq,trkContRec->GetArray()));
1225  if(!constTrackRec) continue;
1226  bRecPrim = kFALSE; //not yet matched to generator level physical primary
1227 
1228  parContGen->ResetCurrentID();
1229  while((constTrackGen = (AliVParticle*)(parContGen->GetNextAcceptParticle()))){
1230  if(!constTrackGen) continue;
1231  if(!IsTrackInAcceptance(constTrackGen, kTRUE)) continue; //gen level physical primary
1232  if(TMath::Abs(constTrackRec->GetLabel()) == TMath::Abs(constTrackGen->GetLabel())){
1233  bRecPrim = kTRUE;
1234  break;
1235  }
1236  }
1237  if(!bRecPrim){ //this is a fake track
1238 
1239  lab = TMath::Abs(constTrackRec->GetLabel());
1240  bStrange = 0;
1241 
1242  if(lab < iNTracksMC){
1243  particleMC = (AliAODMCParticle*) arrayMC->At(lab);
1244  if(particleMC){
1245  iPdgCode = TMath::Abs(particleMC->GetPdgCode());
1246  iIndexMother = TMath::Abs(particleMC->GetMother());
1247  if(IsStrange(iPdgCode)){
1248  sumFakeTrackPtInJetStrange += constTrackRec->Pt();
1249  sumStrangeTrackPtInJet += constTrackRec->Pt();
1250  bStrange = 1;
1251  }
1252 
1253  if(iIndexMother < iNTracksMC && 0 <= iIndexMother && !bStrange){
1254  //check strangeness of mother of the secondary tracks
1255  particleMCMother = (AliAODMCParticle*) arrayMC->At(iIndexMother);
1256  if(particleMCMother){
1257  iPdgCodeMother = TMath::Abs(particleMCMother->GetPdgCode());
1258 
1259  if(IsStrange(iPdgCodeMother)){
1260  sumFakeTrackPtInJetStrange += constTrackRec->Pt();
1261  sumStrangeTrackPtInJet += constTrackRec->Pt();
1262  bStrange = 1;
1263  }
1264  }
1265  }
1266  }
1267  }
1268  if(!bStrange){
1269  sumFakeTrackPtInJetNotStrange += constTrackRec->Pt();
1270  }
1271  }else{
1272  //check strangeness of primary tracks
1273  lab = TMath::Abs(constTrackRec->GetLabel());
1274  if(lab < iNTracksMC){
1275  particleMC = (AliAODMCParticle*) arrayMC->At(lab);
1276  if(particleMC){
1277  iPdgCode = TMath::Abs(particleMC->GetPdgCode());
1278  if(IsStrange(iPdgCode)){
1279  sumStrangeTrackPtInJet += constTrackRec->Pt();
1280  }
1281  }
1282  }
1283  }
1284  }
1285 
1286  for(Int_t iw=0; iw<21; iw++){
1287  enhw = iw*0.05; //enhances weight of strangeness particles from 0 to 100% in steps of 5%
1288  newJetPt = jetRec->Pt() + enhw*sumStrangeTrackPtInJet; //jet contain 100 strg add only enhanc
1289  newSumFakePt = sumFakeTrackPtInJetNotStrange + (1.0+ enhw)*sumFakeTrackPtInJetStrange;
1290  fhPtJetPrimVsPtJetRec[iw]->Fill(newJetPt, newSumFakePt/newJetPt);
1291  }
1292 
1293  //if(fDebug>20) Printf("REC JET phi=%f eta=%f pt=%f",jetRec->Phi(), jetRec->Eta(), jetRec->Pt());
1294 
1295  //Find closest gen level+rec level jets
1296  jetGen = 0x0;
1297  jetGen = jetRec->ClosestJet();
1298 
1299  if(!jetGen){ //did not find matching generator level jet
1300  //if(fDebug>20) Printf("NO MATCH (NO SUCH GEN JET)");
1301 
1302  continue;
1303  }
1304  if(jetGen->Pt()<1e-3){
1305  //if(fDebug>20) Printf("SKIP MATCH WITH GHOST JET");
1306  continue;
1307  }
1308  //corresponding generator level jet found
1309  //if(fDebug>20) Printf("MATCHED WITH phi=%f eta=%f pt=%f",jetGen->Phi(), jetGen->Eta(), jetGen->Pt());
1310 
1311  if(!IsSignalJetInAcceptance(jetGen,kTRUE)) continue; //cuts on eta, pT ,area
1312 
1313  //check fraction of tracks from generator level jet in rec level jet
1314  //if(fDebug>20) Printf("FRACTIONH SHARED = %f ", GetFractionSharedPt(jetRec,jetContRec,jetGen, jetContGen));
1315 
1316  //FK// if(GetFractionSharedPt(jetRec, jetContRec, jetGen, jetContGen) < fMinFractionShared) continue;
1317 
1318  //if(fDebug>20) Printf("PASSED MIN FRACTION CRITERION ");
1319 
1320  for(Int_t ir=0; ir< kRho; ir++){
1321  ptGenCorr = jetGen->Pt() - jetGen->Area()*fRhoMC[kRef][ir]; //perp cone bg correction to pt
1322  ptRecCorr = jetRec->Pt() - jetRec->Area()*fRhoRec[kRef][ir]; //perp cone bg correction to pt
1323 
1324  for(Int_t ic=0; ic<2; ic++){
1325  if(ficb[ic]==-1) continue;
1326 
1327  fhJetPtGenVsJetPtRec[ficb[ic]][ir]->Fill(ptRecCorr, ptGenCorr); //response matrix
1328 
1329  if(ptGenCorr>0){
1330  fhJetPtResolutionVsPtGen[ficb[ic]][ir]->Fill(ptGenCorr,(ptRecCorr-ptGenCorr)/ptGenCorr); //jet pT resolution
1331  }
1332  }
1333  }
1334  }
1335  }//rec jet container exists
1336  }//analyze efficiency mode (response matrix + single particle efficiency)
1337 
1338 
1339  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1340  //if(bStop){
1341  // fHistEvtSelection->Fill(5); //Rejected events
1342  // return kTRUE; // SKIP EVENTS WHERE REF TT CLASS EVENT CONTAINS SIG TT TRACK
1343  // }
1344  fHistEvtSelection->Fill(6); //Accepted events
1345  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1346  Double_t areaJet, pTJet;
1347  Double_t dphi, dfi;
1348  Bool_t bFirstCycle = kTRUE;
1349 
1350 
1351  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1352  //H-JET CORRELATIONS IN MC TRUTH
1353  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1354  if(fTypeOfAnal == kEff || fTypeOfAnal == kKine){
1355  for(Int_t it=0; it< kTT; it++){
1356  if(parContGen && ntriggersGen[it] >0 && ((AliVParticle*)trackTTGen[it])){
1357  bFirstCycle = kTRUE;
1358 
1359  AliVParticle* triggerHadronGen = (AliVParticle*) trackTTGen[it]; //(fTrigTracksGen[it]);
1360  if(!triggerHadronGen) continue;
1361 
1362  for(Int_t ic=0; ic<2; ic++){
1363  if(ficb[ic]==-1) continue;
1364  fh1NtriggersGen[ficb[ic]][it]->Fill((Float_t) triggerHadronGen->Pt()); //trigger pT gen
1365 
1366  if( fTypeOfData == kHijing){ //impact parameter for triggered events
1367  fhImpactParameterTT[ficb[ic]][it]->Fill(fImpParam);
1368  }
1369  }
1370 
1371  //JET LOOP
1372  if(jetContGen){
1373  jetContGen->ResetCurrentID();
1374  while((jetGen = jetContGen->GetNextAcceptJet())) {
1375  if(!jetGen){
1376  AliError(Form("%s: Could not receive gen jet", GetName()));
1377  continue;
1378  }
1379  if(!IsSignalJetInAcceptance(jetGen,kTRUE)) continue;
1380 
1381  areaJet = jetGen->Area();
1382  pTJet = jetGen->Pt();
1383 
1384  if(bFirstCycle){
1385  for(Int_t ic=0; ic<2; ic++){
1386  if(ficb[ic]==-1) continue;
1387 
1388  fhJetPhiGen[ficb[ic]][it]->Fill( pTJet, RelativePhi(jetGen->Phi(),0.0));
1389  fhJetEtaGen[ficb[ic]][it]->Fill( pTJet, jetGen->Eta());
1390  }
1391  }
1392 
1393  dphi = RelativePhi(triggerHadronGen->Phi(), jetGen->Phi());
1394 
1395  dfi = dphi; //-0.5*pi to 1.5*Pi
1396  if(dfi < -0.5*TMath::Pi()) dfi += TMath::TwoPi();
1397  if(dfi > 1.5*TMath::Pi()) dfi -= TMath::TwoPi();
1398 
1399  for(Int_t ir=0; ir< kRho;ir++){
1400  for(Int_t ic=0; ic<2; ic++){
1401  if(ficb[ic]==-1) continue;
1402  fhDphiTriggerJetGen[ficb[ic]][it][ir]->Fill((Float_t) (pTJet - areaJet*fRhoMC[it][ir]), (Float_t) dfi); //Rongrong's analysis
1403  }
1404  }
1405  //-------------------------
1406 
1407  if(TMath::Abs(dphi) < fDphiCut) continue; //Dphi cut between trigger and assoc
1408 
1409  //Centrality, A, pTjet
1410  ftmpArray[0] = areaJet;
1411  for(Int_t ir=0; ir< kRho;ir++){
1412  ftmpArray[1] = pTJet - areaJet*fRhoMC[it][ir];
1413  for(Int_t ic=0; ic<2; ic++){
1414  if(ficb[ic]==-1) continue;
1415  fHJetSpecGen[ficb[ic]][it][ir]->Fill(ftmpArray[0],ftmpArray[1]);
1416 
1417  if(ir==0){
1418  fhJetPhiRecoilGen[ficb[ic]][it]->Fill( pTJet, RelativePhi(jetGen->Phi(),0.0));
1419  fhJetEtaRecoilGen[ficb[ic]][it]->Fill( pTJet, jetGen->Eta());
1420  }
1421 
1422  }
1423  }
1424  }//JET LOOP
1425  bFirstCycle = kFALSE;
1426  }//container exists
1427  }//TT exists
1428  }//TT loop
1429  }
1430 
1431  //++++++++++++++++++++++++++++++++++++++++++++++++++++
1432  if(fTypeOfAnal == kKine) return kTRUE;
1433  //++++++++++++++++++++++++++++++++++++++++++++++++++++
1434 
1435 
1436  for(Int_t it=0; it<kTT; it++){
1437  if((fTypeOfAnal==kEmb || fTypeOfAnal == kEmbSingl) && ((AliVParticle*) trackTT[it])){
1438  //delta pT using embedded pythia events
1439  //delta pT analyzed only in events with REAL EVENT TT present !!!!!!!!!!! (condition above)
1440  // PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskDeltaPtJEmb.cxx
1441  //AliJetResponseMaker
1442  //get deltaPt from embedding
1443  //++++++++++++++++++++++++++++++++++++++++
1444  //ANALYZE DELTA PT FOR ALL EMB MC JETS
1445  //++++++++++++++++++++++++++++++++++++++++
1446 
1447  AliEmcalJet *jetEmb = NULL;
1448  Bool_t bEmbJetCloseToTT = kFALSE;
1449 
1450  if(jetContRec){
1451  jetContRec->ResetCurrentID();
1452  while((jetRec = jetContRec->GetNextAcceptJet())) { //loop over reconstructed jets
1453  if(!jetRec) continue;
1454  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue; //apply cuts on eta, pT ,area
1455 
1456  //skip the jet that contains TT
1457  bEmbJetCloseToTT = kFALSE;
1458 
1459  if(jetRec->Pt() > trackTT[it]->Pt()*0.5){ // if jet contains TT it has to have at least pT of TT
1460 
1461  for(Int_t iq=0; iq < jetRec->GetNumberOfTracks(); iq++) {
1462  constTrackRec = static_cast<AliVParticle*> (jetRec->TrackAt(iq,trkContRec->GetArray())); //matched rec and emb tracks
1463  if(!constTrackRec) continue;
1464  if(constTrackRec != ((AliVParticle*) trackTT[it])) continue;
1465  //if(fDebug>21) Printf("EMB FIND TT COMPARE TRACK PT %f %f", constTrack->Pt(), triggerHadron->Pt());
1466  bEmbJetCloseToTT = kTRUE;
1467  break;
1468  }
1469  }
1470 
1471  if(bEmbJetCloseToTT) continue;//skip the jet that contains TT track
1472 
1473  jetEmb = NULL;
1474  jetEmb = jetRec->ClosestJet();
1475  if(!jetEmb){
1476  //if(fDebug>21) Printf("embedded jet does not exists, returning");
1477  continue;
1478  }
1479  if(jetEmb->Pt()<1e-3){
1480  //if(fDebug>21) Printf("SKIP MATCH WITH EMBEDDED GHOST JET");
1481  continue;
1482  }
1483  //if((fTypeOfAnal==kEmb) && (jetEmb->Pt()<6.0)) continue; // some hard cut on pT of emb jet ???
1484 
1485  if(!IsSignalJetInAcceptance(jetEmb,kTRUE)) continue; //apply cuts on eta, pT ,area on the embedded jet
1486 
1487  //Check fraction of tracks from generator level jet in rec level jet
1488  //At least 50% of embedded jet has to be in the reconstructed jet
1489  if(GetFractionSharedPt(jetRec, jetContRec, jetEmb, jetContGen) < fMinFractionShared) continue;
1490 
1491  for(Int_t ir=0; ir < kRho-1; ir++){
1492  for(Int_t ic=0; ic<2; ic++){
1493  if(ficb[ic]==-1) continue;
1494  //1 Dim distribution
1495  fhDeltaPtEmb[ficb[ic]][it][ir]->Fill(
1496  jetRec->Pt() - jetRec->Area() * fRhoRec[it][ir] - jetEmb->Pt());
1497 
1498  //2 Dim distribution
1499  fhDeltaPtEmb2D[ficb[ic]][it][ir]->Fill(jetEmb->Pt(),
1500  jetRec->Pt() - jetRec->Area() * fRhoRec[it][ir] - jetEmb->Pt());
1501 
1502 
1503  if(trackTT[it]){ //embedded track/jet is perp to TT
1504  dphi = TMath::Abs(RelativePhi(trackTT[it]->Phi(), jetEmb->Phi()));
1505  if(TMath::Pi()/4 < dphi && dphi < 3*TMath::Pi()/4){
1506  //1 Dim distribution
1507  fhDeltaPtEmbPerp[ficb[ic]][it][ir]->Fill(
1508  jetRec->Pt() - jetRec->Area() * fRhoRec[it][ir] - jetEmb->Pt());
1509 
1510  //2 Dim distribution
1511  fhDeltaPtEmbPerp2D[ficb[ic]][it][ir]->Fill(jetEmb->Pt(),
1512  jetRec->Pt() - jetRec->Area() * fRhoRec[it][ir] - jetEmb->Pt());
1513 
1514  }else if(dphi >= 3*TMath::Pi()/4 ){
1515  //embedded track back-to-back w.r.t the TT
1516  //1 Dim distribution
1517  fhDeltaPtEmbBc2Bc[ficb[ic]][it][ir]->Fill(
1518  jetRec->Pt() - jetRec->Area() * fRhoRec[it][ir] - jetEmb->Pt());
1519 
1520  //2 Dim distribution
1521  fhDeltaPtEmbBc2Bc2D[ficb[ic]][it][ir]->Fill(jetEmb->Pt(),
1522  jetRec->Pt() - jetRec->Area() * fRhoRec[it][ir] - jetEmb->Pt());
1523  }
1524  }
1525  }
1526  }
1527  }
1528  }
1529  }
1530  }//end of embedding
1531 
1532 
1533  // RECONSTRUCTED DATA ANALYSIS
1534 
1535  for(Int_t ir=0; ir < kRho-1; ir++){
1536  for(Int_t ic=0; ic<2; ic++){
1537  if(ficb[ic]==-1) continue;
1538 
1539  fhRhoIncl[ficb[ic]][ir]->Fill((Float_t) fRhoRec[kRef][ir]);
1540  }
1541  }
1542  for(Int_t it=0; it<kTT;it++){
1543  if(ntriggers[it]>0){
1544 
1545  for(Int_t ir=0; ir < kRho-1; ir++){
1546  //Estimate UE density in events with TT
1547  //Fill once per event
1548  for(Int_t ic=0; ic<2; ic++){
1549  if(ficb[ic]==-1) continue;
1550 
1551  fhRhoTT[ficb[ic]][it][ir]->Fill( (Float_t) fRhoRec[it][ir]);
1552  }
1553  }
1554  }
1555  }
1556  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1557  if(fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl) return kTRUE;
1558  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1559 
1560 
1561  // CALCULATE DELTA PT IN RECONSTRUCTED DATA WITH RANDOM CONES
1562  Double_t deltapt[kRho-1], phiTT = 0., etaTT = -1000.;
1563  //Double_t ncoll = -1.0;
1564  //if(centralityPercentile>=0.) ncoll = GetNcoll(centralityPercentile);
1565  Bool_t bRecJetCloseToTT = kFALSE;
1566  //Exclude region around TT from delta pt calculation
1567 
1568  for(Int_t it=0; it<kTT;it++){
1569  bRecJetCloseToTT = kFALSE;
1570  phiTT = 0.;
1571  etaTT = -1000.;
1572 
1573  if(trackTT[it]){ //get phi and eta of the TT or the reconstructed jet that contains TT
1574  phiTT = trackTT[it]->Phi(); // TT is proxy for jet
1575  etaTT = trackTT[it]->Eta();
1576 
1577  if(jetContRec){ //find jet that contains TT if it exists
1578  jetContRec->ResetCurrentID();
1579  while((jetRec = jetContRec->GetNextAcceptJet())) { //loop over reconstructed jets
1580  if(!jetRec) continue;
1581  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue;
1582 
1583  //skip the jet that contains TT
1584  bRecJetCloseToTT = kFALSE;
1585 
1586  if(jetRec->Pt() > trackTT[it]->Pt()*0.5){ // if jet contains TT it has to have at leat pT of TT
1587  for(Int_t iq=0; iq < jetRec->GetNumberOfTracks(); iq++){
1588  constTrackRec = (AliVParticle*) (jetRec->TrackAt(iq,trkContRec->GetArray())); //matched rec and emb tracks
1589  if(!constTrackRec) continue;
1590  if(constTrackRec != ((AliVParticle*) trackTT[it])) continue;
1591  phiTT = jetRec->Phi();
1592  etaTT = jetRec->Eta();
1593  bRecJetCloseToTT = kTRUE;
1594  break;
1595  }
1596  }
1597  if(bRecJetCloseToTT) break;
1598  }
1599  }
1600  }
1601 
1602  for(Int_t irc=0; irc<fNofRandomCones; irc++){
1603 
1604  //generate certain number of random cones per event
1605  GetDeltaPt(kRho-1, (fRhoRec[it]), &deltapt[0], phiTT, etaTT, trkContRec, kFALSE);// 1.0/ncoll);//FK//????? prob exlude RC ??? 1/Ncoll
1606 
1607 
1608  for(Int_t ir=0; ir < kRho-1; ir++){
1609  for(Int_t ic=0; ic<2; ic++){
1610  if(ficb[ic]==-1) continue;
1611  //fill delta pt histograms in inclusive events
1612  if(it==kRef) fhDeltaPtIncl[ficb[ic]][ir]->Fill(deltapt[ir]);
1613 
1614  if(ntriggers[it]>0){
1615  //fill delta pt histograms in events with TT (trigger track)
1616  fhDeltaPt[ficb[ic]][it][ir]->Fill( deltapt[ir]);
1617  }
1618  }
1619  }
1620  }
1621  }//trigger loop
1622  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1623  // Track Multiplicity, Track momentum smearing
1624  Double_t xyz[50];
1625  Double_t pxpypz[50];
1626  Double_t cv[21];
1627  Int_t itrkq;
1628  AliAODTrack *atrk=NULL;
1629 
1630  if(trkContRec){
1631 
1632  Int_t mult = 0;
1633  trkContRec->ResetCurrentID();
1634  while((atrk = (AliAODTrack*) (trkContRec->GetNextAcceptParticle()))){
1635  if(!atrk) continue;
1636  if(!IsTrackInAcceptance((AliVParticle*) atrk, kFALSE)) continue; //reconstructed level tracks
1637  mult++;
1638 
1639  dphi = atrk->Phi();
1640  if(dphi < 0) dphi+= 2*TMath::Pi();
1641  if(dphi > 2*TMath::Pi()) dphi-= 2*TMath::Pi();
1642 
1643  if(atrk->IsGlobalConstrained()){
1644  fhTrackPhiCG->Fill(atrk->Pt(), atrk->Phi()); //global constrained
1645  }else{
1646  fhTrackPhiTPCG->Fill(atrk->Pt(), atrk->Phi()); //complementary
1647  }
1648 
1649  itrkq = (atrk->Charge()<0) ? 0 : 1;
1650 
1651  fhInvPtQVsPhi[itrkq]->Fill(atrk->Phi(), 1.0/atrk->Pt());
1652  fhInvPtQVsEta[itrkq]->Fill(atrk->Eta(), 1.0/atrk->Pt());
1653 
1654  if(atrk->Eta()>0){
1655  fhInvPtQVsPhiASide[itrkq]->Fill(atrk->Phi(), 1.0/atrk->Pt());
1656  }else{
1657  fhInvPtQVsPhiCSide[itrkq]->Fill(atrk->Phi(), 1.0/atrk->Pt());
1658  }
1659 
1660  //get sigma pT / pT
1661  //Taken from AliEMCalTriggerExtraCuts::CalculateTPCTrackLength
1662  memset(cv, 0, sizeof(Double_t) * 21); //cleanup arrays
1663  memset(pxpypz, 0, sizeof(Double_t) * 50);
1664  memset(xyz, 0, sizeof(Double_t) * 50);
1665  atrk->GetXYZ(xyz);
1666  atrk->GetPxPyPz(pxpypz);
1667  atrk->GetCovarianceXYZPxPyPz(cv);
1668 
1669  AliExternalTrackParam par(xyz, pxpypz, cv, atrk->Charge());
1670  fhSigmaPtOverPtVsPt[itrkq]->Fill(atrk->Pt(), TMath::Abs(sqrt(par.GetSigma1Pt2())/par.GetSigned1Pt()));
1671 
1672  //XXX DCA distributions
1673  fhDCAinXVsPt->Fill(atrk->Pt(), atrk->XAtDCA());
1674  fhDCAinYVsPt->Fill(atrk->Pt(), atrk->YAtDCA());
1675 
1676  if(fTypeOfAnal == kEff){
1677  lab = TMath::Abs(atrk->GetLabel());
1678 
1679  if(lab < iNTracksMC){
1680  particleMC = (AliAODMCParticle*) arrayMC->At(lab);
1681  if(particleMC){
1682  if(particleMC->IsPhysicalPrimary()){
1683  fhDCAinXVsPtNonStrange->Fill(atrk->Pt(), atrk->XAtDCA());
1684  fhDCAinYVsPtNonStrange->Fill(atrk->Pt(), atrk->YAtDCA());
1685  }else{
1686  fhDCAinXVsPtStrange->Fill(atrk->Pt(), atrk->XAtDCA());
1687  fhDCAinYVsPtStrange->Fill(atrk->Pt(), atrk->YAtDCA());
1688  }
1689  }
1690  }
1691  }
1692  }
1693 
1694  for(Int_t ic=0; ic<2; ic++){
1695  if(ficb[ic]==-1) continue;
1696  //MB or CENT biases
1697  fhTrackMultiplicity[ficb[ic]]->Fill(mult);
1698 
1699  ftmpArrayX[0] = energyZdcNA;
1700  ftmpArrayX[1] = multVzero;
1701  ftmpArrayX[2] = mult;
1702  fhZNAVzeroATrack[ficb[ic]]->Fill(ftmpArrayX);
1703 
1704  for(Int_t it=0; it<kTT;it++){
1705  if(trackTT[it]){//TT biased + TT&&CENT biased distributions
1706  fhTrackMultiplicityTT[ficb[ic]][it]->Fill(mult);
1707  fhVzeroATotMultTT[ficb[ic]][it]->Fill(multVzero);
1708  fhZNAEnergyTT[ficb[ic]][it]->Fill(energyZdcNA);
1709  fhZNAVzeroATrackTT[ficb[ic]][it]->Fill(ftmpArrayX);
1710  }
1711  }
1712  }
1713 
1714 
1715  for(Int_t it=0; it<kTT;it++){
1716  if(trackTT[it]){ //TT biased distributions
1717  fhCentralityTT[it]->Fill((Float_t) centralityPercentile);
1718  }
1719  }
1720  }
1721  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1722  // INCLUSIVE JETS RECONSTRUCTED DATA
1723  if(jetContRec){
1724  jetContRec->ResetCurrentID();
1725  while((jetRec = jetContRec->GetNextAcceptJet())) {
1726  if(!jetRec){
1727  AliError(Form("%s: Could not receive jet", GetName()));
1728  continue;
1729  }
1730  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue;
1731 
1732  pTJet = jetRec->Pt();
1733 
1734  fhJetPhiIncl->Fill( pTJet, RelativePhi(jetRec->Phi(),0.0));
1735  fhJetEtaIncl->Fill( pTJet, jetRec->Eta());
1736  }
1737  }
1738  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1739  // H+JET IN RECONSTRUCTED DATA
1740 
1741  for(Int_t it=0; it<kTT;it++){
1742  if(ntriggers[it]>0 && ((AliVParticle*) trackTT[it])){
1743 
1744  bFirstCycle=kTRUE;
1745  if(trkContRec){
1746 
1747  AliVParticle* triggerHadron = (AliVParticle*) trackTT[it];// (fTrigTracks[it]);
1748  if(!triggerHadron) continue;
1749 
1750  for(Int_t ic=0; ic<2; ic++){
1751  if(ficb[ic]==-1) continue;
1752 
1753  fh1Ntriggers[ficb[ic]][it]->Fill((Float_t) triggerHadron->Pt()); //trigger p
1754  }
1755 
1756  //JET LOOP
1757  if(jetContRec){
1758  jetContRec->ResetCurrentID();
1759  while((jetRec = jetContRec->GetNextAcceptJet())) {
1760  if(!jetRec){
1761  AliError(Form("%s: Could not receive jet", GetName()));
1762  continue;
1763  }
1764  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue;
1765 
1766  areaJet = jetRec->Area();
1767  pTJet = jetRec->Pt();
1768 
1769  if(bFirstCycle){
1770  for(Int_t ic=0; ic<2; ic++){
1771  if(ficb[ic]==-1) continue;
1772  fhJetPhi[ficb[ic]][it]->Fill( pTJet, RelativePhi(jetRec->Phi(),0.0));
1773  fhJetEta[ficb[ic]][it]->Fill( pTJet, jetRec->Eta());
1774  }
1775  }
1776 
1777  dphi = RelativePhi(triggerHadron->Phi(), jetRec->Phi());
1778 
1779  dfi = dphi; //-0.5*pi to 1.5*Pi
1780  if(dfi < -0.5*TMath::Pi()) dfi += TMath::TwoPi();
1781  if(dfi > 1.5*TMath::Pi()) dfi -= TMath::TwoPi();
1782 
1783  for(Int_t ic=0; ic<2; ic++){
1784  if(ficb[ic]==-1) continue;
1785 
1786  for(Int_t ir=0; ir< kRho;ir++){
1787  fhDphiTriggerJet[ficb[ic]][it][ir]->Fill((Float_t) (pTJet - areaJet*fRhoRec[it][ir]), (Float_t) dfi); //Rongrong's analysis
1788  }
1789  }
1790  //-------------------------
1791 
1792  if(TMath::Abs(dphi) < fDphiCut) continue; //Dphi cut between trigger and assoc
1793  if(it==kRef) fhDphiTriggerJetAccept->Fill(dfi); //Accepted
1794 
1795  //Centrality, A, pTjet
1796  ftmpArray[0] = areaJet;
1797  for(Int_t ir=0; ir< kRho;ir++){
1798  ftmpArray[1] = pTJet - areaJet*fRhoRec[it][ir];
1799 
1800  for(Int_t ic=0; ic<2; ic++){
1801  if(ficb[ic]==-1) continue;
1802 
1803  fHJetSpec[ficb[ic]][it][ir]->Fill(ftmpArray[0],ftmpArray[1]);
1804 
1805  if(ir==0){
1806  fhJetPhiRecoil[ficb[ic]][it]->Fill( pTJet, RelativePhi(jetRec->Phi(),0.0));
1807  fhJetEtaRecoil[ficb[ic]][it]->Fill( pTJet, jetRec->Eta());
1808  }
1809 
1810  if(ir<kRho-1){
1811  fARhoTT[ficb[ic]][it][ir]->Fill((Float_t) (areaJet*fRhoRec[it][ir]));
1812  }
1813  }
1814  }
1815  }//JET LOOP
1816 
1817  bFirstCycle = kFALSE;
1818  }//container exists
1819  }
1820  }
1821  } //TT loop
1822 
1823  return kTRUE;
1824 }
1825 
1826 //________________________________________________________________________
1827 /*Double_t AliAnalysisTaskHJetSpectra::GetNcoll(Double_t centr){
1828  //Get Ncoll for given centrality
1829  if(fCollisionSystem == kpPb){
1830  //convert centrality percentle to Ncoll (2014-Sep-26-analysis_note-AnalysisNoteDijetspPb.pdf table 1)
1831  // PLATI JEN PRO pA !!!!!!!!!! CO PP?
1832  if(centr < 0.0 || centr > 100.) return -1;
1833 
1834  if(centr < 5.0) return 14.7; //0-5%
1835  else if(centr < 10) return 13.0; //5-10%
1836  else if(centr < 20) return 11.7; //10-20%
1837  else if(centr < 40) return 9.38; //20-40%
1838  else if(centr < 60) return 6.49; //40-60%
1839  else if(centr < 80) return 3.96; //60-80%
1840  else return 1.52; //80-100%
1841  }
1842 
1843 
1844  return -1.; //pp, pbpb
1845 
1846 }*/
1847 //________________________________________________________________________
1848 void AliAnalysisTaskHJetSpectra::Terminate(Option_t *){
1849  //Treminate
1850  PostData(1, fOutput);
1851 
1852  // Mandatory
1853  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1)); // '1' refers to the output slot
1854  if(!fOutput) {
1855  printf("ERROR: Output list not available\n");
1856  return;
1857  }
1858 }
1859 
1860 //________________________________________________________________________
1861 AliAnalysisTaskHJetSpectra::~AliAnalysisTaskHJetSpectra(){
1862  // Destructor. Clean-up the output list, but not the histograms that are put inside
1863  // (the list is owner and will clean-up these histograms). Protect in PROOF case.
1864  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
1865  delete fOutput;
1866  }
1867  delete fRandom;
1868  delete fHelperClass;
1869 
1870 }
1871 
1872 //________________________________________________________________________
1873 void AliAnalysisTaskHJetSpectra::UserCreateOutputObjects(){
1874  // called once to create user defined output objects like histograms, plots etc.
1875  // and to put it on the output list.
1876  // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
1877  //fOutput TList defined in the mother class
1878  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
1879 
1880  //Label of background
1881  TString bgtype[]={"Perp","CMS","KT","Zero"};
1882 
1883  fRandom = new TRandom3(0);
1884 
1885  Bool_t oldStatus = TH1::AddDirectoryStatus();
1886  TH1::AddDirectory(kFALSE);
1887  TString name;
1888 
1889  Bool_t bHistRec = (fTypeOfAnal != kEmb && fTypeOfAnal != kEmbSingl && fTypeOfAnal != kKine);
1890  Bool_t bNotKine = (fTypeOfAnal != kKine);
1891 
1892  Int_t icmax = (fTypeOfData == kPythia || fTypeOfData == kDmpjet ||fTypeOfAnal == kKine) ? 1 : kCAll;
1893  //__________________________________________________________
1894  // Event statistics
1895  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 7, -0.5, 6.5);
1896  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
1897  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
1898  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"pile up (rejected)");
1899  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
1900  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
1901  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"sig TT (rejected)");
1902  fHistEvtSelection->GetXaxis()->SetBinLabel(7,"sig TT (accepted)");
1903 
1904  fOutput->Add(fHistEvtSelection);
1905  //___________________________________________________________
1906  // Hard trigger counter
1907  for(Int_t it =0; it<kTT; it++){
1908  for(Int_t ic =0; ic<icmax; ic++){
1909  name = (ic==0) ? Form("fh1NtriggersMB") :
1910  Form("fh1Ntriggers%d%d",TMath::Nint(fCentralityBins[ic-1]),TMath::Nint(fCentralityBins[ic]));
1911  name.Append(Form("TT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
1912  fh1Ntriggers[ic][it] = new TH1D(name.Data(),"# of triggers",50,0.0,50.0);
1913  if(bHistRec) fOutput->Add((TH1D*)fh1Ntriggers[ic][it]);
1914 
1915  name = (ic==0) ? Form("fh1TriggerMultMB") :
1916  Form("fh1TriggerMult%d%d",TMath::Nint(fCentralityBins[ic-1]),TMath::Nint(fCentralityBins[ic]));
1917  name.Append(Form("TT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
1918  fh1TriggerMult[ic][it] = new TH1D(name.Data(),"# of triggers",50,0.0,50.0);
1919  if(bHistRec) fOutput->Add((TH1D*)fh1TriggerMult[ic][it]);
1920  }
1921  }
1922  //___________________________________________________________
1923  // trigger associated jet spectra (jet pT not corrected for UE)
1924  Int_t bw = (fUseDoubleBinPrecision==0) ? 1 : 2; //make larger bin width
1925 
1926  //jet associated to given TT
1927  //A, pTjet
1928  const Int_t dimSpec = 2;
1929  const Int_t nBinsSpec[dimSpec] = { 50, bw*160};
1930  const Double_t lowBinSpec[dimSpec] = { 0.0, -20.0};
1931  const Double_t hiBinSpec[dimSpec] = { 2.0, 300.0};
1932 
1933  for(Int_t ic =0; ic<icmax; ic++){
1934  for(Int_t it =0; it<kTT; it++){
1935  for(Int_t ir=0; ir< kRho; ir++){
1936  name = (ic==0) ? "fHJetSpecMB" :
1937  Form("fHJetSpec%d%d",TMath::Nint(fCentralityBins[ic-1]),
1938  TMath::Nint(fCentralityBins[ic]));
1939  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
1940 
1941  fHJetSpec[ic][it][ir] = new TH2D(
1942  name.Data(),
1943  Form("Recoil jet spectrum TT%d%d [A,pTjet-A*rho%s]",
1944  TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()),
1945  nBinsSpec[0], lowBinSpec[0], hiBinSpec[0], nBinsSpec[1], lowBinSpec[1],hiBinSpec[1]);
1946  if(bHistRec) fOutput->Add((TH2D*) fHJetSpec[ic][it][ir]);
1947  }
1948  }
1949  }
1950  //____________________________________________________________________
1951  //UE from cell median [Centrality, rho, pTUe ]
1952 
1953  for(Int_t ic =0; ic<icmax; ic++){
1954  for(Int_t ir=0; ir< kRho-1; ir++){ //Skip Zero bg
1955  for(Int_t it=0; it<kTT;it++){
1956  name = (ic==0) ? "fhRhoMB" :
1957  Form("fhRho%d%d",TMath::Nint(fCentralityBins[ic-1]),
1958  TMath::Nint(fCentralityBins[ic]));
1959 
1960  name.Append(Form("TT%d%d%s",TMath::Nint(fTTlow[it]),TMath::Nint(fTThigh[it]),bgtype[ir].Data()));
1961 
1962  //rho in events with TT
1963  fhRhoTT[ic][it][ir] = new TH1F(name.Data(),
1964  Form("Rho%s",bgtype[ir].Data()),80, 0.0, 40.0);
1965  if(bNotKine) fOutput->Add((TH1F*) fhRhoTT[ic][it][ir]);
1966 
1967 
1968  // rho times area in events with TT
1969  name = (ic==0) ? "fARhoMB" :
1970  Form("fARho%d%d",TMath::Nint(fCentralityBins[ic-1]),
1971  TMath::Nint(fCentralityBins[ic]));
1972 
1973  name.Append(Form("TT%d%d%s",TMath::Nint(fTTlow[it]),TMath::Nint(fTThigh[it]),bgtype[ir].Data()));
1974  fARhoTT[ic][it][ir] = new TH1F(name.Data(),
1975  Form("Area times rho %s",bgtype[ir].Data()),80, 0.0, 40.0);
1976  if(bHistRec) fOutput->Add((TH1F*) fARhoTT[ic][it][ir]);
1977  }
1978  //rho in inclusive events
1979  name = (ic==0) ? Form("fhRhoInclMB%s",bgtype[ir].Data()) :
1980  Form("fhRhoIncl%d%d%s",TMath::Nint(fCentralityBins[ic-1]),
1981  TMath::Nint(fCentralityBins[ic]),bgtype[ir].Data());
1982 
1983  fhRhoIncl[ic][ir] = (TH1F*) fhRhoTT[ic][kRef][ir]->Clone(name.Data());
1984  if(bNotKine) fOutput->Add((TH1F*) fhRhoIncl[ic][ir]);
1985  }
1986  }
1987  //_______________________________________________________________________
1988  // Delta pt distributions
1989  for(Int_t it =0; it< kTT; it++){
1990  for(Int_t ic =0; ic<icmax; ic++){
1991  for(Int_t ir=0; ir< kRho-1; ir++){
1992  //events with TT tracks
1993  name = (ic==0) ? "fhDeltaPtMB" :
1994  Form("fhDeltaPt%d%d",TMath::Nint(fCentralityBins[ic-1]),
1995  TMath::Nint(fCentralityBins[ic]));
1996 
1997  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
1998  fhDeltaPt[ic][it][ir] = new TH1D(name.Data(),
1999  Form("DeltaPt%s",bgtype[ir].Data()), 150, -50, 100);
2000  if(bHistRec) fOutput->Add((TH1D*) fhDeltaPt[ic][it][ir]);
2001 
2002  if(it==kRef){
2003  //inclusive events
2004  name = (ic==0) ? "fhDeltaPtInclMB" :
2005  Form("fhDeltaPtIncl%d%d",TMath::Nint(fCentralityBins[ic-1]),
2006  TMath::Nint(fCentralityBins[ic]));
2007 
2008  name.Append(Form("%s", bgtype[ir].Data()));
2009  fhDeltaPtIncl[ic][ir] = (TH1D*) fhDeltaPt[ic][it][ir]->Clone(name.Data());
2010  if(bHistRec) fOutput->Add((TH1D*) fhDeltaPtIncl[ic][ir]);
2011  }
2012 
2013  if(fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
2014  //Embedded PYTHIA jets
2015  //1D
2016  name = (ic==0) ? "fhDeltaPtEmbMB" :
2017  Form("fhDeltaPtEmb%d%d",TMath::Nint(fCentralityBins[ic-1]),
2018  TMath::Nint(fCentralityBins[ic]));
2019 
2020  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
2021  fhDeltaPtEmb[ic][it][ir] = (TH1D*) fhDeltaPt[ic][it][ir]->Clone(name.Data());
2022  fOutput->Add((TH1D*) fhDeltaPtEmb[ic][it][ir]);
2023 
2024  //1D perp to TT
2025  name = (ic==0) ? "fhDeltaPtEmbPerpMB" :
2026  Form("fhDeltaPtEmbPerp%d%d",TMath::Nint(fCentralityBins[ic-1]),
2027  TMath::Nint(fCentralityBins[ic]));
2028 
2029  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
2030  fhDeltaPtEmbPerp[ic][it][ir] = (TH1D*) fhDeltaPt[ic][it][ir]->Clone(name.Data());
2031  fOutput->Add((TH1D*) fhDeltaPtEmbPerp[ic][it][ir]);
2032 
2033  //1D back-to back to TT
2034  name = (ic==0) ? "fhDeltaPtEmbBc2BcMB" :
2035  Form("fhDeltaPtEmbBc2Bc%d%d",TMath::Nint(fCentralityBins[ic-1]),
2036  TMath::Nint(fCentralityBins[ic]));
2037 
2038  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
2039  fhDeltaPtEmbBc2Bc[ic][it][ir]= (TH1D*) fhDeltaPt[ic][it][ir]->Clone(name.Data());
2040  fOutput->Add((TH1D*) fhDeltaPtEmbBc2Bc[ic][it][ir]);
2041 
2042  //2D
2043  name = (ic==0) ? "fhDeltaPtEmb2DMB" :
2044  Form("fhDeltaPtEmb2D%d%d",TMath::Nint(fCentralityBins[ic-1]),
2045  TMath::Nint(fCentralityBins[ic]));
2046 
2047  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
2048  fhDeltaPtEmb2D[ic][it][ir] = new TH2D(name.Data(),
2049  Form("fhDeltaPtEmb2D%s",bgtype[ir].Data()), 125,0, 250, 150, -50, 100);
2050  fOutput->Add((TH2D*)fhDeltaPtEmb2D[ic][it][ir]);
2051 
2052 
2053  //2D perp to TT
2054  name = (ic==0) ? "fhDeltaPtEmbPerp2DMB" :
2055  Form("fhDeltaPtEmbPerp2D%d%d",TMath::Nint(fCentralityBins[ic-1]),
2056  TMath::Nint(fCentralityBins[ic]));
2057 
2058  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
2059  fhDeltaPtEmbPerp2D[ic][it][ir] = (TH2D*) fhDeltaPtEmb2D[ic][it][ir]->Clone(name.Data());
2060  fOutput->Add((TH2D*)fhDeltaPtEmbPerp2D[ic][it][ir]);
2061 
2062 
2063  //2D back-to-back p to TT
2064  name = (ic==0) ? "fhDeltaPtEmbBc2Bc2DMB" :
2065  Form("fhDeltaPtEmbBc2Bc2D%d%d",TMath::Nint(fCentralityBins[ic-1]),
2066  TMath::Nint(fCentralityBins[ic]));
2067 
2068  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
2069  fhDeltaPtEmbBc2Bc2D[ic][it][ir]=(TH2D*) fhDeltaPtEmb2D[ic][it][ir]->Clone(name.Data());
2070  fOutput->Add((TH2D*)fhDeltaPtEmbBc2Bc2D[ic][it][ir]);
2071 
2072  }
2073  }
2074  }
2075  }
2076 
2077  //_______________________________________________________________________
2078  //KT jets area versus PT
2079  Double_t binsPt [] = {0, 0.05,0.15, 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6.,6.5, 7., 8., 9., 10.,15., 20., 50., 100.};
2080  Int_t nbinsPt = sizeof(binsPt)/sizeof(Double_t)-1;
2081 
2082  fhKTAreaPt = new TH2F("fhKTAreaPt","KT jet Area vs Pt",nbinsPt, binsPt, 50,0,2);
2083  if(fTypeOfAnal == kRec) fOutput->Add(fhKTAreaPt);
2084 
2085  //_______________________________________________________________________
2086  //inclusive azimuthal and pseudorapidity histograms
2087  for(Int_t ic =0; ic<icmax; ic++){
2088  for(Int_t it=0; it<kTT; it++){
2089  name = (ic==0) ? Form("fhJetPhiMBTT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])) :
2090  Form("fhJetPhi%d%dTT%d%d",TMath::Nint(fCentralityBins[ic-1]),
2091  TMath::Nint(fCentralityBins[ic]),
2092  TMath::Nint(fTTlow[it]),
2093  TMath::Nint(fTThigh[it]));
2094 
2095 
2096  fhJetPhi[ic][it] = new TH2F(name.Data(),"Azim dist jets vs pTjet", 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
2097  if(bHistRec) fOutput->Add((TH2F*)fhJetPhi[ic][it]);
2098  }
2099  //-------------------------
2100  for(Int_t it=0; it<kTT; it++){
2101  name = (ic==0) ? Form("fhJetPhiMBTT%d%dRecoil", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])) :
2102  Form("fhJetPhi%d%dTT%d%dRecoil",TMath::Nint(fCentralityBins[ic-1]),
2103  TMath::Nint(fCentralityBins[ic]),
2104  TMath::Nint(fTTlow[it]),
2105  TMath::Nint(fTThigh[it]));
2106 
2107 
2108  fhJetPhiRecoil[ic][it] = new TH2F(name.Data(),"Azim dist jets vs pTjet", 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
2109  if(bHistRec) fOutput->Add((TH2F*)fhJetPhiRecoil[ic][it]);
2110  }
2111  //-------------------------
2112 
2113  name = (ic==0) ? Form("fhTrackPhiMB") :
2114  Form("fhTrackPhi%d%d",TMath::Nint(fCentralityBins[ic-1]),
2115  TMath::Nint(fCentralityBins[ic]));
2116 
2117  fhTrackPhi[ic] = new TH2F(name.Data(),"azim dist trig had vs pT,trk", 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
2118  if(bNotKine) fOutput->Add((TH2F*)fhTrackPhi[ic]);
2119  //-------------------------
2120  for(Int_t it=0; it<kTT; it++){
2121  name = (ic==0) ? Form("fhJetEtaMBTT%d%d",TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]) ) :
2122  Form("fhJetEta%d%dTT%d%d",TMath::Nint(fCentralityBins[ic-1]),
2123  TMath::Nint(fCentralityBins[ic]),
2124  TMath::Nint(fTTlow[it]),
2125  TMath::Nint(fTThigh[it]));
2126 
2127  fhJetEta[ic][it] = new TH2F(name.Data(),"Eta dist jets vs pTjet", 50,0, 100, 40,-0.9,0.9);
2128  if(bHistRec) fOutput->Add((TH2F*)fhJetEta[ic][it]);
2129  }
2130  //-------------------------
2131  for(Int_t it=0; it<kTT; it++){
2132  name = (ic==0) ? Form("fhJetEtaMBTT%d%dRecoil",TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]) ) :
2133  Form("fhJetEta%d%dTT%d%dRecoil",TMath::Nint(fCentralityBins[ic-1]),
2134  TMath::Nint(fCentralityBins[ic]),
2135  TMath::Nint(fTTlow[it]),
2136  TMath::Nint(fTThigh[it]));
2137 
2138  fhJetEtaRecoil[ic][it] = new TH2F(name.Data(),"Eta dist of recoil jets vs pTjet", 50,0, 100, 40,-0.9,0.9);
2139  if(bHistRec) fOutput->Add((TH2F*)fhJetEtaRecoil[ic][it]);
2140  }
2141  //-------------------------
2142  name = (ic==0) ? Form("fhTrackEtaMB") :
2143  Form("fhTrackEta%d%d",TMath::Nint(fCentralityBins[ic-1]),
2144  TMath::Nint(fCentralityBins[ic]));
2145 
2146  fhTrackEta[ic] = new TH2F(name.Data(),"Eta dist trig had vs pT,trk", 50, 0, 50, 40,-0.9,0.9);
2147  if(bNotKine) fOutput->Add((TH2F*) fhTrackEta[ic]);
2148  //-------------------------
2149  name = (ic==0) ? Form("fhTrackPtMB") :
2150  Form("fhTrackPt%d%d",TMath::Nint(fCentralityBins[ic-1]),
2151  TMath::Nint(fCentralityBins[ic]));
2152 
2153  fhTrackPt[ic] = new TH1F(name.Data(),"pT,trk ", 50, 0, 50);
2154  if(bNotKine) fOutput->Add((TH1F*) fhTrackPt[ic]);
2155  }
2156  //-------------------------
2157  fhVertexZ = new TH1F("fhVertexZ","z vertex",40,-20,20);
2158  if(bNotKine) fOutput->Add(fhVertexZ);
2159  //-------------------------
2160  fhVertexXAccept = new TH1F("fhVertexXAccept","vertex after cut",600,-3,3);
2161  if(bNotKine) fOutput->Add(fhVertexXAccept);
2162 
2163  fhVertexYAccept = (TH1F*) fhVertexXAccept->Clone("fhVertexYAccept");
2164  if(bNotKine) fOutput->Add(fhVertexYAccept);
2165 
2166  fhVertexZAccept = new TH1F("fhVertexZAccept","z vertex after cut",40,-20,20);
2167  if(bNotKine) fOutput->Add(fhVertexZAccept);
2168 
2169  fhVertexXAcceptTT = (TH1F*) fhVertexXAccept->Clone("fhVertexXAcceptTT");
2170  if(bNotKine) fOutput->Add(fhVertexXAcceptTT);
2171 
2172  fhVertexYAcceptTT = (TH1F*) fhVertexXAccept->Clone("fhVertexYAcceptTT");
2173  if(bNotKine) fOutput->Add(fhVertexYAcceptTT);
2174 
2175  fhVertexZAcceptTT = (TH1F*) fhVertexZAccept->Clone("fhVertexZAcceptTT");
2176  if(bNotKine) fOutput->Add(fhVertexZAcceptTT);
2177 
2178  //-------------------------
2179 
2180  if(fTypeOfData!=kReal){
2181  fhVertexZMC = new TH1F("fhVertexZMC","z vertex",40,-20,20);
2182  fOutput->Add(fhVertexZMC);
2183  //-------------------------
2184  fhVertexZAcceptMC = new TH1F("fhVertexZAcceptMC","z vertex after cut",40,-20,20);
2185  fOutput->Add(fhVertexZAcceptMC);
2186  }
2187  //-------------------------
2188  for(Int_t ic =0; ic<icmax; ic++){
2189  for(Int_t it =0; it<kTT; it++){
2190  for(Int_t ir=0; ir < kRho; ir++){ //Rongrong's analysis
2191  name = (ic==0) ?
2192  Form("fhDphiTriggerJetMB") :
2193  Form("fhDphiTriggerJet%d%d",TMath::Nint(fCentralityBins[ic-1]), TMath::Nint(fCentralityBins[ic]));
2194 
2195  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]) ,bgtype[ir].Data()));
2196 
2197  fhDphiTriggerJet[ic][it][ir] = new TH2F(name.Data(),"Deltaphi trig-jet",75,-50,100, 100, -0.5*TMath::Pi(),1.5*TMath::Pi());
2198  if(bHistRec) fOutput->Add((TH2F*) fhDphiTriggerJet[ic][it][ir]);
2199  }
2200  }
2201  }
2202  //-------------------------
2203 
2204  fhDphiTriggerJetAccept = new TH1F("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -0.5*TMath::Pi(),1.5*TMath::Pi());
2205  if(bHistRec) fOutput->Add(fhDphiTriggerJetAccept);
2206 
2207  fhDphiTTTT[kRef] = new TH1F("fhDphiTTTT67","Deltaphi between multiple TT",50, -0.5*TMath::Pi(),1.5*TMath::Pi());
2208  if(bHistRec) fOutput->Add((TH1F*) fhDphiTTTT[kRef]);
2209 
2210  fhDphiTTTT[kSig] = (TH1F*) fhDphiTTTT[kRef]->Clone("fhDphiTTTT1250");
2211  if(bHistRec) fOutput->Add((TH1F*) fhDphiTTTT[kSig]);
2212 
2213  fhJetPhiIncl = new TH2F("fhJetPhiIncl","Azim dist jets vs pTjet", 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
2214  if(bHistRec) fOutput->Add((TH2F*) fhJetPhiIncl);
2215 
2216  fhJetEtaIncl = new TH2F("fhJetEtaIncl","Eta dist inclusive jets vs pTjet", 50,0, 100, 40,-0.9,0.9);
2217  if(bHistRec) fOutput->Add((TH2F*) fhJetEtaIncl);
2218 
2219  fhTrackPhiCG = new TH2F("fhTrackPhiCG","azim dist trig had vs pT,trk Glob Const", 50, 0, 50, 50,0,2*TMath::Pi());
2220  if(bNotKine) fOutput->Add((TH2F*) fhTrackPhiCG);
2221 
2222  fhTrackPhiTPCG = new TH2F("fhTrackPhiTPCG","azim dist trig had vs pT,trk TPC Const", 50, 0, 50, 50,0,2*TMath::Pi());
2223  if(bNotKine) fOutput->Add((TH2F*) fhTrackPhiTPCG);
2224 
2225 
2226  for(Int_t i=0; i<2; i++){
2227  name = (i==0) ? "Neg" : "Pos";
2228  fhInvPtQVsPhi[i] = new TH2D(Form("fhInvPtVsPhiQ%s", name.Data()),
2229  Form("%s track 1/pt versus track phi", name.Data()), 36, 0, 2*TMath::Pi(), 40, 0, 0.4);
2230  if(bNotKine) fOutput->Add((TH2D*) fhInvPtQVsPhi[i]);
2231 
2232  fhInvPtQVsEta[i] = new TH2D(Form("fhInvPtVsEtaQ%s", name.Data()),
2233  Form("%s track 1/pt versus track eta", name.Data()), 20, -0.9, 0.9, 40, 0, 0.4);
2234  if(bNotKine) fOutput->Add((TH2D*) fhInvPtQVsEta[i]);
2235 
2236  fhInvPtQVsPhiASide[i] = new TH2D(Form("fhInvPtVsPhiASideQ%s", name.Data()),
2237  Form("%s track 1/pt versus track phi", name.Data()), 36, 0, 2*TMath::Pi(), 40, 0, 0.4);
2238  if(bNotKine) fOutput->Add((TH2D*) fhInvPtQVsPhiASide[i]);
2239 
2240  fhInvPtQVsPhiCSide[i] = new TH2D(Form("fhInvPtVsPhiCSideQ%s", name.Data()),
2241  Form("%s track 1/pt versus track phi", name.Data()), 36, 0, 2*TMath::Pi(), 40, 0, 0.4);
2242  if(bNotKine) fOutput->Add((TH2D*) fhInvPtQVsPhiCSide[i]);
2243 
2244  fhSigmaPtOverPtVsPt[i] = new TH2D(Form("fhSigmaPtOverPtVsPtQ%s", name.Data()),
2245  Form("%s track sigma(1/pt)/ 1/pt vs pt", name.Data()), 100, 0, 100, 250, 0, 1);
2246  if(bNotKine) fOutput->Add((TH2D*) fhSigmaPtOverPtVsPt[i]);
2247  }
2248  //-------------------------
2249  for(Int_t ic =0; ic<icmax; ic++){
2250  name = (ic==0) ? Form("fhCentralityMB") :
2251  Form("fhCentrality%d%d",TMath::Nint(fCentralityBins[ic-1]),
2252  TMath::Nint(fCentralityBins[ic]));
2253 
2254  fhCentrality[ic] = new TH1F(name.Data(),"Centrality",100,0,100);
2255  if(bNotKine) fOutput->Add((TH1F*) fhCentrality[ic]);
2256  }
2257  //-------------------------
2258  fhCentralityV0M = new TH1F("hCentralityV0M","hCentralityV0M",100,0,100);
2259  if(bNotKine)fOutput->Add(fhCentralityV0M);
2260  //-------------------------
2261  fhCentralityV0A = new TH1F("hCentralityV0A","hCentralityV0A",100,0,100);
2262  if(bNotKine) fOutput->Add(fhCentralityV0A);
2263  //-------------------------
2264  fhCentralityV0C = new TH1F("hCentralityV0C","hCentralityV0C",100,0,100);
2265  if(bNotKine) fOutput->Add(fhCentralityV0C);
2266  //-------------------------
2267  fhCentralityZNA = new TH1F("hCentralityZNA","hCentralityZNA",100,0,100);
2268  if(bNotKine) fOutput->Add(fhCentralityZNA);
2269  //-------------------------
2270  for(Int_t it =0; it<kTT;it++){
2271  name = Form("fhCentralityTT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]));
2272  fhCentralityTT[it] = (TH1F*) fhCentralityV0M->Clone(name.Data());
2273  if(bHistRec) fOutput->Add(fhCentralityTT[it]);
2274  }
2275  //-----------------------------------------------------
2276 
2277  for(Int_t ic =0; ic<icmax; ic++){
2278  name = (ic==0) ? Form("fhVzeroATotMultMB") :
2279  Form("fhVzeroATotMult%d%d",TMath::Nint(fCentralityBins[ic-1]),
2280  TMath::Nint(fCentralityBins[ic]));
2281 
2282  // vzero multiplicity
2283  fhVzeroATotMult[ic] = new TH1F(name.Data(),"hVzeroATotMult",1000,0,1000);
2284  if(bHistRec) fOutput->Add((TH1F*) fhVzeroATotMult[ic]);
2285 
2286  for(Int_t it=0; it<kTT; it++){
2287  name = (ic==0) ? "fhVzeroATotMultMB" :
2288  Form("fhVzeroATotMult%d%d",TMath::Nint(fCentralityBins[ic-1]),
2289  TMath::Nint(fCentralityBins[ic]));
2290 
2291  name.Append(Form("TT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
2292  fhVzeroATotMultTT[ic][it] = (TH1F*) fhVzeroATotMult[ic]->Clone(name.Data());
2293  if(bHistRec) fOutput->Add((TH1F*) fhVzeroATotMultTT[ic][it]);
2294  }
2295  }
2296 
2297 
2298  //-----------------------------------------------------
2299  // ZDC ZNA energy
2300 
2301  for(Int_t ic =0; ic<icmax; ic++){
2302  name = (ic==0) ? Form("fhZNAEnergyMB") :
2303  Form("fhZNAEnergy%d%d",TMath::Nint(fCentralityBins[ic-1]),
2304  TMath::Nint(fCentralityBins[ic]));
2305 
2306 
2307  fhZNAEnergy[ic] = new TH1F(name.Data(),"fhZNAEnergy",1000,0,200);
2308  if(bHistRec) fOutput->Add((TH1F*)fhZNAEnergy[ic]);
2309 
2310  for(Int_t it=0; it<kTT; it++){
2311  name = (ic==0) ? Form("fhZNAEnergyMB") :
2312  Form("fhZNAEnergy%d%d",TMath::Nint(fCentralityBins[ic-1]),
2313  TMath::Nint(fCentralityBins[ic]));
2314 
2315  name.Append(Form("TT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
2316  fhZNAEnergyTT[ic][it] = (TH1F*) fhZNAEnergy[ic]->Clone(name.Data());
2317  if(bHistRec) fOutput->Add((TH1F*) fhZNAEnergyTT[ic][it]);
2318  }
2319  }
2320 
2321  //-----------------------------------------------------
2322  // track multiplicity
2323 
2324  for(Int_t ic =0; ic<icmax; ic++){
2325  name = (ic==0) ? Form("fhTrackMultiplicityMB") :
2326  Form("fhTrackMultiplicity%d%d",TMath::Nint(fCentralityBins[ic-1]),
2327  TMath::Nint(fCentralityBins[ic]));
2328 
2329  fhTrackMultiplicity[ic] = new TH1D(name.Data(),"fhTrackMultiplicity",1000,0,1000);
2330  if(bNotKine) fOutput->Add((TH1D*)fhTrackMultiplicity[ic]);
2331 
2332  for(Int_t it=0; it<kTT; it++){
2333  name = (ic==0) ? Form("fhTrackMultiplicityMB") :
2334  Form("fhTrackMultiplicity%d%d",TMath::Nint(fCentralityBins[ic-1]),
2335  TMath::Nint(fCentralityBins[ic]));
2336 
2337  name.Append(Form("TT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
2338  fhTrackMultiplicityTT[ic][it] = (TH1D*)fhTrackMultiplicity[ic]->Clone(name.Data());
2339  if(bNotKine) fOutput->Add(fhTrackMultiplicityTT[ic][it]);
2340  }
2341  }
2342 
2343  //-----------------------------------------------------
2344  // ZNA energy versus Vzero mult. versus track mult. in all events
2345 
2346  const Int_t dimZ = 3;
2347  const Int_t nBinsZ[dimZ] = { 100, 100, 25 };
2348  const Double_t lowBinZ[dimZ] = { 0.0, 0.0, 0.0};
2349  const Double_t hiBinZ[dimZ] = { 200.0, 1000.0, 250};
2350 
2351  for(Int_t ic =0; ic<icmax; ic++){
2352  name = (ic==0) ? Form("fhZNAVzeroATrackMB") :
2353  Form("fhZNAVzeroATrack%d%d",TMath::Nint(fCentralityBins[ic-1]),
2354  TMath::Nint(fCentralityBins[ic]));
2355 
2356  fhZNAVzeroATrack[ic] = new THnSparseF(
2357  name.Data(),
2358  Form("ZNA, V0A mult, track mult"),
2359  dimZ, nBinsZ,lowBinZ,hiBinZ);
2360  if(fTypeOfAnal == kRec) fOutput->Add((THnSparseF*) fhZNAVzeroATrack[ic]);
2361 
2362 
2363  //-------------
2364  for(Int_t it=0; it<kTT; it++){
2365  // ZNA energy versus Vzero mult. versus track mult. in events with TT
2366  name = (ic==0) ? Form("fhZNAVzeroATrackMB") :
2367  Form("fhZNAVzeroATrack%d%d",TMath::Nint(fCentralityBins[ic-1]),
2368  TMath::Nint(fCentralityBins[ic]));
2369 
2370  name.Append(Form("TT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
2371 
2372  fhZNAVzeroATrackTT[ic][it] = (THnSparseF*) fhZNAVzeroATrack[ic]->Clone(name.Data());
2373  if(fTypeOfAnal == kRec) fOutput->Add((THnSparseF*) fhZNAVzeroATrackTT[ic][it]);
2374  }
2375  }
2376 
2377  //-----------------------------------------------------
2378 
2379  /*
2380  if(fTypeOfAnal == kKine){
2381  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
2382  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
2383  fOutput->Add(fh1Xsec);
2384 
2385  fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
2386  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
2387  fOutput->Add(fh1Trials);
2388 
2389  //fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",1000,0,1000);
2390  //fOutput->Add(fh1PtHard);
2391  }
2392 */
2393  if(fTypeOfData == kHijing){
2394  for(Int_t ic =0; ic<kCAll; ic++){
2395  name = (ic==0) ? Form("fhImpactParameterMB") :
2396  Form("fhImpactParameter%d%d",TMath::Nint(fCentralityBins[ic-1]),
2397  TMath::Nint(fCentralityBins[ic]));
2398 
2399  fhImpactParameter[ic] = new TH1D(name.Data(),"impact parameter distribution from HIJING",50,0,10);
2400  fOutput->Add((TH1D*) fhImpactParameter[ic]);
2401 
2402 
2403  for(Int_t it=0; it<kTT;it++){
2404  name = (ic==0) ? Form("fhImpactParameterMB") :
2405  Form("fhImpactParameter%d%d",TMath::Nint(fCentralityBins[ic-1]),
2406  TMath::Nint(fCentralityBins[ic]));
2407 
2408  name.Append(Form("TT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
2409  fhImpactParameterTT[ic][it] = new TH1D(name.Data(),"b versus TT",50,0,10);
2410  fOutput->Add((TH1D*) fhImpactParameterTT[ic][it]);
2411  }
2412  }
2413  }
2414 
2415  Double_t bins [] = {0, 0.2,0.4,0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.5, 3., 3.5, 4., 5., 6., 8., 10., 20., 50.};
2416  Int_t nbins = sizeof(bins)/sizeof(Double_t)-1;
2417 
2418  if(fTypeOfAnal== kEff ){
2419  for(Int_t ic =0; ic<icmax; ic++){
2420  for(Int_t ir=0; ir < kRho; ir++){
2421  name = (ic==0) ? Form("fhJetPtGenMB%s",bgtype[ir].Data()) :
2422  Form("fhJetPtGen%d%d%s",TMath::Nint(fCentralityBins[ic-1]),
2423  TMath::Nint(fCentralityBins[ic]),bgtype[ir].Data());
2424 
2425  fhJetPtGen[ic][ir] = new TH1D(name.Data(),
2426  Form("Jet pT Gen %s",bgtype[ir].Data()),bw*160,-20,300);
2427  fOutput->Add((TH1D*) fhJetPtGen[ic][ir]);
2428 
2429 
2430  name = (ic==0) ? Form("fhJetPtGenVsJetPtRecMB%s",bgtype[ir].Data()) :
2431  Form("fhJetPtGenVsJetPtRec%d%d%s",TMath::Nint(fCentralityBins[ic-1]),
2432  TMath::Nint(fCentralityBins[ic]),bgtype[ir].Data());
2433 
2434  fhJetPtGenVsJetPtRec[ic][ir] = new TH2D(name.Data(),
2435  "", bw*160,-20,300, bw*160,-20,300);
2436  fOutput->Add((TH2D*) fhJetPtGenVsJetPtRec[ic][ir]);
2437 
2438 
2439  name = (ic==0) ? Form("fhJetPtResolutionVsPtGenMB%s",bgtype[ir].Data()) :
2440  Form("fhJetPtResolutionVsPtGen%d%d%s",TMath::Nint(fCentralityBins[ic-1]),
2441  TMath::Nint(fCentralityBins[ic]),bgtype[ir].Data());
2442 
2443  fhJetPtResolutionVsPtGen[ic][ir] = new TH2D(name.Data(),
2444  "Resolution", 20,0,100, 200,-1.,1.);
2445  fOutput->Add((TH2D*) fhJetPtResolutionVsPtGen[ic][ir]);
2446  }
2447  }
2448 
2449 
2450  Double_t binsDiff [] = {
2451  -50, -48., -46., -44, -42.,
2452  -40., -38., -36., -34., -32., -30., -28., -26., -25., -24., -23., -22., -21.,
2453  -20., -19.5, -19., -18.5, -18., -17.5, -17., -16.5, -16., -15.5,
2454  -15., -14.5, -14., -13.5, -13., -12.5, -12., -11.5, -11., -10.5,
2455  -10.0, -9.8, -9.6, -9.4, -9.2,
2456  -9.0, -8.8, -8.6, -8.4, -8.2,
2457  -8.0, -7.8, -7.6, -7.4, -7.2,
2458  -7.0, -6.8, -6.6, -6.4, -6.2,
2459  -6.0, -5.8, -5.6, -5.4, -5.2,
2460  -5.0, -4.8, -4.6, -4.4, -4.2,
2461  -4.0, -3.8, -3.6, -3.4, -3.2,
2462  -3.0, -2.8, -2.6, -2.4, -2.2,
2463  -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1,
2464  -1.0,-0.95, -0.9, -0.85, -0.8, -0.75, -0.7, -0.65, -0.6, -0.55,
2465  -0.5,-0.48, -0.46, -0.44, -0.42,
2466  -0.4, -0.38, -0.36, -0.34, -0.32,
2467  -0.3, -0.28, -0.26, -0.24, -0.22,
2468  -0.2, -0.18, -0.16, -0.14, -0.12,
2469  -0.1, -0.09, -0.08, -0.07, -0.06,
2470  -0.05, -0.045, -0.04,-0.035 ,-0.03, -0.025, -0.02, -0.015, -0.01, -0.005,
2471  0, 0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.045,
2472  0.05, 0.06, 0.07, 0.08, 0.09,
2473  0.1, 0.12, 0.14, 0.16, 0.18,
2474  0.2, 0.22, 0.24, 0.26, 0.28,
2475  0.3, 0.32, 0.34, 0.36, 0.38,
2476  0.4, 0.42, 0.44, 0.46, 0.48,
2477  0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8,0.85, 0.9, 0.95,
2478  1.0,1.1, 1.2,1.3, 1.4,1.5, 1.6, 1.7, 1.8, 1.9,
2479  2.0, 2.2, 2.4, 2.6, 2.8,
2480  3.0, 3.2, 3.4, 3.6, 3.8,
2481  4.0, 4.2, 4.4, 4.6, 4.8,
2482  5.0, 5.2, 5.4, 5.6, 5.8,
2483  6.0, 6.2, 6.4, 6.6, 6.8,
2484  7.0, 7.2, 7.4, 7.6, 7.8,
2485  8.0, 8.2, 8.4, 8.6, 8.8,
2486  9.0, 9.2, 9.4, 9.6, 9.8,
2487  10., 10.5, 11., 11.5, 12., 12.5, 13., 13.5, 14., 14.5,
2488  15., 15.5, 16., 16.5, 17., 17.5, 18., 18.5, 19., 19.5,
2489  20., 21., 22., 23., 24, 25, 26., 28., 30., 32.,34., 36., 38.,
2490  40., 42., 44., 46., 48., 50.};
2491  Int_t nbinsDiff = sizeof(binsDiff)/sizeof(Double_t)-1;
2492 
2493  for(Int_t ic =0; ic<icmax; ic++){
2494  name = (ic==0) ? Form("fhPtTrkTruePrimRecMB") :
2495  Form("fhPtTrkTruePrimRec%d%d",TMath::Nint(fCentralityBins[ic-1]),
2496  TMath::Nint(fCentralityBins[ic]));
2497 
2498  fhPtTrkTruePrimRec[ic] = new TH2D(name.Data(),"",nbins, bins, 18,-0.9,0.9);
2499  fOutput->Add((TH2D*) fhPtTrkTruePrimRec[ic]);
2500 
2501 
2502  name = (ic==0) ? Form("fhPtTrkTruePrimGenMB") :
2503  Form("fhPtTrkTruePrimGen%d%d",TMath::Nint(fCentralityBins[ic-1]),
2504  TMath::Nint(fCentralityBins[ic]));
2505 
2506  fhPtTrkTruePrimGen[ic] = (TH2D*) fhPtTrkTruePrimRec[ic]->Clone(name.Data());
2507  fOutput->Add((TH2D*) fhPtTrkTruePrimGen[ic]);
2508 
2509  name = (ic==0) ? Form("fhPtTrkSecOrFakeRecMB") :
2510  Form("fhPtTrkSecOrFakeRec%d%d",TMath::Nint(fCentralityBins[ic-1]),
2511  TMath::Nint(fCentralityBins[ic]));
2512 
2513  fhPtTrkSecOrFakeRec[ic] = (TH2D*) fhPtTrkTruePrimRec[ic]->Clone(name.Data());
2514  fOutput->Add((TH2D*) fhPtTrkSecOrFakeRec[ic]);
2515  }
2516 
2517  for(Int_t iw=0; iw<21; iw++){
2518  name = Form("fhPtJetPrimVsPtJetRecMB%03d",TMath::Nint(100*iw*0.05));
2519  fhPtJetPrimVsPtJetRec[iw] = new TH2D(name.Data(), name.Data(),50,0,50,210,0,1.05);
2520  fOutput->Add((TH2D*) fhPtJetPrimVsPtJetRec[iw]);
2521  }
2522  name = Form("fhDiffPtVsPtTrackTrueMB");
2523  fhDiffPtVsPtTrackTrue = new TH2D(name.Data(), name.Data(),100,0,100, nbinsDiff,binsDiff);
2524  fOutput->Add((TH2D*) fhDiffPtVsPtTrackTrue);
2525  }
2526 
2527  //DCA DISTRIBUTIONS
2528  fhDCAinXVsPt = new TH2D("fhDCAinXVsPt","fhDCAinXVsPt",nbins, bins, 200, -10.,10);
2529  if(bNotKine) fOutput->Add((TH2D*) fhDCAinXVsPt);
2530  fhDCAinYVsPt = (TH2D*) fhDCAinXVsPt->Clone("fhDCAinYVsPt");
2531  if(bNotKine) fOutput->Add((TH2D*) fhDCAinYVsPt);
2532 
2533  fhDCAinXVsPtStrange = (TH2D*) fhDCAinXVsPt->Clone("fhDCAinXVsPtStrange");
2534  if(fTypeOfAnal == kEff) fOutput->Add((TH2D*) fhDCAinXVsPtStrange);
2535  fhDCAinYVsPtStrange = (TH2D*) fhDCAinXVsPt->Clone("fhDCAinYVsPtStrange");
2536  if(fTypeOfAnal == kEff) fOutput->Add((TH2D*) fhDCAinYVsPtStrange);
2537  fhDCAinXVsPtNonStrange = (TH2D*) fhDCAinXVsPt->Clone("fhDCAinXVsPtNonStrange");
2538  if(fTypeOfAnal == kEff) fOutput->Add((TH2D*) fhDCAinXVsPtNonStrange);
2539  fhDCAinYVsPtNonStrange = (TH2D*) fhDCAinXVsPt->Clone("fhDCAinYVsPtNonStrange");
2540  if(fTypeOfAnal == kEff) fOutput->Add((TH2D*) fhDCAinYVsPtNonStrange);
2541 
2542 
2543 
2544  if(fTypeOfAnal == kEff || fTypeOfAnal == kKine){
2545  for(Int_t ic =0; ic<icmax; ic++){
2546  for(Int_t it=0; it<kTT;it++){
2547  name = (ic==0) ? Form("fh1NtriggersGenMB") :
2548  Form("fh1NtriggersGen%d%d",TMath::Nint(fCentralityBins[ic-1]),
2549  TMath::Nint(fCentralityBins[ic]));
2550  name.Append(Form("TT%d%d",TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
2551 
2552  fh1NtriggersGen[ic][it] = (TH1D*) fh1Ntriggers[ic][it]->Clone(name.Data());
2553  fOutput->Add((TH1D*) fh1NtriggersGen[ic][it]);
2554 
2555 
2556  name = (ic==0) ? Form("fh1TriggerMultGenMB") :
2557  Form("fh1TriggerMultGen%d%d",TMath::Nint(fCentralityBins[ic-1]),
2558  TMath::Nint(fCentralityBins[ic]));
2559 
2560  name.Append(Form("TT%d%d",TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])));
2561  fh1TriggerMultGen[ic][it] = (TH1D*) fh1TriggerMult[ic][it]->Clone();
2562  fOutput->Add((TH1D*) fh1TriggerMultGen[ic][it]);
2563 
2564  //-------------------------
2565  for(Int_t ir=0; ir< kRho; ir++){
2566  name = (ic==0) ? "fHJetSpecGenMB" :
2567  Form("fHJetSpecGen%d%d",TMath::Nint(fCentralityBins[ic-1]),
2568  TMath::Nint(fCentralityBins[ic]));
2569 
2570  name.Append(Form("TT%d%d%s", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it]), bgtype[ir].Data()));
2571  fHJetSpecGen[ic][it][ir] = (TH2D*)fHJetSpec[ic][it][ir]->Clone(name.Data());
2572  fOutput->Add((TH2D*) fHJetSpecGen[ic][it][ir]);
2573  }
2574  }
2575  //-------------------------
2576  for(Int_t it=0; it<kTT; it++){
2577  name = (ic==0) ? Form("fhJetPhiGenMBTT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])) :
2578  Form("fhJetPhiGen%d%dTT%d%d",TMath::Nint(fCentralityBins[ic-1]),
2579  TMath::Nint(fCentralityBins[ic]),
2580  TMath::Nint(fTTlow[it]),
2581  TMath::Nint(fTThigh[it]));
2582 
2583  fhJetPhiGen[ic][it] = (TH2F*) fhJetPhi[ic][it]->Clone(name.Data());
2584  fOutput->Add((TH2F*) fhJetPhiGen[ic][it]);
2585  }
2586  //-------------------------
2587  for(Int_t it=0; it<kTT; it++){
2588  name = (ic==0) ? Form("fhJetEtaGenMBTT%d%d", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])) :
2589  Form("fhJetEtaGen%d%dTT%d%d",TMath::Nint(fCentralityBins[ic-1]),
2590  TMath::Nint(fCentralityBins[ic]),
2591  TMath::Nint(fTTlow[it]),
2592  TMath::Nint(fTThigh[it])
2593  );
2594 
2595 
2596  fhJetEtaGen[ic][it] = (TH2F*) fhJetEta[ic][it]->Clone(name.Data());
2597  fOutput->Add((TH2F*) fhJetEtaGen[ic][it]);
2598  }
2599  //-------------------------
2600  for(Int_t it=0; it<kTT; it++){
2601  name = (ic==0) ? Form("fhJetPhiGenMBTT%d%dRecoil", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])) :
2602  Form("fhJetPhiGen%d%dTT%d%dRecoil",TMath::Nint(fCentralityBins[ic-1]),
2603  TMath::Nint(fCentralityBins[ic]),
2604  TMath::Nint(fTTlow[it]),
2605  TMath::Nint(fTThigh[it]));
2606 
2607  fhJetPhiRecoilGen[ic][it] = (TH2F*) fhJetPhi[ic][it]->Clone(name.Data());
2608  fOutput->Add((TH2F*) fhJetPhiRecoilGen[ic][it]);
2609  }
2610  //-------------------------
2611  for(Int_t it=0; it<kTT; it++){
2612  name = (ic==0) ? Form("fhJetEtaGenMBTT%d%dRecoil", TMath::Nint(fTTlow[it]), TMath::Nint(fTThigh[it])) :
2613  Form("fhJetEtaGen%d%dTT%d%dRecoil",TMath::Nint(fCentralityBins[ic-1]),
2614  TMath::Nint(fCentralityBins[ic]),
2615  TMath::Nint(fTTlow[it]),
2616  TMath::Nint(fTThigh[it])
2617  );
2618 
2619 
2620  fhJetEtaRecoilGen[ic][it] = (TH2F*) fhJetEta[ic][it]->Clone(name.Data());
2621  fOutput->Add((TH2F*) fhJetEtaRecoilGen[ic][it]);
2622  }
2623 
2624  //-------------------------
2625  name = (ic==0) ? Form("fhTrackPtGenMB") :
2626  Form("fhTrackPtGen%d%d",TMath::Nint(fCentralityBins[ic-1]),
2627  TMath::Nint(fCentralityBins[ic]));
2628 
2629  fhTrackPtGen[ic] = (TH1F*) fhTrackPt[ic]->Clone(name.Data());
2630  fOutput->Add((TH1F*) fhTrackPtGen[ic]);
2631  //------------------------- Rongrong's analysis
2632  for(Int_t it=0; it< kTT; it++){
2633  for(Int_t ir=0; ir< kRho; ir++){
2634 
2635  name = Form("%sGen",fhDphiTriggerJet[ic][it][ir]->GetName());
2636  fhDphiTriggerJetGen[ic][it][ir] = (TH2F*) fhDphiTriggerJet[ic][it][ir]->Clone(name.Data());
2637  fOutput->Add((TH2F*) fhDphiTriggerJetGen[ic][it][ir]);
2638  }
2639  }
2640  }
2641  }
2642 
2643  // =========== Switch on Sumw2 for all histos ===========
2644  for(Int_t i=0; i<fOutput->GetEntries(); i++){
2645  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
2646  if(h1){
2647  h1->Sumw2();
2648  continue;
2649  }
2650  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
2651  if(hn){
2652  hn->Sumw2();
2653  }
2654  }
2655  TH1::AddDirectory(oldStatus);
2656 
2657 
2658  PostData(1, fOutput);
2659 }
2660 //________________________________________________________________________
2661 Bool_t AliAnalysisTaskHJetSpectra::RetrieveEventObjects() {
2662  //
2663  // retrieve event objects
2664  //
2665  if(!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) return kFALSE;
2666 
2667  return kTRUE;
2668 }
2669 //________________________________________________________________________
2670 Bool_t AliAnalysisTaskHJetSpectra::Run()
2671 {
2672  // Run analysis code here, if needed. It will be executed before FillHistograms().
2673 
2674  return kTRUE;
2675 }
2676 
2677 //________________________________________________________________________
2678 
2679 Double_t AliAnalysisTaskHJetSpectra::RelativePhi(Double_t mphi,Double_t vphi){
2680  //Get relative azimuthal angle of two particles -pi to pi
2681  if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
2682  else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
2683 
2684  if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2685  else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
2686 
2687  Double_t dphi = mphi - vphi;
2688  if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2689  else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
2690 
2691  return dphi;//dphi in [-Pi, Pi]
2692 }
2693 
2694 
2695 //________________________________________________________________________
2696 Double_t AliAnalysisTaskHJetSpectra::EstimateBgCone(AliJetContainer *jetCont, AliParticleContainer *trkCont, AliVParticle* triggerHadron, Bool_t isGen){
2697  //Estimate background rho by means of integrating track pT outside TT jet + recoil jet region
2698  //if TT exists find jet that containts TT and exclude range +- phiCut around the TT/TTjet in azimuth
2699 
2700  if(!trkCont) return 0.0;
2701 
2702  AliEmcalJet* jet = NULL;
2703  AliVParticle* track = NULL;
2704  Double_t phiTT = fRandom->Rndm()*TMath::TwoPi(); //in case of no TT make random dice
2705  Double_t etaTT = -fTrackEtaWindow + fRandom->Rndm()*2*fTrackEtaWindow;
2706  Bool_t bTTJetFound = kFALSE;
2707 
2708  if(triggerHadron){
2709 
2710  phiTT = triggerHadron->Phi();
2711  etaTT = triggerHadron->Eta();
2712 
2713  if(jetCont){
2714  //find ANY jet that contains TT if it exists
2715  jetCont->ResetCurrentID();
2716  while((jet = jetCont->GetNextAcceptJet())){ //loop over reconstructed jets
2717  if(!jet) continue;
2718  if(jet->Pt() < triggerHadron->Pt()*0.5) continue;
2719  //CUT ON JET ACCEPTANCE IS NOT NEEDED, WE SEARCH FOR ANY JET THAT CONTAINS TT
2720 
2721  for(Int_t iq=0; iq < jet->GetNumberOfTracks(); iq++){
2722  track = (AliVParticle*) (jet->TrackAt(iq,trkCont->GetArray())); //matched rec and emb tracks
2723  if(!track) continue;
2724  if(track != triggerHadron) continue;
2725 
2726  phiTT = jet->Phi(); // used phi,eta coordinates of the jet to exclude the TT jet
2727  etaTT = jet->Eta();
2728  bTTJetFound = kTRUE;
2729  break;
2730  }
2731  if(bTTJetFound) break; //skip the rest of jets when the jet with TT is found
2732  }
2733  }
2734  }
2735 
2736  phiTT = RelativePhi(phiTT,0.); //convert phi TT to (-pi,pi)
2737 
2738  if(TMath::Abs(etaTT) > fTrackEtaWindow){
2739  etaTT = (etaTT<0) ? -fTrackEtaWindow : fTrackEtaWindow;
2740  }
2741  //Sum pT outside TT+recoil jet region
2742  Double_t sumPt = 0.;
2743 
2744  trkCont->ResetCurrentID();
2745  while((track = trkCont->GetNextAcceptParticle())){
2746  if(!track) continue;
2747  if(!IsTrackInAcceptance(track, isGen)) continue;
2748 
2749  if(TMath::Abs(RelativePhi(phiTT+TMath::Pi(),track->Phi())) < fCutPhi) continue; //exclude recoil region of TT
2750  if(GetDeltaR(phiTT, track->Phi(), etaTT, track->Eta()) < fCutPhi) continue; //exclude region around TT
2751 
2752  if(fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
2753  if(TMath::Abs(track->GetLabel()) == 99999) continue;//reject embedded stuff 9999 set in AddTaskHJetSpectra.C
2754  }
2755 
2756  sumPt += track->Pt();
2757  }
2758  //Calculate area
2759  Double_t area = 2*fTrackEtaWindow*2*TMath::Pi();
2760  Double_t alpha;
2761 
2762  area -= 2*fTrackEtaWindow*2*fCutPhi; // subtract area of the recoil region
2763 
2764  if(TMath::Abs(etaTT) < fTrackEtaWindow - fCutPhi){ //TT cicle fully inside acceptance
2765  area -= fCutPhi*fCutPhi*TMath::Pi(); // subtract area of the trigger region
2766  }else{ //TT circle partly around acceptance
2767  alpha = TMath::ACos((fTrackEtaWindow - TMath::Abs(etaTT))/fCutPhi);
2768  area -= (fCutPhi*fCutPhi*(TMath::Pi()-alpha) + fCutPhi*TMath::Sin(alpha)*(fTrackEtaWindow - TMath::Abs(etaTT)));
2769  }
2770 
2771  return sumPt/area;
2772 
2773 }
2774 //________________________________________________________________________
2775 Double_t AliAnalysisTaskHJetSpectra::EstimateBgKT(AliJetContainer *jetCont, AliParticleContainer *trkCont, AliVParticle* triggerHadron){
2776  //Estimate rho from KT jet median. Ignore jet that contains TT
2777  Double_t rhoKT = 0.0;
2778 
2779  if(!jetCont) return rhoKT;
2780 
2781  AliEmcalJet* jet = NULL;
2782  AliVParticle* constTrack = NULL;
2783  Bool_t bKTJetCloseToTT = kFALSE;
2784  Int_t nJetAcc = 0;
2785  Double_t jetpt;
2786  Double_t sumEmbPt;
2787 
2788  jetCont->ResetCurrentID();
2789  while((jet = jetCont->GetNextAcceptJet())){ //loop over KT jets
2790  if(!jet) continue;
2791  if(!IsSignalJetInAcceptance(jet,kFALSE)) continue;
2792 
2793  bKTJetCloseToTT = kFALSE;
2794 
2795  if(triggerHadron){ //identify the KT jet which contains TT
2796  if(jet->Pt() > triggerHadron->Pt()*0.5){ //jet containing TT has pT larger than pT of TT
2797  for(Int_t iq=0; iq < jet->GetNumberOfTracks(); iq++) {
2798  constTrack = (AliVParticle*) (jet->TrackAt(iq,trkCont->GetArray())); //matched rec and emb tracks
2799  if(!constTrack) continue;
2800  if(constTrack == triggerHadron){
2801  bKTJetCloseToTT = kTRUE;
2802  break;
2803  }
2804  }
2805  }
2806  }
2807  if(bKTJetCloseToTT) continue; //skip the jet that contains TT
2808 
2809  sumEmbPt = 0.;//sum pt of embedded tracks in jet which is to be subtracted
2810  if(fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
2811  for(Int_t iq=0; iq < jet->GetNumberOfTracks(); iq++) {
2812  constTrack = (AliVParticle*) (jet->TrackAt(iq,trkCont->GetArray()));
2813  if(!constTrack) continue;
2814  if(TMath::Abs(constTrack->GetLabel()) == 99999){
2815  sumEmbPt += constTrack->Pt();
2816  }
2817  }
2818  }
2819 
2820  jetpt = jet->Pt()- sumEmbPt; //subtract embedded pt
2821  if(triggerHadron) fhKTAreaPt->Fill(jetpt,jet->Area());
2822 
2823  if(jetpt <0.005) jetpt = 0.; //set pt of ghost jets identical to zero
2824  frhovec[nJetAcc] = jetpt/jet->Area();
2825  nJetAcc++;
2826  }
2827 
2828  if(nJetAcc>0){
2829  rhoKT = TMath::Median(nJetAcc, frhovec);
2830  }
2831 
2832  return rhoKT;
2833 }
2834 //________________________________________________________________________
2835 Double_t AliAnalysisTaskHJetSpectra::EstimateBgKTcms(AliJetContainer *jetCont, AliParticleContainer *trkCont, AliVParticle* triggerHadron){
2836  //Estimate rho from KT jet median ala CMS. Ignore jet that contains TT
2837  Double_t rhoKTcms = 0.0;
2838 
2839  if(!jetCont) return rhoKTcms;
2840 
2841  AliEmcalJet* jet = NULL;
2842  AliVParticle* constTrack = NULL;
2843  Bool_t bKTJetCloseToTT = kFALSE;
2844  Int_t nJetAcc = 0;
2845  Double_t areaPhysJets = 0.0;
2846  Double_t areaAllJets = 0.0;
2847  Double_t jetpt;
2848  Double_t sumEmbPt = 0.;
2849 
2850  jetCont->ResetCurrentID();
2851  while((jet = jetCont->GetNextAcceptJet())){ //loop over KT jets
2852  if(!jet) continue;
2853  if(!IsSignalJetInAcceptance(jet,kFALSE)) continue;
2854 
2855  bKTJetCloseToTT = kFALSE;
2856 
2857  if(triggerHadron){ //identify the KT jet which contains TT
2858  if(jet->Pt() > triggerHadron->Pt()*0.5){ //jet containing TT has pT larger than pT of TT
2859  for(Int_t iq=0; iq < jet->GetNumberOfTracks(); iq++) {
2860  constTrack = (AliVParticle*) (jet->TrackAt(iq,trkCont->GetArray())); //matched rec and emb tracks
2861  if(!constTrack) continue;
2862  if(constTrack != triggerHadron) continue;
2863  bKTJetCloseToTT = kTRUE;
2864  break;
2865  }
2866  }
2867  }
2868  if(bKTJetCloseToTT) continue; //skip the jet that contains TT
2869 
2870  sumEmbPt = 0.;
2871  if(fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
2872  for(Int_t iq=0; iq < jet->GetNumberOfTracks(); iq++) {
2873  constTrack = (AliVParticle*) (jet->TrackAt(iq,trkCont->GetArray())); //matched rec and emb tracks
2874  if(!constTrack) continue;
2875  if(TMath::Abs(constTrack->GetLabel()) == 99999){
2876  sumEmbPt += constTrack->Pt();
2877  }
2878  }
2879  }
2880 
2881 
2882  areaAllJets += jet->Area();
2883 
2884  jetpt = jet->Pt()- sumEmbPt; //subtract pt of embedded tracks
2885 
2886  if(jetpt > 0.1){
2887  areaPhysJets += jet->Area();
2888  frhovec[nJetAcc] = jetpt/jet->Area();
2889  nJetAcc++;
2890  }
2891  }
2892 
2893  if(nJetAcc>0){
2894  rhoKTcms = TMath::Median(nJetAcc, frhovec)*(areaPhysJets/areaAllJets);
2895  }
2896 
2897  return rhoKTcms;
2898 }
2899 
2900 //________________________________________________________________________
2901 Double_t AliAnalysisTaskHJetSpectra::GetDeltaR(Double_t phi1, Double_t phi2, Double_t eta1, Double_t eta2){
2902  //angular distance between two jets
2903  Double_t dphi = RelativePhi(phi1,phi2);
2904  Double_t deta = eta1 - eta2;
2905  return sqrt(dphi*dphi + deta*deta);
2906 
2907 }
2908 
2909 //________________________________________________________________________
2910 Double_t AliAnalysisTaskHJetSpectra::GetFractionSharedPt(AliEmcalJet *jRec, AliJetContainer *jconRec, AliEmcalJet *jGen, AliJetContainer *jconGen){
2911 
2912  //get fraction of pT shared by reconstructed and generated level jet
2913  if(!jRec) return -1.0;
2914  if(!jconRec) return -1.0;
2915  if(!jGen) return -1.0;
2916  if(!jconGen) return -1.0;
2917 
2918  Double_t fraction = 0., sumPt = 0.;
2919  Double_t jetPt2 = jGen->Pt();
2920  //Int_t idxGen, idxRec;
2921  AliVParticle *pgen, *prec;
2922  if(jetPt2>0){
2923 
2924  for(Int_t ig=0; ig< jGen->GetNumberOfTracks(); ig++) {
2925  pgen = (AliVParticle*) (jGen->TrackAt(ig, jconGen->GetParticleContainer()->GetArray()));
2926  if(!pgen) continue;
2927 
2928  for(Int_t ir=0; ir< jRec->GetNumberOfTracks(); ir++){
2929  prec = (AliVParticle*) (jRec->TrackAt(ir, jconRec->GetParticleContainer()->GetArray()));
2930  if(!prec) continue;
2931 
2932  if(TMath::Abs(prec->GetLabel()) == TMath::Abs(pgen->GetLabel())){
2933 
2934  if(fTypeOfAnal == kEmb || fTypeOfAnal == kEmbSingl){
2935  //All embedded tracks have the same label check also spatial coordinates
2936  if(TMath::Abs(prec->Eta() - pgen->Eta()) > 1e-4) continue;
2937  if(TMath::Abs(RelativePhi(prec->Phi(), pgen->Phi())) > 1e-4) continue;
2938  if(TMath::Abs(prec->Pt() - pgen->Pt()) > 1e-4) continue;
2939  //if(fDebug>20){
2940  //Printf("fraction TRACK REC eta = %f; phi = %f; pt = %f", prec->Eta(), prec->Phi(), prec->Pt());
2941  //Printf("fraction TRACK GEN eta = %f; phi = %f; pt = %f", pgen->Eta(), pgen->Phi(), pgen->Pt());
2942  //}
2943 
2944  }
2945 
2946  sumPt += pgen->Pt();
2947  break;
2948  }
2949  }
2950  }
2951 
2952 
2953  fraction = sumPt/jetPt2;
2954  } else{
2955  fraction = -1;
2956  }
2957 
2958  //if(fDebug>20) Printf("fraction return = %f ",fraction);
2959 
2960  return fraction;
2961 
2962 }
2963 
2964 //________________________________________________________________________
2965 Bool_t AliAnalysisTaskHJetSpectra::IsStrange(Int_t ip){
2966  //Check whtether the particle is strange
2967  if(ip==321) return kTRUE; //K+
2968  if(ip==130) return kTRUE; //KOL
2969  if(ip==310) return kTRUE; //K0s
2970  if(ip==3112) return kTRUE; //Sigma -
2971  if(ip==3122) return kTRUE; //Lambda 0
2972  if(ip==3222) return kTRUE; //Sigma+
2973  if(ip==3312) return kTRUE; //Xi-
2974  if(ip==3322) return kTRUE; //Xi0
2975  if(ip==3334) return kTRUE; //Omega-
2976 
2977  return kFALSE;
2978 }
2979 
AliAnalysisTaskHJetSpectra::GetDeltaR
Double_t GetDeltaR(Double_t phi1, Double_t phi2, Double_t eta1, Double_t eta2)
Definition: AliAnalysisTaskHJetSpectra.cxx:2901
AliAnalysisTaskHJetSpectra::fhPtTrkSecOrFakeRec
TH2D * fhPtTrkSecOrFakeRec[kCAll]
pt spectrum of true generated primary track
Definition: AliAnalysisTaskHJetSpectra.h:302
AliAnalysisTaskHJetSpectra::fhVertexZAcceptTT
TH1F * fhVertexZAcceptTT
gc vertexZ accepted after vtx cut
Definition: AliAnalysisTaskHJetSpectra.h:260
AliAnalysisTaskHJetSpectra::SetTT
void SetTT(Double_t tlr, Double_t thr, Double_t tls, Double_t ths)
Definition: AliAnalysisTaskHJetSpectra.cxx:402
AliAnalysisTaskHJetSpectra::fhJetPhiRecoil
TH2F * fhJetPhiRecoil[kCAll][kTT]
minimum bias eta for recoil jets gen
Definition: AliAnalysisTaskHJetSpectra.h:249
AliAnalysisTaskHJetSpectra::fhVertexZ
TH1F * fhVertexZ
gc X=centrality; Y= track pT
Definition: AliAnalysisTaskHJetSpectra.h:254
AliAnalysisTaskHJetSpectra::fTTlow
Double_t fTTlow[kTT]
gc trigger if tracks/jets are loaded initiates calling ExecOnce
Definition: AliAnalysisTaskHJetSpectra.h:209
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:130
AliParticleContainer::GetNextAcceptParticle
virtual AliVParticle * GetNextAcceptParticle()
Definition: AliParticleContainer.cxx:132
AliAnalysisTaskHJetSpectra::fhVzeroATotMultTT
TH1F * fhVzeroATotMultTT[kCAll][kTT]
V0A multiplicity for given V0A centrality selection.
Definition: AliAnalysisTaskHJetSpectra.h:279
AliPicoTrack.h
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskHJetSpectra.h
AliAnalysisTaskHJetSpectra::fhDeltaPtEmbBc2Bc
TH1D * fhDeltaPtEmbBc2Bc[kCAll][kTT][kRho]
embedded delta pT versus pT of the embedded jet (emb track is perp to TT)
Definition: AliAnalysisTaskHJetSpectra.h:235
AliAnalysisTaskHJetSpectra::fTrigTracksGen
AliVParticle * fTrigTracksGen[kTT][999]
Definition: AliAnalysisTaskHJetSpectra.h:332
AliAnalysisTaskHJetSpectra::fhDCAinXVsPt
TH2D * fhDCAinXVsPt
hybrid TPC constrained track phi vs track pT
Definition: AliAnalysisTaskHJetSpectra.h:314
AliAnalysisTaskEmcalJet.h
TH2F
Definition: External.C:236
AliAnalysisTaskHJetSpectra::fhCentralityV0C
TH1F * fhCentralityV0C
centrality from V0A
Definition: AliAnalysisTaskHJetSpectra.h:275
AliAnalysisTaskHJetSpectra::fhVertexZAccept
TH1F * fhVertexZAccept
gc vertexZ accepted after vtx cut
Definition: AliAnalysisTaskHJetSpectra.h:257
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliAnalysisTaskHJetSpectra::fhTrackPhi
TH2F * fhTrackPhi[kCAll]
gc jet phi vs jet pT
Definition: AliAnalysisTaskHJetSpectra.h:244
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:377
AliAnalysisTaskHJetSpectra::EstimateBgKTcms
Double_t EstimateBgKTcms(AliJetContainer *jetCont, AliParticleContainer *trkArray, AliVParticle *triggerHadron)
Definition: AliAnalysisTaskHJetSpectra.cxx:2835
AliAnalysisTaskHJetSpectra::fhJetPtResolutionVsPtGen
TH2D * fhJetPtResolutionVsPtGen[kCAll][kRho]
pt jet gen level vs pT jet rec level
Definition: AliAnalysisTaskHJetSpectra.h:299
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:121
AliAnalysisTaskHJetSpectra::EstimateBgCone
Double_t EstimateBgCone(AliJetContainer *jetCont, AliParticleContainer *trkArray, AliVParticle *triggerHadron, Bool_t isGen=kFALSE)
Definition: AliAnalysisTaskHJetSpectra.cxx:2696
AliAnalysisTaskHJetSpectra::fhPtTrkTruePrimGen
TH2D * fhPtTrkTruePrimGen[kCAll]
pt spectrum of true reconstructed primary tracks
Definition: AliAnalysisTaskHJetSpectra.h:301
AliEmcalJet::Phi
Double_t Phi() const
Definition: AliEmcalJet.h:117
AliAnalysisTaskHJetSpectra::FillHistograms
Bool_t FillHistograms()
Function filling histograms.
Definition: AliAnalysisTaskHJetSpectra.cxx:723
AliAnalysisTaskHJetSpectra::fhKTAreaPt
TH2F * fhKTAreaPt
delta pT from RndCone using rho from perp cone inclusive event
Definition: AliAnalysisTaskHJetSpectra.h:240
AliAnalysisTaskHJetSpectra::IsEventInAcceptance
Bool_t IsEventInAcceptance(AliVEvent *event)
Definition: AliAnalysisTaskHJetSpectra.cxx:546
rl
AliRunLoader * rl
Definition: AliMCTreeTools.cxx:73
AliAnalysisTaskHJetSpectra::fhDeltaPt
TH1D * fhDeltaPt[kCAll][kTT][kRho]
jet area times rho from perp cone
Definition: AliAnalysisTaskHJetSpectra.h:230
AliAnalysisTaskHJetSpectra::fhDeltaPtEmbPerp
TH1D * fhDeltaPtEmbPerp[kCAll][kTT][kRho]
embedded delta pT versus pT of the embedded jet
Definition: AliAnalysisTaskHJetSpectra.h:233
AliAnalysisTaskHJetSpectra::IsSignalJetInAcceptance
Bool_t IsSignalJetInAcceptance(AliEmcalJet *jet, Bool_t suppressGhost=1)
Definition: AliAnalysisTaskHJetSpectra.cxx:615
AliAnalysisTaskHJetSpectra::fhDeltaPtEmb2D
TH2D * fhDeltaPtEmb2D[kCAll][kTT][kRho]
embedded delta pT
Definition: AliAnalysisTaskHJetSpectra.h:232
AliAnalysisTaskHJetSpectra::fh1NtriggersGen
TH1D * fh1NtriggersGen[kCAll][kTT]
tirgger multiplicity in event
Definition: AliAnalysisTaskHJetSpectra.h:219
AliAnalysisTaskHJetSpectra::fhJetEta
TH2F * fhJetEta[kCAll][kTT]
gc track phi vs track pT
Definition: AliAnalysisTaskHJetSpectra.h:245
AliAnalysisTaskHJetSpectra::fhRhoIncl
TH1F * fhRhoIncl[kCAll][kRho]
gc X=rho from perp cone, Y=centrality
Definition: AliAnalysisTaskHJetSpectra.h:225
fCollisionSystem
Int_t fCollisionSystem
Definition: makeTFile4CutsLctoV0bachelor.C:26
AliAnalysisTaskHJetSpectra::fRhoRec
TArrayD fRhoRec[kTT]
Y DCA versus pT of non strange tracks.
Definition: AliAnalysisTaskHJetSpectra.h:322
AliAnalysisTaskHJetSpectra::fhVertexXAccept
TH1F * fhVertexXAccept
gc vertexZ inclusive
Definition: AliAnalysisTaskHJetSpectra.h:255
AliAnalysisTaskHJetSpectra::fhDeltaPtEmb
TH1D * fhDeltaPtEmb[kCAll][kTT][kRho]
delta pT
Definition: AliAnalysisTaskHJetSpectra.h:231
AliAnalysisTaskHJetSpectra::fh1TriggerMultGen
TH1D * fh1TriggerMultGen[kCAll][kTT]
trigger counter
Definition: AliAnalysisTaskHJetSpectra.h:220
AliAnalysisHelperJetTasks.h
AliAnalysisTaskHJetSpectra::Run
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Definition: AliAnalysisTaskHJetSpectra.cxx:2670
AliAnalysisTaskHJetSpectra::fhTrackPt
TH1F * fhTrackPt[kCAll]
track eta vs track pT
Definition: AliAnalysisTaskHJetSpectra.h:252
AliAnalysisTaskHJetSpectra::fhDeltaPtIncl
TH1D * fhDeltaPtIncl[kCAll][kRho]
embedded delta pT versus pT of the embedded jet (emb track is backtoback in azimtuh w...
Definition: AliAnalysisTaskHJetSpectra.h:238
AliAnalysisTaskHJetSpectra::fhInvPtQVsPhi
TH2D * fhInvPtQVsPhi[2]
track Y= rec pt - true pt X= true track pT
Definition: AliAnalysisTaskHJetSpectra.h:306
AliAnalysisTaskHJetSpectra::fRandom
TRandom3 * fRandom
impact parameter from hijing
Definition: AliAnalysisTaskHJetSpectra.h:204
AliAnalysisTaskHJetSpectra::fh1TriggerMult
TH1D * fh1TriggerMult[kCAll][kTT]
trigger counter
Definition: AliAnalysisTaskHJetSpectra.h:218
AliAnalysisTaskHJetSpectra::fhSigmaPtOverPtVsPt
TH2D * fhSigmaPtOverPtVsPt[2]
q*1/pT versus eta
Definition: AliAnalysisTaskHJetSpectra.h:310
AliAnalysisTaskHJetSpectra::fhVertexXAcceptTT
TH1F * fhVertexXAcceptTT
gc vertexZ accepted after vtx cut
Definition: AliAnalysisTaskHJetSpectra.h:258
AliParticleContainer
Container for particles within the EMCAL framework.
Definition: AliParticleContainer.h:31
AliEmcalJet::TrackAt
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
AliEmcalJet::GetNumberOfTracks
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
AliAnalysisTaskHJetSpectra::GetFractionSharedPt
Double_t GetFractionSharedPt(AliEmcalJet *jRec, AliJetContainer *jconRec, AliEmcalJet *jGen, AliJetContainer *jconGen)
Definition: AliAnalysisTaskHJetSpectra.cxx:2910
AliAnalysisTaskHJetSpectra::fhZNAVzeroATrackTT
THnSparse * fhZNAVzeroATrackTT[kCAll][kTT]
ZNA energy versus Vzero A mult versus track mult.
Definition: AliAnalysisTaskHJetSpectra.h:288
AliAnalysisTaskHJetSpectra::fhTrackMultiplicity
TH1D * fhTrackMultiplicity[kCAll]
ZDC A neutral energy.
Definition: AliAnalysisTaskHJetSpectra.h:284
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1639
AliAnalysisTaskHJetSpectra::fhDphiTriggerJet
TH2F * fhDphiTriggerJet[kCAll][kTT][kRho]
gc vertexZ accepted after vtx cut in MC
Definition: AliAnalysisTaskHJetSpectra.h:263
AliAnalysisTaskHJetSpectra::fhCentrality
TH1F * fhCentrality[kCAll]
minimum bias eta inclusive
Definition: AliAnalysisTaskHJetSpectra.h:271
AliAnalysisTaskHJetSpectra::IsMCEventInAcceptance
Bool_t IsMCEventInAcceptance(AliVEvent *event)
Definition: AliAnalysisTaskHJetSpectra.cxx:523
AliJetContainer::GetParticleContainer
AliParticleContainer * GetParticleContainer() const
Definition: AliJetContainer.h:180
AliAnalysisTaskHJetSpectra::fhPtJetPrimVsPtJetRec
TH2D * fhPtJetPrimVsPtJetRec[21]
pt spectrum of reconstructed fake or secondary tracks
Definition: AliAnalysisTaskHJetSpectra.h:303
AliAnalysisTaskHJetSpectra::fhTrackPtGen
TH1F * fhTrackPtGen[kCAll]
gc X=centrality; Y= track pT
Definition: AliAnalysisTaskHJetSpectra.h:253
AliAnalysisTaskHJetSpectra::fhZNAVzeroATrack
THnSparse * fhZNAVzeroATrack[kCAll]
multiplicity of tracks in event with TT track
Definition: AliAnalysisTaskHJetSpectra.h:287
AliAnalysisTaskHJetSpectra::fh1Ntriggers
TH1D * fh1Ntriggers[kCAll][kTT]
gc event statistics
Definition: AliAnalysisTaskHJetSpectra.h:217
AliAnalysisTaskHJetSpectra::fhJetPtGen
TH1D * fhJetPtGen[kCAll][kRho]
impact parameter distribution hijing versus TT
Definition: AliAnalysisTaskHJetSpectra.h:297
AliAnalysisTaskHJetSpectra::fhDCAinYVsPtNonStrange
TH2D * fhDCAinYVsPtNonStrange
X DCA versus pT of non strange tracks.
Definition: AliAnalysisTaskHJetSpectra.h:319
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskHJetSpectra::fhCentralityV0A
TH1F * fhCentralityV0A
centrality V0 multiplicity A+C
Definition: AliAnalysisTaskHJetSpectra.h:274
AliAnalysisTaskHJetSpectra::fhJetEtaGen
TH2F * fhJetEtaGen[kCAll][kTT]
jet eta vs jet pT
Definition: AliAnalysisTaskHJetSpectra.h:246
TH1I
Definition: External.C:204
UInt_t
unsigned int UInt_t
Definition: External.C:33
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskHJetSpectra::IsTrackInAcceptance
Bool_t IsTrackInAcceptance(AliVParticle *track, Bool_t isGen=0)
Definition: AliAnalysisTaskHJetSpectra.cxx:594
AliAnalysisTaskHJetSpectra::fhZNAEnergy
TH1F * fhZNAEnergy[kCAll]
V0A multiplicity.
Definition: AliAnalysisTaskHJetSpectra.h:281
AliAnalysisTaskHJetSpectra::fhDeltaPtEmbPerp2D
TH2D * fhDeltaPtEmbPerp2D[kCAll][kTT][kRho]
embedded delta pT (emb track is perp to TT)
Definition: AliAnalysisTaskHJetSpectra.h:234
AliAnalysisTaskHJetSpectra::fHJetSpecGen
TH2D * fHJetSpecGen[kCAll][kTT][kRho]
TT associated spectrum of jets.
Definition: AliAnalysisTaskHJetSpectra.h:222
AliAnalysisTaskHJetSpectra::EstimateBgKT
Double_t EstimateBgKT(AliJetContainer *jetCont, AliParticleContainer *trkArray, AliVParticle *trackTT)
Definition: AliAnalysisTaskHJetSpectra.cxx:2775
AliAnalysisTaskHJetSpectra::fhJetEtaRecoilGen
TH2F * fhJetEtaRecoilGen[kCAll][kTT]
minimum bias eta fore recoil jets
Definition: AliAnalysisTaskHJetSpectra.h:248
AliAnalysisTaskHJetSpectra::fhVertexZAcceptMC
TH1F * fhVertexZAcceptMC
gc vertexZ inclusive in MC
Definition: AliAnalysisTaskHJetSpectra.h:262
AliAnalysisTaskHJetSpectra::fhInvPtQVsPhiCSide
TH2D * fhInvPtQVsPhiCSide[2]
q*1/pT versus eta
Definition: AliAnalysisTaskHJetSpectra.h:309
AliParticleContainer.h
AliAnalysisTaskHJetSpectra::fhJetPhiRecoilGen
TH2F * fhJetPhiRecoilGen[kCAll][kTT]
minimum bias phi for recoil jets
Definition: AliAnalysisTaskHJetSpectra.h:250
AliAnalysisTaskHJetSpectra::fHelperClass
AliAnalysisUtils * fHelperClass
A random number.
Definition: AliAnalysisTaskHJetSpectra.h:205
TH2D
Definition: External.C:228
TH1D
Definition: External.C:212
AliJetContainer.h
AliAnalysisTaskHJetSpectra::GetConePt
Double_t GetConePt(Double_t eta, Double_t phi, Double_t radius, AliParticleContainer *trkArray, Bool_t isGen)
Definition: AliAnalysisTaskHJetSpectra.cxx:410
AliAnalysisTaskHJetSpectra::fhDphiTTTT
TH1F * fhDphiTTTT[kTT]
Dphi of accepted jets after dphi cut.
Definition: AliAnalysisTaskHJetSpectra.h:266
AliAnalysisTaskHJetSpectra::fhJetPhi
TH2F * fhJetPhi[kCAll][kTT]
KT jets area versus PT.
Definition: AliAnalysisTaskHJetSpectra.h:242
AliAnalysisTaskHJetSpectra::fhTrackMultiplicityTT
TH1D * fhTrackMultiplicityTT[kCAll][kTT]
multiplicity of tracks
Definition: AliAnalysisTaskHJetSpectra.h:285
AliRhoParameter.h
AliAnalysisTaskHJetSpectra::fhJetPtGenVsJetPtRec
TH2D * fhJetPtGenVsJetPtRec[kCAll][kRho]
pt distribution of generator level jets
Definition: AliAnalysisTaskHJetSpectra.h:298
AliAnalysisTaskHJetSpectra::fhDphiTriggerJetGen
TH2F * fhDphiTriggerJetGen[kCAll][kTT][kRho]
gc Delta phi versus jet pT
Definition: AliAnalysisTaskHJetSpectra.h:264
AliJetContainer::GetNextAcceptJet
AliEmcalJet * GetNextAcceptJet()
Definition: AliJetContainer.cxx:299
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskHJetSpectra::RelativePhi
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Definition: AliAnalysisTaskHJetSpectra.cxx:2679
AliAnalysisTaskHJetSpectra::fhVertexYAccept
TH1F * fhVertexYAccept
gc vertexZ accepted after vtx cut
Definition: AliAnalysisTaskHJetSpectra.h:256
AliAnalysisTaskHJetSpectra::fhZNAEnergyTT
TH1F * fhZNAEnergyTT[kCAll][kTT]
ZDC A neutral energy for given V0A centrality selection.
Definition: AliAnalysisTaskHJetSpectra.h:282
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliAnalysisTaskHJetSpectra::fhDCAinXVsPtStrange
TH2D * fhDCAinXVsPtStrange
Y DCA versus pT.
Definition: AliAnalysisTaskHJetSpectra.h:316
AliEmcalList
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
AliAnalysisTaskHJetSpectra::fhJetEtaIncl
TH2F * fhJetEtaIncl
minimum bias phi inclusive
Definition: AliAnalysisTaskHJetSpectra.h:269
AliAnalysisTaskHJetSpectra::fhRhoTT
TH1F * fhRhoTT[kCAll][kTT][kRho]
TT associated spectrum of jets.
Definition: AliAnalysisTaskHJetSpectra.h:224
AliAnalysisTaskHJetSpectra::fhDphiTriggerJetAccept
TH1F * fhDphiTriggerJetAccept
gc Delta phi versus jet pT
Definition: AliAnalysisTaskHJetSpectra.h:265
AliAnalysisTaskHJetSpectra::fhJetPhiIncl
TH2F * fhJetPhiIncl
Dphi between multiple trigger tracks.
Definition: AliAnalysisTaskHJetSpectra.h:268
AliAnalysisTaskHJetSpectra::fHJetSpec
TH2D * fHJetSpec[kCAll][kTT][kRho]
trigger multiplicity in event
Definition: AliAnalysisTaskHJetSpectra.h:221
AliAnalysisTaskHJetSpectra::fImpParam
Double_t fImpParam
gc value is filled, if pythia header is accessible
Definition: AliAnalysisTaskHJetSpectra.h:201
AliAnalysisTaskHJetSpectra::fhJetEtaRecoil
TH2F * fhJetEtaRecoil[kCAll][kTT]
jet eta vs jet pT
Definition: AliAnalysisTaskHJetSpectra.h:247
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:924
AliAnalysisTaskEmcal.h
AliAnalysisTaskHJetSpectra::fhDiffPtVsPtTrackTrue
TH2D * fhDiffPtVsPtTrackTrue
pt spectrum of reconstructed jets without fake track pT vs reconstructed jet pT
Definition: AliAnalysisTaskHJetSpectra.h:304
AliAnalysisTaskHJetSpectra::fARhoTT
TH1F * fARhoTT[kCAll][kTT][kRho]
gc X=rho from perp cone, Y=centrality
Definition: AliAnalysisTaskHJetSpectra.h:227
AliAnalysisTaskHJetSpectra::fhTrackEta
TH2F * fhTrackEta[kCAll]
minimum bias phi for recoil jets gen
Definition: AliAnalysisTaskHJetSpectra.h:251
AliAnalysisTaskHJetSpectra::GetDeltaPt
void GetDeltaPt(Int_t nrho, TArrayD &rho, Double_t *dpt, Double_t ttPhi, Double_t ttEta, AliParticleContainer *trkArray, Bool_t isGen)
Definition: AliAnalysisTaskHJetSpectra.cxx:646
AliAnalysisTaskHJetSpectra::fhVertexYAcceptTT
TH1F * fhVertexYAcceptTT
gc vertexZ accepted after vtx cut
Definition: AliAnalysisTaskHJetSpectra.h:259
AliAnalysisTaskEmcalJet
Base task in the EMCAL jet framework.
Definition: AliAnalysisTaskEmcalJet.h:30
AliAnalysisTaskHJetSpectra::fhInvPtQVsPhiASide
TH2D * fhInvPtQVsPhiASide[2]
q*1/pT versus eta
Definition: AliAnalysisTaskHJetSpectra.h:308
AliEmcalJet
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskEmcal::UserCreateOutputObjects
void UserCreateOutputObjects()
Main initialization function on the worker.
Definition: AliAnalysisTaskEmcal.cxx:342
AliAnalysisTaskHJetSpectra::fhDCAinYVsPtStrange
TH2D * fhDCAinYVsPtStrange
X DCA versus pT of strange tracks.
Definition: AliAnalysisTaskHJetSpectra.h:317
nbins
const Int_t nbins
Definition: AverageDmesonRaa.C:63
AliAnalysisTaskHJetSpectra::fhVzeroATotMult
TH1F * fhVzeroATotMult[kCAll]
centrality from ZNA
Definition: AliAnalysisTaskHJetSpectra.h:278
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskHJetSpectra::fhVertexZMC
TH1F * fhVertexZMC
gc vertexZ accepted after vtx cut
Definition: AliAnalysisTaskHJetSpectra.h:261
AliAnalysisTaskHJetSpectra::fhInvPtQVsEta
TH2D * fhInvPtQVsEta[2]
q*1/pT versus phi
Definition: AliAnalysisTaskHJetSpectra.h:307
AliAnalysisTaskHJetSpectra::fTrigTracks
AliVParticle * fTrigTracks[kTT][999]
Definition: AliAnalysisTaskHJetSpectra.h:333
AliAnalysisTaskHJetSpectra::fhImpactParameter
TH1D * fhImpactParameter[kCAll]
ZNA energy versus Vzero mult. versus track mult. in events with TT.
Definition: AliAnalysisTaskHJetSpectra.h:294
AliAnalysisTaskHJetSpectra::fhTrackPhiTPCG
TH2F * fhTrackPhiTPCG
hybrid constrained global track phi vs track pT
Definition: AliAnalysisTaskHJetSpectra.h:312
AliAnalysisTaskHJetSpectra::RetrieveEventObjects
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskHJetSpectra.cxx:2661
AliAnalysisTaskHJetSpectra::fhDCAinXVsPtNonStrange
TH2D * fhDCAinXVsPtNonStrange
Y DCA versus pT of strange tracks.
Definition: AliAnalysisTaskHJetSpectra.h:318
AliAnalysisTaskHJetSpectra::fhDeltaPtEmbBc2Bc2D
TH2D * fhDeltaPtEmbBc2Bc2D[kCAll][kTT][kRho]
embedded delta pT (emb track is back-to-back in azimuth to TT)
Definition: AliAnalysisTaskHJetSpectra.h:236
AliAnalysisTaskHJetSpectra::fhCentralityZNA
TH1F * fhCentralityZNA
centrality from V0C
Definition: AliAnalysisTaskHJetSpectra.h:276
AliAnalysisTaskHJetSpectra::fhDCAinYVsPt
TH2D * fhDCAinYVsPt
X DCA versus pT.
Definition: AliAnalysisTaskHJetSpectra.h:315
AliAnalysisTaskHJetSpectra::fInitializedLocal
Bool_t fInitializedLocal
gc Vertex selection helper
Definition: AliAnalysisTaskHJetSpectra.h:206
AliEmcalJet.h
AliAnalysisTaskHJetSpectra::fhPtTrkTruePrimRec
TH2D * fhPtTrkTruePrimRec[kCAll]
pt jet resolution
Definition: AliAnalysisTaskHJetSpectra.h:300
AliJetContainer
Container for jet within the EMCAL jet framework.
Definition: AliJetContainer.h:37
AliAnalysisTaskHJetSpectra::fhJetPhiGen
TH2F * fhJetPhiGen[kCAll][kTT]
gc jet phi vs jet pT
Definition: AliAnalysisTaskHJetSpectra.h:243
TH1
Definition: External.C:196
AliAnalysisTaskHJetSpectra::fhImpactParameterTT
TH1D * fhImpactParameterTT[kCAll][kTT]
impact parameter distribution hijing
Definition: AliAnalysisTaskHJetSpectra.h:295
AliAnalysisTaskHJetSpectra::fhCentralityV0M
TH1F * fhCentralityV0M
centrality V0 multiplicity A+C when TT is present
Definition: AliAnalysisTaskHJetSpectra.h:273