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