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