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