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