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