AliPhysics  32e057f (32e057f)
AliAnalysisTaskFullpAJets.cxx
Go to the documentation of this file.
2 
3 #include <Riostream.h>
4 #include <ctime>
5 #include <TString.h>
6 #include <TChain.h>
7 #include <TTree.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <THnSparse.h>
12 #include <TCanvas.h>
13 #include <TList.h>
14 #include <TLorentzVector.h>
15 #include <TProfile.h>
16 #include <TProfile2D.h>
17 #include <TProfile3D.h>
18 #include <TRandom.h>
19 #include <TRandom3.h>
20 #include <TClonesArray.h>
21 #include <TObjArray.h>
22 
23 #include "AliAnalysisTaskEmcal.h"
25 #include "AliAnalysisManager.h"
26 #include "AliStack.h"
27 #include "AliESDtrackCuts.h"
28 #include "AliESDEvent.h"
29 #include "AliESDInputHandler.h"
30 #include "AliAODEvent.h"
31 #include "AliVEvent.h"
32 #include "AliMCEvent.h"
33 #include "AliVTrack.h"
34 #include "AliVCluster.h"
35 #include "AliEmcalJet.h"
36 #include "AliEMCALGeometry.h"
37 #include "AliEMCALRecoUtils.h"
38 #include "AliVCaloCells.h"
39 #include "AliPicoTrack.h"
40 #include "AliAnalysisTaskEmcal.h"
42 #include "Rtypes.h"
43 
45 
46 //________________________________________________________________________
49 
50  fOutput(0),
51  flTrack(0),
52  flCluster(0),
53  fhTrackPt(0),
54  fhTrackEta(0),
55  fhTrackPhi(0),
56  fhGlobalTrackPt(0),
57  fhGlobalTrackEta(0),
58  fhGlobalTrackPhi(0),
59  fhComplementaryTrackPt(0),
60  fhComplementaryTrackEta(0),
61  fhComplementaryTrackPhi(0),
62  fhClusterPt(0),
63  fhClusterEta(0),
64  fhClusterPhi(0),
65  fhCentrality(0),
66  fhEMCalCellCounts(0),
67 
68  fhChargeAndNeutralEvents(0),
69  fhChargeOnlyEvents(0),
70  fhNeutralOnlyEvents(0),
71  fhNothingEvents(0),
72  fhEMCalChargeAndNeutralEvents(0),
73  fhEMCalChargeOnlyEvents(0),
74  fhEMCalNeutralOnlyEvents(0),
75  fhEMCalNothingEvents(0),
76 
77  fhTrackEtaPhi(0),
78  fhTrackPhiPt(0),
79  fhTrackEtaPt(0),
80  fhGlobalTrackEtaPhi(0),
81  fhGlobalTrackPhiPt(0),
82  fhGlobalTrackEtaPt(0),
83  fhComplementaryTrackEtaPhi(0),
84  fhComplementaryTrackPhiPt(0),
85  fhComplementaryTrackEtaPt(0),
86  fhClusterEtaPhi(0),
87  fhClusterPhiPt(0),
88  fhClusterEtaPt(0),
89  fhEMCalEventMult(0),
90  fhTPCEventMult(0),
91  fhEMCalTrackEventMult(0),
92 
93  fpEMCalEventMult(0),
94  fpTPCEventMult(0),
95 
96  fpTrackPtProfile(0),
97  fpClusterPtProfile(0),
98 
99  fpFullJetEDProfile(0),
100  fpChargedJetEDProfile(0),
101  fpChargedJetEDProfileScaled(0),
102 
103  fTPCRawJets(0),
104  fEMCalRawJets(0),
105  fRhoChargedCMSScale(0),
106  fRhoChargedScale(0),
107  fRhoFull0(0),
108  fRhoFull1(0),
109  fRhoFull2(0),
110  fRhoFullN(0),
111  fRhoFullDijet(0),
112  fRhoFullkT(0),
113  fRhoFullCMS(0),
114  fRhoCharged0(0),
115  fRhoCharged1(0),
116  fRhoCharged2(0),
117  fRhoChargedN(0),
118  fRhoChargedkT(0),
119  fRhoChargedkTScale(0),
120  fRhoChargedCMS(0),
121 
122  fTPCJet(0),
123  fTPCFullJet(0),
124  fTPCOnlyJet(0),
125  fTPCJetUnbiased(0),
126  fTPCkTFullJet(0),
127  fEMCalJet(0),
128  fEMCalFullJet(0),
129  fEMCalPartJet(0),
130  fEMCalPartJetUnbiased(0),
131  fEMCalkTFullJet(0),
132 
133  fIsInitialized(0),
134  fRJET(4),
135  fnEvents(0),
136  fnEventsCharged(0),
137  fnDiJetEvents(0),
138  fEvent(0),
139  fRecoUtil(0),
140  fEMCALGeometry(0),
141  fCells(0),
142  fDoNEF(0),
143  fDoNEFSignalOnly(1),
144  fSignalTrackBias(0),
145  fTrackQA(0),
146  fClusterQA(0),
147  fCalculateRhoJet(0),
148  fDoVertexRCut(0),
149  fMCPartLevel(0),
150  fDoTHnSparse(0),
151  fDoJetRhoDensity(0),
152  fDo3DHistos(0),
153  fEMCalPhiMin(1.39626),
154  fEMCalPhiMax(3.26377),
155  fEMCalPhiTotal(1.86750),
156  fEMCalEtaMin(-0.7),
157  fEMCalEtaMax(0.7),
158  fEMCalEtaTotal(1.4),
159  fEMCalArea(2.61450),
160  fTPCPhiMin(0),
161  fTPCPhiMax(6.28319),
162  fTPCPhiTotal(6.28319),
163  fTPCEtaMin(-0.9),
164  fTPCEtaMax(0.9),
165  fTPCEtaTotal(1.8),
166  fTPCArea(11.30973),
167  fParticlePtLow(0.0),
168  fParticlePtUp(200.0),
169  fParticlePtBins(200),
170  fJetR(0.4),
171  fJetRAccept(0.4),
172  fFullEDJetR(0.6),
173  fChargedEDJetR(0.8),
174  fJetRForRho(0.5),
175  fJetAreaCutFrac(0.6),
176  fJetAreaThreshold(0.30159),
177  fnEMCalCells(12288),
178  fScaleFactor(1.28),
179  fNColl(7),
180  fTrackMinPt(0.15),
181  fClusterMinPt(0.3),
182  fNEFSignalJetCut(0.9),
183  fCentralityTag("V0A"),
184  fCentralityBins(10),
185  fCentralityLow(0),
186  fCentralityUp(100),
187  fEventCentrality(0),
188  fRhoFull(0),
189  fRhoCharged(0),
190  fnTracks(0),
191  fnEMCalTracks(0),
192  fnClusters(0),
193  fnCaloClusters(0),
194  fnAKTFullJets(0),
195  fnAKTChargedJets(0),
196  fnKTFullJets(0),
197  fnKTChargedJets(0),
198  fnBckgClusters(0),
199  fTPCJetThreshold(5),
200  fEMCalJetThreshold(5),
201  fVertexWindow(10),
202  fVertexMaxR(1),
203  fTrackName(0),
204  fClusName(0),
205  fkTChargedName(0),
206  fAkTChargedName(0),
207  fkTFullName(0),
208  fAkTFullName(0),
209  fOrgTracks(0),
210  fOrgClusters(0),
211  fmyAKTFullJets(0),
212  fmyAKTChargedJets(0),
213  fmyKTFullJets(0),
214  fmyKTChargedJets(0),
215  fmyTracks(0),
216  fmyClusters(0),
217  fEMCalRCBckgFluc(0),
218  fTPCRCBckgFluc(0),
219  fEMCalRCBckgFlucSignal(0),
220  fTPCRCBckgFlucSignal(0),
221  fEMCalRCBckgFlucNColl(0),
222  fTPCRCBckgFlucNColl(0)
223 {
224  // Dummy constructor ALWAYS needed for I/O.
225  fVertex[0]=0.0,fVertex[1]=0.0,fVertex[2]=0.0;
226 }
227 
228 //________________________________________________________________________
231 
232  fOutput(0),
233  flTrack(0),
234  flCluster(0),
235  fhTrackPt(0),
236  fhTrackEta(0),
237  fhTrackPhi(0),
238  fhGlobalTrackPt(0),
239  fhGlobalTrackEta(0),
240  fhGlobalTrackPhi(0),
241  fhComplementaryTrackPt(0),
242  fhComplementaryTrackEta(0),
243  fhComplementaryTrackPhi(0),
244  fhClusterPt(0),
245  fhClusterEta(0),
246  fhClusterPhi(0),
247  fhCentrality(0),
248  fhEMCalCellCounts(0),
249 
250  fhChargeAndNeutralEvents(0),
251  fhChargeOnlyEvents(0),
252  fhNeutralOnlyEvents(0),
253  fhNothingEvents(0),
254  fhEMCalChargeAndNeutralEvents(0),
255  fhEMCalChargeOnlyEvents(0),
256  fhEMCalNeutralOnlyEvents(0),
257  fhEMCalNothingEvents(0),
258 
259  fhTrackEtaPhi(0),
260  fhTrackPhiPt(0),
261  fhTrackEtaPt(0),
262  fhGlobalTrackEtaPhi(0),
263  fhGlobalTrackPhiPt(0),
264  fhGlobalTrackEtaPt(0),
265  fhComplementaryTrackEtaPhi(0),
266  fhComplementaryTrackPhiPt(0),
267  fhComplementaryTrackEtaPt(0),
268  fhClusterEtaPhi(0),
269  fhClusterPhiPt(0),
270  fhClusterEtaPt(0),
271  fhEMCalEventMult(0),
272  fhTPCEventMult(0),
273  fhEMCalTrackEventMult(0),
274 
275  fpEMCalEventMult(0),
276  fpTPCEventMult(0),
277 
278  fpTrackPtProfile(0),
279  fpClusterPtProfile(0),
280 
281  fpFullJetEDProfile(0),
282  fpChargedJetEDProfile(0),
283  fpChargedJetEDProfileScaled(0),
284 
285  fTPCRawJets(0),
286  fEMCalRawJets(0),
287  fRhoChargedCMSScale(0),
288  fRhoChargedScale(0),
289  fRhoFull0(0),
290  fRhoFull1(0),
291  fRhoFull2(0),
292  fRhoFullN(0),
293  fRhoFullDijet(0),
294  fRhoFullkT(0),
295  fRhoFullCMS(0),
296  fRhoCharged0(0),
297  fRhoCharged1(0),
298  fRhoCharged2(0),
299  fRhoChargedN(0),
300  fRhoChargedkT(0),
301  fRhoChargedkTScale(0),
302  fRhoChargedCMS(0),
303 
304  fTPCJet(0),
305  fTPCFullJet(0),
306  fTPCOnlyJet(0),
307  fTPCJetUnbiased(0),
308  fTPCkTFullJet(0),
309  fEMCalJet(0),
310  fEMCalFullJet(0),
311  fEMCalPartJet(0),
312  fEMCalPartJetUnbiased(0),
313  fEMCalkTFullJet(0),
314 
315  fIsInitialized(0),
316  fRJET(4),
317  fnEvents(0),
318  fnEventsCharged(0),
319  fnDiJetEvents(0),
320  fEvent(0),
321  fRecoUtil(0),
322  fEMCALGeometry(0),
323  fCells(0),
324  fDoNEF(0),
325  fDoNEFSignalOnly(1),
326  fSignalTrackBias(0),
327  fTrackQA(0),
328  fClusterQA(0),
329  fCalculateRhoJet(0),
330  fDoVertexRCut(0),
331  fMCPartLevel(0),
332  fDoTHnSparse(0),
333  fDoJetRhoDensity(0),
334  fDo3DHistos(0),
335  fEMCalPhiMin(1.39626),
336  fEMCalPhiMax(3.26377),
337  fEMCalPhiTotal(1.86750),
338  fEMCalEtaMin(-0.7),
339  fEMCalEtaMax(0.7),
340  fEMCalEtaTotal(1.4),
341  fEMCalArea(2.61450),
342  fTPCPhiMin(0),
343  fTPCPhiMax(6.28319),
344  fTPCPhiTotal(6.28319),
345  fTPCEtaMin(-0.9),
346  fTPCEtaMax(0.9),
347  fTPCEtaTotal(1.8),
348  fTPCArea(11.30973),
349  fParticlePtLow(0.0),
350  fParticlePtUp(200.0),
351  fParticlePtBins(2000),
352  fJetR(0.4),
353  fJetRAccept(0.4),
354  fFullEDJetR(0.6),
355  fChargedEDJetR(0.8),
356  fJetRForRho(0.5),
357  fJetAreaCutFrac(0.6),
358  fJetAreaThreshold(0.30159),
359  fnEMCalCells(12288),
360  fScaleFactor(1.28),
361  fNColl(7),
362  fTrackMinPt(0.15),
363  fClusterMinPt(0.3),
364  fNEFSignalJetCut(0.9),
365  fCentralityTag("V0A"),
366  fCentralityBins(10),
367  fCentralityLow(0),
368  fCentralityUp(100),
369  fEventCentrality(0),
370  fRhoFull(0),
371  fRhoCharged(0),
372  fnTracks(0),
373  fnEMCalTracks(0),
374  fnClusters(0),
375  fnCaloClusters(0),
376  fnAKTFullJets(0),
377  fnAKTChargedJets(0),
378  fnKTFullJets(0),
379  fnKTChargedJets(0),
380  fnBckgClusters(0),
381  fTPCJetThreshold(5),
382  fEMCalJetThreshold(5),
383  fVertexWindow(10),
384  fVertexMaxR(1),
385  fTrackName(0),
386  fClusName(0),
387  fkTChargedName(0),
388  fAkTChargedName(0),
389  fkTFullName(0),
390  fAkTFullName(0),
391  fOrgTracks(0),
392  fOrgClusters(0),
393  fmyAKTFullJets(0),
394  fmyAKTChargedJets(0),
395  fmyKTFullJets(0),
396  fmyKTChargedJets(0),
397  fmyTracks(0),
398  fmyClusters(0),
399  fEMCalRCBckgFluc(0),
400  fTPCRCBckgFluc(0),
401  fEMCalRCBckgFlucSignal(0),
402  fTPCRCBckgFlucSignal(0),
403  fEMCalRCBckgFlucNColl(0),
404  fTPCRCBckgFlucNColl(0)
405 {
406  // Constructor
407  // Define input and output slots here (never in the dummy constructor)
408  // Input slot #0 works with a TChain - it is connected to the default input container
409  // Output slot #1 writes into a TH1 container
410  fVertex[0]=0.0,fVertex[1]=0.0,fVertex[2]=0.0;
411 
412  DefineOutput(1,TList::Class()); // for output list
413 }
414 
415 //________________________________________________________________________
417 {
418  // Destructor. Clean-up the output list, but not the histograms that are put inside
419  // (the list is owner and will clean-up these histograms). Protect in PROOF case.
420  if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
421  {
422  delete fOutput;
423  }
424 }
425 
426 //________________________________________________________________________
428 {
429  // Create histograms
430  // Called once (on the worker node)
431  fIsInitialized=kFALSE;
432  fOutput = new TList();
433  fOutput->SetOwner(); // IMPORTANT!
434 
435  // Initialize Global Variables
436  fnEvents=0;
437  fnEventsCharged=0;
438  fnDiJetEvents=0;
439 
440  // fRJET=4 -> fJetR=0.4 && fRJET=25 -> fJetR=0.25, but for writing files, should be 4 and 25 respectively
441  if (fRJET>10)
442  {
443  fJetR=(Double_t)fRJET/100.0;
444  }
445  else
446  {
447  fJetR=(Double_t)fRJET/10.0;
448  }
449  fJetRForRho=0.5;
450 
451  fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
452  fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
454  fEMCalEtaMin=-0.7;
455  fEMCalEtaMax=0.7;
458 
459  fTPCPhiMin=(0/(double)360)*2*TMath::Pi();
460  fTPCPhiMax=(360/(double)360)*2*TMath::Pi();
462  fTPCEtaMin=-0.9;
463  fTPCEtaMax=0.9;
466 
467  fParticlePtLow=0.0;
468  fParticlePtUp=200.0;
470 
471  fCentralityBins=10;
472  fCentralityLow=0.0;
473  fCentralityUp=100.0;
474  Int_t CentralityBinMult=10;
475 
476  fJetAreaCutFrac =0.6; // Fudge factor for selecting on jets with threshold Area or higher
478  fTPCJetThreshold=5.0; // Threshold required for an Anti-kT Charged jet to be considered a "true" jet in TPC
479  fEMCalJetThreshold=5.0; // Threshold required for an Anti-kT Charged+Neutral jet to be considered a "true" jet in EMCal
480  fVertexWindow=10.0;
481  fVertexMaxR=1.0;
482 
483  fnBckgClusters=1;
490  for (Int_t i=0;i<fnBckgClusters;i++)
491  {
492  fEMCalRCBckgFluc[i]=0.0;
493  fTPCRCBckgFluc[i]=0.0;
494  fEMCalRCBckgFlucSignal[i]=0.0;
495  fTPCRCBckgFlucSignal[i]=0.0;
496  fEMCalRCBckgFlucNColl[i]=0.0;
497  fTPCRCBckgFlucNColl[i]=0.0;
498  }
499 
500  fnEMCalCells=12288; // sMods 1-10 have 24x48 cells, sMods 11&12 have 8x48 cells...
501 
502  Int_t TCBins=100;
503  Int_t multBins=200;
504  Double_t multLow=0;
505  Double_t multUp=200;
506 
507  // Track QA Plots
508  if (fTrackQA==kTRUE)
509  {
510  flTrack = new TList();
511  flTrack->SetName("TrackQA");
512 
513  // Hybrid Tracks
514  fhTrackPt = new TH1F("fhTrackPt","p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
515  fhTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
516  fhTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
517  fhTrackPt->Sumw2();
518 
519  fhTrackPhi = new TH1F("fhTrackPhi","#varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
520  fhTrackPhi->GetXaxis()->SetTitle("#varphi");
521  fhTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
522  fhTrackPhi->Sumw2();
523 
524  fhTrackEta = new TH1F("fhTrackEta","#eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
525  fhTrackEta->GetXaxis()->SetTitle("#eta");
526  fhTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
527  fhTrackEta->Sumw2();
528 
529  fhTrackEtaPhi = new TH2F("fhTrackEtaPhi","#eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
530  fhTrackEtaPhi->GetXaxis()->SetTitle("#eta");
531  fhTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
532  fhTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
533  fhTrackEtaPhi->Sumw2();
534 
535  fhTrackPhiPt = new TH2F("fhTrackPhiPt","#varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
536  fhTrackPhiPt->GetXaxis()->SetTitle("#varphi");
537  fhTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
538  fhTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
539  fhTrackPhiPt->Sumw2();
540 
541  fhTrackEtaPt = new TH2F("fhTrackEtaPt","#eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
542  fhTrackEtaPt->GetXaxis()->SetTitle("#varphi");
543  fhTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
544  fhTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
545  fhTrackEtaPt->Sumw2();
546 
547  // Global Tracks
548  fhGlobalTrackPt = new TH1F("fhGlobalTrackPt","Global p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
549  fhGlobalTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
550  fhGlobalTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
551  fhGlobalTrackPt->Sumw2();
552 
553  fhGlobalTrackPhi = new TH1F("fhGlobalTrackPhi","Global #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
554  fhGlobalTrackPhi->GetXaxis()->SetTitle("#varphi");
555  fhGlobalTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
556  fhGlobalTrackPhi->Sumw2();
557 
558  fhGlobalTrackEta = new TH1F("fhGlobalTrackEta","Global #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
559  fhGlobalTrackEta->GetXaxis()->SetTitle("#eta");
560  fhGlobalTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
561  fhGlobalTrackEta->Sumw2();
562 
563  fhGlobalTrackEtaPhi = new TH2F("fhGlobalTrackEtaPhi","Global #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
564  fhGlobalTrackEtaPhi->GetXaxis()->SetTitle("#eta");
565  fhGlobalTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
566  fhGlobalTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
567  fhGlobalTrackEtaPhi->Sumw2();
568 
569  fhGlobalTrackPhiPt = new TH2F("fhGlobalTrackPhiPt","Global #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
570  fhGlobalTrackPhiPt->GetXaxis()->SetTitle("#varphi");
571  fhGlobalTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
572  fhGlobalTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
573  fhGlobalTrackPhiPt->Sumw2();
574 
575  fhGlobalTrackEtaPt = new TH2F("fhGlobalTrackEtaPt","Global #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
576  fhGlobalTrackEtaPt->GetXaxis()->SetTitle("#varphi");
577  fhGlobalTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
578  fhGlobalTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
579  fhGlobalTrackEtaPt->Sumw2();
580 
581  // Complementary Tracks
582  fhComplementaryTrackPt = new TH1F("fhComplementaryTrackPt","Complementary p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
583  fhComplementaryTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
584  fhComplementaryTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
585  fhComplementaryTrackPt->Sumw2();
586 
587  fhComplementaryTrackPhi = new TH1F("fhComplementaryTrackPhi","Complementary #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
588  fhComplementaryTrackPhi->GetXaxis()->SetTitle("#varphi");
589  fhComplementaryTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
590  fhComplementaryTrackPhi->Sumw2();
591 
592  fhComplementaryTrackEta = new TH1F("fhComplementaryTrackEta","Complementary #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
593  fhComplementaryTrackEta->GetXaxis()->SetTitle("#eta");
594  fhComplementaryTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
595  fhComplementaryTrackEta->Sumw2();
596 
597  fhComplementaryTrackEtaPhi = new TH2F("fhComplementaryTrackEtaPhi","Complementary #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
598  fhComplementaryTrackEtaPhi->GetXaxis()->SetTitle("#eta");
599  fhComplementaryTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
600  fhComplementaryTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
602 
603  fhComplementaryTrackPhiPt = new TH2F("fhComplementaryTrackPhiPt","Complementary #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
604  fhComplementaryTrackPhiPt->GetXaxis()->SetTitle("#varphi");
605  fhComplementaryTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
606  fhComplementaryTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
607  fhComplementaryTrackPhiPt->Sumw2();
608 
609  fhComplementaryTrackEtaPt = new TH2F("fhComplementaryTrackEtaPt","Complementary #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
610  fhComplementaryTrackEtaPt->GetXaxis()->SetTitle("#varphi");
611  fhComplementaryTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
612  fhComplementaryTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
613  fhComplementaryTrackEtaPt->Sumw2();
614 
615  fhTPCEventMult = new TH2F("fhTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
616  fhTPCEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
617  fhTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
618  fhTPCEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Charged}");
619  fhTPCEventMult->Sumw2();
620 
621  fhEMCalTrackEventMult = new TH2F("fhEMCalTrackEventMult","EMCal Track Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
622  fhEMCalTrackEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
623  fhEMCalTrackEventMult->GetYaxis()->SetTitle("Multiplicity");
624  fhEMCalTrackEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Neutral}");
625  fhEMCalTrackEventMult->Sumw2();
626 
627  fpTPCEventMult = new TProfile("fpTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
628  fpTPCEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
629  fpTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
630 
631  // QA::2D Energy Density Profiles for Tracks and Clusters
632  fpTrackPtProfile = new TProfile2D("fpTrackPtProfile","2D Profile of track pT density throughout the TPC",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
633  fpTrackPtProfile->GetXaxis()->SetTitle("#eta");
634  fpTrackPtProfile->GetYaxis()->SetTitle("#varphi");
635  fpTrackPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
636 
637  flTrack->Add(fhTrackPt);
638  flTrack->Add(fhTrackEta);
639  flTrack->Add(fhTrackPhi);
640  flTrack->Add(fhTrackEtaPhi);
641  flTrack->Add(fhTrackPhiPt);
642  flTrack->Add(fhTrackEtaPt);
643  flTrack->Add(fhGlobalTrackPt);
655  flTrack->Add(fhTPCEventMult);
657  flTrack->Add(fpTPCEventMult);
659  fOutput->Add(flTrack);
660  }
661 
662  if (fClusterQA==kTRUE)
663  {
664  flCluster = new TList();
665  flCluster->SetName("ClusterQA");
666 
667  fhClusterPt = new TH1F("fhClusterPt","p_{T} distribution of clusters in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
668  fhClusterPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
669  fhClusterPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
670  fhClusterPt->Sumw2();
671 
672  fhClusterPhi = new TH1F("fhClusterPhi","#varphi distribution of clusters in event",TCBins,fTPCPhiMin,fTPCPhiMax);
673  fhClusterPhi->GetXaxis()->SetTitle("#varphi");
674  fhClusterPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
675  fhClusterPhi->Sumw2();
676 
677  fhClusterEta = new TH1F("fhClusterEta","#eta distribution of clusters in event",TCBins,fTPCEtaMin,fTPCEtaMax);
678  fhClusterEta->GetXaxis()->SetTitle("#eta");
679  fhClusterEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
680  fhClusterEta->Sumw2();
681 
682  fhClusterEtaPhi = new TH2F("fhClusterEtaPhi","#eta-#varphi distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
683  fhClusterEtaPhi->GetXaxis()->SetTitle("#eta");
684  fhClusterEtaPhi->GetYaxis()->SetTitle("#varphi");
685  fhClusterEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
686  fhClusterEtaPhi->Sumw2();
687 
688  fhClusterPhiPt = new TH2F("fhClusterPhiPt","#varphi-p_{T} distribution of clusters in event",TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
689  fhClusterPhiPt->GetXaxis()->SetTitle("#varphi");
690  fhClusterPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
691  fhClusterPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
692  fhClusterPhiPt->Sumw2();
693 
694  fhClusterEtaPt = new TH2F("fhClusterEtaPt","#eta-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
695  fhClusterEtaPt->GetXaxis()->SetTitle("#varphi");
696  fhClusterEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
697  fhClusterEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
698  fhClusterEtaPt->Sumw2();
699 
700  fhEMCalEventMult = new TH2F("fhEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
701  fhEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
702  fhEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
703  fhEMCalEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Neutral}");
704  fhEMCalEventMult->Sumw2();
705 
706  fpEMCalEventMult = new TProfile("fpEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
707  fpEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
708  fpEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
709 
710  fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
711  fpClusterPtProfile->GetXaxis()->SetTitle("#eta");
712  fpClusterPtProfile->GetYaxis()->SetTitle("#varphi");
713  fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
714 
715  fhEMCalCellCounts = new TH1F("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
716  fhEMCalCellCounts->GetXaxis()->SetTitle("Absoulute Cell Id");
717  fhEMCalCellCounts->GetYaxis()->SetTitle("Counts per Event");
718  fhEMCalCellCounts->Sumw2();
719 
720  flCluster->Add(fhClusterPt);
721  flCluster->Add(fhClusterEta);
722  flCluster->Add(fhClusterPhi);
730  fOutput->Add(flCluster);
731  }
732 
733  if (fCalculateRhoJet>=0 && fMCPartLevel==kFALSE) // Default Rho & Raw Jet Spectra
734  {
736 
738  fRhoChargedCMSScale->SetSignalTrackPtBias(fSignalTrackBias);
740  fOutput->Add(fRhoChargedCMSScale->GetOutputHistos());
741 
742  }
743  if (fCalculateRhoJet>=1) // Basic Rho & Raw Jet Spectra
744  {
747 
748  if (fMCPartLevel==kTRUE)
749  {
752  fRhoChargedCMS->SetSignalTrackPtBias(fSignalTrackBias);
753 
755  fOutput->Add(fRhoChargedCMS->GetOutputHistos());
756  }
757 
758  if (fDoJetRhoDensity == kTRUE && fDo3DHistos == kTRUE)
759  {
760  Double_t fullEDR = fFullEDJetR *100;
761  Double_t chargedEDR = fChargedEDJetR *100;
762  const Int_t fullEDBins = (Int_t) fullEDR;
763  const Int_t chargedEDBins = (Int_t) chargedEDR;
764 
765  fpFullJetEDProfile = new TProfile3D("fpFullJetEDProfile","Jet ED Profile for #varphi_{jet}#in(#varphi_{EMCal min} + R,#varphi_{EMCal max} - R) & #eta_{jet}#in(#eta_{EMCal min} + R,#eta_{EMCal max} - R)",fParticlePtBins,fParticlePtLow,fParticlePtUp,fCentralityBins,fCentralityLow,fCentralityUp,fullEDBins,0,fFullEDJetR);
766  fpFullJetEDProfile->GetXaxis()->SetTitle("p_{T,jet}^{ch+em} (GeV)");
767  fpFullJetEDProfile->GetYaxis()->SetTitle(fCentralityTag.Data());
768  fpFullJetEDProfile->GetZaxis()->SetTitle("R");
769 
770  fpChargedJetEDProfile = new TProfile3D("fpChargedJetEDProfile","Charged Jet ED Profile for #eta_{jet}#in(#eta_{TPC min} + R,#eta_{TPC max} - R)",fParticlePtBins,fParticlePtLow,fParticlePtUp,fCentralityBins,fCentralityLow,fCentralityUp,chargedEDBins,0,fChargedEDJetR);
771  fpChargedJetEDProfile->GetXaxis()->SetTitle("p_{T,jet}^{ch} (GeV)");
772  fpChargedJetEDProfile->GetYaxis()->SetTitle(fCentralityTag.Data());
773  fpChargedJetEDProfile->GetZaxis()->SetTitle("R");
774 
775  fpChargedJetEDProfileScaled = new TProfile3D("fpChargedJetEDProfileScaled","Charged Jet ED Profile x SF for #eta_{jet}#in(#eta_{TPC min} + R,#eta_{TPC max} - R)",fParticlePtBins,fParticlePtLow,fParticlePtUp,fCentralityBins,fCentralityLow,fCentralityUp,chargedEDBins,0,fChargedEDJetR);
776  fpChargedJetEDProfileScaled->GetXaxis()->SetTitle("p_{T,jet}^{ch} (GeV)");
777  fpChargedJetEDProfileScaled->GetYaxis()->SetTitle(fCentralityTag.Data());
778  fpChargedJetEDProfileScaled->GetZaxis()->SetTitle("R");
779 
783  }
785  }
786  if (fCalculateRhoJet>=2 && fMCPartLevel==kFALSE) // Charged Rho & Raw Jet Spectra
787  {
790  fRhoChargedCMS->SetSignalTrackPtBias(fSignalTrackBias);
791 
793  fOutput->Add(fRhoChargedCMS->GetOutputHistos());
794  }
795  if (fCalculateRhoJet>=3 && fMCPartLevel==kFALSE) // Basic Rho & Raw Jet Spectra
796  {
797  fRhoFull0 = new AlipAJetHistos("RhoFull0",fCentralityTag.Data());
799  fRhoFull1 = new AlipAJetHistos("RhoFull1",fCentralityTag.Data());
801  fRhoFull2 = new AlipAJetHistos("RhoFull2",fCentralityTag.Data());
803  fRhoFullN = new AlipAJetHistos("RhoFullN",fCentralityTag.Data());
805  fRhoFullDijet = new AlipAJetHistos("RhoFullDijet",fCentralityTag.Data());
807  fRhoFullkT = new AlipAJetHistos("RhoFullkT",fCentralityTag.Data());
809  fRhoFullCMS = new AlipAJetHistos("RhoFullCMS",fCentralityTag.Data());
811  fRhoCharged0 = new AlipAJetHistos("RhoCharged0",fCentralityTag.Data());
813  fRhoCharged1 = new AlipAJetHistos("RhoCharged1",fCentralityTag.Data());
815  fRhoCharged2 = new AlipAJetHistos("RhoCharged2",fCentralityTag.Data());
817  fRhoChargedN = new AlipAJetHistos("RhoChargedN",fCentralityTag.Data());
819  fRhoChargedkT = new AlipAJetHistos("RhoChargedkT",fCentralityTag.Data());
821  fRhoChargedkTScale = new AlipAJetHistos("RhoChargedkTScale",fCentralityTag.Data());
823 
837  }
838 
839  fhCentrality = new TH1F("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
840  fhCentrality->GetXaxis()->SetTitle(fCentralityTag.Data());
841  fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
842  fhCentrality->Sumw2();
843 
844  fOutput->Add(fhCentrality);
845 
846  if (fMCPartLevel == kFALSE)
847  {
848  fhChargeAndNeutralEvents = new TH1F("fhChargeAndNeutralEvents","Events with n_{tracks}>0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
849  fhChargeAndNeutralEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
850  fhChargeAndNeutralEvents->GetYaxis()->SetTitle("1/N_{Events}");
851  fhChargeAndNeutralEvents->Sumw2();
852 
853  fhChargeOnlyEvents = new TH1F("fhChargeOnlyEvents","Events with n_{tracks}>0 & n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
854  fhChargeOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
855  fhChargeOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
856  fhChargeOnlyEvents->Sumw2();
857 
858  fhNeutralOnlyEvents = new TH1F("fhNeutralOnlyEvents","Events with n_{tracks}=0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
859  fhNeutralOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
860  fhNeutralOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
861  fhNeutralOnlyEvents->Sumw2();
862 
863  fhNothingEvents = new TH1F("fhNothingEvents","Events with n_{tracks}=n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
864  fhNothingEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
865  fhNothingEvents->GetYaxis()->SetTitle("1/N_{Events}");
866  fhNothingEvents->Sumw2();
867 
868  fhEMCalChargeAndNeutralEvents = new TH1F("fhEMCalChargeAndNeutralEvents","Events with n_{emcal tracks}>0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
869  fhEMCalChargeAndNeutralEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
870  fhEMCalChargeAndNeutralEvents->GetYaxis()->SetTitle("1/N_{Events}");
872 
873  fhEMCalChargeOnlyEvents = new TH1F("fhEMCalChargeOnlyEvents","Events with n_{emcal tracks}>0 & n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
874  fhEMCalChargeOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
875  fhEMCalChargeOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
876  fhEMCalChargeOnlyEvents->Sumw2();
877 
878  fhEMCalNeutralOnlyEvents = new TH1F("fhEMCalNeutralOnlyEvents","Events with n_{tracks}>0 & n_{emcal tracks}=0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
879  fhEMCalNeutralOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
880  fhEMCalNeutralOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
881  fhEMCalNeutralOnlyEvents->Sumw2();
882 
883  fhEMCalNothingEvents = new TH1F("fhEMCalNothingEvents","Events with n_{emcal tracks}=n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
884  fhEMCalNothingEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
885  fhEMCalNothingEvents->GetYaxis()->SetTitle("1/N_{Events}");
886  fhEMCalNothingEvents->Sumw2();
887 
891  fOutput->Add(fhNothingEvents);
896  }
897 
898  // Post data for ALL output slots >0 here,
899  // To get at least an empty histogram
900  // 1 is the outputnumber of a certain weg of task 1
901  PostData(1, fOutput);
902 }
903 
905 {
906  // Get the event tracks
907  fOrgTracks = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fTrackName));
908 
909  // Get the event caloclusters
910  fOrgClusters = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fClusName));
911 
912  // Get charged jets
913  fmyKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fkTChargedName));
914  fmyAKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fAkTChargedName));
915 
916  // Get the full jets
917  fmyKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fkTFullName));
918  fmyAKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fAkTFullName));
919 
920  fIsInitialized=kTRUE;
921 }
922 //________________________________________________________________________
924 {
925  if (fIsInitialized==kFALSE)
926  {
927  UserExecOnce();
928  }
929 
930  // Get pointer to reconstructed event
931  fEvent = InputEvent();
932  if (!fEvent)
933  {
934  AliError("Pointer == 0, this can not happen!");
935  return 0;
936  }
937 
938  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent);
939  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fEvent);
940 
941  fRecoUtil = new AliEMCALRecoUtils();
942  fEMCALGeometry = AliEMCALGeometry::GetInstance();
943 
944  if (esd)
945  {
946  fEventCentrality=esd->GetCentrality()->GetCentralityPercentile(fCentralityTag.Data());
947 
948  if (fDoVertexRCut == kTRUE)
949  {
950  // Vertex R cut not done in AliAnalysisTaskEmcal
951  if (TMath::Sqrt(TMath::Power(esd->GetPrimaryVertex()->GetX(),2)+TMath::Power(esd->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
952  {
953  DeleteJetData(0);
954  return 0;
955  }
956  }
957 
958  esd->GetPrimaryVertex()->GetXYZ(fVertex);
959  fCells = (AliVCaloCells*) esd->GetEMCALCells();
960  fnCaloClusters = esd->GetNumberOfCaloClusters();
961  }
962  else if (aod)
963  {
964  fEventCentrality=aod->GetCentrality()->GetCentralityPercentile(fCentralityTag.Data());
965 
966  if (fDoVertexRCut == kTRUE)
967  {
968  // Vertex R cut not done in AliAnalysisTaskEmcal
969  if (TMath::Sqrt(TMath::Power(aod->GetPrimaryVertex()->GetX(),2)+TMath::Power(aod->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
970  {
971  DeleteJetData(0);
972  return 0;
973  }
974  }
975 
976  aod->GetPrimaryVertex()->GetXYZ(fVertex);
977  fCells = (AliVCaloCells*) aod->GetEMCALCells();
978  fnCaloClusters = aod->GetNumberOfCaloClusters();
979  }
980  else
981  {
982  AliError("Cannot get AOD/ESD event. Rejecting Event");
983  DeleteJetData(0);
984  return 0;
985  }
986 
987  // Make sure Centrality isn't exactly 100% (to avoid bin filling errors in profile plots. Make it 99.99
988  if (fEventCentrality>99.99)
989  {
990  fEventCentrality=99.99;
991  }
992  fhCentrality->Fill(fEventCentrality); // Counts total events
993 
994 
995  // Count Events
996  EventCounts();
997 
998  // Reject any event that doesn't have any tracks, i.e. TPC is off
999  if (fnTracks<1)
1000  {
1001  if (fTrackQA==kTRUE)
1002  {
1003  fhTPCEventMult->Fill(fEventCentrality,0.0);
1004  fpTPCEventMult->Fill(fEventCentrality,0.0);
1006  }
1007  AliWarning("No PicoTracks, Rejecting Event");
1008  DeleteJetData(1);
1009  PostData(1,fOutput);
1010  return 0;
1011  }
1012 
1013  if (fnClusters<1)
1014  {
1015  //AliInfo("No Corrected CaloClusters, using only charged jets");
1016 
1017  if (fTrackQA==kTRUE)
1018  {
1019  TrackHisto();
1020  }
1021  InitChargedJets();
1023 
1024  if (fClusterQA==kTRUE)
1025  {
1028  }
1029 
1030  // Rho's
1031  if (fCalculateRhoJet>=1)
1032  {
1033  if (fDoJetRhoDensity == kTRUE && fDo3DHistos == kTRUE)
1034  {
1036  }
1037  if (fMCPartLevel==kTRUE)
1038  {
1040  }
1041  }
1042  if (fCalculateRhoJet>=2 && fMCPartLevel==kFALSE)
1043  {
1045  }
1046  if (fCalculateRhoJet>=3 && fMCPartLevel==kFALSE)
1047  {
1053  }
1054 
1055  DeleteJetData(2);
1056  fnEventsCharged++;
1057  PostData(1,fOutput);
1058  return 0;
1059  }
1060 
1061  if (fTrackQA==kTRUE)
1062  {
1063  TrackHisto();
1064  }
1065  if (fClusterQA==kTRUE)
1066  {
1067  ClusterHisto();
1068  }
1069 
1070  // Prep the jets
1071  InitChargedJets();
1073 
1074  InitFullJets();
1076 
1077  if (fCalculateRhoJet>=0)
1078  {
1080  }
1081  if (fCalculateRhoJet>=1)
1082  {
1084  if (fDoJetRhoDensity == kTRUE && fDo3DHistos == kTRUE)
1085  {
1088  }
1089  }
1090  if (fCalculateRhoJet>=2)
1091  {
1093  }
1094  if (fCalculateRhoJet>=3)
1095  {
1101 
1102  EstimateFullRho0();
1103  EstimateFullRho1();
1104  EstimateFullRho2();
1105  EstimateFullRhoN();
1108 
1110 
1111  // Dijet
1112  if (IsDiJetEvent()==kTRUE)
1113  {
1115  }
1116  }
1117 
1118  // Delete Dynamic Arrays
1119  DeleteJetData(3);
1120  fnEvents++;
1121  PostData(1,fOutput);
1122  return 1;
1123 }
1124 
1125 //________________________________________________________________________
1126 void AliAnalysisTaskFullpAJets::Terminate(Option_t *) //specify what you want to have done
1127 {
1128  // Called once at the end of the query. Done nothing here.
1129 }
1130 
1134 
1136 {
1137  // Fill a TObjArray with the tracks from a TClonesArray which grabs the picotracks.
1138  Int_t i;
1139  Int_t j=0;
1140 
1141  fmyTracks = new TObjArray();
1142  for (i=0;i<fOrgTracks->GetEntries();i++)
1143  {
1144  AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(i);
1145  if (vtrack->Pt()>=fTrackMinPt)
1146  {
1147  fmyTracks->Add(vtrack);
1148  if (IsInEMCal(vtrack->Phi(),vtrack->Eta()) == kTRUE)
1149  {
1150  j++;
1151  }
1152  }
1153  }
1154  fnTracks = fmyTracks->GetEntries();
1155  fnEMCalTracks = j;
1156 }
1157 
1159 {
1160  // Fill a TObjArray with the clusters from a TClonesArray which grabs the caloclusterscorr.
1161  Int_t i;
1162 
1163  fmyClusters = new TObjArray();
1164  TLorentzVector *cluster_vec = new TLorentzVector;
1165  for (i=0;i<fOrgClusters->GetEntries();i++)
1166  {
1167  AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(i);
1168  vcluster->GetMomentum(*cluster_vec,fVertex);
1169 
1170  if (cluster_vec->Pt()>=fClusterMinPt && vcluster->IsEMCAL()==kTRUE)
1171  {
1172  fmyClusters->Add(vcluster);
1173  }
1174  }
1175  fnClusters = fmyClusters->GetEntries();
1176  delete cluster_vec;
1177 }
1178 
1180 {
1181  TrackCuts();
1182  if (fMCPartLevel == kTRUE)
1183  {
1184  return;
1185  }
1186 
1187  ClusterCuts();
1188  if (fnTracks>0)
1189  {
1190  if (fnClusters>0)
1191  {
1193  }
1194  else
1195  {
1197  }
1198  if (fnEMCalTracks>0)
1199  {
1200  if (fnClusters>0)
1201  {
1203  }
1204  else
1205  {
1207  }
1208  }
1209  else
1210  {
1211  if (fnClusters>0)
1212  {
1214  }
1215  else
1216  {
1218  }
1219  }
1220  }
1221  else
1222  {
1223  if (fnClusters>0)
1224  {
1226  }
1227  else
1228  {
1230  }
1231  }
1232 }
1233 
1235 {
1236  // Fill track histograms: Phi,Eta,Pt
1237  Int_t i,j;
1238  Int_t TCBins=100;
1239  TH2F *hdummypT= new TH2F("hdummypT","",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
1240 
1241  for (i=0;i<fnTracks;i++)
1242  {
1243  AliPicoTrack* vtrack = (AliPicoTrack*) fmyTracks->At(i);
1244  fhTrackPt->Fill(vtrack->Pt());
1245  fhTrackEta->Fill(vtrack->Eta());
1246  fhTrackPhi->Fill(vtrack->Phi());
1247  fhTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1248  fhTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1249  fhTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1250 
1251  // Fill Associated Track Distributions for AOD QA Productions
1252  // Global Tracks
1253  if (vtrack->GetTrackType()==0)
1254  {
1255  fhGlobalTrackPt->Fill(vtrack->Pt());
1256  fhGlobalTrackEta->Fill(vtrack->Eta());
1257  fhGlobalTrackPhi->Fill(vtrack->Phi());
1258  fhGlobalTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1259  fhGlobalTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1260  fhGlobalTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1261  }
1262  // Complementary Tracks
1263  else if (vtrack->GetTrackType()==1)
1264  {
1265  fhComplementaryTrackPt->Fill(vtrack->Pt());
1266  fhComplementaryTrackEta->Fill(vtrack->Eta());
1267  fhComplementaryTrackPhi->Fill(vtrack->Phi());
1268  fhComplementaryTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1269  fhComplementaryTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1270  fhComplementaryTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1271  }
1272  hdummypT->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
1273  }
1274  for (i=1;i<=TCBins;i++)
1275  {
1276  for (j=1;j<=TCBins;j++)
1277  {
1278  fpTrackPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fTPCArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1279  }
1280  }
1281  delete hdummypT;
1282 }
1283 
1285 {
1286  // Fill cluster histograms: Phi,Eta,Pt
1287  Int_t i,j;
1288  Int_t TCBins=100;
1289  TH2F *hdummypT= new TH2F("hdummypT","",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
1290  Int_t myCellID=-2;
1291  TLorentzVector *cluster_vec = new TLorentzVector;
1292 
1293  for (i=0;i<fnClusters;i++)
1294  {
1295  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1296  vcluster->GetMomentum(*cluster_vec,fVertex);
1297 
1298  fhClusterPt->Fill(cluster_vec->Pt());
1299  fhClusterEta->Fill(cluster_vec->Eta());
1300  fhClusterPhi->Fill(cluster_vec->Phi());
1301  fhClusterEtaPhi->Fill(cluster_vec->Eta(),cluster_vec->Phi());
1302  fhClusterPhiPt->Fill(cluster_vec->Phi(),cluster_vec->Pt());
1303  fhClusterEtaPt->Fill(cluster_vec->Eta(),cluster_vec->Pt());
1304  hdummypT->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
1305  fEMCALGeometry->GetAbsCellIdFromEtaPhi(cluster_vec->Eta(),cluster_vec->Phi(),myCellID);
1306  fhEMCalCellCounts->Fill(myCellID);
1307  }
1308  for (i=1;i<=TCBins;i++)
1309  {
1310  for (j=1;j<=TCBins;j++)
1311  {
1312  fpClusterPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fEMCalArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1313  }
1314  }
1315  delete hdummypT;
1316  delete cluster_vec;
1317 }
1318 
1320 {
1321  // Preliminary Jet Placement and Selection Cuts
1322  Int_t i;
1323 
1324  fnAKTChargedJets = fmyAKTChargedJets->GetEntries();
1325  fnKTChargedJets = fmyKTChargedJets->GetEntries();
1326 
1327  fTPCJet = new AlipAJetData("fTPCJet",kFALSE,fnAKTChargedJets);
1328  fTPCFullJet = new AlipAJetData("fTPCFullJet",kFALSE,fnAKTChargedJets);
1329  fTPCOnlyJet = new AlipAJetData("fTPCOnlyJet",kFALSE,fnAKTChargedJets);
1330  fTPCJetUnbiased = new AlipAJetData("fTPCJetUnbiased",kFALSE,fnAKTChargedJets);
1331 
1334  fTPCJet->SetJetR(fJetR);
1348 
1349  // Initialize Jet Data
1350  for (i=0;i<fnAKTChargedJets;i++)
1351  {
1352  AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(i);
1353 
1354  fTPCJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
1355  fTPCFullJet->SetIsJetInArray(IsInTPC(fJetRAccept,myJet->Phi(),myJet->Eta(),kTRUE),i);
1356  fTPCOnlyJet->SetIsJetInArray(IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1357  fTPCJetUnbiased->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
1358  }
1359 
1360  fTPCJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1361  fTPCFullJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1362  fTPCOnlyJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1364 
1365  // kT Jets
1366  fTPCkTFullJet = new AlipAJetData("fTPCkTFullJet",kFALSE,fnKTChargedJets);
1370 
1371  for (i=0;i<fnKTChargedJets;i++)
1372  {
1373  AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(i);
1374  fTPCkTFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
1375  }
1377 
1378  // Raw Charged Jet Spectra
1379  if (fCalculateRhoJet>=2 || (fCalculateRhoJet>=1 && fMCPartLevel==kTRUE))
1380  {
1382  }
1383 }
1384 
1386 {
1387  // Preliminary Jet Placement and Selection Cuts
1388  Int_t i;
1389 
1390  fnAKTFullJets = fmyAKTFullJets->GetEntries();
1391  fnKTFullJets = fmyKTFullJets->GetEntries();
1392 
1393  fEMCalJet = new AlipAJetData("fEMCalJet",kTRUE,fnAKTFullJets);
1394  fEMCalFullJet = new AlipAJetData("fEMCalFullJet",kTRUE,fnAKTFullJets);
1395  fEMCalPartJet = new AlipAJetData("fEMCalPartJet",kTRUE,fnAKTFullJets);
1396  fEMCalPartJetUnbiased = new AlipAJetData("fEMCalPartJetUnbiased",kTRUE,fnAKTFullJets);
1397 
1418 
1419  // Initialize Jet Data
1420  for (i=0;i<fnAKTFullJets;i++)
1421  {
1422  AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
1423 
1424  fEMCalJet->SetIsJetInArray(IsInEMCal(myJet->Phi(),myJet->Eta()),i);
1426  fEMCalPartJet->SetIsJetInArray(IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta()),i);
1428  }
1429  fEMCalJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
1433 
1434  // kT Jets
1435  fEMCalkTFullJet = new AlipAJetData("fEMCalkTFullJet",kTRUE,fnKTFullJets);
1441 
1442  for (i=0;i<fnKTFullJets;i++)
1443  {
1444  AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(i);
1445  fEMCalkTFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1446  }
1448 
1449  // Raw Full Jet Spectra
1450  if (fMCPartLevel==kFALSE)
1451  {
1453  }
1454 }
1455 
1457 {
1458  Int_t i,j;
1459  Double_t E_tracks_total=0.;
1460  TRandom3 u(time(NULL));
1461 
1462  Double_t Eta_Center=0.5*(fTPCEtaMin+fTPCEtaMax);
1463  Double_t Phi_Center=0.5*(fTPCPhiMin+fTPCPhiMax);
1464  Int_t event_mult=0;
1465  Int_t event_track_mult=0;
1466  Int_t clus_mult=0;
1467 
1468  for (i=0;i<fnBckgClusters;i++)
1469  {
1470  fTPCRCBckgFluc[i]=0.0;
1471  fTPCRCBckgFlucSignal[i]=0.0;
1472  fTPCRCBckgFlucNColl[i]=0.0;
1473  }
1474 
1475  TLorentzVector *dummy= new TLorentzVector;
1476  TLorentzVector *temp_jet= new TLorentzVector;
1477  TLorentzVector *track_vec = new TLorentzVector;
1478 
1479  // First, consider the RC with no spatial restrictions
1480  for (j=0;j<fnBckgClusters;j++)
1481  {
1482  E_tracks_total=0.;
1483 
1484  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1485  // Loop over all tracks
1486  for (i=0;i<fnTracks;i++)
1487  {
1488  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1489  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1490  {
1491  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1492  if (dummy->DeltaR(*track_vec)<fJetR)
1493  {
1494  E_tracks_total+=vtrack->Pt();
1495  }
1496  }
1497  }
1498  fTPCRCBckgFlucSignal[j]=E_tracks_total;
1499  }
1500 
1501  // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1502  E_tracks_total=0.0;
1503  if (fTPCJetUnbiased->GetLeadingPt()<0.0)
1504  {
1505  temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1506  }
1507  else
1508  {
1510  myJet->GetMomentum(*temp_jet);
1511  }
1512 
1513  for (j=0;j<fnBckgClusters;j++)
1514  {
1515  event_mult=0;
1516  event_track_mult=0;
1517  clus_mult=0;
1518  E_tracks_total=0.;
1519 
1520  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1521  while (temp_jet->DeltaR(*dummy)<fJetR)
1522  {
1523  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1524  }
1525  // Loop over all tracks
1526  for (i=0;i<fnTracks;i++)
1527  {
1528  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1529  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE) == kTRUE)
1530  {
1531  event_mult++;
1532  if (IsInEMCal(vtrack->Phi(),vtrack->Eta()) == kTRUE)
1533  {
1534  event_track_mult++;
1535  }
1536  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1537  if (dummy->DeltaR(*track_vec)<fJetR)
1538  {
1539  clus_mult++;
1540  E_tracks_total+=vtrack->Pt();
1541  }
1542  }
1543  }
1544  fTPCRCBckgFluc[j]=E_tracks_total;
1545  }
1546  if (fTrackQA==kTRUE)
1547  {
1548  fhTPCEventMult->Fill(fEventCentrality,event_mult);
1549  fpTPCEventMult->Fill(fEventCentrality,event_mult);
1550  fhEMCalTrackEventMult->Fill(fEventCentrality,event_track_mult);
1551  }
1552  if (fCalculateRhoJet>=2 || (fCalculateRhoJet>=1 && fMCPartLevel==kTRUE))
1553  {
1555  }
1556 
1557  // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
1558  Double_t exclusion_prob;
1559  for (j=0;j<fnBckgClusters;j++)
1560  {
1561  exclusion_prob = u.Uniform(0,1);
1562  if (exclusion_prob<(1/fNColl))
1563  {
1565  }
1566  else
1567  {
1569  }
1570  }
1571 
1572  delete dummy;
1573  delete temp_jet;
1574  delete track_vec;
1575 }
1576 
1578 {
1579  Int_t i,j;
1580  Double_t E_tracks_total=0.;
1581  Double_t E_caloclusters_total=0.;
1582  TRandom3 u(time(NULL));
1583 
1584  Double_t Eta_Center=0.5*(fEMCalEtaMin+fEMCalEtaMax);
1585  Double_t Phi_Center=0.5*(fEMCalPhiMin+fEMCalPhiMax);
1586  Int_t event_mult=0;
1587  Int_t clus_mult=0;
1588 
1589  for (i=0;i<fnBckgClusters;i++)
1590  {
1591  fEMCalRCBckgFluc[i]=0.0;
1592  fEMCalRCBckgFlucSignal[i]=0.0;
1593  fEMCalRCBckgFlucNColl[i]=0.0;
1594  }
1595 
1596  TLorentzVector *dummy= new TLorentzVector;
1597  TLorentzVector *temp_jet= new TLorentzVector;
1598  TLorentzVector *track_vec = new TLorentzVector;
1599  TLorentzVector *cluster_vec = new TLorentzVector;
1600 
1601  // First, consider the RC with no spatial restrictions
1602  for (j=0;j<fnBckgClusters;j++)
1603  {
1604  E_tracks_total=0.;
1605  E_caloclusters_total=0.;
1606 
1607  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1608  // Loop over all tracks
1609  for (i=0;i<fnTracks;i++)
1610  {
1611  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1612  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1613  {
1614  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1615  if (dummy->DeltaR(*track_vec)<fJetR)
1616  {
1617  E_tracks_total+=vtrack->Pt();
1618  }
1619  }
1620  }
1621 
1622  // Loop over all caloclusters
1623  for (i=0;i<fnClusters;i++)
1624  {
1625  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1626  vcluster->GetMomentum(*cluster_vec,fVertex);
1627  if (dummy->DeltaR(*cluster_vec)<fJetR)
1628  {
1629  clus_mult++;
1630  E_caloclusters_total+=cluster_vec->Pt();
1631  }
1632  }
1633  fEMCalRCBckgFlucSignal[j]=E_tracks_total+E_caloclusters_total;
1634  }
1635 
1636  // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1637  E_tracks_total=0.;
1638  E_caloclusters_total=0.;
1640  {
1641  temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1642  }
1643  else
1644  {
1646  myJet->GetMomentum(*temp_jet);
1647  }
1648 
1649  for (j=0;j<fnBckgClusters;j++)
1650  {
1651  event_mult=0;
1652  clus_mult=0;
1653  E_tracks_total=0.;
1654  E_caloclusters_total=0.;
1655 
1656  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1657  while (temp_jet->DeltaR(*dummy)<fJetR)
1658  {
1659  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1660  }
1661  // Loop over all tracks
1662  for (i=0;i<fnTracks;i++)
1663  {
1664  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1665  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1666  {
1667  event_mult++;
1668  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1669  if (dummy->DeltaR(*track_vec)<fJetR)
1670  {
1671  clus_mult++;
1672  E_tracks_total+=vtrack->Pt();
1673  }
1674  }
1675  }
1676 
1677  // Loop over all caloclusters
1678  for (i=0;i<fnClusters;i++)
1679  {
1680  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1681  vcluster->GetMomentum(*cluster_vec,fVertex);
1682  event_mult++;
1683  if (dummy->DeltaR(*cluster_vec)<fJetR)
1684  {
1685  clus_mult++;
1686  E_caloclusters_total+=cluster_vec->Pt();
1687  }
1688  }
1689  fEMCalRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
1690  }
1691  if (fClusterQA==kTRUE)
1692  {
1693  fhEMCalEventMult->Fill(fEventCentrality,event_mult);
1694  fpEMCalEventMult->Fill(fEventCentrality,event_mult);
1695  }
1697 
1698  // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
1699  Double_t exclusion_prob;
1700  for (j=0;j<fnBckgClusters;j++)
1701  {
1702  exclusion_prob = u.Uniform(0,1);
1703  if (exclusion_prob<(1/fNColl))
1704  {
1706  }
1707  else
1708  {
1710  }
1711  }
1712 
1713  delete dummy;
1714  delete temp_jet;
1715  delete track_vec;
1716  delete cluster_vec;
1717 }
1718 
1719 // Charged Rho's
1721 {
1722  Int_t i;
1723  Double_t E_tracks_total=0.0;
1724  Double_t TPC_rho=0.0;
1725 
1726  // Loop over all tracks
1727  for (i=0;i<fnTracks;i++)
1728  {
1729  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1730  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1731  {
1732  E_tracks_total+=vtrack->Pt();
1733  }
1734  }
1735 
1736  // Calculate the mean Background density
1737  TPC_rho=E_tracks_total/fTPCArea;
1738 
1739  fRhoCharged = TPC_rho;
1740 
1741  // Fill Histograms
1749 
1750 }
1751 
1753 {
1754  Int_t i;
1755  Double_t E_tracks_total=0.0;
1756  Double_t TPC_rho=0.;
1757 
1758  TLorentzVector *temp_jet= new TLorentzVector;
1759  TLorentzVector *track_vec = new TLorentzVector;
1760 
1762  {
1764  myJet->GetMomentum(*temp_jet);
1765 
1766  // Loop over all tracks
1767  for (i=0;i<fnTracks;i++)
1768  {
1769  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1770  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1771  {
1772  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1773  if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
1774  {
1775  E_tracks_total+=vtrack->Pt();
1776  }
1777  }
1778  }
1779 
1780  // Calculate the mean Background density
1781  TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1782  }
1783  else // i.e. No signal jets -> same as total background density
1784  {
1785  // Loop over all tracks
1786  for (i=0;i<fnTracks;i++)
1787  {
1788  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1789  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1790  {
1791  E_tracks_total+=vtrack->Pt();
1792  }
1793  }
1794  // Calculate the mean Background density
1795  TPC_rho=E_tracks_total/fTPCArea;
1796  }
1797  delete track_vec;
1798  delete temp_jet;
1799 
1800  // Fill histograms
1808 }
1809 
1811 {
1812  Int_t i;
1813  Double_t E_tracks_total=0.0;
1814  Double_t TPC_rho=0.;
1815 
1816  TLorentzVector *temp_jet1= new TLorentzVector;
1817  TLorentzVector *temp_jet2= new TLorentzVector;
1818  TLorentzVector *track_vec = new TLorentzVector;
1819 
1821  {
1823  myhJet->GetMomentum(*temp_jet1);
1824 
1826  myJet->GetMomentum(*temp_jet2);
1827 
1828  // Loop over all tracks
1829  for (i=0;i<fnTracks;i++)
1830  {
1831  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1832  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1833  {
1834  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1835  if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
1836  {
1837  E_tracks_total+=vtrack->Pt();
1838  }
1839  }
1840  }
1841 
1842  // Calculate the mean Background density
1843  TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myhJet->Eta())-AreaWithinTPC(fJetR,myJet->Eta()));
1844  }
1846  {
1848  myJet->GetMomentum(*temp_jet1);
1849 
1850  // Loop over all tracks
1851  for (i=0;i<fnTracks;i++)
1852  {
1853  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1854  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1855  {
1856  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1857  if (temp_jet1->DeltaR(*track_vec)>fJetRForRho)
1858  {
1859  E_tracks_total+=vtrack->Pt();
1860  }
1861  }
1862  }
1863 
1864  // Calculate the mean Background density
1865  TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1866  }
1867  else // i.e. No signal jets -> same as total background density
1868  {
1869  // Loop over all tracks
1870  for (i=0;i<fnTracks;i++)
1871  {
1872  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1873  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1874  {
1875  E_tracks_total+=vtrack->Pt();
1876  }
1877  }
1878 
1879  // Calculate the mean Background density
1880  TPC_rho=E_tracks_total/fTPCArea;
1881  }
1882  delete temp_jet1;
1883  delete temp_jet2;
1884  delete track_vec;
1885 
1886  // Fill histograms
1894 }
1895 
1897 {
1898  Int_t i,j;
1899  Bool_t track_away_from_jet;
1900  Double_t E_tracks_total=0.0;
1901  Double_t TPC_rho=0.0;
1902  Double_t jet_area_total=0.0;
1903 
1904  TLorentzVector *jet_vec= new TLorentzVector;
1905  TLorentzVector *track_vec = new TLorentzVector;
1906 
1907  // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1908  for (i=0;i<fnTracks;i++)
1909  {
1910  // First, check if track is in the EMCal!!
1911  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1912  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1913  {
1915  {
1916  E_tracks_total+=vtrack->Pt();
1917  }
1918  else
1919  {
1920  track_away_from_jet=kTRUE;
1921  j=0;
1922  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1923  while (track_away_from_jet==kTRUE && j<fTPCJetUnbiased->GetTotalSignalJets())
1924  {
1926  myJet->GetMomentum(*jet_vec);
1927  if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
1928  {
1929  track_away_from_jet=kFALSE;
1930  }
1931  j++;
1932  }
1933  if (track_away_from_jet==kTRUE)
1934  {
1935  E_tracks_total+=vtrack->Pt();
1936  }
1937  }
1938  }
1939  }
1940 
1941  // Determine area of all Jets that are within the EMCal
1943  {
1944  jet_area_total=0.0;
1945  }
1946  else
1947  {
1948  for (i=0;i<fTPCJetUnbiased->GetTotalSignalJets();i++)
1949  {
1951  jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
1952  }
1953  }
1954  delete jet_vec;
1955  delete track_vec;
1956 
1957  // Calculate Rho
1958  TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
1959 
1960  // Fill Histogram
1968 
1969 }
1970 
1972 {
1973  Int_t i,j;
1974  Bool_t track_away_from_jet;
1975  Double_t E_tracks_total=0.0;
1976  Double_t TPC_rho=0.0;
1977  Double_t jet_area_total=0.0;
1978 
1979  TLorentzVector *jet_vec= new TLorentzVector;
1980  TLorentzVector *track_vec = new TLorentzVector;
1981 
1982  // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1983  for (i=0;i<fnTracks;i++)
1984  {
1985  // First, check if track is in the EMCal!!
1986  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1987  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1988  {
1990  {
1991  E_tracks_total+=vtrack->Pt();
1992  }
1993  else
1994  {
1995  track_away_from_jet=kTRUE;
1996  j=0;
1997  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1998  while (track_away_from_jet==kTRUE && j<fTPCJetUnbiased->GetTotalSignalJets())
1999  {
2001  myJet->GetMomentum(*jet_vec);
2002  if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
2003  {
2004  track_away_from_jet=kFALSE;
2005  }
2006  j++;
2007  }
2008  if (track_away_from_jet==kTRUE)
2009  {
2010  E_tracks_total+=vtrack->Pt();
2011  }
2012  }
2013  }
2014  }
2015 
2016  // Determine area of all Jets that are within the TPC
2018  {
2019  jet_area_total=0.0;
2020  }
2021  else
2022  {
2023  for (i=0;i<fTPCJetUnbiased->GetTotalSignalJets();i++)
2024  {
2026  jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
2027  }
2028  }
2029  delete jet_vec;
2030  delete track_vec;
2031 
2032  // Calculate Rho
2033  TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
2034  TPC_rho*=fScaleFactor;
2035 
2036  // Fill Histogram
2045 }
2046 
2048 {
2049  Int_t i;
2050  Double_t kTRho = 0.0;
2051  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2052  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2053 
2054  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2055  {
2057  pTArray[i]=myJet->Pt();
2058  RhoArray[i]=myJet->Pt()/myJet->Area();
2059  }
2060 
2061  if (fTPCkTFullJet->GetTotalJets()>=2)
2062  {
2063  kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
2064 
2072  }
2073  delete [] RhoArray;
2074  delete [] pTArray;
2075 }
2076 
2078 {
2079  Int_t i;
2080  Double_t kTRho = 0.0;
2081  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2082  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2083 
2084  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2085  {
2087  pTArray[i]=myJet->Pt();
2088  RhoArray[i]=myJet->Pt()/myJet->Area();
2089  }
2090 
2091  if (fTPCkTFullJet->GetTotalJets()>=2)
2092  {
2093  kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
2094  kTRho*=fScaleFactor;
2095 
2103  }
2104  delete [] RhoArray;
2105  delete [] pTArray;
2106 }
2107 
2109 {
2110  Int_t i,k;
2111  Double_t kTRho = 0.0;
2112  Double_t CMSTotalkTArea = 0.0;
2113  Double_t CMSTrackArea = 0.0;
2114  Double_t CMSCorrectionFactor = 1.0;
2115  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2116  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2117 
2118  k=0;
2120  {
2123 
2124  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2125  {
2127 
2128  CMSTotalkTArea+=myJet->Area();
2129  if (myJet->GetNumberOfTracks()>0)
2130  {
2131  CMSTrackArea+=myJet->Area();
2132  }
2133  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
2134  {
2135  pTArray[k]=myJet->Pt();
2136  RhoArray[k]=myJet->Pt()/myJet->Area();
2137  k++;
2138  }
2139  }
2140  if (k>0)
2141  {
2142  kTRho=MedianRhokT(pTArray,RhoArray,k);
2143  }
2144  else
2145  {
2146  kTRho=0.0;
2147  }
2148  }
2149  else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
2150  {
2152 
2153  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2154  {
2156 
2157  CMSTotalkTArea+=myJet->Area();
2158  if (myJet->GetNumberOfTracks()>0)
2159  {
2160  CMSTrackArea+=myJet->Area();
2161  }
2162  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
2163  {
2164  pTArray[k]=myJet->Pt();
2165  RhoArray[k]=myJet->Pt()/myJet->Area();
2166  k++;
2167  }
2168  }
2169  if (k>0)
2170  {
2171  kTRho=MedianRhokT(pTArray,RhoArray,k);
2172  }
2173  else
2174  {
2175  kTRho=0.0;
2176  }
2177  }
2178  else
2179  {
2180  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2181  {
2183 
2184  CMSTotalkTArea+=myJet->Area();
2185  if (myJet->GetNumberOfTracks()>0)
2186  {
2187  CMSTrackArea+=myJet->Area();
2188  }
2189  pTArray[k]=myJet->Pt();
2190  RhoArray[k]=myJet->Pt()/myJet->Area();
2191  k++;
2192  }
2193  if (k>0)
2194  {
2195  kTRho=MedianRhokT(pTArray,RhoArray,k);
2196  }
2197  else
2198  {
2199  kTRho=0.0;
2200  }
2201  }
2202  // Scale CMS Rho by Correction factor
2203  if (CMSTotalkTArea==0.0)
2204  {
2205  CMSCorrectionFactor = 1.0;
2206  }
2207  else
2208  {
2209  //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2210  CMSCorrectionFactor = CMSTrackArea/(fTPCPhiTotal*(fTPCEtaTotal-2*fJetR)); // The total physical area should be reduced by the eta cut due to looping over only fully contained kT jets within the TPC
2211  }
2212  kTRho*=CMSCorrectionFactor;
2220  delete [] RhoArray;
2221  delete [] pTArray;
2222 }
2223 
2225 {
2226  Int_t i,k;
2227  Double_t kTRho = 0.0;
2228  Double_t CMSTotalkTArea = 0.0;
2229  Double_t CMSTrackArea = 0.0;
2230  Double_t CMSCorrectionFactor = 1.0;
2231  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2232  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2233 
2234  k=0;
2236  {
2239 
2240  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2241  {
2243 
2244  CMSTotalkTArea+=myJet->Area();
2245  if (myJet->GetNumberOfTracks()>0)
2246  {
2247  CMSTrackArea+=myJet->Area();
2248  }
2249  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
2250  {
2251  pTArray[k]=myJet->Pt();
2252  RhoArray[k]=myJet->Pt()/myJet->Area();
2253  k++;
2254  }
2255  }
2256  if (k>0)
2257  {
2258  kTRho=MedianRhokT(pTArray,RhoArray,k);
2259  }
2260  else
2261  {
2262  kTRho=0.0;
2263  }
2264  }
2265  else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
2266  {
2268 
2269  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2270  {
2272 
2273  CMSTotalkTArea+=myJet->Area();
2274  if (myJet->GetNumberOfTracks()>0)
2275  {
2276  CMSTrackArea+=myJet->Area();
2277  }
2278  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
2279  {
2280  pTArray[k]=myJet->Pt();
2281  RhoArray[k]=myJet->Pt()/myJet->Area();
2282  k++;
2283  }
2284  }
2285  if (k>0)
2286  {
2287  kTRho=MedianRhokT(pTArray,RhoArray,k);
2288  }
2289  else
2290  {
2291  kTRho=0.0;
2292  }
2293  }
2294  else
2295  {
2296  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2297  {
2299 
2300  CMSTotalkTArea+=myJet->Area();
2301  if (myJet->GetNumberOfTracks()>0)
2302  {
2303  CMSTrackArea+=myJet->Area();
2304  }
2305  pTArray[k]=myJet->Pt();
2306  RhoArray[k]=myJet->Pt()/myJet->Area();
2307  k++;
2308  }
2309  if (k>0)
2310  {
2311  kTRho=MedianRhokT(pTArray,RhoArray,k);
2312  }
2313  else
2314  {
2315  kTRho=0.0;
2316  }
2317  }
2318  kTRho*=fScaleFactor;
2319  // Scale CMS Rho by Correction factor
2320  if (CMSTotalkTArea==0.0)
2321  {
2322  CMSCorrectionFactor = 1.0;
2323  }
2324  else
2325  {
2326  CMSCorrectionFactor = CMSTrackArea/(fTPCPhiTotal*(fTPCEtaTotal-2*fJetR)); // The total physical area should be reduced by the eta cut due to looping over only fully contained kT jets within the TPC
2327  }
2328  kTRho*=CMSCorrectionFactor;
2329 
2340 
2341  delete [] RhoArray;
2342  delete [] pTArray;
2343 }
2344 
2345 // Full Rho's
2347 {
2348  Int_t i;
2349  Double_t E_tracks_total=0.0;
2350  Double_t E_caloclusters_total=0.0;
2351  Double_t EMCal_rho=0.0;
2352 
2353  TLorentzVector *cluster_vec = new TLorentzVector;
2354 
2355  // Loop over all tracks
2356  for (i=0;i<fnTracks;i++)
2357  {
2358  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2359  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2360  {
2361  E_tracks_total+=vtrack->Pt();
2362  }
2363  }
2364 
2365  // Loop over all caloclusters
2366  for (i=0;i<fnClusters;i++)
2367  {
2368  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2369  vcluster->GetMomentum(*cluster_vec,fVertex);
2370  E_caloclusters_total+=cluster_vec->Pt();
2371  }
2372  delete cluster_vec;
2373 
2374  // Calculate the mean Background density
2375  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2376  fRhoFull=EMCal_rho;
2377 
2378  // Fill Histograms
2379  fRhoFull0->FillRho(fEventCentrality,EMCal_rho);
2386 }
2387 
2389 {
2390  Int_t i;
2391  Double_t E_tracks_total=0.0;
2392  Double_t E_caloclusters_total=0.0;
2393  Double_t EMCal_rho=0.0;
2394 
2395  TLorentzVector *temp_jet= new TLorentzVector;
2396  TLorentzVector *track_vec = new TLorentzVector;
2397  TLorentzVector *cluster_vec = new TLorentzVector;
2398 
2400  {
2402  myJet->GetMomentum(*temp_jet);
2403 
2404  // Loop over all tracks
2405  for (i=0;i<fnTracks;i++)
2406  {
2407  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2408  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2409  {
2410  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2411  if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
2412  {
2413  E_tracks_total+=vtrack->Pt();
2414  }
2415  }
2416  }
2417 
2418  // Loop over all caloclusters
2419  for (i=0;i<fnClusters;i++)
2420  {
2421  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2422  vcluster->GetMomentum(*cluster_vec,fVertex);
2423  if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
2424  {
2425  E_caloclusters_total+=cluster_vec->Pt();
2426  }
2427  }
2428  // Calculate the mean Background density
2429  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2430  }
2431  else // i.e. No signal jets -> same as total background density
2432  {
2433  // Loop over all tracks
2434  for (i=0;i<fnTracks;i++)
2435  {
2436  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2437  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2438  {
2439  E_tracks_total+=vtrack->Pt();
2440  }
2441  }
2442 
2443  // Loop over all caloclusters
2444  for (i=0;i<fnClusters;i++)
2445  {
2446  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2447  vcluster->GetMomentum(*cluster_vec,fVertex);
2448  E_caloclusters_total+=cluster_vec->Pt();
2449  }
2450  // Calculate the mean Background density
2451  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2452  }
2453  delete temp_jet;
2454  delete track_vec;
2455  delete cluster_vec;
2456 
2457  // Fill histograms
2458  fRhoFull1->FillRho(fEventCentrality,EMCal_rho);
2465 }
2466 
2468 {
2469  Int_t i;
2470  Double_t E_tracks_total=0.0;
2471  Double_t E_caloclusters_total=0.0;
2472  Double_t EMCal_rho=0.0;
2473 
2474  TLorentzVector *temp_jet1 = new TLorentzVector;
2475  TLorentzVector *temp_jet2 = new TLorentzVector;
2476  TLorentzVector *track_vec = new TLorentzVector;
2477  TLorentzVector *cluster_vec = new TLorentzVector;
2478 
2480  {
2482  myhJet->GetMomentum(*temp_jet1);
2483 
2485  myJet->GetMomentum(*temp_jet2);
2486 
2487  // Loop over all tracks
2488  for (i=0;i<fnTracks;i++)
2489  {
2490  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2491  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2492  {
2493  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2494  if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
2495  {
2496  E_tracks_total+=vtrack->Pt();
2497  }
2498  }
2499  }
2500 
2501  // Loop over all caloclusters
2502  for (i=0;i<fnClusters;i++)
2503  {
2504  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2505  vcluster->GetMomentum(*cluster_vec,fVertex);
2506  if ((temp_jet1->DeltaR(*cluster_vec)>fJetRForRho) && (temp_jet2->DeltaR(*cluster_vec)>fJetRForRho))
2507  {
2508  E_caloclusters_total+=cluster_vec->Pt();
2509  }
2510  }
2511 
2512  // Calculate the mean Background density
2513  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myhJet->Phi(),myhJet->Eta())-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2514  }
2516  {
2518  myJet->GetMomentum(*temp_jet1);
2519 
2520  // Loop over all tracks
2521  for (i=0;i<fnTracks;i++)
2522  {
2523  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2524  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2525  {
2526  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2527  if (temp_jet1->DeltaR(*track_vec)>fJetRForRho)
2528  {
2529  E_tracks_total+=vtrack->Pt();
2530  }
2531  }
2532  }
2533 
2534  // Loop over all caloclusters
2535  for (i=0;i<fnClusters;i++)
2536  {
2537  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2538  vcluster->GetMomentum(*cluster_vec,fVertex);
2539  if (temp_jet1->DeltaR(*cluster_vec)>fJetRForRho)
2540  {
2541  E_caloclusters_total+=cluster_vec->Pt();
2542  }
2543  }
2544  // Calculate the mean Background density
2545  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2546  }
2547  else // i.e. No signal jets -> same as total background density
2548  {
2549  // Loop over all tracks
2550  for (i=0;i<fnTracks;i++)
2551  {
2552  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2553  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2554  {
2555  E_tracks_total+=vtrack->Pt();
2556  }
2557  }
2558 
2559  // Loop over all caloclusters
2560  for (i=0;i<fnClusters;i++)
2561  {
2562  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2563  vcluster->GetMomentum(*cluster_vec,fVertex);
2564  E_caloclusters_total+=cluster_vec->Pt();
2565  }
2566  // Calculate the mean Background density
2567  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2568  }
2569  delete temp_jet1;
2570  delete temp_jet2;
2571  delete track_vec;
2572  delete cluster_vec;
2573 
2574  // Fill histograms
2575  fRhoFull2->FillRho(fEventCentrality,EMCal_rho);
2582 }
2583 
2585 {
2586  Int_t i,j;
2587  Bool_t track_away_from_jet;
2588  Bool_t cluster_away_from_jet;
2589  Double_t E_tracks_total=0.0;
2590  Double_t E_caloclusters_total=0.0;
2591  Double_t EMCal_rho=0.0;
2592  Double_t jet_area_total=0.0;
2593 
2594  TLorentzVector *jet_vec= new TLorentzVector;
2595  TLorentzVector *track_vec = new TLorentzVector;
2596  TLorentzVector *cluster_vec = new TLorentzVector;
2597 
2598  // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
2599  for (i=0;i<fnTracks;i++)
2600  {
2601  // First, check if track is in the EMCal!!
2602  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
2603  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2604  {
2606  {
2607  E_tracks_total+=vtrack->Pt();
2608  }
2609  else
2610  {
2611  track_away_from_jet=kTRUE;
2612  j=0;
2613  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2614  while (track_away_from_jet==kTRUE && j<fEMCalPartJetUnbiased->GetTotalSignalJets())
2615  {
2617  myJet->GetMomentum(*jet_vec);
2618  if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
2619  {
2620  track_away_from_jet=kFALSE;
2621  }
2622  j++;
2623  }
2624  if (track_away_from_jet==kTRUE)
2625  {
2626  E_tracks_total+=vtrack->Pt();
2627  }
2628  }
2629  }
2630  }
2631 
2632  // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
2633  for (i=0;i<fnClusters;i++)
2634  {
2635  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2636  vcluster->GetMomentum(*cluster_vec,fVertex);
2638  {
2639  E_caloclusters_total+=cluster_vec->Pt();
2640  }
2641  else
2642  {
2643  cluster_away_from_jet=kTRUE;
2644  j=0;
2645 
2646  while (cluster_away_from_jet==kTRUE && j<fEMCalPartJetUnbiased->GetTotalSignalJets())
2647  {
2649  myJet->GetMomentum(*jet_vec);
2650  if (cluster_vec->DeltaR(*jet_vec)<=fJetRForRho)
2651  {
2652  cluster_away_from_jet=kFALSE;
2653  }
2654  j++;
2655  }
2656  if (cluster_away_from_jet==kTRUE)
2657  {
2658  E_caloclusters_total+=cluster_vec->Pt();
2659  }
2660  }
2661  }
2662 
2663  // Determine area of all Jets that are within the EMCal
2665  {
2666  jet_area_total=0.0;
2667  }
2668  else
2669  {
2670  for (i=0;i<fEMCalPartJetUnbiased->GetTotalSignalJets();i++)
2671  {
2673  jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
2674  }
2675  }
2676  delete jet_vec;
2677  delete track_vec;
2678  delete cluster_vec;
2679 
2680  // Calculate Rho
2681  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
2682 
2683  // Fill Histogram
2684  fRhoFullN->FillRho(fEventCentrality,EMCal_rho);
2691 }
2692 
2694 {
2695  Int_t i;
2696  Double_t E_tracks_total=0.0;
2697  Double_t E_caloclusters_total=0.0;
2698  Double_t EMCal_rho=0.0;
2699 
2700  TLorentzVector *cluster_vec = new TLorentzVector;
2701 
2702  // Loop over all tracks
2703  for (i=0;i<fnTracks;i++)
2704  {
2705  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2706  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2707  {
2708  E_tracks_total+=vtrack->Pt();
2709  }
2710  }
2711 
2712  // Loop over all caloclusters
2713  for (i=0;i<fnClusters;i++)
2714  {
2715  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2716  vcluster->GetMomentum(*cluster_vec,fVertex);
2717  E_caloclusters_total+=cluster_vec->Pt();
2718  }
2719 
2720  delete cluster_vec;
2721 
2722  // Calculate the mean Background density
2723  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2724 
2725  // Fill Histograms
2733 }
2734 
2736 {
2737  Int_t i;
2738  Double_t kTRho = 0.0;
2739  Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2740  Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2741 
2742  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2743  {
2745  pTArray[i]=myJet->Pt();
2746  RhoArray[i]=myJet->Pt()/myJet->Area();
2747  }
2748 
2749  if (fEMCalkTFullJet->GetTotalJets()>0)
2750  {
2751  kTRho=MedianRhokT(pTArray,RhoArray,fEMCalkTFullJet->GetTotalJets());
2752  }
2753  else
2754  {
2755  kTRho=0.0;
2756  }
2764  delete [] RhoArray;
2765  delete [] pTArray;
2766 }
2767 
2769 {
2770  Int_t i,k;
2771  Double_t kTRho = 0.0;
2772  Double_t CMSTotalkTArea = 0.0;
2773  Double_t CMSParticleArea = 0.0;
2774  Double_t CMSCorrectionFactor = 1.0;
2775  Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2776  Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2777 
2778  k=0;
2780  {
2783 
2784  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2785  {
2787 
2788  CMSTotalkTArea+=myJet->Area();
2789  if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2790  {
2791  CMSParticleArea+=myJet->Area();
2792  }
2793  if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kTRUE)
2794  {
2795  pTArray[k]=myJet->Pt();
2796  RhoArray[k]=myJet->Pt()/myJet->Area();
2797  k++;
2798  }
2799  }
2800  if (k>0)
2801  {
2802  kTRho=MedianRhokT(pTArray,RhoArray,k);
2803  }
2804  else
2805  {
2806  kTRho=0.0;
2807  }
2808  }
2810  {
2812 
2813  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2814  {
2816 
2817  CMSTotalkTArea+=myJet->Area();
2818  if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2819  {
2820  CMSParticleArea+=myJet->Area();
2821  }
2822  if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE)
2823  {
2824  pTArray[k]=myJet->Pt();
2825  RhoArray[k]=myJet->Pt()/myJet->Area();
2826  k++;
2827  }
2828  }
2829  if (k>0)
2830  {
2831  kTRho=MedianRhokT(pTArray,RhoArray,k);
2832  }
2833  else
2834  {
2835  kTRho=0.0;
2836  }
2837  }
2838  else
2839  {
2840  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2841  {
2843 
2844  CMSTotalkTArea+=myJet->Area();
2845  if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2846  {
2847  CMSParticleArea+=myJet->Area();
2848  }
2849  pTArray[k]=myJet->Pt();
2850  RhoArray[k]=myJet->Pt()/myJet->Area();
2851  k++;
2852  }
2853  if (k>0)
2854  {
2855  kTRho=MedianRhokT(pTArray,RhoArray,k);
2856  }
2857  else
2858  {
2859  kTRho=0.0;
2860  }
2861  }
2862  // Scale CMS Rho by Correction factor
2863  if (CMSTotalkTArea==0.0)
2864  {
2865  CMSCorrectionFactor = 1.0;
2866  }
2867  else
2868  {
2869  //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2870  CMSCorrectionFactor = CMSParticleArea/((fEMCalPhiTotal-2*fJetR)*(fEMCalEtaTotal-2*fJetR)); // The total physical area should be reduced by the eta & phi cuts due to looping over only fully contained kT jets within the EMCal
2871  }
2872  kTRho*=CMSCorrectionFactor;
2873 
2881  delete [] RhoArray;
2882  delete [] pTArray;
2883 }
2884 
2886 {
2887  Int_t i,j;
2888  Double_t delta_R;
2889  Double_t dR=0.0;
2890  Double_t fullEDR = fFullEDJetR*100;
2891  const Int_t fullEDBins = (Int_t) fullEDR;
2892  Double_t ED_pT[fullEDBins];
2893  dR = fFullEDJetR/Double_t(fullEDBins); // Should be 0.01 be default
2894 
2895  TLorentzVector *jet_vec= new TLorentzVector;
2896  TLorentzVector *track_vec = new TLorentzVector;
2897  TLorentzVector *cluster_vec = new TLorentzVector;
2898 
2899  for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
2900  {
2902 
2903  if (IsInEMCalFull(fFullEDJetR,myJet->Phi(),myJet->Eta())==kTRUE)
2904  {
2905  for (j=0;j<fullEDBins;j++)
2906  {
2907  ED_pT[j]=0;
2908  }
2909  myJet->GetMomentum(*jet_vec);
2910 
2911  // Sum all tracks in concentric rings around jet vertex
2912  for (j=0;j<fnTracks;j++)
2913  {
2914  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2915  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2916  delta_R=jet_vec->DeltaR(*track_vec);
2917  if (delta_R<fFullEDJetR)
2918  {
2919  ED_pT[TMath::FloorNint((delta_R/fFullEDJetR)*fullEDBins)]+=vtrack->Pt();
2920  }
2921  }
2922 
2923  // Sum all clusters in concentric rings around jet vertex
2924  for (j=0;j<fnClusters;j++)
2925  {
2926  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2927  vcluster->GetMomentum(*cluster_vec,fVertex);
2928  delta_R=jet_vec->DeltaR(*cluster_vec);
2929  if (delta_R<fFullEDJetR)
2930  {
2931  ED_pT[TMath::FloorNint((delta_R/fFullEDJetR)*fullEDBins)]+=cluster_vec->Pt();
2932  }
2933  }
2934 
2935  for (j=0;j<fullEDBins;j++)
2936  {
2937  ED_pT[j] /= TMath::Pi()*dR*dR*(2*j+1);
2938  fpFullJetEDProfile->Fill(myJet->Pt(),fEventCentrality,j*dR,ED_pT[j]);
2939  }
2940  }
2941  }
2942  delete cluster_vec;
2943  delete track_vec;
2944  delete jet_vec;
2945 }
2946 
2947 // Although this calculates the charged pT density in concentric rings around a charged jet contained fudically within the TPC (R=0.8), we actually scale the density by the scale factor (fScaleFactor) in order to compare to 'Full Charge Density' and 'Rho Median Occupancy Approach' (RhoChargedCMSScaled). We still keep the charged value for completeness.
2949 {
2950  Int_t i,j;
2951  Double_t delta_R;
2952  Double_t dR=0.0;
2953  Double_t chargedEDR = fChargedEDJetR*100;
2954  const Int_t chargedEDBins = (Int_t) chargedEDR;
2955  Double_t ED_pT[chargedEDBins];
2956  dR = fChargedEDJetR/Double_t(chargedEDBins); // Should be 0.01 be default
2957 
2958  TLorentzVector *jet_vec= new TLorentzVector;
2959  TLorentzVector *track_vec = new TLorentzVector;
2960 
2961  for (i=0;i<fTPCFullJet->GetTotalSignalJets();i++)
2962  {
2964  if (IsInTPC(fChargedEDJetR,myJet->Phi(),myJet->Eta(),kTRUE)==kTRUE)
2965  {
2966  for (j=0;j<chargedEDBins;j++)
2967  {
2968  ED_pT[j]=0;
2969  }
2970  myJet->GetMomentum(*jet_vec);
2971 
2972  // Sum all tracks in concentric rings around jet vertex
2973  for (j=0;j<fnTracks;j++)
2974  {
2975  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2976  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2977  delta_R=jet_vec->DeltaR(*track_vec);
2978  if (delta_R<fChargedEDJetR)
2979  {
2980  ED_pT[TMath::FloorNint((delta_R/fChargedEDJetR)*chargedEDBins)]+=vtrack->Pt();
2981  }
2982  }
2983  for (j=0;j<chargedEDBins;j++)
2984  {
2985  ED_pT[j] /= TMath::Pi()*dR*dR*(2*j+1);
2986  fpChargedJetEDProfile->Fill(myJet->Pt(),fEventCentrality,j*dR,ED_pT[j]);
2987  fpChargedJetEDProfileScaled->Fill(myJet->Pt(),fEventCentrality,j*dR,fScaleFactor*ED_pT[j]);
2988  }
2989  }
2990  }
2991  delete track_vec;
2992  delete jet_vec;
2993 }
2994 
2996 {
2997  if (delOption ==0)
2998  {
2999  delete fRecoUtil;
3000 
3001  fRecoUtil = NULL;
3002  }
3003  else if (delOption==1)
3004  {
3005  delete fRecoUtil;
3006  delete fmyTracks;
3007  delete fmyClusters;
3008 
3009  fRecoUtil = NULL;
3010  fmyTracks = NULL;
3011  fmyClusters = NULL;
3012  }
3013  else if (delOption==2)
3014  {
3015  delete fRecoUtil;
3016  delete fmyTracks;
3017  delete fmyClusters;
3018  delete fTPCJet;
3019  delete fTPCFullJet;
3020  delete fTPCOnlyJet;
3021  delete fTPCJetUnbiased;
3022  delete fTPCkTFullJet;
3023 
3024  fRecoUtil = NULL;
3025  fmyTracks = NULL;
3026  fmyClusters = NULL;
3027  fTPCJet = NULL;
3028  fTPCFullJet = NULL;
3029  fTPCOnlyJet = NULL;
3030  fTPCJetUnbiased = NULL;
3031  fTPCkTFullJet = NULL;
3032  }
3033  if (delOption==3)
3034  {
3035  delete fRecoUtil;
3036  delete fmyTracks;
3037  delete fmyClusters;
3038  delete fTPCJet;
3039  delete fTPCFullJet;
3040  delete fTPCOnlyJet;
3041  delete fTPCJetUnbiased;
3042  delete fTPCkTFullJet;
3043  delete fEMCalJet;
3044  delete fEMCalFullJet;
3045  delete fEMCalPartJet;
3046  delete fEMCalPartJetUnbiased;
3047  delete fEMCalkTFullJet;
3048 
3049  fRecoUtil = NULL;
3050  fmyTracks = NULL;
3051  fmyClusters = NULL;
3052  fTPCJet = NULL;
3053  fTPCFullJet = NULL;
3054  fTPCOnlyJet = NULL;
3055  fTPCJetUnbiased = NULL;
3056  fEMCalJet = NULL;
3057  fEMCalFullJet = NULL;
3058  fEMCalPartJet = NULL;
3059  fEMCalPartJetUnbiased = NULL;
3060  fEMCalkTFullJet = NULL;
3061  }
3062 }
3063 
3067 
3069 {
3070  // Determine if event contains a di-jet within the detector. Uses charged jets.
3071  // Requires the delta phi of the jets to be 180 +/- 15 degrees.
3072  // Requires both jets to be outside of the EMCal
3073  // Requires both jets to be signal jets
3074 
3075  const Double_t dijet_delta_phi=(180/360.)*2*TMath::Pi();
3076  const Double_t dijet_phi_acceptance=0.5*(30/360.)*2*TMath::Pi(); //Input the total acceptance within the paraenthesis to be +/- dijet_phi_acceptance
3077  Double_t dummy_phi=0.0;
3078 
3080  {
3083  dummy_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myJet->Phi()));
3084  if (dummy_phi>(dijet_delta_phi-dijet_phi_acceptance))
3085  {
3086  fnDiJetEvents++;
3087  return kTRUE;
3088  }
3089  }
3090  return kFALSE;
3091 }
3092 
3094 {
3095  if (phi>phi_min && phi<phi_max)
3096  {
3097  if (eta>eta_min && eta<eta_max)
3098  {
3099  return kTRUE;
3100  }
3101  }
3102  return kFALSE;
3103 }
3104 
3106 {
3108 }
3109 
3111 {
3113 }
3114 
3116 {
3118 }
3119 
3121 {
3124 
3125  if (in_EMCal==kFALSE && in_TPC==kTRUE)
3126  {
3127  return kTRUE;
3128  }
3129  return kFALSE;
3130 }
3131 
3133 {
3134  if (Complete==kTRUE)
3135  {
3137  }
3139 }
3140 
3142 {
3143  Int_t i,j;
3144  Int_t jetTrack1=0;
3145  Int_t jetTrack2=0;
3146  Int_t jetCluster1=0;
3147  Int_t jetCluster2=0;
3148 
3149  for (i=0;i<jet1->GetNumberOfTracks();i++)
3150  {
3151  jetTrack1=jet1->TrackAt(i);
3152  for (j=0;j<jet2->GetNumberOfTracks();j++)
3153  {
3154  jetTrack2=jet2->TrackAt(j);
3155  if (jetTrack1 == jetTrack2)
3156  {
3157  return kTRUE;
3158  }
3159  }
3160  }
3161  if (EMCalOn == kTRUE)
3162  {
3163  for (i=0;i<jet1->GetNumberOfClusters();i++)
3164  {
3165  jetCluster1=jet1->ClusterAt(i);
3166  for (j=0;j<jet2->GetNumberOfClusters();j++)
3167  {
3168  jetCluster2=jet2->ClusterAt(j);
3169  if (jetCluster1 == jetCluster2)
3170  {
3171  return kTRUE;
3172  }
3173  }
3174  }
3175  }
3176  return kFALSE;
3177 }
3178 
3180 {
3181  Double_t z;
3182  if (eta<(fTPCEtaMin+r))
3183  {
3184  z=eta-fTPCEtaMin;
3185  }
3186  else if(eta>(fTPCEtaMax-r))
3187  {
3188  z=fTPCEtaMax-eta;
3189  }
3190  else
3191  {
3192  z=r;
3193  }
3194  return r*r*TMath::Pi()-AreaEdge(r,z);
3195 }
3196 
3198 {
3199  Double_t x,y;
3200 
3201  if (phi<(fEMCalPhiMin-r) || phi>(fEMCalPhiMax+r))
3202  {
3203  x=-r;
3204  }
3205  else if (phi<(fEMCalPhiMin+r))
3206  {
3207  x=phi-fEMCalPhiMin;
3208  }
3209  else if (phi>(fEMCalPhiMin+r) && phi<(fEMCalPhiMax-r))
3210  {
3211  x=r;
3212  }
3213  else
3214  {
3215  x=fEMCalPhiMax-phi;
3216  }
3217 
3218  if (eta<(fEMCalEtaMin-r) || eta>(fEMCalEtaMax+r))
3219  {
3220  y=-r;
3221  }
3222  else if (eta<(fEMCalEtaMin+r))
3223  {
3224  y=eta-fEMCalEtaMin;
3225  }
3226  else if (eta>(fEMCalEtaMin+r) && eta<(fEMCalEtaMax-r))
3227  {
3228  y=r;
3229  }
3230  else
3231  {
3232  y=fEMCalEtaMax-eta;
3233  }
3234 
3235  if (x>=0 && y>=0)
3236  {
3237  if (TMath::Sqrt(x*x+y*y)>=r)
3238  {
3239  return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3240  }
3241  return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y)+AreaOverlap(r,x,y);
3242  }
3243  else if ((x>=r && y<0) || (y>=r && x<0))
3244  {
3245  return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3246  }
3247  else if (x>0 && x<r && y<0)
3248  {
3249  Double_t a=TMath::Sqrt(r*r-x*x);
3250  Double_t b=TMath::Sqrt(r*r-y*y);
3251  if ((x-b)>0)
3252  {
3253  return r*r*TMath::ASin(b/r)+y*b;
3254  }
3255  else
3256  {
3257  return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
3258  }
3259  }
3260  else if (y>0 && y<r && x<0)
3261  {
3262  Double_t a=TMath::Sqrt(r*r-x*x);
3263  Double_t b=TMath::Sqrt(r*r-y*y);
3264  if ((y-a)>0)
3265  {
3266  return r*r*TMath::ASin(a/r)+x*a;
3267  }
3268  else
3269  {
3270  return 0.5*y*b+0.5*r*r*TMath::ASin(y/r)+0.5*x*a+x*y+0.5*r*r*TMath::ASin(a/r);
3271  }
3272  }
3273  else
3274  {
3275  Double_t a=TMath::Sqrt(r*r-x*x);
3276  Double_t b=TMath::Sqrt(r*r-y*y);
3277  if ((x+b)<0)
3278  {
3279  return 0;
3280  }
3281  else
3282  {
3283  return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
3284  }
3285  }
3286 }
3287 
3289 {
3290  Double_t a=TMath::Sqrt(r*r-z*z);
3291  return r*r*TMath::ASin(a/r)-a*z;
3292 }
3293 
3295 {
3296  Double_t a=TMath::Sqrt(r*r-x*x);
3297  Double_t b=TMath::Sqrt(r*r-y*y);
3298  return x*y-0.5*(x*a+y*b)+0.5*r*r*(TMath::ASin(b/r)-TMath::ASin(x/r));
3299 }
3300 
3302 {
3303  Double_t area_left=0;
3304  Double_t area_right=0;
3305  Double_t eta_a=0;
3306  Double_t eta_b=0;
3307  Double_t eta_up=0;
3308  Double_t eta_down=0;
3309 
3310  Double_t u=eta-fEMCalEtaMin;
3311  Double_t v=fEMCalEtaMax-eta;
3312 
3313  Double_t phi1=phi+u*TMath::Tan(psi0);
3314  Double_t phi2=phi-u*TMath::Tan(psi0);
3315  Double_t phi3=phi+v*TMath::Tan(psi0);
3316  Double_t phi4=phi-v*TMath::Tan(psi0);
3317 
3318  //Calculate the Left side area
3319  if (phi1>=fEMCalPhiMax)
3320  {
3321  eta_a=eta-u*((fEMCalPhiMax-phi)/(phi1-phi));
3322  }
3323  if (phi2<=fEMCalPhiMin)
3324  {
3325  eta_b=eta-u*((phi-fEMCalPhiMin)/(phi-phi2));
3326  }
3327 
3328  if ((phi1>=fEMCalPhiMax) && (phi2<=fEMCalPhiMin))
3329  {
3330  eta_up=TMath::Max(eta_a,eta_b);
3331  eta_down=TMath::Min(eta_a,eta_b);
3332 
3333  area_left=(eta_down-fEMCalEtaMin)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta-eta_up)*TMath::Tan(psi0))*(eta_up-eta_down) + (eta-eta_up+r)*TMath::Tan(psi0)*(eta-eta_up-r);
3334  }
3335  else if (phi1>=fEMCalPhiMax)
3336  {
3337  area_left=0.5*(fEMCalPhiMax-phi2+2*(eta-eta_a)*TMath::Tan(psi0))*(eta_a-fEMCalEtaMin) + (eta-eta_a+r)*TMath::Tan(psi0)*(eta-eta_a-r);
3338  }
3339  else if (phi2<=fEMCalPhiMin)
3340  {
3341  area_left=0.5*(phi1-fEMCalPhiMin+2*(eta-eta_b)*TMath::Tan(psi0))*(eta_b-fEMCalEtaMin) + (eta-eta_b+r)*TMath::Tan(psi0)*(eta-eta_b-r);
3342  }
3343  else
3344  {
3345  area_left=0.5*(phi1-phi2+2*r*TMath::Tan(psi0))*(u-r);
3346  }
3347 
3348  // Calculate the Right side area
3349  if (phi3>=fEMCalPhiMax)
3350  {
3351  eta_a=eta+v*((fEMCalPhiMax-phi)/(phi3-phi));
3352  }
3353  if (phi4<=fEMCalPhiMin)
3354  {
3355  eta_b=eta+v*((phi-fEMCalPhiMin)/(phi-phi4));
3356  }
3357 
3358  if ((phi3>=fEMCalPhiMax) && (phi4<=fEMCalPhiMin))
3359  {
3360  eta_up=TMath::Max(eta_a,eta_b);
3361  eta_down=TMath::Min(eta_a,eta_b);
3362 
3363  area_right=(fEMCalEtaMax-eta_up)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta_down-eta)*TMath::Tan(psi0))*(eta_down-eta_up) + (eta_down-eta+r)*TMath::Tan(psi0)*(eta_up-eta-r);
3364  }
3365  else if (phi3>=fEMCalPhiMax)
3366  {
3367  area_right=0.5*(fEMCalPhiMax-phi4+2*(eta_a-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_a) + (eta_a-eta+r)*TMath::Tan(psi0)*(eta_a-eta-r);
3368  }
3369  else if (phi4<=fEMCalPhiMin)
3370  {
3371  area_right=0.5*(phi3-fEMCalPhiMin+2*(eta_b-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_b) + (eta_b-eta+r)*TMath::Tan(psi0)*(eta_b-eta-r);
3372  }
3373  else
3374  {
3375  area_right=0.5*(phi3-phi4+2*r*TMath::Tan(psi0))*(v-r);
3376  }
3377  return area_left+area_right;
3378 }
3379 
3381 {
3382  // This function is used to calculate the median Rho kT value. The procedure is:
3383  // - Order the kT cluster array from highest rho value to lowest
3384  // - Exclude highest rho kT cluster
3385  // - Return the median rho value of the remaining subset
3386 
3387  // Sort Array
3388  const Double_t rho_min=-9.9999E+99;
3389  Int_t j,k;
3390  Double_t w[nEntries]; // Used for sorting
3391  Double_t smax=rho_min;
3392  Int_t sindex=-1;
3393 
3394  Double_t pTmax=0.0;
3395  Int_t pTmaxID=-1;
3396 
3397  for (j=0;j<nEntries;j++)
3398  {
3399  w[j]=0.0;
3400  }
3401 
3402  for (j=0;j<nEntries;j++)
3403  {
3404  if (pTkTEntries[j]>pTmax)
3405  {
3406  pTmax=pTkTEntries[j];
3407  pTmaxID=j;
3408  }
3409  }
3410 
3411  for (j=0;j<nEntries;j++)
3412  {
3413  for (k=0;k<nEntries;k++)
3414  {
3415  if (RhokTEntries[k]>smax)
3416  {
3417  smax=RhokTEntries[k];
3418  sindex=k;
3419  }
3420  }
3421  w[j]=smax;
3422  RhokTEntries[sindex]=rho_min;
3423  smax=rho_min;
3424  sindex=-1;
3425  }
3426  return w[nEntries/2];
3427 }
3428 
3429 
3430 // AlipAJetData Class Member Defs
3431 // Constructors
3433 
3434  fName(0),
3435  fIsJetsFull(0),
3436  fnTotal(0),
3437  fnJets(0),
3438  fnJetsSC(0),
3439  fJetR(0),
3440  fSignalPt(0),
3441  fAreaCutFrac(0.6),
3442  fNEF(1.0),
3443  fSignalTrackBias(0),
3444  fPtMaxIndex(0),
3445  fPtMax(0),
3446  fPtSubLeadingIndex(0),
3447  fPtSubLeading(0),
3448  fJetsIndex(0),
3449  fJetsSCIndex(0),
3450  fIsJetInArray(0),
3451  fJetMaxChargedPt(0)
3452 {
3453  fnTotal=0;
3454  // Dummy constructor ALWAYS needed for I/O.
3455 }
3456 
3458 
3459  fName(0),
3460  fIsJetsFull(0),
3461  fnTotal(0),
3462  fnJets(0),
3463  fnJetsSC(0),
3464  fJetR(0),
3465  fSignalPt(0),
3466  fAreaCutFrac(0.6),
3467  fNEF(1.0),
3468  fSignalTrackBias(0),
3469  fPtMaxIndex(0),
3470  fPtMax(0),
3471  fPtSubLeadingIndex(0),
3472  fPtSubLeading(0),
3473  fJetsIndex(0),
3474  fJetsSCIndex(0),
3475  fIsJetInArray(0),
3476  fJetMaxChargedPt(0)
3477 {
3478  SetName(name);
3479  SetIsJetsFull(isFull);
3480  SetTotalEntries(nEntries);
3481  SetLeading(0,-9.99E+099);
3482  SetSubLeading(0,-9.99E+099);
3483  SetSignalCut(0);
3484  SetAreaCutFraction(0.6);
3485  SetJetR(fJetR);
3487 }
3488 
3489 // Destructor
3491 {
3492  SetName("");
3493  SetIsJetsFull(kFALSE);
3494  SetTotalJets(0);
3495  SetTotalSignalJets(0);
3496  SetLeading(0,0);
3497  SetSubLeading(0,0);
3498  SetSignalCut(0);
3499  SetAreaCutFraction(0);
3500  SetJetR(0);
3501  SetNEF(0);
3502  SetSignalTrackPtBias(kFALSE);
3503 
3504  delete [] fJetsIndex;
3505  delete [] fJetsSCIndex;
3506  delete [] fIsJetInArray;
3507  delete [] fJetMaxChargedPt;
3508 
3509  fnTotal = 0;
3510 }
3511 
3512 // User Defined Sub-Routines
3514 {
3515  Int_t i=0;
3516  Int_t k=0;
3517  Int_t l=0;
3518  Double_t AreaThreshold = fAreaCutFrac*TMath::Pi()*TMath::Power(fJetR,2);
3519 
3520  // Initialize Jet Data
3521  for (i=0;i<nEntries;i++)
3522  {
3523  AliEmcalJet *myJet =(AliEmcalJet*) jetList->At(i);
3524 
3525  if (fIsJetInArray[i]==kTRUE && myJet->Area()>AreaThreshold)
3526  {
3527  SetJetIndex(i,k);
3528  if (myJet->Pt()>fPtMax)
3529  {
3531  SetLeading(i,myJet->Pt());
3532  }
3533  else if (myJet->Pt()>fPtSubLeading)
3534  {
3535  SetSubLeading(i,myJet->Pt());
3536  }
3537  // require leading charged constituent to have a pT greater then the signal threshold & Jet NEF to be less then the Signal Jet NEF cut
3538  fJetMaxChargedPt[i] = myJet->MaxTrackPt();
3539  if (fSignalTrackBias==kTRUE)
3540  {
3541  if (fJetMaxChargedPt[i]>=fSignalPt && myJet->NEF()<=fNEF)
3542  {
3543  SetSignalJetIndex(i,l);
3544  l++;
3545  }
3546  }
3547  else
3548  {
3549  if (myJet->Pt()>=fSignalPt && myJet->NEF()<=fNEF)
3550  {
3551  SetSignalJetIndex(i,l);
3552  l++;
3553  }
3554  }
3555  k++;
3556  }
3557  }
3558  SetTotalJets(k);
3559  SetTotalSignalJets(l);
3560 }
3561 
3562 // Setters
3564 {
3565  fName = name;
3566 }
3567 
3569 {
3570  fIsJetsFull = isFull;
3571 }
3572 
3574 {
3575  fnTotal = nEntries;
3576  fJetsIndex = new Int_t[fnTotal];
3577  fJetsSCIndex = new Int_t[fnTotal];
3578  fIsJetInArray = new Bool_t[fnTotal];
3580 }
3581 
3583 {
3584  fnJets = nJets;
3585 }
3586 
3588 {
3589  fnJetsSC = nSignalJets;
3590 }
3591 
3593 {
3594  fSignalPt = Pt;
3595 }
3596 
3598 {
3599  fPtMaxIndex = index;
3600  fPtMax = Pt;
3601 }
3602 
3604 {
3605  fPtSubLeadingIndex = index;
3606  fPtSubLeading = Pt;
3607 }
3608 
3610 {
3611  fJetsIndex[At] = index;
3612 }
3613 
3615 {
3616  fJetsSCIndex[At] = index;
3617 }
3618 
3620 {
3621  fIsJetInArray[At] = isInArray;
3622 }
3623 
3625 {
3626  fAreaCutFrac = areaFraction;
3627 }
3628 
3630 {
3631  fJetR = jetR;
3632 }
3633 
3635 {
3636  fNEF = nef;
3637 }
3638 
3640 {
3641  fSignalTrackBias = chargedBias;
3642 }
3643 
3644 // Getters
3646 {
3647  return fnTotal;
3648 }
3649 
3651 {
3652  return fnJets;
3653 }
3654 
3656 {
3657  return fnJetsSC;
3658 }
3659 
3661 {
3662  return fSignalPt;
3663 }
3664 
3666 {
3667  return fPtMaxIndex;
3668 }
3669 
3671 {
3672  return fPtMax;
3673 }
3674 
3676 {
3677  return fPtSubLeadingIndex;
3678 }
3679 
3681 {
3682  return fPtSubLeading;
3683 }
3684 
3686 {
3687  return fJetsIndex[At];
3688 }
3689 
3691 {
3692  return fJetsSCIndex[At];
3693 }
3694 
3696 {
3697  return fIsJetInArray[At];
3698 }
3699 
3701 {
3702  return fJetMaxChargedPt[At];
3703 }
3704 
3706 {
3707  return fNEF;
3708 }
3709 
3710 // AlipAJetHistos Class Member Defs
3711 // Constructors
3713 
3714  fOutput(0),
3715 
3716  fh020Rho(0),
3717  fh80100Rho(0),
3718  fhRho(0),
3719  fhRhoCen(0),
3720  fh020BSPt(0),
3721  fh80100BSPt(0),
3722  fhBSPt(0),
3723  fhBSPtCen(0),
3724  fh020BSPtSignal(0),
3725  fh80100BSPtSignal(0),
3726  fhBSPtSignal(0),
3727  fhBSPtCenSignal(0),
3728  fh020DeltaPt(0),
3729  fh80100DeltaPt(0),
3730  fhDeltaPt(0),
3731  fhDeltaPtCen(0),
3732  fh020DeltaPtSignal(0),
3733  fh80100DeltaPtSignal(0),
3734  fhDeltaPtSignal(0),
3735  fhDeltaPtCenSignal(0),
3736  fh020DeltaPtNColl(0),
3737  fh80100DeltaPtNColl(0),
3738  fhDeltaPtNColl(0),
3739  fhDeltaPtCenNColl(0),
3740  fh020BckgFlucPt(0),
3741  fh80100BckgFlucPt(0),
3742  fhBckgFlucPt(0),
3743  fhBckgFlucPtCen(0),
3744 
3745  fpRho(0),
3746  fpLJetRho(0),
3747 
3748  fhJetPtEtaPhi(0),
3749  fhJetPtArea(0),
3750  fhJetConstituentPt(0),
3751  fhJetTracksPt(0),
3752  fhJetClustersPt(0),
3753  fhJetConstituentCounts(0),
3754  fhJetTracksCounts(0),
3755  fhJetClustersCounts(0),
3756  fhJetPtZConstituent(0),
3757  fhJetPtZTrack(0),
3758  fhJetPtZCluster(0),
3759  fhJetPtZLeadingConstituent(0),
3760  fhJetPtZLeadingTrack(0),
3761  fhJetPtZLeadingCluster(0),
3762 
3763  fhEventCentralityVsZNA(0),
3764  fhEventCentralityVsZNAPt(0),
3765 
3766  fNEFOutput(0),
3767  fhJetPtNEF(0),
3768  fhJetNEFInfo(0),
3769  fhJetNEFSignalInfo(0),
3770  fhClusterNEFInfo(0),
3771  fhClusterNEFSignalInfo(0),
3772  fhClusterShapeAll(0),
3773  fhClusterPtCellAll(0),
3774 
3775  fName(0),
3776  fCentralityTag(0),
3777  fCentralityBins(0),
3778  fCentralityLow(0),
3779  fCentralityUp(0),
3780  fPtBins(0),
3781  fPtLow(0),
3782  fPtUp(0),
3783  fRhoPtBins(0),
3784  fRhoPtLow(0),
3785  fRhoPtUp(0),
3786  fDeltaPtBins(0),
3787  fDeltaPtLow(0),
3788  fDeltaPtUp(0),
3789  fBckgFlucPtBins(0),
3790  fBckgFlucPtLow(0),
3791  fBckgFlucPtUp(0),
3792  fLJetPtBins(0),
3793  fLJetPtLow(0),
3794  fLJetPtUp(0),
3795  fRhoValue(0),
3796  fLChargedTrackPtBins(0),
3797  fLChargedTrackPtLow(0),
3798  fLChargedTrackPtUp(0),
3799  fDoNEFQAPlots(0),
3800  fDoNEFSignalOnly(1),
3801  fSignalTrackBias(0),
3802  fDoTHnSparse(0),
3803  fDo3DHistos(0),
3804  fNEFBins(0),
3805  fNEFLow(0),
3806  fNEFUp(0),
3807  fnDimJet(0),
3808  fnDimCluster(0),
3809  fEMCalPhiMin(1.39626),
3810  fEMCalPhiMax(3.26377),
3811  fEMCalEtaMin(-0.7),
3812  fEMCalEtaMax(0.7)
3813 
3814 {
3815  // Dummy constructor ALWAYS needed for I/O.
3816 }
3817 
3819 
3820  fOutput(0),
3821 
3822  fh020Rho(0),
3823  fh80100Rho(0),
3824  fhRho(0),
3825  fhRhoCen(0),
3826  fh020BSPt(0),
3827  fh80100BSPt(0),
3828  fhBSPt(0),
3829  fhBSPtCen(0),
3830  fh020BSPtSignal(0),
3831  fh80100BSPtSignal(0),
3832  fhBSPtSignal(0),
3833  fhBSPtCenSignal(0),
3834  fh020DeltaPt(0),
3835  fh80100DeltaPt(0),
3836  fhDeltaPt(0),
3837  fhDeltaPtCen(0),
3838  fh020DeltaPtSignal(0),
3840  fhDeltaPtSignal(0),
3841  fhDeltaPtCenSignal(0),
3842  fh020DeltaPtNColl(0),
3844  fhDeltaPtNColl(0),
3845  fhDeltaPtCenNColl(0),
3846  fh020BckgFlucPt(0),
3847  fh80100BckgFlucPt(0),
3848  fhBckgFlucPt(0),
3849  fhBckgFlucPtCen(0),
3850 
3851  fpRho(0),
3852  fpLJetRho(0),
3853 
3854  fhJetPtEtaPhi(0),
3855  fhJetPtArea(0),
3856  fhJetConstituentPt(0),
3857  fhJetTracksPt(0),
3858  fhJetClustersPt(0),
3860  fhJetTracksCounts(0),
3863  fhJetPtZTrack(0),
3864  fhJetPtZCluster(0),
3868 
3871 
3872  fNEFOutput(0),
3873  fhJetPtNEF(0),
3874  fhJetNEFInfo(0),
3875  fhJetNEFSignalInfo(0),
3876  fhClusterNEFInfo(0),
3878  fhClusterShapeAll(0),
3879  fhClusterPtCellAll(0),
3880 
3881  fName(0),
3882  fCentralityTag(0),
3883  fCentralityBins(0),
3884  fCentralityLow(0),
3885  fCentralityUp(0),
3886  fPtBins(0),
3887  fPtLow(0),
3888  fPtUp(0),
3889  fRhoPtBins(0),
3890  fRhoPtLow(0),
3891  fRhoPtUp(0),
3892  fDeltaPtBins(0),
3893  fDeltaPtLow(0),
3894  fDeltaPtUp(0),
3895  fBckgFlucPtBins(0),
3896  fBckgFlucPtLow(0),
3897  fBckgFlucPtUp(0),
3898  fLJetPtBins(0),
3899  fLJetPtLow(0),
3900  fLJetPtUp(0),
3901  fRhoValue(0),
3904  fLChargedTrackPtUp(0),
3905  fDoNEFQAPlots(0),
3906  fDoNEFSignalOnly(1),
3907  fSignalTrackBias(0),
3908  fDoTHnSparse(0),
3909  fDo3DHistos(0),
3910  fNEFBins(0),
3911  fNEFLow(0),
3912  fNEFUp(0),
3913  fnDimJet(0),
3914  fnDimCluster(0),
3915  fEMCalPhiMin(1.39626),
3916  fEMCalPhiMax(3.26377),
3917  fEMCalEtaMin(-0.7),
3918  fEMCalEtaMax(0.7)
3919 
3920 {
3921  SetName(name);
3922  SetCentralityTag("ZNA");
3923  SetCentralityRange(100,0,100);
3924  SetPtRange(250,-50,200);
3925  SetRhoPtRange(500,0,50);
3926  SetDeltaPtRange(200,-100,100);
3928  SetLeadingJetPtRange(200,0,200);
3929  SetLeadingChargedTrackPtRange(100,0,100);
3930  SetNEFRange(100,0,1);
3931  DoNEFQAPlots(kFALSE);
3932 
3933  Init();
3934 }
3935 
3937 
3938  fOutput(0),
3939 
3940  fh020Rho(0),
3941  fh80100Rho(0),
3942  fhRho(0),
3943  fhRhoCen(0),
3944  fh020BSPt(0),
3945  fh80100BSPt(0),
3946  fhBSPt(0),
3947  fhBSPtCen(0),
3948  fh020BSPtSignal(0),
3949  fh80100BSPtSignal(0),
3950  fhBSPtSignal(0),
3951  fhBSPtCenSignal(0),
3952  fh020DeltaPt(0),
3953  fh80100DeltaPt(0),
3954  fhDeltaPt(0),
3955  fhDeltaPtCen(0),
3956  fh020DeltaPtSignal(0),
3958  fhDeltaPtSignal(0),
3959  fhDeltaPtCenSignal(0),
3960  fh020DeltaPtNColl(0),
3962  fhDeltaPtNColl(0),
3963  fhDeltaPtCenNColl(0),
3964  fh020BckgFlucPt(0),
3965  fh80100BckgFlucPt(0),
3966  fhBckgFlucPt(0),
3967  fhBckgFlucPtCen(0),
3968 
3969  fpRho(0),
3970  fpLJetRho(0),
3971 
3972  fhJetPtEtaPhi(0),
3973  fhJetPtArea(0),
3974  fhJetConstituentPt(0),
3975  fhJetTracksPt(0),
3976  fhJetClustersPt(0),
3978  fhJetTracksCounts(0),
3981  fhJetPtZTrack(0),
3982  fhJetPtZCluster(0),
3986 
3989 
3990  fNEFOutput(0),
3991  fhJetPtNEF(0),
3992  fhJetNEFInfo(0),
3993  fhJetNEFSignalInfo(0),
3994  fhClusterNEFInfo(0),
3996  fhClusterShapeAll(0),
3997  fhClusterPtCellAll(0),
3998 
3999  fName(0),
4000  fCentralityTag(0),
4001  fCentralityBins(0),
4002  fCentralityLow(0),
4003  fCentralityUp(0),
4004  fPtBins(0),
4005  fPtLow(0),
4006  fPtUp(0),
4007  fRhoPtBins(0),
4008  fRhoPtLow(0),
4009  fRhoPtUp(0),
4010  fDeltaPtBins(0),
4011  fDeltaPtLow(0),
4012  fDeltaPtUp(0),
4013  fBckgFlucPtBins(0),
4014  fBckgFlucPtLow(0),
4015  fBckgFlucPtUp(0),
4016  fLJetPtBins(0),
4017  fLJetPtLow(0),
4018  fLJetPtUp(0),
4019  fRhoValue(0),
4022  fLChargedTrackPtUp(0),
4023  fDoNEFQAPlots(0),
4024  fDoNEFSignalOnly(1),
4025  fSignalTrackBias(0),
4026  fDoTHnSparse(0),
4027  fDo3DHistos(0),
4028  fNEFBins(0),
4029  fNEFLow(0),
4030  fNEFUp(0),
4031  fnDimJet(0),
4032  fnDimCluster(0),
4033  fEMCalPhiMin(1.39626),
4034  fEMCalPhiMax(3.26377),
4035  fEMCalEtaMin(-0.7),
4036  fEMCalEtaMax(0.7)
4037 
4038 {
4039  SetName(name);
4040  SetCentralityTag(centag.Data());
4041  SetCentralityRange(100,0,100);
4042  SetPtRange(250,-50,200);
4043  SetRhoPtRange(500,0,50);
4044  SetDeltaPtRange(200,-100,100);
4046  SetLeadingJetPtRange(200,0,200);
4047  SetLeadingChargedTrackPtRange(100,0,100);
4048  SetNEFRange(100,0,1);
4049  DoNEFQAPlots(doNEF);
4050 
4051  Init();
4052 }
4053 
4054 AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, TString centag, Bool_t doNEF, Bool_t doNEFSignalOnly, Bool_t doTHnSparse, Bool_t do3DPlotting) :
4055 
4056  fOutput(0),
4057 
4058  fh020Rho(0),
4059  fh80100Rho(0),
4060  fhRho(0),
4061  fhRhoCen(0),
4062  fh020BSPt(0),
4063  fh80100BSPt(0),
4064  fhBSPt(0),
4065  fhBSPtCen(0),
4066  fh020BSPtSignal(0),
4067  fh80100BSPtSignal(0),
4068  fhBSPtSignal(0),
4069  fhBSPtCenSignal(0),
4070  fh020DeltaPt(0),
4071  fh80100DeltaPt(0),
4072  fhDeltaPt(0),
4073  fhDeltaPtCen(0),
4074  fh020DeltaPtSignal(0),
4076  fhDeltaPtSignal(0),
4077  fhDeltaPtCenSignal(0),
4078  fh020DeltaPtNColl(0),
4080  fhDeltaPtNColl(0),
4081  fhDeltaPtCenNColl(0),
4082  fh020BckgFlucPt(0),
4083  fh80100BckgFlucPt(0),
4084  fhBckgFlucPt(0),
4085  fhBckgFlucPtCen(0),
4086 
4087  fpRho(0),
4088  fpLJetRho(0),
4089 
4090  fhJetPtEtaPhi(0),
4091  fhJetPtArea(0),
4092  fhJetConstituentPt(0),
4093  fhJetTracksPt(0),
4094  fhJetClustersPt(0),
4096  fhJetTracksCounts(0),
4099  fhJetPtZTrack(0),
4100  fhJetPtZCluster(0),
4104 
4107 
4108  fNEFOutput(0),
4109  fhJetPtNEF(0),
4110  fhJetNEFInfo(0),
4111  fhJetNEFSignalInfo(0),
4112  fhClusterNEFInfo(0),
4114  fhClusterShapeAll(0),
4115  fhClusterPtCellAll(0),
4116 
4117  fName(0),
4118  fCentralityTag(0),
4119  fCentralityBins(0),
4120  fCentralityLow(0),
4121  fCentralityUp(0),
4122  fPtBins(0),
4123  fPtLow(0),
4124  fPtUp(0),
4125  fRhoPtBins(0),
4126  fRhoPtLow(0),
4127  fRhoPtUp(0),
4128  fDeltaPtBins(0),
4129  fDeltaPtLow(0),
4130  fDeltaPtUp(0),
4131  fBckgFlucPtBins(0),
4132  fBckgFlucPtLow(0),
4133  fBckgFlucPtUp(0),
4134  fLJetPtBins(0),
4135  fLJetPtLow(0),
4136  fLJetPtUp(0),
4137  fRhoValue(0),
4140  fLChargedTrackPtUp(0),
4141  fDoNEFQAPlots(0),
4142  fDoNEFSignalOnly(1),
4143  fSignalTrackBias(0),
4144  fDoTHnSparse(0),
4145  fDo3DHistos(0),
4146  fNEFBins(0),
4147  fNEFLow(0),
4148  fNEFUp(0),
4149  fnDimJet(0),
4150  fnDimCluster(0),
4151  fEMCalPhiMin(1.39626),
4152  fEMCalPhiMax(3.26377),
4153  fEMCalEtaMin(-0.7),
4154  fEMCalEtaMax(0.7)
4155 
4156 {
4157  SetName(name);
4158  SetCentralityTag(centag.Data());
4159  SetCentralityRange(100,0,100);
4160  SetPtRange(250,-50,200);
4161  SetRhoPtRange(500,0,50);
4162  SetDeltaPtRange(200,-100,100);
4164  SetLeadingJetPtRange(200,0,200);
4165  SetLeadingChargedTrackPtRange(100,0,100);
4166  SetNEFRange(100,0,1);
4167  DoNEFQAPlots(doNEF);
4168  DoNEFSignalOnly(doNEFSignalOnly);
4169  DoTHnSparse(doTHnSparse);
4170  Do3DPlotting(do3DPlotting);
4171 
4172  Init();
4173 }
4174 
4175 // Destructor
4177 {
4178  if (fOutput)
4179  {
4180  delete fOutput;
4181  }
4182 }
4183 
4185 {
4186  Int_t i;
4187 
4188  // Initialize Private Variables
4189  Int_t TCBins = 100;
4190  fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
4191  fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
4192  fEMCalEtaMin=-0.7;
4193  fEMCalEtaMax=0.7;
4194 
4195  fOutput = new TList();
4196  fOutput->SetOwner();
4197  fOutput->SetName(fName);
4198 
4199  TString RhoString="";
4200  TString PtString="";
4201  TString DeltaPtString="";
4202  TString BckgFlucPtString="";
4203  TString CentralityString;
4204  TString EventCentralityString;
4205  CentralityString = Form("Centrality (%s) ",fCentralityTag.Data());
4206  EventCentralityString = Form("%s vs ZNA Centrality ",fCentralityTag.Data());
4207 
4208  // Rho Spectral Plots
4209  RhoString = Form("%d-%d Centrality, Rho Spectrum",0,20);
4210  fh020Rho = new TH1F("fh020Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
4211  fh020Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4212  fh020Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
4213  fh020Rho->Sumw2();
4214 
4215  RhoString = Form("%d-%d Centrality, Rho Spectrum",80,100);
4216  fh80100Rho = new TH1F("fh80100Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
4217  fh80100Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4218  fh80100Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
4219  fh80100Rho->Sumw2();
4220 
4221  RhoString = Form("%d-%d Centrality, Rho Spectrum",0,100);
4222  fhRho = new TH1F("fhRho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
4223  fhRho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4224  fhRho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
4225  fhRho->Sumw2();
4226 
4227  RhoString = "Rho Spectrum vs Centrality";
4229  fhRhoCen->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4230  fhRhoCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4231  fhRhoCen->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
4232  fhRhoCen->Sumw2();
4233 
4234  // Background Subtracted Plots
4235  PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,20);
4236  fh020BSPt = new TH1F("fh020BSPt",PtString,fPtBins,fPtLow,fPtUp);
4237  fh020BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4238  fh020BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4239  fh020BSPt->Sumw2();
4240 
4241  PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",80,100);
4242  fh80100BSPt = new TH1F("fh80100BSPt",PtString,fPtBins,fPtLow,fPtUp);
4243  fh80100BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4244  fh80100BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4245  fh80100BSPt->Sumw2();
4246 
4247  PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,100);
4248  fhBSPt = new TH1F("fhBSPt",PtString,fPtBins,fPtLow,fPtUp);
4249  fhBSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4250  fhBSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4251  fhBSPt->Sumw2();
4252 
4253  PtString = "Background Subtracted Jet Spectrum vs Centrality";
4255  fhBSPtCen->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4256  fhBSPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4257  fhBSPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4258  fhBSPtCen->Sumw2();
4259 
4260  PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,20);
4261  fh020BSPtSignal = new TH1F("fh020BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
4262  fh020BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4263  fh020BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4264  fh020BSPtSignal->Sumw2();
4265 
4266  PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",80,100);
4267  fh80100BSPtSignal = new TH1F("fh80100BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
4268  fh80100BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4269  fh80100BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4270  fh80100BSPtSignal->Sumw2();
4271 
4272  PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,100);
4273  fhBSPtSignal = new TH1F("fhBSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
4274  fhBSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4275  fhBSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4276  fhBSPtSignal->Sumw2();
4277 
4278  PtString = "Background Subtracted Signal Jet Spectrum vs Centrality";
4280  fhBSPtCenSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4281  fhBSPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4282  fhBSPtCenSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4283  fhBSPtCenSignal->Sumw2();
4284 
4285  // Delta Pt Plots with RC at least 2R away from Leading Signal
4286  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
4287  fh020DeltaPt = new TH1F("fh020DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4288  fh020DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4289  fh020DeltaPt->GetYaxis()->SetTitle("Probability Density");
4290  fh020DeltaPt->Sumw2();
4291 
4292  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
4293  fh80100DeltaPt = new TH1F("fh80100DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4294  fh80100DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4295  fh80100DeltaPt->GetYaxis()->SetTitle("Probability Density");
4296  fh80100DeltaPt->Sumw2();
4297 
4298  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
4299  fhDeltaPt = new TH1F("fhDeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4300  fhDeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4301  fhDeltaPt->GetYaxis()->SetTitle("Probability Density");
4302  fhDeltaPt->Sumw2();
4303 
4304  DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
4306  fhDeltaPtCen->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4307  fhDeltaPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4308  fhDeltaPtCen->GetZaxis()->SetTitle("Probability Density");
4309  fhDeltaPtCen->Sumw2();
4310 
4311  // Delta Pt Plots with no spatial restrictions on RC
4312  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
4313  fh020DeltaPtSignal = new TH1F("fh020DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4314  fh020DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4315  fh020DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
4316  fh020DeltaPtSignal->Sumw2();
4317 
4318  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
4319  fh80100DeltaPtSignal = new TH1F("fh80100DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4320  fh80100DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4321  fh80100DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
4322  fh80100DeltaPtSignal->Sumw2();
4323 
4324  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
4325  fhDeltaPtSignal = new TH1F("fhDeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4326  fhDeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4327  fhDeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
4328  fhDeltaPtSignal->Sumw2();
4329 
4330  DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
4332  fhDeltaPtCenSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4333  fhDeltaPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4334  fhDeltaPtCenSignal->GetZaxis()->SetTitle("Probability Density");
4335  fhDeltaPtCenSignal->Sumw2();
4336 
4337  // Delta Pt Plots with NColl restrictions on RC
4338  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
4339  fh020DeltaPtNColl = new TH1F("fh020DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4340  fh020DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4341  fh020DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
4342  fh020DeltaPtNColl->Sumw2();
4343 
4344  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
4345  fh80100DeltaPtNColl = new TH1F("fh80100DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4346  fh80100DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4347  fh80100DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
4348  fh80100DeltaPtNColl->Sumw2();
4349 
4350  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
4351  fhDeltaPtNColl = new TH1F("fhDeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4352  fhDeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4353  fhDeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
4354  fhDeltaPtNColl->Sumw2();
4355 
4356  DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
4358  fhDeltaPtCenNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4359  fhDeltaPtCenNColl->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4360  fhDeltaPtCenNColl->GetZaxis()->SetTitle("Probability Density");
4361  fhDeltaPtCenNColl->Sumw2();
4362 
4363  // Background Fluctuations Pt Plots
4364  BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,20);
4365  fh020BckgFlucPt = new TH1F("fh020BckgFlucPt",PtString,fPtBins,fPtLow,fPtUp);
4366  fh020BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4367  fh020BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4368  fh020BckgFlucPt->Sumw2();
4369 
4370  BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",80,100);
4371  fh80100BckgFlucPt = new TH1F("fh80100BckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
4372  fh80100BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4373  fh80100BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4374  fh80100BckgFlucPt->Sumw2();
4375 
4376  BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,100);
4377  fhBckgFlucPt = new TH1F("fhBckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
4378  fhBckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4379  fhBckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4380  fhBckgFlucPt->Sumw2();
4381 
4382  BckgFlucPtString = "Background Fluctuation p_{T} Spectrum vs Centrality";
4384  fhBckgFlucPtCen->GetXaxis()->SetTitle("#p_{T} (GeV/c)");
4385  fhBckgFlucPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4386  fhBckgFlucPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4387  fhBckgFlucPtCen->Sumw2();
4388 
4389  // Background Density vs Centrality Profile
4390  RhoString = "Background Density vs Centrality";
4391  fpRho = new TProfile("fpRho",RhoString,fCentralityBins,fCentralityLow,fCentralityUp);
4392  fpRho->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
4393  fpRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
4394 
4395  // Background Density vs Leading Jet Profile
4396  fpLJetRho = new TProfile("fpLJetRho","#rho vs Leading Jet p_{T}",fLJetPtBins,fLJetPtLow,fLJetPtUp);
4397  fpLJetRho->GetXaxis()->SetTitle("Leading Jet p_{T}");
4398  fpLJetRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
4399 
4400  // Jet pT vs Area
4401  Int_t JetPtAreaBins=200;
4402  Double_t JetPtAreaLow=0.0;
4403  Double_t JetPtAreaUp=2.0;
4404 
4405  fhJetPtArea = new TH2F("fhJetPtArea","Jet Area Distribution",fPtBins,fPtLow,fPtUp,JetPtAreaBins,JetPtAreaLow,JetPtAreaUp);
4406  fhJetPtArea->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4407  fhJetPtArea->GetYaxis()->SetTitle("A_{jet}");
4408  fhJetPtArea->GetZaxis()->SetTitle("1/N_{Events} dN/dA_{jet}dp_{T}");
4409  fhJetPtArea->Sumw2();
4410 
4411  // Jet pT vs Constituent pT
4412  fhJetConstituentPt = new TH2F("fhJetConstituentPt","Jet constituents p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
4413  fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4414  fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
4415  fhJetConstituentPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,track}");
4416  fhJetConstituentPt->Sumw2();
4417 
4418  fhJetTracksPt = new TH2F("fhJetTracksPt","Jet constituents Tracks p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
4419  fhJetTracksPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4420  fhJetTracksPt->GetYaxis()->SetTitle("Constituent Track p_{T} (GeV/c)");
4421  fhJetTracksPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,track}");
4422  fhJetTracksPt->Sumw2();
4423 
4424  fhJetClustersPt = new TH2F("fhJetClustersPt","Jet constituents Clusters p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
4425  fhJetClustersPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4426  fhJetClustersPt->GetYaxis()->SetTitle("Constituent Cluster p_{T} (GeV/c)");
4427  fhJetClustersPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,cluster}");
4428  fhJetClustersPt->Sumw2();
4429 
4430  // Jet pT vs Constituent Counts
4431  fhJetConstituentCounts = new TH2F("fhJetConstituentCounts","Jet constituents distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
4432  fhJetConstituentCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4433  fhJetConstituentCounts->GetYaxis()->SetTitle("Constituent Count");
4434  fhJetConstituentCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{constituent}");
4435  fhJetConstituentCounts->Sumw2();
4436 
4437  fhJetTracksCounts = new TH2F("fhJetTracksCounts","Jet constituents Tracks distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
4438  fhJetTracksCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4439  fhJetTracksCounts->GetYaxis()->SetTitle("Constituent Track Count");
4440  fhJetTracksCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{track}");
4441  fhJetTracksCounts->Sumw2();
4442 
4443  fhJetClustersCounts = new TH2F("fhJetClustersCounts","Jet constituents Clusters distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
4444  fhJetClustersCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4445  fhJetClustersCounts->GetYaxis()->SetTitle("Constituent Cluster Count");
4446  fhJetClustersCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{cluster}");
4447  fhJetClustersCounts->Sumw2();
4448 
4449  // Jet pT vs z_{constituent/track/cluster}
4450  fhJetPtZConstituent = new TH2F("fhJetPtZConstituent","Jet z_{constituent} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4451  fhJetPtZConstituent->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4452  fhJetPtZConstituent->GetYaxis()->SetTitle("z_{constituent}");
4453  fhJetPtZConstituent->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{constituent}");
4454  fhJetPtZConstituent->Sumw2();
4455 
4456  fhJetPtZTrack = new TH2F("fhJetPtZTrack","Jet z_{track} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4457  fhJetPtZTrack->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4458  fhJetPtZTrack->GetYaxis()->SetTitle("z_{track}");
4459  fhJetPtZTrack->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{track}");
4460  fhJetPtZTrack->Sumw2();
4461 
4462  fhJetPtZCluster = new TH2F("fhJetPtZCluster","Jet z_{cluster} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4463  fhJetPtZCluster->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4464  fhJetPtZCluster->GetYaxis()->SetTitle("z_{cluster}");
4465  fhJetPtZCluster->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{cluster}");
4466  fhJetPtZCluster->Sumw2();
4467 
4468  // Jet pT vs z_Leading{constituent/track/cluster}
4469  fhJetPtZLeadingConstituent = new TH2F("fhJetPtZLeadingConstituent","Jet z_{Leading,constituent} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4470  fhJetPtZLeadingConstituent->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4471  fhJetPtZLeadingConstituent->GetYaxis()->SetTitle("z_{Leading,constituent}");
4472  fhJetPtZLeadingConstituent->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{constituent}");
4473  fhJetPtZLeadingConstituent->Sumw2();
4474 
4475  fhJetPtZLeadingTrack = new TH2F("fhJetPtZLeadingTrack","Jet z_{Leading,track} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4476  fhJetPtZLeadingTrack->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4477  fhJetPtZLeadingTrack->GetYaxis()->SetTitle("z_{Leading,track}");
4478  fhJetPtZLeadingTrack->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{track}");
4479  fhJetPtZLeadingTrack->Sumw2();
4480 
4481  fhJetPtZLeadingCluster = new TH2F("fhJetPtZLeadingCluster","Jet z_{Leading,cluster} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4482  fhJetPtZLeadingCluster->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4483  fhJetPtZLeadingCluster->GetYaxis()->SetTitle("z_{Leading,cluster}");
4484  fhJetPtZLeadingCluster->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{cluster}");
4485  fhJetPtZLeadingCluster->Sumw2();
4486 
4487  // Event Centralities vs Leading Jet Pt
4489  fhEventCentralityVsZNA->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
4490  fhEventCentralityVsZNA->GetYaxis()->SetTitle("Centrality (ZNA)");
4491  fhEventCentralityVsZNA->GetZaxis()->SetTitle("Probability Density");
4492  fhEventCentralityVsZNA->Sumw2();
4493 
4494  // 3D Histos
4495  if (fDo3DHistos == kTRUE)
4496  {
4497  // Jet pT, Eta, Phi
4498  fhJetPtEtaPhi = new TH3F("fhJetPtEtaPhi","Jet p_{T} vs #eta-#varphi",fPtBins,fPtLow,fPtUp,TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
4499  fhJetPtEtaPhi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4500  fhJetPtEtaPhi->GetYaxis()->SetTitle("#eta");
4501  fhJetPtEtaPhi->GetZaxis()->SetTitle("#varphi");
4502  fhJetPtEtaPhi->Sumw2();
4503 
4504  EventCentralityString = Form("%s vs ZNA Centrality vs Leading Jet p_{T} ",fCentralityTag.Data());
4506  fhEventCentralityVsZNAPt->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
4507  fhEventCentralityVsZNAPt->GetYaxis()->SetTitle("Centrality (ZNA)");
4508  fhEventCentralityVsZNAPt->GetZaxis()->SetTitle("Leading Jet p_{T} (GeV/c)");
4509  fhEventCentralityVsZNAPt->Sumw2();
4510 
4511  fOutput->Add(fhJetPtEtaPhi);
4512  fOutput->Add(fhEventCentralityVsZNAPt);
4513  }
4514 
4515  // Neutral Energy Fraction Histograms & QA
4516  if (fDoNEFQAPlots==kTRUE)
4517  {
4518  fNEFOutput = new TList();
4519  fNEFOutput->SetOwner();
4520  fNEFOutput->SetName("ListNEFQAPlots");
4521 
4522  fhJetPtNEF = new TH2F("fhJetPtNEF","Jet p_{T} vs NEF",fPtBins,fPtLow,fPtUp,fNEFBins,fNEFLow,fNEFUp);
4523  fhJetPtNEF->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4524  fhJetPtNEF->GetYaxis()->SetTitle("Neutral Energy Fraction");
4525  fhJetPtNEF->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T}dNEF");
4526  fhJetPtNEF->Sumw2();
4527 
4528  SetNEFJetDimensions(10); // Order: nef,Jet Pt,Eta,Phi,Centrality,Constituent mult,Charged mult, Neutral mult, z_leading
4529  SetNEFClusterDimensions(11); // Order: nef,Jet Pt,Eta,Phi,Centrality,F_Cross,z_leading,t_cell,deltat_cell,Cell Count, E_Cluster
4530 
4531  Int_t dimJetBins[fnDimJet];
4532  Double_t dimJetLow[fnDimJet];
4533  Double_t dimJetUp[fnDimJet];
4534 
4535  Int_t dimClusterBins[fnDimCluster];
4536  Double_t dimClusterLow[fnDimCluster];
4537  Double_t dimClusterUp[fnDimCluster];
4538 
4539  // Establish dimensinality bin counts and bounds
4540  // NEF
4541  dimJetBins[0] = dimClusterBins[0] = fNEFBins;
4542  dimJetLow[0] = dimClusterLow[0] = fNEFLow;
4543  dimJetUp[0] = dimClusterUp[0] = fNEFUp;
4544 
4545  // Jet Pt
4546  dimJetBins[1] = dimClusterBins[1] = fPtBins;
4547  dimJetLow[1] = dimClusterLow[1] = fPtLow;
4548  dimJetUp[1] = dimClusterUp[1] = fPtUp;
4549 
4550  // Eta-Phi
4551  dimJetBins[2] = dimJetBins[3] = dimClusterBins[2] = dimClusterBins[3] = TCBins;