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