AliPhysics  c69cd47 (c69cd47)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskFullpAJets_Eli_Mod.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 "AliClusterContainer.h"
24 #include "AliAnalysisTaskEmcal.h"
26 #include "AliAnalysisManager.h"
27 #include "AliStack.h"
28 #include "AliESDtrackCuts.h"
29 #include "AliESDEvent.h"
30 #include "AliESDInputHandler.h"
31 #include "AliAODEvent.h"
32 #include "AliVEvent.h"
33 #include "AliMCEvent.h"
34 #include "AliVTrack.h"
35 #include "AliVCluster.h"
36 #include "AliEmcalJet.h"
37 #include "AliEMCALGeometry.h"
38 #include "AliEMCALRecoUtils.h"
39 #include "AliVCaloCells.h"
40 #include "AliPicoTrack.h"
41 #include "AliAnalysisTaskEmcal.h"
43 #include "Rtypes.h"
44 
45 using std::cout;
46 using std::endl;
47 
49 
50 //________________________________________________________________________
53 
54 fOutput(0),
55 flTrack(0),
56 flCluster(0),
57 fhTrackPt(0),
58 fhTrackEta(0),
59 fhTrackPhi(0),
60 fhGlobalTrackPt(0),
61 fhGlobalTrackEta(0),
62 fhGlobalTrackPhi(0),
63 fhComplementaryTrackPt(0),
64 fhComplementaryTrackEta(0),
65 fhComplementaryTrackPhi(0),
66 fhClusterPt(0),
67 fhClusterEta(0),
68 fhClusterPhi(0),
69 fhCentrality(0),
70 fhEMCalCellCounts(0),
71 
72 fhChargeAndNeutralEvents(0),
73 fhChargeOnlyEvents(0),
74 fhNeutralOnlyEvents(0),
75 fhNothingEvents(0),
76 fhEMCalChargeAndNeutralEvents(0),
77 fhEMCalChargeOnlyEvents(0),
78 fhEMCalNeutralOnlyEvents(0),
79 fhEMCalNothingEvents(0),
80 
81 fhTrackEtaPhi(0),
82 fhTrackPhiPt(0),
83 fhTrackEtaPt(0),
84 fhGlobalTrackEtaPhi(0),
85 fhGlobalTrackPhiPt(0),
86 fhGlobalTrackEtaPt(0),
87 fhComplementaryTrackEtaPhi(0),
88 fhComplementaryTrackPhiPt(0),
89 fhComplementaryTrackEtaPt(0),
90 fhClusterEtaPhi(0),
91 fhClusterPhiPt(0),
92 fhClusterEtaPt(0),
93 fhEMCalEventMult(0),
94 fhTPCEventMult(0),
95 fhEMCalTrackEventMult(0),
96 
97 fpEMCalEventMult(0),
98 fpTPCEventMult(0),
99 
100 fpTrackPtProfile(0),
101 fpClusterPtProfile(0),
102 
103 fpFullJetEDProfile(0),
104 fpChargedJetEDProfile(0),
105 fpChargedJetEDProfileScaled(0),
106 
107 fTPCRawJets(0),
108 fEMCalRawJets(0),
109 fRhoChargedCMSScale(0),
110 fRhoChargedScale(0),
111 fRhoFull0(0),
112 fRhoFull1(0),
113 fRhoFull2(0),
114 fRhoFullN(0),
115 fRhoFullDijet(0),
116 fRhoFullkT(0),
117 fRhoFullCMS(0),
118 fRhoCharged0(0),
119 fRhoCharged1(0),
120 fRhoCharged2(0),
121 fRhoChargedN(0),
122 fRhoChargedkT(0),
123 fRhoChargedkTScale(0),
124 fRhoChargedCMS(0),
125 
126 fTPCJet(0),
127 fTPCFullJet(0),
128 fTPCOnlyJet(0),
129 fTPCJetUnbiased(0),
130 fTPCkTFullJet(0),
131 fEMCalJet(0),
132 fEMCalFullJet(0),
133 fEMCalPartJet(0),
134 fEMCalPartJetUnbiased(0),
135 fEMCalkTFullJet(0),
136 
137 fIsInitialized(0),
138 fRJET(4),
139 fnEvents(0),
140 fnEventsCharged(0),
141 fnDiJetEvents(0),
142 fEvent(0),
143 fRecoUtil(0),
144 fEMCALGeometry(0),
145 fCells(0),
146 fDoNEF(0),
147 fDoNEFSignalOnly(1),
148 fSignalTrackBias(0),
149 fTrackQA(0),
150 fClusterQA(0),
151 fCalculateRhoJet(0),
152 fDoVertexRCut(0),
153 fMCPartLevel(0),
154 fDoTHnSparse(0),
155 fDoJetRhoDensity(0),
156 fDo3DHistos(0),
157 fEMCalPhiMin(1.39626),
158 fEMCalPhiMax(3.26377),
159 fEMCalPhiTotal(1.86750),
160 fEMCalEtaMin(-0.7),
161 fEMCalEtaMax(0.7),
162 fEMCalEtaTotal(1.4),
163 fEMCalArea(2.61450),
164 fTPCPhiMin(0),
165 fTPCPhiMax(6.28319),
166 fTPCPhiTotal(6.28319),
167 fTPCEtaMin(-0.9),
168 fTPCEtaMax(0.9),
169 fTPCEtaTotal(1.8),
170 fTPCArea(11.30973),
171 fParticlePtLow(0.0),
172 fParticlePtUp(200.0),
173 fParticlePtBins(200),
174 fJetR(0.4),
175 fJetRAccept(0.4),
176 fFullEDJetR(0.6),
177 fChargedEDJetR(0.8),
178 fJetRForRho(0.5),
179 fJetAreaCutFrac(0.6),
180 fJetAreaThreshold(0.30159),
181 fnEMCalCells(12288),
182 fScaleFactor(1.28),
183 fNColl(7),
184 fTrackMinPt(0.15),
185 fClusterMinPt(0.3),
186 fNEFSignalJetCut(0.9),
187 fCentralityTag("V0A"),
188 fCentralityBins(10),
189 fCentralityLow(0),
190 fCentralityUp(100),
191 fEventCentrality(0),
192 fRhoFull(0),
193 fRhoCharged(0),
194 fnTracks(0),
195 fnEMCalTracks(0),
196 fnClusters(0),
197 fnCaloClusters(0),
198 fnAKTFullJets(0),
199 fnAKTChargedJets(0),
200 fnKTFullJets(0),
201 fnKTChargedJets(0),
202 fnBckgClusters(0),
203 fTPCJetThreshold(5),
204 fEMCalJetThreshold(5),
205 fVertexWindow(10),
206 fVertexMaxR(1),
207 fTrackName(0),
208 fClusName(0),
209 fkTChargedName(0),
210 fAkTChargedName(0),
211 fkTFullName(0),
212 fAkTFullName(0),
213 fOrgTracks(0),
214 fOrgClusters(0),
215 fmyAKTFullJets(0),
216 fmyAKTChargedJets(0),
217 fmyKTFullJets(0),
218 fmyKTChargedJets(0),
219 fmyTracks(0),
220 fmyClusters(0),
221 fEMCalRCBckgFluc(0),
222 fTPCRCBckgFluc(0),
223 fEMCalRCBckgFlucSignal(0),
224 fTPCRCBckgFlucSignal(0),
225 fEMCalRCBckgFlucNColl(0),
226 fTPCRCBckgFlucNColl(0)
227 {
228  // Dummy constructor ALWAYS needed for I/O.
229  fVertex[0]=0.0,fVertex[1]=0.0,fVertex[2]=0.0;
230 }
231 
232 //________________________________________________________________________
235 
236  fOutput(0),
237  flTrack(0),
238  flCluster(0),
239  fhTrackPt(0),
240  fhTrackEta(0),
241  fhTrackPhi(0),
242  fhGlobalTrackPt(0),
243  fhGlobalTrackEta(0),
244  fhGlobalTrackPhi(0),
245  fhComplementaryTrackPt(0),
246  fhComplementaryTrackEta(0),
247  fhComplementaryTrackPhi(0),
248  fhClusterPt(0),
249  fhClusterEta(0),
250  fhClusterPhi(0),
251  fhCentrality(0),
252  fhEMCalCellCounts(0),
253 
254  fhChargeAndNeutralEvents(0),
255  fhChargeOnlyEvents(0),
256  fhNeutralOnlyEvents(0),
257  fhNothingEvents(0),
258  fhEMCalChargeAndNeutralEvents(0),
259  fhEMCalChargeOnlyEvents(0),
260  fhEMCalNeutralOnlyEvents(0),
261  fhEMCalNothingEvents(0),
262 
263  fhTrackEtaPhi(0),
264  fhTrackPhiPt(0),
265  fhTrackEtaPt(0),
266  fhGlobalTrackEtaPhi(0),
267  fhGlobalTrackPhiPt(0),
268  fhGlobalTrackEtaPt(0),
269  fhComplementaryTrackEtaPhi(0),
270  fhComplementaryTrackPhiPt(0),
271  fhComplementaryTrackEtaPt(0),
272  fhClusterEtaPhi(0),
273  fhClusterPhiPt(0),
274  fhClusterEtaPt(0),
275  fhEMCalEventMult(0),
276  fhTPCEventMult(0),
277  fhEMCalTrackEventMult(0),
278 
279  fpEMCalEventMult(0),
280  fpTPCEventMult(0),
281 
282  fpTrackPtProfile(0),
283  fpClusterPtProfile(0),
284 
285  fpFullJetEDProfile(0),
286  fpChargedJetEDProfile(0),
287  fpChargedJetEDProfileScaled(0),
288 
289  fTPCRawJets(0),
290  fEMCalRawJets(0),
291  fRhoChargedCMSScale(0),
292  fRhoChargedScale(0),
293  fRhoFull0(0),
294  fRhoFull1(0),
295  fRhoFull2(0),
296  fRhoFullN(0),
297  fRhoFullDijet(0),
298  fRhoFullkT(0),
299  fRhoFullCMS(0),
300  fRhoCharged0(0),
301  fRhoCharged1(0),
302  fRhoCharged2(0),
303  fRhoChargedN(0),
304  fRhoChargedkT(0),
305  fRhoChargedkTScale(0),
306  fRhoChargedCMS(0),
307 
308  fTPCJet(0),
309  fTPCFullJet(0),
310  fTPCOnlyJet(0),
311  fTPCJetUnbiased(0),
312  fTPCkTFullJet(0),
313  fEMCalJet(0),
314  fEMCalFullJet(0),
315  fEMCalPartJet(0),
316  fEMCalPartJetUnbiased(0),
317  fEMCalkTFullJet(0),
318 
319  fIsInitialized(0),
320  fRJET(4),
321  fnEvents(0),
322  fnEventsCharged(0),
323  fnDiJetEvents(0),
324  fEvent(0),
325  fRecoUtil(0),
326  fEMCALGeometry(0),
327  fCells(0),
328  fDoNEF(0),
329  fDoNEFSignalOnly(1),
330  fSignalTrackBias(0),
331  fTrackQA(0),
332  fClusterQA(0),
333  fCalculateRhoJet(0),
334  fDoVertexRCut(0),
335  fMCPartLevel(0),
336  fDoTHnSparse(0),
337  fDoJetRhoDensity(0),
338  fDo3DHistos(0),
339  fEMCalPhiMin(1.39626),
340  fEMCalPhiMax(3.26377),
341  fEMCalPhiTotal(1.86750),
342  fEMCalEtaMin(-0.7),
343  fEMCalEtaMax(0.7),
344  fEMCalEtaTotal(1.4),
345  fEMCalArea(2.61450),
346  fTPCPhiMin(0),
347  fTPCPhiMax(6.28319),
348  fTPCPhiTotal(6.28319),
349  fTPCEtaMin(-0.9),
350  fTPCEtaMax(0.9),
351  fTPCEtaTotal(1.8),
352  fTPCArea(11.30973),
353  fParticlePtLow(0.0),
354  fParticlePtUp(200.0),
355  fParticlePtBins(2000),
356  fJetR(0.4),
357  fJetRAccept(0.4),
358  fFullEDJetR(0.6),
359  fChargedEDJetR(0.8),
360  fJetRForRho(0.5),
361  fJetAreaCutFrac(0.6),
362  fJetAreaThreshold(0.30159),
363  fnEMCalCells(12288),
364  fScaleFactor(1.28),
365  fNColl(7),
366  fTrackMinPt(0.15),
367  fClusterMinPt(0.3),
368  fNEFSignalJetCut(0.9),
369  fCentralityTag("V0A"),
370  fCentralityBins(10),
371  fCentralityLow(0),
372  fCentralityUp(100),
373  fEventCentrality(0),
374  fRhoFull(0),
375  fRhoCharged(0),
376  fnTracks(0),
377  fnEMCalTracks(0),
378  fnClusters(0),
379  fnCaloClusters(0),
380  fnAKTFullJets(0),
381  fnAKTChargedJets(0),
382  fnKTFullJets(0),
383  fnKTChargedJets(0),
384  fnBckgClusters(0),
385  fTPCJetThreshold(5),
386  fEMCalJetThreshold(5),
387  fVertexWindow(10),
388  fVertexMaxR(1),
389  fTrackName(0),
390  fClusName(0),
391  fkTChargedName(0),
392  fAkTChargedName(0),
393  fkTFullName(0),
394  fAkTFullName(0),
395  fOrgTracks(0),
396  fOrgClusters(0),
397  fmyAKTFullJets(0),
398  fmyAKTChargedJets(0),
399  fmyKTFullJets(0),
400  fmyKTChargedJets(0),
401  fmyTracks(0),
402  fmyClusters(0),
403  fEMCalRCBckgFluc(0),
404  fTPCRCBckgFluc(0),
405  fEMCalRCBckgFlucSignal(0),
406  fTPCRCBckgFlucSignal(0),
407  fEMCalRCBckgFlucNColl(0),
408  fTPCRCBckgFlucNColl(0)
409 {
410  // Constructor
411  // Define input and output slots here (never in the dummy constructor)
412  // Input slot #0 works with a TChain - it is connected to the default input container
413  // Output slot #1 writes into a TH1 container
414  fVertex[0]=0.0,fVertex[1]=0.0,fVertex[2]=0.0;
415 
416  DefineOutput(1,TList::Class()); // for output list
417 }
418 
419 //________________________________________________________________________
421 {
422  // Destructor. Clean-up the output list, but not the histograms that are put inside
423  // (the list is owner and will clean-up these histograms). Protect in PROOF case.
424  if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
425  {
426  delete fOutput;
427  }
428 }
429 
430 //________________________________________________________________________
432 {
433  // Create histograms
434  // Called once (on the worker node)
435  fIsInitialized=kFALSE;
436  fOutput = new TList();
437  fOutput->SetOwner(); // IMPORTANT!
438 
439  // Initialize Global Variables
440  fnEvents=0;
441  fnEventsCharged=0;
442  fnDiJetEvents=0;
443 
444  // fRJET=4 -> fJetR=0.4 && fRJET=25 -> fJetR=0.25, but for writing files, should be 4 and 25 respectively
445  if (fRJET>10)
446  {
447  fJetR=(Double_t)fRJET/100.0;
448  }
449  else
450  {
451  fJetR=(Double_t)fRJET/10.0;
452  }
453  fJetRForRho=0.5;
454 
455  fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
456  fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
458  fEMCalEtaMin=-0.7;
459  fEMCalEtaMax=0.7;
462 
463  fTPCPhiMin=(0/(double)360)*2*TMath::Pi();
464  fTPCPhiMax=(360/(double)360)*2*TMath::Pi();
466  fTPCEtaMin=-0.9;
467  fTPCEtaMax=0.9;
470 
471  fParticlePtLow=0.0;
472  fParticlePtUp=200.0;
474 
475  fCentralityBins=10;
476  fCentralityLow=0.0;
477  fCentralityUp=100.0;
478  Int_t CentralityBinMult=10;
479 
480  fJetAreaCutFrac =0.6; // Fudge factor for selecting on jets with threshold Area or higher
482  fTPCJetThreshold=5.0; // Threshold required for an Anti-kT Charged jet to be considered a "true" jet in TPC
483  fEMCalJetThreshold=5.0; // Threshold required for an Anti-kT Charged+Neutral jet to be considered a "true" jet in EMCal
484  fVertexWindow=10.0;
485  fVertexMaxR=1.0;
486 
487  fnBckgClusters=1;
488  fEMCalRCBckgFluc = new Double_t[fnBckgClusters];
489  fTPCRCBckgFluc = new Double_t[fnBckgClusters];
490  fEMCalRCBckgFlucSignal = new Double_t[fnBckgClusters];
491  fTPCRCBckgFlucSignal = new Double_t[fnBckgClusters];
492  fEMCalRCBckgFlucNColl = new Double_t[fnBckgClusters];
493  fTPCRCBckgFlucNColl = new Double_t[fnBckgClusters];
494  for (Int_t i=0;i<fnBckgClusters;i++)
495  {
496  fEMCalRCBckgFluc[i]=0.0;
497  fTPCRCBckgFluc[i]=0.0;
498  fEMCalRCBckgFlucSignal[i]=0.0;
499  fTPCRCBckgFlucSignal[i]=0.0;
500  fEMCalRCBckgFlucNColl[i]=0.0;
501  fTPCRCBckgFlucNColl[i]=0.0;
502  }
503 
504  fnEMCalCells=12288; // sMods 1-10 have 24x48 cells, sMods 11&12 have 8x48 cells...
505 
506  Int_t TCBins=100;
507  Int_t multBins=200;
508  Double_t multLow=0;
509  Double_t multUp=200;
510 
511  // Track QA Plots
512  if (fTrackQA==kTRUE)
513  {
514  flTrack = new TList();
515  flTrack->SetName("TrackQA");
516 
517  // Hybrid Tracks
518  fhTrackPt = new TH1F("fhTrackPt","p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
519  fhTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
520  fhTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
521  fhTrackPt->Sumw2();
522 
523  fhTrackPhi = new TH1F("fhTrackPhi","#varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
524  fhTrackPhi->GetXaxis()->SetTitle("#varphi");
525  fhTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
526  fhTrackPhi->Sumw2();
527 
528  fhTrackEta = new TH1F("fhTrackEta","#eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
529  fhTrackEta->GetXaxis()->SetTitle("#eta");
530  fhTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
531  fhTrackEta->Sumw2();
532 
533  fhTrackEtaPhi = new TH2F("fhTrackEtaPhi","#eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
534  fhTrackEtaPhi->GetXaxis()->SetTitle("#eta");
535  fhTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
536  fhTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
537  fhTrackEtaPhi->Sumw2();
538 
539  fhTrackPhiPt = new TH2F("fhTrackPhiPt","#varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
540  fhTrackPhiPt->GetXaxis()->SetTitle("#varphi");
541  fhTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
542  fhTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
543  fhTrackPhiPt->Sumw2();
544 
545  fhTrackEtaPt = new TH2F("fhTrackEtaPt","#eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
546  fhTrackEtaPt->GetXaxis()->SetTitle("#varphi");
547  fhTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
548  fhTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
549  fhTrackEtaPt->Sumw2();
550 
551  // Global Tracks
552  fhGlobalTrackPt = new TH1F("fhGlobalTrackPt","Global p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
553  fhGlobalTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
554  fhGlobalTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
555  fhGlobalTrackPt->Sumw2();
556 
557  fhGlobalTrackPhi = new TH1F("fhGlobalTrackPhi","Global #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
558  fhGlobalTrackPhi->GetXaxis()->SetTitle("#varphi");
559  fhGlobalTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
560  fhGlobalTrackPhi->Sumw2();
561 
562  fhGlobalTrackEta = new TH1F("fhGlobalTrackEta","Global #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
563  fhGlobalTrackEta->GetXaxis()->SetTitle("#eta");
564  fhGlobalTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
565  fhGlobalTrackEta->Sumw2();
566 
567  fhGlobalTrackEtaPhi = new TH2F("fhGlobalTrackEtaPhi","Global #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
568  fhGlobalTrackEtaPhi->GetXaxis()->SetTitle("#eta");
569  fhGlobalTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
570  fhGlobalTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
571  fhGlobalTrackEtaPhi->Sumw2();
572 
573  fhGlobalTrackPhiPt = new TH2F("fhGlobalTrackPhiPt","Global #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
574  fhGlobalTrackPhiPt->GetXaxis()->SetTitle("#varphi");
575  fhGlobalTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
576  fhGlobalTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
577  fhGlobalTrackPhiPt->Sumw2();
578 
579  fhGlobalTrackEtaPt = new TH2F("fhGlobalTrackEtaPt","Global #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
580  fhGlobalTrackEtaPt->GetXaxis()->SetTitle("#varphi");
581  fhGlobalTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
582  fhGlobalTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
583  fhGlobalTrackEtaPt->Sumw2();
584 
585  // Complementary Tracks
586  fhComplementaryTrackPt = new TH1F("fhComplementaryTrackPt","Complementary p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
587  fhComplementaryTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
588  fhComplementaryTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
589  fhComplementaryTrackPt->Sumw2();
590 
591  fhComplementaryTrackPhi = new TH1F("fhComplementaryTrackPhi","Complementary #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
592  fhComplementaryTrackPhi->GetXaxis()->SetTitle("#varphi");
593  fhComplementaryTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
594  fhComplementaryTrackPhi->Sumw2();
595 
596  fhComplementaryTrackEta = new TH1F("fhComplementaryTrackEta","Complementary #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
597  fhComplementaryTrackEta->GetXaxis()->SetTitle("#eta");
598  fhComplementaryTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
599  fhComplementaryTrackEta->Sumw2();
600 
601  fhComplementaryTrackEtaPhi = new TH2F("fhComplementaryTrackEtaPhi","Complementary #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
602  fhComplementaryTrackEtaPhi->GetXaxis()->SetTitle("#eta");
603  fhComplementaryTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
604  fhComplementaryTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
606 
607  fhComplementaryTrackPhiPt = new TH2F("fhComplementaryTrackPhiPt","Complementary #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
608  fhComplementaryTrackPhiPt->GetXaxis()->SetTitle("#varphi");
609  fhComplementaryTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
610  fhComplementaryTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
611  fhComplementaryTrackPhiPt->Sumw2();
612 
613  fhComplementaryTrackEtaPt = new TH2F("fhComplementaryTrackEtaPt","Complementary #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
614  fhComplementaryTrackEtaPt->GetXaxis()->SetTitle("#varphi");
615  fhComplementaryTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
616  fhComplementaryTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
617  fhComplementaryTrackEtaPt->Sumw2();
618 
619  fhTPCEventMult = new TH2F("fhTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
620  fhTPCEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
621  fhTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
622  fhTPCEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Charged}");
623  fhTPCEventMult->Sumw2();
624 
625  fhEMCalTrackEventMult = new TH2F("fhEMCalTrackEventMult","EMCal Track Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
626  fhEMCalTrackEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
627  fhEMCalTrackEventMult->GetYaxis()->SetTitle("Multiplicity");
628  fhEMCalTrackEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Neutral}");
629  fhEMCalTrackEventMult->Sumw2();
630 
631  fpTPCEventMult = new TProfile("fpTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
632  fpTPCEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
633  fpTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
634 
635  // QA::2D Energy Density Profiles for Tracks and Clusters
636  fpTrackPtProfile = new TProfile2D("fpTrackPtProfile","2D Profile of track pT density throughout the TPC",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
637  fpTrackPtProfile->GetXaxis()->SetTitle("#eta");
638  fpTrackPtProfile->GetYaxis()->SetTitle("#varphi");
639  fpTrackPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
640 
641  flTrack->Add(fhTrackPt);
642  flTrack->Add(fhTrackEta);
643  flTrack->Add(fhTrackPhi);
644  flTrack->Add(fhTrackEtaPhi);
645  flTrack->Add(fhTrackPhiPt);
646  flTrack->Add(fhTrackEtaPt);
647  flTrack->Add(fhGlobalTrackPt);
659  flTrack->Add(fhTPCEventMult);
661  flTrack->Add(fpTPCEventMult);
663  fOutput->Add(flTrack);
664  }
665 
666  if (fClusterQA==kTRUE)
667  {
668  flCluster = new TList();
669  flCluster->SetName("ClusterQA");
670 
671  fhClusterPt = new TH1F("fhClusterPt","p_{T} distribution of clusters in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
672  fhClusterPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
673  fhClusterPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
674  fhClusterPt->Sumw2();
675 
676  fhClusterPhi = new TH1F("fhClusterPhi","#varphi distribution of clusters in event",TCBins,fTPCPhiMin,fTPCPhiMax);
677  fhClusterPhi->GetXaxis()->SetTitle("#varphi");
678  fhClusterPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
679  fhClusterPhi->Sumw2();
680 
681  fhClusterEta = new TH1F("fhClusterEta","#eta distribution of clusters in event",TCBins,fTPCEtaMin,fTPCEtaMax);
682  fhClusterEta->GetXaxis()->SetTitle("#eta");
683  fhClusterEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
684  fhClusterEta->Sumw2();
685 
686  fhClusterEtaPhi = new TH2F("fhClusterEtaPhi","#eta-#varphi distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
687  fhClusterEtaPhi->GetXaxis()->SetTitle("#eta");
688  fhClusterEtaPhi->GetYaxis()->SetTitle("#varphi");
689  fhClusterEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
690  fhClusterEtaPhi->Sumw2();
691 
692  fhClusterPhiPt = new TH2F("fhClusterPhiPt","#varphi-p_{T} distribution of clusters in event",TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
693  fhClusterPhiPt->GetXaxis()->SetTitle("#varphi");
694  fhClusterPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
695  fhClusterPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
696  fhClusterPhiPt->Sumw2();
697 
698  fhClusterEtaPt = new TH2F("fhClusterEtaPt","#eta-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
699  fhClusterEtaPt->GetXaxis()->SetTitle("#varphi");
700  fhClusterEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
701  fhClusterEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
702  fhClusterEtaPt->Sumw2();
703 
704  fhEMCalEventMult = new TH2F("fhEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
705  fhEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
706  fhEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
707  fhEMCalEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Neutral}");
708  fhEMCalEventMult->Sumw2();
709 
710  fpEMCalEventMult = new TProfile("fpEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
711  fpEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
712  fpEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
713 
714  fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
715  fpClusterPtProfile->GetXaxis()->SetTitle("#eta");
716  fpClusterPtProfile->GetYaxis()->SetTitle("#varphi");
717  fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
718 
719  fhEMCalCellCounts = new TH1F("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
720  fhEMCalCellCounts->GetXaxis()->SetTitle("Absoulute Cell Id");
721  fhEMCalCellCounts->GetYaxis()->SetTitle("Counts per Event");
722  fhEMCalCellCounts->Sumw2();
723 
724  flCluster->Add(fhClusterPt);
725  flCluster->Add(fhClusterEta);
726  flCluster->Add(fhClusterPhi);
734  fOutput->Add(flCluster);
735  }
736 
737  if (fCalculateRhoJet>=0 && fMCPartLevel==kFALSE) // Default Rho & Raw Jet Spectra
738  {
740 
742  fRhoChargedCMSScale->SetSignalTrackPtBias(fSignalTrackBias);
744  fOutput->Add(fRhoChargedCMSScale->GetOutputHistos());
745 
746  }
747  if (fCalculateRhoJet>=1) // Basic Rho & Raw Jet Spectra
748  {
751 
752  if (fMCPartLevel==kTRUE)
753  {
756  fRhoChargedCMS->SetSignalTrackPtBias(fSignalTrackBias);
757 
759  fOutput->Add(fRhoChargedCMS->GetOutputHistos());
760  }
761 
762  if (fDoJetRhoDensity == kTRUE && fDo3DHistos == kTRUE)
763  {
764  Double_t fullEDR = fFullEDJetR *100;
765  Double_t chargedEDR = fChargedEDJetR *100;
766  const Int_t fullEDBins = (Int_t) fullEDR;
767  const Int_t chargedEDBins = (Int_t) chargedEDR;
768 
769  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);
770  fpFullJetEDProfile->GetXaxis()->SetTitle("p_{T,jet}^{ch+em} (GeV)");
771  fpFullJetEDProfile->GetYaxis()->SetTitle(fCentralityTag.Data());
772  fpFullJetEDProfile->GetZaxis()->SetTitle("R");
773 
774  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);
775  fpChargedJetEDProfile->GetXaxis()->SetTitle("p_{T,jet}^{ch} (GeV)");
776  fpChargedJetEDProfile->GetYaxis()->SetTitle(fCentralityTag.Data());
777  fpChargedJetEDProfile->GetZaxis()->SetTitle("R");
778 
779  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);
780  fpChargedJetEDProfileScaled->GetXaxis()->SetTitle("p_{T,jet}^{ch} (GeV)");
781  fpChargedJetEDProfileScaled->GetYaxis()->SetTitle(fCentralityTag.Data());
782  fpChargedJetEDProfileScaled->GetZaxis()->SetTitle("R");
783 
787  }
789  }
790  if (fCalculateRhoJet>=2 && fMCPartLevel==kFALSE) // Charged Rho & Raw Jet Spectra
791  {
794  fRhoChargedCMS->SetSignalTrackPtBias(fSignalTrackBias);
795 
797  fOutput->Add(fRhoChargedCMS->GetOutputHistos());
798  }
799  if (fCalculateRhoJet>=3 && fMCPartLevel==kFALSE) // Basic Rho & Raw Jet Spectra
800  {
801  fRhoFull0 = new AlipAJetHistos("RhoFull0",fCentralityTag.Data());
803  fRhoFull1 = new AlipAJetHistos("RhoFull1",fCentralityTag.Data());
805  fRhoFull2 = new AlipAJetHistos("RhoFull2",fCentralityTag.Data());
807  fRhoFullN = new AlipAJetHistos("RhoFullN",fCentralityTag.Data());
809  fRhoFullDijet = new AlipAJetHistos("RhoFullDijet",fCentralityTag.Data());
811  fRhoFullkT = new AlipAJetHistos("RhoFullkT",fCentralityTag.Data());
813  fRhoFullCMS = new AlipAJetHistos("RhoFullCMS",fCentralityTag.Data());
815  fRhoCharged0 = new AlipAJetHistos("RhoCharged0",fCentralityTag.Data());
817  fRhoCharged1 = new AlipAJetHistos("RhoCharged1",fCentralityTag.Data());
819  fRhoCharged2 = new AlipAJetHistos("RhoCharged2",fCentralityTag.Data());
821  fRhoChargedN = new AlipAJetHistos("RhoChargedN",fCentralityTag.Data());
823  fRhoChargedkT = new AlipAJetHistos("RhoChargedkT",fCentralityTag.Data());
825  fRhoChargedkTScale = new AlipAJetHistos("RhoChargedkTScale",fCentralityTag.Data());
827 
841  }
842 
843  fhCentrality = new TH1F("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
844  fhCentrality->GetXaxis()->SetTitle(fCentralityTag.Data());
845  fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
846  fhCentrality->Sumw2();
847 
848  fOutput->Add(fhCentrality);
849 
850  if (fMCPartLevel == kFALSE)
851  {
852  fhChargeAndNeutralEvents = new TH1F("fhChargeAndNeutralEvents","Events with n_{tracks}>0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
853  fhChargeAndNeutralEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
854  fhChargeAndNeutralEvents->GetYaxis()->SetTitle("1/N_{Events}");
855  fhChargeAndNeutralEvents->Sumw2();
856 
857  fhChargeOnlyEvents = new TH1F("fhChargeOnlyEvents","Events with n_{tracks}>0 & n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
858  fhChargeOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
859  fhChargeOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
860  fhChargeOnlyEvents->Sumw2();
861 
862  fhNeutralOnlyEvents = new TH1F("fhNeutralOnlyEvents","Events with n_{tracks}=0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
863  fhNeutralOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
864  fhNeutralOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
865  fhNeutralOnlyEvents->Sumw2();
866 
867  fhNothingEvents = new TH1F("fhNothingEvents","Events with n_{tracks}=n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
868  fhNothingEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
869  fhNothingEvents->GetYaxis()->SetTitle("1/N_{Events}");
870  fhNothingEvents->Sumw2();
871 
872  fhEMCalChargeAndNeutralEvents = new TH1F("fhEMCalChargeAndNeutralEvents","Events with n_{emcal tracks}>0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
873  fhEMCalChargeAndNeutralEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
874  fhEMCalChargeAndNeutralEvents->GetYaxis()->SetTitle("1/N_{Events}");
876 
877  fhEMCalChargeOnlyEvents = new TH1F("fhEMCalChargeOnlyEvents","Events with n_{emcal tracks}>0 & n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
878  fhEMCalChargeOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
879  fhEMCalChargeOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
880  fhEMCalChargeOnlyEvents->Sumw2();
881 
882  fhEMCalNeutralOnlyEvents = new TH1F("fhEMCalNeutralOnlyEvents","Events with n_{tracks}>0 & n_{emcal tracks}=0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
883  fhEMCalNeutralOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
884  fhEMCalNeutralOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
885  fhEMCalNeutralOnlyEvents->Sumw2();
886 
887  fhEMCalNothingEvents = new TH1F("fhEMCalNothingEvents","Events with n_{emcal tracks}=n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
888  fhEMCalNothingEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
889  fhEMCalNothingEvents->GetYaxis()->SetTitle("1/N_{Events}");
890  fhEMCalNothingEvents->Sumw2();
891 
895  fOutput->Add(fhNothingEvents);
900  }
901 
902  // Post data for ALL output slots >0 here,
903  // To get at least an empty histogram
904  // 1 is the outputnumber of a certain weg of task 1
905  PostData(1, fOutput);
906 }
907 
909 {
911 // cout<<"Start of ...FullpAJets_Eli_Mod::UserExecOnce()"<<endl;
912  // Get the event tracks
913  fOrgTracks = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fTrackName));
914 
915  // Get the event caloclusters
916  fOrgClusters = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fClusName));
917 
918  if(fAkTFullName!="")fmyAKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fAkTFullName));
919 
920  fIsInitialized=kTRUE;
921 // cout<<"End of ...FullpAJets_Eli_Mod::UserExecOnce()"<<endl;
922 }
923 //________________________________________________________________________
925 {
926 // cout<<"start of ...FullpAJets_Eli_Mod::UserExec(Option_t *)"<<endl;
927 
928  if (fIsInitialized==kFALSE)
929  {
930  UserExecOnce();
931  }
932  // Get pointer to reconstructed event
933  fEvent = InputEvent();
934  if (!fEvent)
935  {
936  AliError("Pointer == 0, this can not happen!");
937  return;
938  }
939  // Check if Event is selected (for triggers)
940  if (IsEventSelected() == kFALSE)
941  {
942  return;
943  }
944  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent);
945  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fEvent);
946 
947  fRecoUtil = new AliEMCALRecoUtils();
948  fEMCALGeometry = AliEMCALGeometry::GetInstance();
949 
950  if (esd)
951  {
952  fEventCentrality=esd->GetCentrality()->GetCentralityPercentile(fCentralityTag.Data());
953 
954  if (fDoVertexRCut == kTRUE)
955  {
956  // Vertex R cut not done in AliAnalysisTaskEmcal
957  if (TMath::Sqrt(TMath::Power(esd->GetPrimaryVertex()->GetX(),2)+TMath::Power(esd->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
958  {
959  DeleteJetData(0);
960  return;
961  }
962  }
963 
964  esd->GetPrimaryVertex()->GetXYZ(fVertex);
965  fCells = (AliVCaloCells*) esd->GetEMCALCells();
966  fnCaloClusters = esd->GetNumberOfCaloClusters();
967  }
968  else if (aod)
969  {
970  fEventCentrality=aod->GetCentrality()->GetCentralityPercentile(fCentralityTag.Data());
971 
972  if (fDoVertexRCut == kTRUE)
973  {
974  // Vertex R cut not done in AliAnalysisTaskEmcal
975  if (TMath::Sqrt(TMath::Power(aod->GetPrimaryVertex()->GetX(),2)+TMath::Power(aod->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
976  {
977  DeleteJetData(0);
978  return;
979  }
980  }
981 
982  aod->GetPrimaryVertex()->GetXYZ(fVertex);
983  fCells = (AliVCaloCells*) aod->GetEMCALCells();
984  fnCaloClusters = aod->GetNumberOfCaloClusters();
985  }
986  else
987  {
988  AliError("Cannot get AOD/ESD event. Rejecting Event");
989  DeleteJetData(0);
990  return;
991  }
992 
993  // Make sure Centrality isn't exactly 100% (to avoid bin filling errors in profile plots. Make it 99.99
994  if (fEventCentrality>99.99)
995  {
996  fEventCentrality=99.99;
997  }
998  fhCentrality->Fill(fEventCentrality); // Counts total events
999 
1000  // Count Events
1001  EventCounts();
1002 
1003  // Reject any event that doesn't have any tracks, i.e. TPC is off
1004  if (fnTracks<1)
1005  {
1006  if (fTrackQA==kTRUE)
1007  {
1008  fhTPCEventMult->Fill(fEventCentrality,0.0);
1009  fpTPCEventMult->Fill(fEventCentrality,0.0);
1011  }
1012  AliWarning("No PicoTracks, Rejecting Event");
1013  DeleteJetData(1);
1014  PostData(1,fOutput);
1015  return;
1016  }
1017  if (fnClusters<1)
1018  {
1019  //AliInfo("No Corrected CaloClusters, using only charged jets");
1020  if (fTrackQA==kTRUE)
1021  {
1022  TrackHisto();
1023  }
1024  if (fClusterQA==kTRUE)
1025  {
1028  }
1029  DeleteJetData(2);
1030  fnEventsCharged++;
1031  PostData(1,fOutput);
1032  return;
1033  }
1034 
1035  if (fTrackQA==kTRUE)
1036  {
1037  TrackHisto();
1038  }
1039  if (fClusterQA==kTRUE)
1040  {
1041  ClusterHisto();
1042  }
1043 
1044  InitFullJets();
1045 //ELI check that GenerateEMCalRandomConesPt();
1046 
1047  if (fCalculateRhoJet>=1)
1048  {
1049  if (fDoJetRhoDensity == kTRUE && fDo3DHistos == kTRUE)
1050  {
1052  }
1053  }
1054 
1055  // Delete Dynamic Arrays
1056  DeleteJetData(3);
1057  fnEvents++;
1058  PostData(1,fOutput);
1059 // cout<<"end of ...FullpAJets_Eli_Mod::UserExec(Option_t *)"<<endl;
1060 }
1061 //________________________________________________________________________
1062 void AliAnalysisTaskFullpAJets_Eli_Mod::Terminate(Option_t *) //specify what you want to have done
1063 {
1064  // Called once at the end of the query. Done nothing here.
1065 }
1066 
1070 
1072 {
1073  // Fill a TObjArray with the tracks from a TClonesArray which grabs the picotracks.
1074  Int_t i;
1075  Int_t j=0;
1076 
1077  fmyTracks = new TObjArray();
1078  for (i=0;i<fOrgTracks->GetEntries();i++)
1079  {
1080  AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(i);
1081  if (vtrack->Pt()>=fTrackMinPt)
1082  {
1083  fmyTracks->Add(vtrack);
1084  if (IsInEMCal(vtrack->Phi(),vtrack->Eta()) == kTRUE)
1085  {
1086  j++;
1087  }
1088  }
1089  }
1090  fnTracks = fmyTracks->GetEntries();
1091  fnEMCalTracks = j;
1092 }
1093 
1095 {
1096  // Fill a TObjArray with the clusters from a TClonesArray which grabs the caloclusterscorr.
1097  Int_t i;
1098 
1099  fmyClusters = new TObjArray();
1100  TLorentzVector *cluster_vec = new TLorentzVector;
1101  for (i=0;i<fOrgClusters->GetEntries();i++)
1102  {
1103  AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(i);
1104  vcluster->GetMomentum(*cluster_vec,fVertex);
1105 
1106  if (cluster_vec->Pt()>=fClusterMinPt && vcluster->IsEMCAL()==kTRUE)
1107  {
1108  fmyClusters->Add(vcluster);
1109  }
1110  }
1111  fnClusters = fmyClusters->GetEntries();
1112  delete cluster_vec;
1113 }
1114 
1116 {
1117  TrackCuts();
1118  if (fMCPartLevel == kTRUE)
1119  {
1120  return;
1121  }
1122 
1123  ClusterCuts();
1124  if (fnTracks>0)
1125  {
1126  if (fnClusters>0)
1127  {
1129  }
1130  else
1131  {
1133  }
1134  if (fnEMCalTracks>0)
1135  {
1136  if (fnClusters>0)
1137  {
1139  }
1140  else
1141  {
1143  }
1144  }
1145  else
1146  {
1147  if (fnClusters>0)
1148  {
1150  }
1151  else
1152  {
1154  }
1155  }
1156  }
1157  else
1158  {
1159  if (fnClusters>0)
1160  {
1162  }
1163  else
1164  {
1166  }
1167  }
1168 }
1169 
1171 {
1172  // Fill track histograms: Phi,Eta,Pt
1173  Int_t i,j;
1174  Int_t TCBins=100;
1175  TH2F *hdummypT= new TH2F("hdummypT","",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
1176 
1177  for (i=0;i<fnTracks;i++)
1178  {
1179  AliPicoTrack* vtrack = (AliPicoTrack*) fmyTracks->At(i);
1180  fhTrackPt->Fill(vtrack->Pt());
1181  fhTrackEta->Fill(vtrack->Eta());
1182  fhTrackPhi->Fill(vtrack->Phi());
1183  fhTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1184  fhTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1185  fhTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1186 
1187  // Fill Associated Track Distributions for AOD QA Productions
1188  // Global Tracks
1189  if (vtrack->GetTrackType()==0)
1190  {
1191  fhGlobalTrackPt->Fill(vtrack->Pt());
1192  fhGlobalTrackEta->Fill(vtrack->Eta());
1193  fhGlobalTrackPhi->Fill(vtrack->Phi());
1194  fhGlobalTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1195  fhGlobalTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1196  fhGlobalTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1197  }
1198  // Complementary Tracks
1199  else if (vtrack->GetTrackType()==1)
1200  {
1201  fhComplementaryTrackPt->Fill(vtrack->Pt());
1202  fhComplementaryTrackEta->Fill(vtrack->Eta());
1203  fhComplementaryTrackPhi->Fill(vtrack->Phi());
1204  fhComplementaryTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
1205  fhComplementaryTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
1206  fhComplementaryTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
1207  }
1208  hdummypT->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
1209  }
1210  for (i=1;i<=TCBins;i++)
1211  {
1212  for (j=1;j<=TCBins;j++)
1213  {
1214  fpTrackPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fTPCArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1215  }
1216  }
1217  delete hdummypT;
1218 }
1219 
1221 {
1222  // Fill cluster histograms: Phi,Eta,Pt
1223  Int_t i,j;
1224  Int_t TCBins=100;
1225  TH2F *hdummypT= new TH2F("hdummypT","",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
1226  Int_t myCellID=-2;
1227  TLorentzVector *cluster_vec = new TLorentzVector;
1228 
1229  for (i=0;i<fnClusters;i++)
1230  {
1231  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1232 
1233  if(!vcluster)continue;
1234  vcluster->GetMomentum(*cluster_vec,fVertex);
1235 
1236 //orig ELI fhClusterPt->Fill(cluster_vec->Pt());
1237  fhClusterPt->Fill(cluster_vec->E());
1238  fhClusterEta->Fill(cluster_vec->Eta());
1239  fhClusterPhi->Fill(cluster_vec->Phi());
1240  fhClusterEtaPhi->Fill(cluster_vec->Eta(),cluster_vec->Phi());
1241  fhClusterPhiPt->Fill(cluster_vec->Phi(),cluster_vec->Pt());
1242  fhClusterEtaPt->Fill(cluster_vec->Eta(),cluster_vec->Pt());
1243  hdummypT->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
1244  fEMCALGeometry->GetAbsCellIdFromEtaPhi(cluster_vec->Eta(),cluster_vec->Phi(),myCellID);
1245  fhEMCalCellCounts->Fill(myCellID);
1246  }
1247  for (i=1;i<=TCBins;i++)
1248  {
1249  for (j=1;j<=TCBins;j++)
1250  {
1251  fpClusterPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fEMCalArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
1252  }
1253  }
1254 
1255  delete hdummypT;
1256  delete cluster_vec;
1257 }
1258 
1260 {
1261  // Preliminary Jet Placement and Selection Cuts
1262  Int_t i;
1263 
1264  fnAKTChargedJets = fmyAKTChargedJets->GetEntries();
1265  fnKTChargedJets = fmyKTChargedJets->GetEntries();
1266 
1267  fTPCJet = new AlipAJetData("fTPCJet",kFALSE,fnAKTChargedJets);
1268  fTPCFullJet = new AlipAJetData("fTPCFullJet",kFALSE,fnAKTChargedJets);
1269  fTPCOnlyJet = new AlipAJetData("fTPCOnlyJet",kFALSE,fnAKTChargedJets);
1270  fTPCJetUnbiased = new AlipAJetData("fTPCJetUnbiased",kFALSE,fnAKTChargedJets);
1271 
1274  fTPCJet->SetJetR(fJetR);
1288 
1289  // Initialize Jet Data
1290  for (i=0;i<fnAKTChargedJets;i++)
1291  {
1292  AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(i);
1293 
1294  fTPCJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
1295  fTPCFullJet->SetIsJetInArray(IsInTPC(fJetRAccept,myJet->Phi(),myJet->Eta(),kTRUE),i);
1296  fTPCOnlyJet->SetIsJetInArray(IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta()),i);
1297  fTPCJetUnbiased->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
1298  }
1299 
1300  fTPCJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1301  fTPCFullJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1302  fTPCOnlyJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
1304 
1305  // kT Jets
1306  fTPCkTFullJet = new AlipAJetData("fTPCkTFullJet",kFALSE,fnKTChargedJets);
1310 
1311  for (i=0;i<fnKTChargedJets;i++)
1312  {
1313  AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(i);
1314  fTPCkTFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
1315  }
1317 
1318  // Raw Charged Jet Spectra
1319  if (fCalculateRhoJet>=2 || (fCalculateRhoJet>=1 && fMCPartLevel==kTRUE))
1320  {
1322  }
1323 }
1324 
1326 {
1327 // cout<<"start of ...FullpAJets_Eli_Mod::InitFullJets()"<<endl;
1328  // Preliminary Jet Placement and Selection Cuts
1329  Int_t i;
1330 
1331  if(fmyAKTFullJets)fnAKTFullJets = fmyAKTFullJets->GetEntries();
1332  else fnAKTFullJets = 0;
1333 
1334  if(fmyKTFullJets)fnKTFullJets = fmyKTFullJets->GetEntries();
1335  else fnKTFullJets=0;
1336 
1337  fEMCalFullJet = new AlipAJetData("fEMCalFullJet",kTRUE,fnAKTFullJets);
1338  //ELI thats the only one that should matter- whats the rest for??
1339  fEMCalFullJet->SetSignalCut(fEMCalJetThreshold); //ELI - fSignalPt = 5
1341  fEMCalFullJet->SetJetR(fJetR); //ELI - 0.2
1344 
1345  // Initialize Jet Data
1346  for (i=0;i<fnAKTFullJets;i++)
1347  {
1348  AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
1349  //fEMCalFullJet->SetIsJetInArray(IsInEMCalFull(fJetRAccept,myJet->Phi(),myJet->Eta()),i);
1350  /*No acceptance cuts*/fEMCalFullJet->SetIsJetInArray(1,i);
1351  }
1353 
1354  // Raw Full Jet Spectra
1355  if (fMCPartLevel==kFALSE)
1356  {
1357  //ELI fill some histograms
1359  //ELI FillBSJS(Double_t eventCentrality, Double_t rho, Double_t signalCut, TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList)
1360  }
1361 // cout<<"end of ...FullpAJets_Eli_Mod::InitFullJets()"<<endl;
1362 }
1363 
1365 {
1366  Int_t i,j;
1367  Double_t E_tracks_total=0.;
1368  TRandom3 u(time(NULL));
1369 
1370  Double_t Eta_Center=0.5*(fTPCEtaMin+fTPCEtaMax);
1371  Double_t Phi_Center=0.5*(fTPCPhiMin+fTPCPhiMax);
1372  Int_t event_mult=0;
1373  Int_t event_track_mult=0;
1374  Int_t clus_mult=0;
1375 
1376  for (i=0;i<fnBckgClusters;i++)
1377  {
1378  fTPCRCBckgFluc[i]=0.0;
1379  fTPCRCBckgFlucSignal[i]=0.0;
1380  fTPCRCBckgFlucNColl[i]=0.0;
1381  }
1382 
1383  TLorentzVector *dummy= new TLorentzVector;
1384  TLorentzVector *temp_jet= new TLorentzVector;
1385  TLorentzVector *track_vec = new TLorentzVector;
1386 
1387  // First, consider the RC with no spatial restrictions
1388  for (j=0;j<fnBckgClusters;j++)
1389  {
1390  E_tracks_total=0.;
1391 
1392  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1393  // Loop over all tracks
1394  for (i=0;i<fnTracks;i++)
1395  {
1396  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1397  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1398  {
1399  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1400  if (dummy->DeltaR(*track_vec)<fJetR)
1401  {
1402  E_tracks_total+=vtrack->Pt();
1403  }
1404  }
1405  }
1406  fTPCRCBckgFlucSignal[j]=E_tracks_total;
1407  }
1408 
1409  // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1410  E_tracks_total=0.0;
1411  if (fTPCJetUnbiased->GetLeadingPt()<0.0)
1412  {
1413  temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1414  }
1415  else
1416  {
1418  myJet->GetMomentum(*temp_jet);
1419  }
1420 
1421  for (j=0;j<fnBckgClusters;j++)
1422  {
1423  event_mult=0;
1424  event_track_mult=0;
1425  clus_mult=0;
1426  E_tracks_total=0.;
1427 
1428  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1429  while (temp_jet->DeltaR(*dummy)<fJetR)
1430  {
1431  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1432  }
1433  // Loop over all tracks
1434  for (i=0;i<fnTracks;i++)
1435  {
1436  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1437  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE) == kTRUE)
1438  {
1439  event_mult++;
1440  if (IsInEMCal(vtrack->Phi(),vtrack->Eta()) == kTRUE)
1441  {
1442  event_track_mult++;
1443  }
1444  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1445  if (dummy->DeltaR(*track_vec)<fJetR)
1446  {
1447  clus_mult++;
1448  E_tracks_total+=vtrack->Pt();
1449  }
1450  }
1451  }
1452  fTPCRCBckgFluc[j]=E_tracks_total;
1453  }
1454  if (fTrackQA==kTRUE)
1455  {
1456  fhTPCEventMult->Fill(fEventCentrality,event_mult);
1457  fpTPCEventMult->Fill(fEventCentrality,event_mult);
1458  fhEMCalTrackEventMult->Fill(fEventCentrality,event_track_mult);
1459  }
1460  if (fCalculateRhoJet>=2 || (fCalculateRhoJet>=1 && fMCPartLevel==kTRUE))
1461  {
1463  }
1464 
1465  // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
1466  Double_t exclusion_prob;
1467  for (j=0;j<fnBckgClusters;j++)
1468  {
1469  exclusion_prob = u.Uniform(0,1);
1470  if (exclusion_prob<(1/fNColl))
1471  {
1473  }
1474  else
1475  {
1477  }
1478  }
1479 
1480  delete dummy;
1481  delete temp_jet;
1482  delete track_vec;
1483 }
1484 
1486 {
1487  Int_t i,j;
1488  Double_t E_tracks_total=0.;
1489  Double_t E_caloclusters_total=0.;
1490  TRandom3 u(time(NULL));
1491 
1492  Double_t Eta_Center=0.5*(fEMCalEtaMin+fEMCalEtaMax);
1493  Double_t Phi_Center=0.5*(fEMCalPhiMin+fEMCalPhiMax);
1494  Int_t event_mult=0;
1495  Int_t clus_mult=0;
1496 
1497  for (i=0;i<fnBckgClusters;i++)
1498  {
1499  fEMCalRCBckgFluc[i]=0.0;
1500  fEMCalRCBckgFlucSignal[i]=0.0;
1501  fEMCalRCBckgFlucNColl[i]=0.0;
1502  }
1503 
1504  TLorentzVector *dummy= new TLorentzVector;
1505  TLorentzVector *temp_jet= new TLorentzVector;
1506  TLorentzVector *track_vec = new TLorentzVector;
1507  TLorentzVector *cluster_vec = new TLorentzVector;
1508 
1509  // First, consider the RC with no spatial restrictions
1510  for (j=0;j<fnBckgClusters;j++)
1511  {
1512  E_tracks_total=0.;
1513  E_caloclusters_total=0.;
1514 
1515  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1516  // Loop over all tracks
1517  for (i=0;i<fnTracks;i++)
1518  {
1519  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1520  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1521  {
1522  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1523  if (dummy->DeltaR(*track_vec)<fJetR)
1524  {
1525  E_tracks_total+=vtrack->Pt();
1526  }
1527  }
1528  }
1529 
1530  // Loop over all caloclusters
1531  for (i=0;i<fnClusters;i++)
1532  {
1533  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1534  vcluster->GetMomentum(*cluster_vec,fVertex);
1535  if (dummy->DeltaR(*cluster_vec)<fJetR)
1536  {
1537  clus_mult++;
1538  E_caloclusters_total+=cluster_vec->Pt();
1539  }
1540  }
1541  fEMCalRCBckgFlucSignal[j]=E_tracks_total+E_caloclusters_total;
1542  }
1543 
1544  // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
1545  E_tracks_total=0.;
1546  E_caloclusters_total=0.;
1548  {
1549  temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1550  }
1551  else
1552  {
1554  myJet->GetMomentum(*temp_jet);
1555  }
1556 
1557  for (j=0;j<fnBckgClusters;j++)
1558  {
1559  event_mult=0;
1560  clus_mult=0;
1561  E_tracks_total=0.;
1562  E_caloclusters_total=0.;
1563 
1564  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1565  while (temp_jet->DeltaR(*dummy)<fJetR)
1566  {
1567  dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
1568  }
1569  // Loop over all tracks
1570  for (i=0;i<fnTracks;i++)
1571  {
1572  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1573  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
1574  {
1575  event_mult++;
1576  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1577  if (dummy->DeltaR(*track_vec)<fJetR)
1578  {
1579  clus_mult++;
1580  E_tracks_total+=vtrack->Pt();
1581  }
1582  }
1583  }
1584 
1585  // Loop over all caloclusters
1586  for (i=0;i<fnClusters;i++)
1587  {
1588  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
1589  vcluster->GetMomentum(*cluster_vec,fVertex);
1590  event_mult++;
1591  if (dummy->DeltaR(*cluster_vec)<fJetR)
1592  {
1593  clus_mult++;
1594  E_caloclusters_total+=cluster_vec->Pt();
1595  }
1596  }
1597  fEMCalRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
1598  }
1599  if (fClusterQA==kTRUE)
1600  {
1601  fhEMCalEventMult->Fill(fEventCentrality,event_mult);
1602  fpEMCalEventMult->Fill(fEventCentrality,event_mult);
1603  }
1605 
1606  // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
1607  Double_t exclusion_prob;
1608  for (j=0;j<fnBckgClusters;j++)
1609  {
1610  exclusion_prob = u.Uniform(0,1);
1611  if (exclusion_prob<(1/fNColl))
1612  {
1614  }
1615  else
1616  {
1618  }
1619  }
1620 
1621  delete dummy;
1622  delete temp_jet;
1623  delete track_vec;
1624  delete cluster_vec;
1625 }
1626 
1627 // Charged Rho's
1629 {
1630  Int_t i;
1631  Double_t E_tracks_total=0.0;
1632  Double_t TPC_rho=0.0;
1633 
1634  // Loop over all tracks
1635  for (i=0;i<fnTracks;i++)
1636  {
1637  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1638  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1639  {
1640  E_tracks_total+=vtrack->Pt();
1641  }
1642  }
1643 
1644  // Calculate the mean Background density
1645  TPC_rho=E_tracks_total/fTPCArea;
1646 
1647  fRhoCharged = TPC_rho;
1648 
1649  // Fill Histograms
1657 
1658 }
1659 
1661 {
1662  Int_t i;
1663  Double_t E_tracks_total=0.0;
1664  Double_t TPC_rho=0.;
1665 
1666  TLorentzVector *temp_jet= new TLorentzVector;
1667  TLorentzVector *track_vec = new TLorentzVector;
1668 
1670  {
1672  myJet->GetMomentum(*temp_jet);
1673 
1674  // Loop over all tracks
1675  for (i=0;i<fnTracks;i++)
1676  {
1677  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1678  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1679  {
1680  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1681  if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
1682  {
1683  E_tracks_total+=vtrack->Pt();
1684  }
1685  }
1686  }
1687 
1688  // Calculate the mean Background density
1689  TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1690  }
1691  else // i.e. No signal jets -> same as total background density
1692  {
1693  // Loop over all tracks
1694  for (i=0;i<fnTracks;i++)
1695  {
1696  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1697  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1698  {
1699  E_tracks_total+=vtrack->Pt();
1700  }
1701  }
1702  // Calculate the mean Background density
1703  TPC_rho=E_tracks_total/fTPCArea;
1704  }
1705  delete track_vec;
1706  delete temp_jet;
1707 
1708  // Fill histograms
1716 }
1717 
1719 {
1720  Int_t i;
1721  Double_t E_tracks_total=0.0;
1722  Double_t TPC_rho=0.;
1723 
1724  TLorentzVector *temp_jet1= new TLorentzVector;
1725  TLorentzVector *temp_jet2= new TLorentzVector;
1726  TLorentzVector *track_vec = new TLorentzVector;
1727 
1729  {
1731  myhJet->GetMomentum(*temp_jet1);
1732 
1734  myJet->GetMomentum(*temp_jet2);
1735 
1736  // Loop over all tracks
1737  for (i=0;i<fnTracks;i++)
1738  {
1739  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1740  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1741  {
1742  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1743  if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
1744  {
1745  E_tracks_total+=vtrack->Pt();
1746  }
1747  }
1748  }
1749 
1750  // Calculate the mean Background density
1751  TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myhJet->Eta())-AreaWithinTPC(fJetR,myJet->Eta()));
1752  }
1754  {
1756  myJet->GetMomentum(*temp_jet1);
1757 
1758  // Loop over all tracks
1759  for (i=0;i<fnTracks;i++)
1760  {
1761  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1762  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1763  {
1764  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1765  if (temp_jet1->DeltaR(*track_vec)>fJetRForRho)
1766  {
1767  E_tracks_total+=vtrack->Pt();
1768  }
1769  }
1770  }
1771 
1772  // Calculate the mean Background density
1773  TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
1774  }
1775  else // i.e. No signal jets -> same as total background density
1776  {
1777  // Loop over all tracks
1778  for (i=0;i<fnTracks;i++)
1779  {
1780  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
1781  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1782  {
1783  E_tracks_total+=vtrack->Pt();
1784  }
1785  }
1786 
1787  // Calculate the mean Background density
1788  TPC_rho=E_tracks_total/fTPCArea;
1789  }
1790  delete temp_jet1;
1791  delete temp_jet2;
1792  delete track_vec;
1793 
1794  // Fill histograms
1802 }
1803 
1805 {
1806  Int_t i,j;
1807  Bool_t track_away_from_jet;
1808  Double_t E_tracks_total=0.0;
1809  Double_t TPC_rho=0.0;
1810  Double_t jet_area_total=0.0;
1811 
1812  TLorentzVector *jet_vec= new TLorentzVector;
1813  TLorentzVector *track_vec = new TLorentzVector;
1814 
1815  // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1816  for (i=0;i<fnTracks;i++)
1817  {
1818  // First, check if track is in the EMCal!!
1819  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1820  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1821  {
1823  {
1824  E_tracks_total+=vtrack->Pt();
1825  }
1826  else
1827  {
1828  track_away_from_jet=kTRUE;
1829  j=0;
1830  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1831  while (track_away_from_jet==kTRUE && j<fTPCJetUnbiased->GetTotalSignalJets())
1832  {
1834  myJet->GetMomentum(*jet_vec);
1835  if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
1836  {
1837  track_away_from_jet=kFALSE;
1838  }
1839  j++;
1840  }
1841  if (track_away_from_jet==kTRUE)
1842  {
1843  E_tracks_total+=vtrack->Pt();
1844  }
1845  }
1846  }
1847  }
1848 
1849  // Determine area of all Jets that are within the EMCal
1851  {
1852  jet_area_total=0.0;
1853  }
1854  else
1855  {
1856  for (i=0;i<fTPCJetUnbiased->GetTotalSignalJets();i++)
1857  {
1859  jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
1860  }
1861  }
1862  delete jet_vec;
1863  delete track_vec;
1864 
1865  // Calculate Rho
1866  TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
1867 
1868  // Fill Histogram
1876 
1877 }
1878 
1880 {
1881  Int_t i,j;
1882  Bool_t track_away_from_jet;
1883  Double_t E_tracks_total=0.0;
1884  Double_t TPC_rho=0.0;
1885  Double_t jet_area_total=0.0;
1886 
1887  TLorentzVector *jet_vec= new TLorentzVector;
1888  TLorentzVector *track_vec = new TLorentzVector;
1889 
1890  // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
1891  for (i=0;i<fnTracks;i++)
1892  {
1893  // First, check if track is in the EMCal!!
1894  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
1895  if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
1896  {
1898  {
1899  E_tracks_total+=vtrack->Pt();
1900  }
1901  else
1902  {
1903  track_away_from_jet=kTRUE;
1904  j=0;
1905  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
1906  while (track_away_from_jet==kTRUE && j<fTPCJetUnbiased->GetTotalSignalJets())
1907  {
1909  myJet->GetMomentum(*jet_vec);
1910  if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
1911  {
1912  track_away_from_jet=kFALSE;
1913  }
1914  j++;
1915  }
1916  if (track_away_from_jet==kTRUE)
1917  {
1918  E_tracks_total+=vtrack->Pt();
1919  }
1920  }
1921  }
1922  }
1923 
1924  // Determine area of all Jets that are within the TPC
1926  {
1927  jet_area_total=0.0;
1928  }
1929  else
1930  {
1931  for (i=0;i<fTPCJetUnbiased->GetTotalSignalJets();i++)
1932  {
1934  jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
1935  }
1936  }
1937  delete jet_vec;
1938  delete track_vec;
1939 
1940  // Calculate Rho
1941  TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
1942  TPC_rho*=fScaleFactor;
1943 
1944  // Fill Histogram
1953 }
1954 
1956 {
1957  Int_t i;
1958  Double_t kTRho = 0.0;
1959  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1960  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1961 
1962  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1963  {
1965  pTArray[i]=myJet->Pt();
1966  RhoArray[i]=myJet->Pt()/myJet->Area();
1967  }
1968 
1969  if (fTPCkTFullJet->GetTotalJets()>=2)
1970  {
1971  kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
1972 
1980  }
1981  delete [] RhoArray;
1982  delete [] pTArray;
1983 }
1984 
1986 {
1987  Int_t i;
1988  Double_t kTRho = 0.0;
1989  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1990  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
1991 
1992  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
1993  {
1995  pTArray[i]=myJet->Pt();
1996  RhoArray[i]=myJet->Pt()/myJet->Area();
1997  }
1998 
1999  if (fTPCkTFullJet->GetTotalJets()>=2)
2000  {
2001  kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
2002  kTRho*=fScaleFactor;
2003 
2011  }
2012  delete [] RhoArray;
2013  delete [] pTArray;
2014 }
2015 
2017 {
2018  Int_t i,k;
2019  Double_t kTRho = 0.0;
2020  Double_t CMSTotalkTArea = 0.0;
2021  Double_t CMSTrackArea = 0.0;
2022  Double_t CMSCorrectionFactor = 1.0;
2023  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2024  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2025 
2026  k=0;
2028  {
2031 
2032  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2033  {
2035 
2036  CMSTotalkTArea+=myJet->Area();
2037  if (myJet->GetNumberOfTracks()>0)
2038  {
2039  CMSTrackArea+=myJet->Area();
2040  }
2041  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
2042  {
2043  pTArray[k]=myJet->Pt();
2044  RhoArray[k]=myJet->Pt()/myJet->Area();
2045  k++;
2046  }
2047  }
2048  if (k>0)
2049  {
2050  kTRho=MedianRhokT(pTArray,RhoArray,k);
2051  }
2052  else
2053  {
2054  kTRho=0.0;
2055  }
2056  }
2057  else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
2058  {
2060 
2061  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2062  {
2064 
2065  CMSTotalkTArea+=myJet->Area();
2066  if (myJet->GetNumberOfTracks()>0)
2067  {
2068  CMSTrackArea+=myJet->Area();
2069  }
2070  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
2071  {
2072  pTArray[k]=myJet->Pt();
2073  RhoArray[k]=myJet->Pt()/myJet->Area();
2074  k++;
2075  }
2076  }
2077  if (k>0)
2078  {
2079  kTRho=MedianRhokT(pTArray,RhoArray,k);
2080  }
2081  else
2082  {
2083  kTRho=0.0;
2084  }
2085  }
2086  else
2087  {
2088  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2089  {
2091 
2092  CMSTotalkTArea+=myJet->Area();
2093  if (myJet->GetNumberOfTracks()>0)
2094  {
2095  CMSTrackArea+=myJet->Area();
2096  }
2097  pTArray[k]=myJet->Pt();
2098  RhoArray[k]=myJet->Pt()/myJet->Area();
2099  k++;
2100  }
2101  if (k>0)
2102  {
2103  kTRho=MedianRhokT(pTArray,RhoArray,k);
2104  }
2105  else
2106  {
2107  kTRho=0.0;
2108  }
2109  }
2110  // Scale CMS Rho by Correction factor
2111  if (CMSTotalkTArea==0.0)
2112  {
2113  CMSCorrectionFactor = 1.0;
2114  }
2115  else
2116  {
2117  //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2118  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
2119  }
2120  kTRho*=CMSCorrectionFactor;
2128  delete [] RhoArray;
2129  delete [] pTArray;
2130 }
2131 
2133 {
2134  Int_t i,k;
2135  Double_t kTRho = 0.0;
2136  Double_t CMSTotalkTArea = 0.0;
2137  Double_t CMSTrackArea = 0.0;
2138  Double_t CMSCorrectionFactor = 1.0;
2139  Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2140  Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
2141 
2142  k=0;
2144  {
2147 
2148  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2149  {
2151 
2152  CMSTotalkTArea+=myJet->Area();
2153  if (myJet->GetNumberOfTracks()>0)
2154  {
2155  CMSTrackArea+=myJet->Area();
2156  }
2157  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
2158  {
2159  pTArray[k]=myJet->Pt();
2160  RhoArray[k]=myJet->Pt()/myJet->Area();
2161  k++;
2162  }
2163  }
2164  if (k>0)
2165  {
2166  kTRho=MedianRhokT(pTArray,RhoArray,k);
2167  }
2168  else
2169  {
2170  kTRho=0.0;
2171  }
2172  }
2173  else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
2174  {
2176 
2177  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2178  {
2180 
2181  CMSTotalkTArea+=myJet->Area();
2182  if (myJet->GetNumberOfTracks()>0)
2183  {
2184  CMSTrackArea+=myJet->Area();
2185  }
2186  if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
2187  {
2188  pTArray[k]=myJet->Pt();
2189  RhoArray[k]=myJet->Pt()/myJet->Area();
2190  k++;
2191  }
2192  }
2193  if (k>0)
2194  {
2195  kTRho=MedianRhokT(pTArray,RhoArray,k);
2196  }
2197  else
2198  {
2199  kTRho=0.0;
2200  }
2201  }
2202  else
2203  {
2204  for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
2205  {
2207 
2208  CMSTotalkTArea+=myJet->Area();
2209  if (myJet->GetNumberOfTracks()>0)
2210  {
2211  CMSTrackArea+=myJet->Area();
2212  }
2213  pTArray[k]=myJet->Pt();
2214  RhoArray[k]=myJet->Pt()/myJet->Area();
2215  k++;
2216  }
2217  if (k>0)
2218  {
2219  kTRho=MedianRhokT(pTArray,RhoArray,k);
2220  }
2221  else
2222  {
2223  kTRho=0.0;
2224  }
2225  }
2226  kTRho*=fScaleFactor;
2227  // Scale CMS Rho by Correction factor
2228  if (CMSTotalkTArea==0.0)
2229  {
2230  CMSCorrectionFactor = 1.0;
2231  }
2232  else
2233  {
2234  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
2235  }
2236  kTRho*=CMSCorrectionFactor;
2237 
2248 
2249  delete [] RhoArray;
2250  delete [] pTArray;
2251 }
2252 
2253 // Full Rho's
2255 {
2256  Int_t i;
2257  Double_t E_tracks_total=0.0;
2258  Double_t E_caloclusters_total=0.0;
2259  Double_t EMCal_rho=0.0;
2260 
2261  TLorentzVector *cluster_vec = new TLorentzVector;
2262 
2263  // Loop over all tracks
2264  for (i=0;i<fnTracks;i++)
2265  {
2266  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2267  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2268  {
2269  E_tracks_total+=vtrack->Pt();
2270  }
2271  }
2272 
2273  // Loop over all caloclusters
2274  for (i=0;i<fnClusters;i++)
2275  {
2276  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2277  vcluster->GetMomentum(*cluster_vec,fVertex);
2278  E_caloclusters_total+=cluster_vec->Pt();
2279  }
2280  delete cluster_vec;
2281 
2282  // Calculate the mean Background density
2283  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2284  fRhoFull=EMCal_rho;
2285 
2286  // Fill Histograms
2287  fRhoFull0->FillRho(fEventCentrality,EMCal_rho);
2294 }
2295 
2297 {
2298  Int_t i;
2299  Double_t E_tracks_total=0.0;
2300  Double_t E_caloclusters_total=0.0;
2301  Double_t EMCal_rho=0.0;
2302 
2303  TLorentzVector *temp_jet= new TLorentzVector;
2304  TLorentzVector *track_vec = new TLorentzVector;
2305  TLorentzVector *cluster_vec = new TLorentzVector;
2306 
2308  {
2310  myJet->GetMomentum(*temp_jet);
2311 
2312  // Loop over all tracks
2313  for (i=0;i<fnTracks;i++)
2314  {
2315  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2316  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2317  {
2318  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2319  if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
2320  {
2321  E_tracks_total+=vtrack->Pt();
2322  }
2323  }
2324  }
2325 
2326  // Loop over all caloclusters
2327  for (i=0;i<fnClusters;i++)
2328  {
2329  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2330  vcluster->GetMomentum(*cluster_vec,fVertex);
2331  if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
2332  {
2333  E_caloclusters_total+=cluster_vec->Pt();
2334  }
2335  }
2336  // Calculate the mean Background density
2337  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2338  }
2339  else // i.e. No signal jets -> same as total background density
2340  {
2341  // Loop over all tracks
2342  for (i=0;i<fnTracks;i++)
2343  {
2344  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2345  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2346  {
2347  E_tracks_total+=vtrack->Pt();
2348  }
2349  }
2350 
2351  // Loop over all caloclusters
2352  for (i=0;i<fnClusters;i++)
2353  {
2354  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2355  vcluster->GetMomentum(*cluster_vec,fVertex);
2356  E_caloclusters_total+=cluster_vec->Pt();
2357  }
2358  // Calculate the mean Background density
2359  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2360  }
2361  delete temp_jet;
2362  delete track_vec;
2363  delete cluster_vec;
2364 
2365  // Fill histograms
2366  fRhoFull1->FillRho(fEventCentrality,EMCal_rho);
2373 }
2374 
2376 {
2377  Int_t i;
2378  Double_t E_tracks_total=0.0;
2379  Double_t E_caloclusters_total=0.0;
2380  Double_t EMCal_rho=0.0;
2381 
2382  TLorentzVector *temp_jet1 = new TLorentzVector;
2383  TLorentzVector *temp_jet2 = new TLorentzVector;
2384  TLorentzVector *track_vec = new TLorentzVector;
2385  TLorentzVector *cluster_vec = new TLorentzVector;
2386 
2388  {
2390  myhJet->GetMomentum(*temp_jet1);
2391 
2393  myJet->GetMomentum(*temp_jet2);
2394 
2395  // Loop over all tracks
2396  for (i=0;i<fnTracks;i++)
2397  {
2398  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2399  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2400  {
2401  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2402  if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
2403  {
2404  E_tracks_total+=vtrack->Pt();
2405  }
2406  }
2407  }
2408 
2409  // Loop over all caloclusters
2410  for (i=0;i<fnClusters;i++)
2411  {
2412  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2413  vcluster->GetMomentum(*cluster_vec,fVertex);
2414  if ((temp_jet1->DeltaR(*cluster_vec)>fJetRForRho) && (temp_jet2->DeltaR(*cluster_vec)>fJetRForRho))
2415  {
2416  E_caloclusters_total+=cluster_vec->Pt();
2417  }
2418  }
2419 
2420  // Calculate the mean Background density
2421  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myhJet->Phi(),myhJet->Eta())-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2422  }
2424  {
2426  myJet->GetMomentum(*temp_jet1);
2427 
2428  // Loop over all tracks
2429  for (i=0;i<fnTracks;i++)
2430  {
2431  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2432  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2433  {
2434  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2435  if (temp_jet1->DeltaR(*track_vec)>fJetRForRho)
2436  {
2437  E_tracks_total+=vtrack->Pt();
2438  }
2439  }
2440  }
2441 
2442  // Loop over all caloclusters
2443  for (i=0;i<fnClusters;i++)
2444  {
2445  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2446  vcluster->GetMomentum(*cluster_vec,fVertex);
2447  if (temp_jet1->DeltaR(*cluster_vec)>fJetRForRho)
2448  {
2449  E_caloclusters_total+=cluster_vec->Pt();
2450  }
2451  }
2452  // Calculate the mean Background density
2453  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
2454  }
2455  else // i.e. No signal jets -> same as total background density
2456  {
2457  // Loop over all tracks
2458  for (i=0;i<fnTracks;i++)
2459  {
2460  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2461  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2462  {
2463  E_tracks_total+=vtrack->Pt();
2464  }
2465  }
2466 
2467  // Loop over all caloclusters
2468  for (i=0;i<fnClusters;i++)
2469  {
2470  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2471  vcluster->GetMomentum(*cluster_vec,fVertex);
2472  E_caloclusters_total+=cluster_vec->Pt();
2473  }
2474  // Calculate the mean Background density
2475  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2476  }
2477  delete temp_jet1;
2478  delete temp_jet2;
2479  delete track_vec;
2480  delete cluster_vec;
2481 
2482  // Fill histograms
2483  fRhoFull2->FillRho(fEventCentrality,EMCal_rho);
2490 }
2491 
2493 {
2494  Int_t i,j;
2495  Bool_t track_away_from_jet;
2496  Bool_t cluster_away_from_jet;
2497  Double_t E_tracks_total=0.0;
2498  Double_t E_caloclusters_total=0.0;
2499  Double_t EMCal_rho=0.0;
2500  Double_t jet_area_total=0.0;
2501 
2502  TLorentzVector *jet_vec= new TLorentzVector;
2503  TLorentzVector *track_vec = new TLorentzVector;
2504  TLorentzVector *cluster_vec = new TLorentzVector;
2505 
2506  // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
2507  for (i=0;i<fnTracks;i++)
2508  {
2509  // First, check if track is in the EMCal!!
2510  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
2511  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2512  {
2514  {
2515  E_tracks_total+=vtrack->Pt();
2516  }
2517  else
2518  {
2519  track_away_from_jet=kTRUE;
2520  j=0;
2521  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2522  while (track_away_from_jet==kTRUE && j<fEMCalPartJetUnbiased->GetTotalSignalJets())
2523  {
2525  myJet->GetMomentum(*jet_vec);
2526  if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
2527  {
2528  track_away_from_jet=kFALSE;
2529  }
2530  j++;
2531  }
2532  if (track_away_from_jet==kTRUE)
2533  {
2534  E_tracks_total+=vtrack->Pt();
2535  }
2536  }
2537  }
2538  }
2539 
2540  // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
2541  for (i=0;i<fnClusters;i++)
2542  {
2543  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2544  vcluster->GetMomentum(*cluster_vec,fVertex);
2546  {
2547  E_caloclusters_total+=cluster_vec->Pt();
2548  }
2549  else
2550  {
2551  cluster_away_from_jet=kTRUE;
2552  j=0;
2553 
2554  while (cluster_away_from_jet==kTRUE && j<fEMCalPartJetUnbiased->GetTotalSignalJets())
2555  {
2557  myJet->GetMomentum(*jet_vec);
2558  if (cluster_vec->DeltaR(*jet_vec)<=fJetRForRho)
2559  {
2560  cluster_away_from_jet=kFALSE;
2561  }
2562  j++;
2563  }
2564  if (cluster_away_from_jet==kTRUE)
2565  {
2566  E_caloclusters_total+=cluster_vec->Pt();
2567  }
2568  }
2569  }
2570 
2571  // Determine area of all Jets that are within the EMCal
2573  {
2574  jet_area_total=0.0;
2575  }
2576  else
2577  {
2578  for (i=0;i<fEMCalPartJetUnbiased->GetTotalSignalJets();i++)
2579  {
2581  jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
2582  }
2583  }
2584  delete jet_vec;
2585  delete track_vec;
2586  delete cluster_vec;
2587 
2588  // Calculate Rho
2589  EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
2590 
2591  // Fill Histogram
2592  fRhoFullN->FillRho(fEventCentrality,EMCal_rho);
2599 }
2600 
2602 {
2603  Int_t i;
2604  Double_t E_tracks_total=0.0;
2605  Double_t E_caloclusters_total=0.0;
2606  Double_t EMCal_rho=0.0;
2607 
2608  TLorentzVector *cluster_vec = new TLorentzVector;
2609 
2610  // Loop over all tracks
2611  for (i=0;i<fnTracks;i++)
2612  {
2613  AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
2614  if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
2615  {
2616  E_tracks_total+=vtrack->Pt();
2617  }
2618  }
2619 
2620  // Loop over all caloclusters
2621  for (i=0;i<fnClusters;i++)
2622  {
2623  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
2624  vcluster->GetMomentum(*cluster_vec,fVertex);
2625  E_caloclusters_total+=cluster_vec->Pt();
2626  }
2627 
2628  delete cluster_vec;
2629 
2630  // Calculate the mean Background density
2631  EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
2632 
2633  // Fill Histograms
2641 }
2642 
2644 {
2645  Int_t i;
2646  Double_t kTRho = 0.0;
2647  Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2648  Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2649 
2650  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2651  {
2653  pTArray[i]=myJet->Pt();
2654  RhoArray[i]=myJet->Pt()/myJet->Area();
2655  }
2656 
2657  if (fEMCalkTFullJet->GetTotalJets()>0)
2658  {
2659  kTRho=MedianRhokT(pTArray,RhoArray,fEMCalkTFullJet->GetTotalJets());
2660  }
2661  else
2662  {
2663  kTRho=0.0;
2664  }
2672  delete [] RhoArray;
2673  delete [] pTArray;
2674 }
2675 
2677 {
2678  Int_t i,k;
2679  Double_t kTRho = 0.0;
2680  Double_t CMSTotalkTArea = 0.0;
2681  Double_t CMSParticleArea = 0.0;
2682  Double_t CMSCorrectionFactor = 1.0;
2683  Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2684  Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
2685 
2686  k=0;
2688  {
2691 
2692  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2693  {
2695 
2696  CMSTotalkTArea+=myJet->Area();
2697  if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2698  {
2699  CMSParticleArea+=myJet->Area();
2700  }
2701  if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kTRUE)
2702  {
2703  pTArray[k]=myJet->Pt();
2704  RhoArray[k]=myJet->Pt()/myJet->Area();
2705  k++;
2706  }
2707  }
2708  if (k>0)
2709  {
2710  kTRho=MedianRhokT(pTArray,RhoArray,k);
2711  }
2712  else
2713  {
2714  kTRho=0.0;
2715  }
2716  }
2718  {
2720 
2721  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2722  {
2724 
2725  CMSTotalkTArea+=myJet->Area();
2726  if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2727  {
2728  CMSParticleArea+=myJet->Area();
2729  }
2730  if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE)
2731  {
2732  pTArray[k]=myJet->Pt();
2733  RhoArray[k]=myJet->Pt()/myJet->Area();
2734  k++;
2735  }
2736  }
2737  if (k>0)
2738  {
2739  kTRho=MedianRhokT(pTArray,RhoArray,k);
2740  }
2741  else
2742  {
2743  kTRho=0.0;
2744  }
2745  }
2746  else
2747  {
2748  for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
2749  {
2751 
2752  CMSTotalkTArea+=myJet->Area();
2753  if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
2754  {
2755  CMSParticleArea+=myJet->Area();
2756  }
2757  pTArray[k]=myJet->Pt();
2758  RhoArray[k]=myJet->Pt()/myJet->Area();
2759  k++;
2760  }
2761  if (k>0)
2762  {
2763  kTRho=MedianRhokT(pTArray,RhoArray,k);
2764  }
2765  else
2766  {
2767  kTRho=0.0;
2768  }
2769  }
2770  // Scale CMS Rho by Correction factor
2771  if (CMSTotalkTArea==0.0)
2772  {
2773  CMSCorrectionFactor = 1.0;
2774  }
2775  else
2776  {
2777  //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
2778  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
2779  }
2780  kTRho*=CMSCorrectionFactor;
2781 
2789  delete [] RhoArray;
2790  delete [] pTArray;
2791 }
2792 
2794 {
2795  Int_t i,j;
2796  Double_t delta_R;
2797  Double_t dR=0.0;
2798  Double_t fullEDR = fFullEDJetR*100;
2799  const Int_t fullEDBins = (Int_t) fullEDR;
2800  Double_t ED_pT[fullEDBins];
2801  dR = fFullEDJetR/Double_t(fullEDBins); // Should be 0.01 be default
2802 
2803  TLorentzVector *jet_vec= new TLorentzVector;
2804  TLorentzVector *track_vec = new TLorentzVector;
2805  TLorentzVector *cluster_vec = new TLorentzVector;
2806 
2807  for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
2808  {
2810 
2811  if (IsInEMCalFull(fFullEDJetR,myJet->Phi(),myJet->Eta())==kTRUE)
2812  {
2813  for (j=0;j<fullEDBins;j++)
2814  {
2815  ED_pT[j]=0;
2816  }
2817  myJet->GetMomentum(*jet_vec);
2818 
2819  // Sum all tracks in concentric rings around jet vertex
2820  for (j=0;j<fnTracks;j++)
2821  {
2822  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2823  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2824  delta_R=jet_vec->DeltaR(*track_vec);
2825  if (delta_R<fFullEDJetR)
2826  {
2827  ED_pT[TMath::FloorNint((delta_R/fFullEDJetR)*fullEDBins)]+=vtrack->Pt();
2828  }
2829  }
2830 
2831  // Sum all clusters in concentric rings around jet vertex
2832  for (j=0;j<fnClusters;j++)
2833  {
2834  AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
2835  vcluster->GetMomentum(*cluster_vec,fVertex);
2836  delta_R=jet_vec->DeltaR(*cluster_vec);
2837  if (delta_R<fFullEDJetR)
2838  {
2839  ED_pT[TMath::FloorNint((delta_R/fFullEDJetR)*fullEDBins)]+=cluster_vec->Pt();
2840  }
2841  }
2842 
2843  for (j=0;j<fullEDBins;j++)
2844  {
2845  ED_pT[j] /= TMath::Pi()*dR*dR*(2*j+1);
2846  fpFullJetEDProfile->Fill(myJet->Pt(),fEventCentrality,j*dR,ED_pT[j]);
2847  }
2848  }
2849  }
2850  delete cluster_vec;
2851  delete track_vec;
2852  delete jet_vec;
2853 }
2854 
2855 // 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.
2857 {
2858  Int_t i,j;
2859  Double_t delta_R;
2860  Double_t dR=0.0;
2861  Double_t chargedEDR = fChargedEDJetR*100;
2862  const Int_t chargedEDBins = (Int_t) chargedEDR;
2863  Double_t ED_pT[chargedEDBins];
2864  dR = fChargedEDJetR/Double_t(chargedEDBins); // Should be 0.01 be default
2865 
2866  TLorentzVector *jet_vec= new TLorentzVector;
2867  TLorentzVector *track_vec = new TLorentzVector;
2868 
2869  for (i=0;i<fTPCFullJet->GetTotalSignalJets();i++)
2870  {
2872  if (IsInTPC(fChargedEDJetR,myJet->Phi(),myJet->Eta(),kTRUE)==kTRUE)
2873  {
2874  for (j=0;j<chargedEDBins;j++)
2875  {
2876  ED_pT[j]=0;
2877  }
2878  myJet->GetMomentum(*jet_vec);
2879 
2880  // Sum all tracks in concentric rings around jet vertex
2881  for (j=0;j<fnTracks;j++)
2882  {
2883  AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
2884  track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
2885  delta_R=jet_vec->DeltaR(*track_vec);
2886  if (delta_R<fChargedEDJetR)
2887  {
2888  ED_pT[TMath::FloorNint((delta_R/fChargedEDJetR)*chargedEDBins)]+=vtrack->Pt();
2889  }
2890  }
2891  for (j=0;j<chargedEDBins;j++)
2892  {
2893  ED_pT[j] /= TMath::Pi()*dR*dR*(2*j+1);
2894  fpChargedJetEDProfile->Fill(myJet->Pt(),fEventCentrality,j*dR,ED_pT[j]);
2895  fpChargedJetEDProfileScaled->Fill(myJet->Pt(),fEventCentrality,j*dR,fScaleFactor*ED_pT[j]);
2896  }
2897  }
2898  }
2899  delete track_vec;
2900  delete jet_vec;
2901 }
2902 
2904 {
2905  if (delOption ==0)
2906  {
2907  delete fRecoUtil;
2908 
2909  fRecoUtil = NULL;
2910  }
2911  else if (delOption==1)
2912  {
2913  delete fRecoUtil;
2914  delete fmyTracks;
2915  delete fmyClusters;
2916 
2917  fRecoUtil = NULL;
2918  fmyTracks = NULL;
2919  fmyClusters = NULL;
2920  }
2921  else if (delOption==2)
2922  {
2923  delete fRecoUtil;
2924  delete fmyTracks;
2925  delete fmyClusters;
2926  delete fTPCJet;
2927  delete fTPCFullJet;
2928  delete fTPCOnlyJet;
2929  delete fTPCJetUnbiased;
2930  delete fTPCkTFullJet;
2931 
2932  fRecoUtil = NULL;
2933  fmyTracks = NULL;
2934  fmyClusters = NULL;
2935  fTPCJet = NULL;
2936  fTPCFullJet = NULL;
2937  fTPCOnlyJet = NULL;
2938  fTPCJetUnbiased = NULL;
2939  fTPCkTFullJet = NULL;
2940  }
2941  if (delOption==3)
2942  {
2943  delete fRecoUtil;
2944  delete fmyTracks;
2945  delete fmyClusters;
2946  delete fTPCJet;
2947  delete fTPCFullJet;
2948  delete fTPCOnlyJet;
2949  delete fTPCJetUnbiased;
2950  delete fTPCkTFullJet;
2951  delete fEMCalJet;
2952  delete fEMCalFullJet;
2953  delete fEMCalPartJet;
2954  delete fEMCalPartJetUnbiased;
2955  delete fEMCalkTFullJet;
2956 
2957  fRecoUtil = NULL;
2958  fmyTracks = NULL;
2959  fmyClusters = NULL;
2960  fTPCJet = NULL;
2961  fTPCFullJet = NULL;
2962  fTPCOnlyJet = NULL;
2963  fTPCJetUnbiased = NULL;
2964  fEMCalJet = NULL;
2965  fEMCalFullJet = NULL;
2966  fEMCalPartJet = NULL;
2967  fEMCalPartJetUnbiased = NULL;
2968  fEMCalkTFullJet = NULL;
2969  }
2970 }
2971 
2975 
2977 {
2978  // Determine if event contains a di-jet within the detector. Uses charged jets.
2979  // Requires the delta phi of the jets to be 180 +/- 15 degrees.
2980  // Requires both jets to be outside of the EMCal
2981  // Requires both jets to be signal jets
2982 
2983  const Double_t dijet_delta_phi=(180/360.)*2*TMath::Pi();
2984  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
2985  Double_t dummy_phi=0.0;
2986 
2988  {
2991  dummy_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myJet->Phi()));
2992  if (dummy_phi>(dijet_delta_phi-dijet_phi_acceptance))
2993  {
2994  fnDiJetEvents++;
2995  return kTRUE;
2996  }
2997  }
2998  return kFALSE;
2999 }
3000 
3001 Bool_t AliAnalysisTaskFullpAJets_Eli_Mod::InsideRect(Double_t phi,Double_t phi_min,Double_t phi_max,Double_t eta,Double_t eta_min,Double_t eta_max)
3002 {
3003  if (phi>phi_min && phi<phi_max)
3004  {
3005  if (eta>eta_min && eta<eta_max)
3006  {
3007  return kTRUE;
3008  }
3009  }
3010  return kFALSE;
3011 }
3012 
3013 Bool_t AliAnalysisTaskFullpAJets_Eli_Mod::IsInEMCal(Double_t phi,Double_t eta)
3014 {
3016 }
3017 
3018 Bool_t AliAnalysisTaskFullpAJets_Eli_Mod::IsInEMCalFull(Double_t r,Double_t phi,Double_t eta)
3019 {
3021 }
3022 
3023 Bool_t AliAnalysisTaskFullpAJets_Eli_Mod::IsInEMCalPart(Double_t r,Double_t phi,Double_t eta)
3024 {
3026 }
3027 
3028 Bool_t AliAnalysisTaskFullpAJets_Eli_Mod::IsInTPCFull(Double_t r,Double_t phi,Double_t eta)
3029 {
3030  Bool_t in_EMCal= InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
3031  Bool_t in_TPC= InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
3032 
3033  if (in_EMCal==kFALSE && in_TPC==kTRUE)
3034  {
3035  return kTRUE;
3036  }
3037  return kFALSE;
3038 }
3039 
3040 Bool_t AliAnalysisTaskFullpAJets_Eli_Mod::IsInTPC(Double_t r,Double_t phi,Double_t eta,Bool_t Complete)
3041 {
3042  if (Complete==kTRUE)
3043  {
3045  }
3047 }
3048 
3050 {
3051  Int_t i,j;
3052  Int_t jetTrack1=0;
3053  Int_t jetTrack2=0;
3054  Int_t jetCluster1=0;
3055  Int_t jetCluster2=0;
3056 
3057  for (i=0;i<jet1->GetNumberOfTracks();i++)
3058  {
3059  jetTrack1=jet1->TrackAt(i);
3060  for (j=0;j<jet2->GetNumberOfTracks();j++)
3061  {
3062  jetTrack2=jet2->TrackAt(j);
3063  if (jetTrack1 == jetTrack2)
3064  {
3065  return kTRUE;
3066  }
3067  }
3068  }
3069  if (EMCalOn == kTRUE)
3070  {
3071  for (i=0;i<jet1->GetNumberOfClusters();i++)
3072  {
3073  jetCluster1=jet1->ClusterAt(i);
3074  for (j=0;j<jet2->GetNumberOfClusters();j++)
3075  {
3076  jetCluster2=jet2->ClusterAt(j);
3077  if (jetCluster1 == jetCluster2)
3078  {
3079  return kTRUE;
3080  }
3081  }
3082  }
3083  }
3084  return kFALSE;
3085 }
3086 
3087 Double_t AliAnalysisTaskFullpAJets_Eli_Mod::AreaWithinTPC(Double_t r,Double_t eta)
3088 {
3089  Double_t z;
3090  if (eta<(fTPCEtaMin+r))
3091  {
3092  z=eta-fTPCEtaMin;
3093  }
3094  else if(eta>(fTPCEtaMax-r))
3095  {
3096  z=fTPCEtaMax-eta;
3097  }
3098  else
3099  {
3100  z=r;
3101  }
3102  return r*r*TMath::Pi()-AreaEdge(r,z);
3103 }
3104 
3105 Double_t AliAnalysisTaskFullpAJets_Eli_Mod::AreaWithinEMCal(Double_t r,Double_t phi,Double_t eta)
3106 {
3107  Double_t x,y;
3108 
3109  if (phi<(fEMCalPhiMin-r) || phi>(fEMCalPhiMax+r))
3110  {
3111  x=-r;
3112  }
3113  else if (phi<(fEMCalPhiMin+r))
3114  {
3115  x=phi-fEMCalPhiMin;
3116  }
3117  else if (phi>(fEMCalPhiMin+r) && phi<(fEMCalPhiMax-r))
3118  {
3119  x=r;
3120  }
3121  else
3122  {
3123  x=fEMCalPhiMax-phi;
3124  }
3125 
3126  if (eta<(fEMCalEtaMin-r) || eta>(fEMCalEtaMax+r))
3127  {
3128  y=-r;
3129  }
3130  else if (eta<(fEMCalEtaMin+r))
3131  {
3132  y=eta-fEMCalEtaMin;
3133  }
3134  else if (eta>(fEMCalEtaMin+r) && eta<(fEMCalEtaMax-r))
3135  {
3136  y=r;
3137  }
3138  else
3139  {
3140  y=fEMCalEtaMax-eta;
3141  }
3142 
3143  if (x>=0 && y>=0)
3144  {
3145  if (TMath::Sqrt(x*x+y*y)>=r)
3146  {
3147  return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3148  }
3149  return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y)+AreaOverlap(r,x,y);
3150  }
3151  else if ((x>=r && y<0) || (y>=r && x<0))
3152  {
3153  return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
3154  }
3155  else if (x>0 && x<r && y<0)
3156  {
3157  Double_t a=TMath::Sqrt(r*r-x*x);
3158  Double_t b=TMath::Sqrt(r*r-y*y);
3159  if ((x-b)>0)
3160  {
3161  return r*r*TMath::ASin(b/r)+y*b;
3162  }
3163  else
3164  {
3165  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);
3166  }
3167  }
3168  else if (y>0 && y<r && x<0)
3169  {
3170  Double_t a=TMath::Sqrt(r*r-x*x);
3171  Double_t b=TMath::Sqrt(r*r-y*y);
3172  if ((y-a)>0)
3173  {
3174  return r*r*TMath::ASin(a/r)+x*a;
3175  }
3176  else
3177  {
3178  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);
3179  }
3180  }
3181  else
3182  {
3183  Double_t a=TMath::Sqrt(r*r-x*x);
3184  Double_t b=TMath::Sqrt(r*r-y*y);
3185  if ((x+b)<0)
3186  {
3187  return 0;
3188  }
3189  else
3190  {
3191  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);
3192  }
3193  }
3194 }
3195 
3196 Double_t AliAnalysisTaskFullpAJets_Eli_Mod::AreaEdge(Double_t r,Double_t z)
3197 {
3198  Double_t a=TMath::Sqrt(r*r-z*z);
3199  return r*r*TMath::ASin(a/r)-a*z;
3200 }
3201 
3202 Double_t AliAnalysisTaskFullpAJets_Eli_Mod::AreaOverlap(Double_t r,Double_t x,Double_t y)
3203 {
3204  Double_t a=TMath::Sqrt(r*r-x*x);
3205  Double_t b=TMath::Sqrt(r*r-y*y);
3206  return x*y-0.5*(x*a+y*b)+0.5*r*r*(TMath::ASin(b/r)-TMath::ASin(x/r));
3207 }
3208 
3209 Double_t AliAnalysisTaskFullpAJets_Eli_Mod::TransverseArea(Double_t r,Double_t psi0,Double_t phi,Double_t eta)
3210 {
3211  Double_t area_left=0;
3212  Double_t area_right=0;
3213  Double_t eta_a=0;
3214  Double_t eta_b=0;
3215  Double_t eta_up=0;
3216  Double_t eta_down=0;
3217 
3218  Double_t u=eta-fEMCalEtaMin;
3219  Double_t v=fEMCalEtaMax-eta;
3220 
3221  Double_t phi1=phi+u*TMath::Tan(psi0);
3222  Double_t phi2=phi-u*TMath::Tan(psi0);
3223  Double_t phi3=phi+v*TMath::Tan(psi0);
3224  Double_t phi4=phi-v*TMath::Tan(psi0);
3225 
3226  //Calculate the Left side area
3227  if (phi1>=fEMCalPhiMax)
3228  {
3229  eta_a=eta-u*((fEMCalPhiMax-phi)/(phi1-phi));
3230  }
3231  if (phi2<=fEMCalPhiMin)
3232  {
3233  eta_b=eta-u*((phi-fEMCalPhiMin)/(phi-phi2));
3234  }
3235 
3236  if ((phi1>=fEMCalPhiMax) && (phi2<=fEMCalPhiMin))
3237  {
3238  eta_up=TMath::Max(eta_a,eta_b);
3239  eta_down=TMath::Min(eta_a,eta_b);
3240 
3241  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);
3242  }
3243  else if (phi1>=fEMCalPhiMax)
3244  {
3245  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);
3246  }
3247  else if (phi2<=fEMCalPhiMin)
3248  {
3249  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);
3250  }
3251  else
3252  {
3253  area_left=0.5*(phi1-phi2+2*r*TMath::Tan(psi0))*(u-r);
3254  }
3255 
3256  // Calculate the Right side area
3257  if (phi3>=fEMCalPhiMax)
3258  {
3259  eta_a=eta+v*((fEMCalPhiMax-phi)/(phi3-phi));
3260  }
3261  if (phi4<=fEMCalPhiMin)
3262  {
3263  eta_b=eta+v*((phi-fEMCalPhiMin)/(phi-phi4));
3264  }
3265 
3266  if ((phi3>=fEMCalPhiMax) && (phi4<=fEMCalPhiMin))
3267  {
3268  eta_up=TMath::Max(eta_a,eta_b);
3269  eta_down=TMath::Min(eta_a,eta_b);
3270 
3271  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);
3272  }
3273  else if (phi3>=fEMCalPhiMax)
3274  {
3275  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);
3276  }
3277  else if (phi4<=fEMCalPhiMin)
3278  {
3279  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);
3280  }
3281  else
3282  {
3283  area_right=0.5*(phi3-phi4+2*r*TMath::Tan(psi0))*(v-r);
3284  }
3285  return area_left+area_right;
3286 }
3287 
3288 Double_t AliAnalysisTaskFullpAJets_Eli_Mod::MedianRhokT(Double_t *pTkTEntries, Double_t *RhokTEntries, Int_t nEntries)
3289 {
3290  // This function is used to calculate the median Rho kT value. The procedure is:
3291  // - Order the kT cluster array from highest rho value to lowest
3292  // - Exclude highest rho kT cluster
3293  // - Return the median rho value of the remaining subset
3294 
3295  // Sort Array
3296  const Double_t rho_min=-9.9999E+99;
3297  Int_t j,k;
3298  Double_t w[nEntries]; // Used for sorting
3299  Double_t smax=rho_min;
3300  Int_t sindex=-1;
3301 
3302  Double_t pTmax=0.0;
3303  Int_t pTmaxID=-1;
3304 
3305  for (j=0;j<nEntries;j++)
3306  {
3307  w[j]=0.0;
3308  }
3309 
3310  for (j=0;j<nEntries;j++)
3311  {
3312  if (pTkTEntries[j]>pTmax)
3313  {
3314  pTmax=pTkTEntries[j];
3315  pTmaxID=j;
3316  }
3317  }
3318 
3319  for (j=0;j<nEntries;j++)
3320  {
3321  for (k=0;k<nEntries;k++)
3322  {
3323  if (RhokTEntries[k]>smax)
3324  {
3325  smax=RhokTEntries[k];
3326  sindex=k;
3327  }
3328  }
3329  w[j]=smax;
3330  RhokTEntries[sindex]=rho_min;
3331  smax=rho_min;
3332  sindex=-1;
3333  }
3334  return w[nEntries/2];
3335 }
3336 
3337 
3338 //
3339 //------------------------------------------------------------------------------------
3340 //------------------------------------------------------------------------------------
3341 //------------------------------------------------------------------------------------
3342 //
3343 
3344 
3345 // AlipAJetData Class Member Defs
3346 // Constructors
3348 
3349  fName(0),
3350  fIsJetsFull(0),
3351  fnTotal(0),
3352  fnJets(0),
3353  fnJetsSC(0),
3354  fJetR(0),
3355  fSignalPt(0),
3356  fAreaCutFrac(0.6),
3357  fNEF(1.0),
3358  fSignalTrackBias(0),
3359  fPtMaxIndex(0),
3360  fPtMax(0),
3361  fPtSubLeadingIndex(0),
3362  fPtSubLeading(0),
3363  fJetsIndex(0),
3364  fJetsSCIndex(0),
3365  fIsJetInArray(0),
3366  fJetMaxChargedPt(0)
3367 {
3368  fnTotal=0;
3369  // Dummy constructor ALWAYS needed for I/O.
3370 }
3371 
3372 AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetData::AlipAJetData(const char *name, Bool_t isFull, Int_t nEntries) :
3373 
3374  fName(0),
3375  fIsJetsFull(0),
3376  fnTotal(0),
3377  fnJets(0),
3378  fnJetsSC(0),
3379  fJetR(0),
3380  fSignalPt(0),
3381  fAreaCutFrac(0.6),
3382  fNEF(1.0),
3383  fSignalTrackBias(0),
3384  fPtMaxIndex(0),
3385  fPtMax(0),
3386  fPtSubLeadingIndex(0),
3387  fPtSubLeading(0),
3388  fJetsIndex(0),
3389  fJetsSCIndex(0),
3390  fIsJetInArray(0),
3391  fJetMaxChargedPt(0)
3392 {
3393  SetName(name);
3394  SetIsJetsFull(isFull);
3395  SetTotalEntries(nEntries);
3396  SetLeading(0,-9.99E+099);
3397  SetSubLeading(0,-9.99E+099);
3398  SetSignalCut(0);
3399  SetAreaCutFraction(0.6);
3400  SetJetR(fJetR);
3402 }
3403 
3404 // Destructor
3406 {
3407  SetName("");
3408  SetIsJetsFull(kFALSE);
3409  SetTotalJets(0);
3410  SetTotalSignalJets(0);
3411  SetLeading(0,0);
3412  SetSubLeading(0,0);
3413  SetSignalCut(0);
3414  SetAreaCutFraction(0);
3415  SetJetR(0);
3416  SetNEF(0);
3417  SetSignalTrackPtBias(kFALSE);
3418 
3419  delete [] fJetsIndex;
3420  delete [] fJetsSCIndex;
3421  delete [] fIsJetInArray;
3422  delete [] fJetMaxChargedPt;
3423 
3424  fnTotal = 0;
3425 }
3426 
3427 // User Defined Sub-Routines
3429 {
3430  Int_t i=0;
3431  Int_t k=0;
3432  Int_t l=0;
3433  Double_t AreaThreshold = fAreaCutFrac*TMath::Pi()*TMath::Power(fJetR,2);
3434 
3435  // Initialize Jet Data
3436  for (i=0;i<nEntries;i++)
3437  {
3438  AliEmcalJet *myJet =(AliEmcalJet*) jetList->At(i);
3439 
3440  /*no cuts*///if (fIsJetInArray[i]==kTRUE && myJet->Area()>AreaThreshold)
3441  {
3442  SetJetIndex(i,k);
3443  if (myJet->Pt()>fPtMax)
3444  {
3445  SetSubLeading(fPtMaxIndex,fPtMax);
3446  SetLeading(i,myJet->Pt());
3447  }
3448  else if (myJet->Pt()>fPtSubLeading)
3449  {
3450  SetSubLeading(i,myJet->Pt());
3451  }
3452  // require leading charged constituent to have a pT greater then the signal threshold & Jet NEF to be less then the Signal Jet NEF cut
3453  fJetMaxChargedPt[i] = myJet->MaxTrackPt();
3454  if (fSignalTrackBias==kTRUE)
3455  {
3456  if (fJetMaxChargedPt[i]>=fSignalPt && myJet->NEF()<=fNEF)
3457  {
3458  SetSignalJetIndex(i,l);
3459  l++;
3460  }
3461  }
3462  else
3463  {
3464  if (myJet->Pt()>=fSignalPt && myJet->NEF()<=fNEF)
3465  {
3466  SetSignalJetIndex(i,l);
3467  l++;
3468  }
3469  }
3470  k++;
3471  }
3472  }
3473  SetTotalJets(k);
3474  SetTotalSignalJets(l);
3475 }
3476 
3477 // Setters
3479 {
3480  fName = name;
3481 }
3482 
3484 {
3485  fIsJetsFull = isFull;
3486 }
3487 
3489 {
3490  fnTotal = nEntries;
3491  fJetsIndex = new Int_t[fnTotal];
3492  fJetsSCIndex = new Int_t[fnTotal];
3493  fIsJetInArray = new Bool_t[fnTotal];
3494  fJetMaxChargedPt = new Double_t[fnTotal];
3495 }
3496 
3498 {
3499  fnJets = nJets;
3500 }
3501 
3503 {
3504  fnJetsSC = nSignalJets;
3505 }
3506 
3508 {
3509  fSignalPt = Pt;
3510 }
3511 
3513 {
3514  fPtMaxIndex = index;
3515  fPtMax = Pt;
3516 }
3517 
3519 {
3520  fPtSubLeadingIndex = index;
3521  fPtSubLeading = Pt;
3522 }
3523 
3525 {
3526  fJetsIndex[At] = index;
3527 }
3528 
3530 {
3531  fJetsSCIndex[At] = index;
3532 }
3533 
3535 {
3536  fIsJetInArray[At] = isInArray;
3537 }
3538 
3540 {
3541  fAreaCutFrac = areaFraction;
3542 }
3543 
3545 {
3546  fJetR = jetR;
3547 }
3548 
3550 {
3551  fNEF = nef;
3552 }
3553 
3555 {
3556  fSignalTrackBias = chargedBias;
3557 }
3558 
3559 // Getters
3561 {
3562  return fnTotal;
3563 }
3564 
3566 {
3567  return fnJets;
3568 }
3569 
3571 {
3572  return fnJetsSC;
3573 }
3574 
3576 {
3577  return fSignalPt;
3578 }
3579 
3581 {
3582  return fPtMaxIndex;
3583 }
3584 
3586 {
3587  return fPtMax;
3588 }
3589 
3591 {
3592  return fPtSubLeadingIndex;
3593 }
3594 
3596 {
3597  return fPtSubLeading;
3598 }
3599 
3601 {
3602  return fJetsIndex[At];
3603 }
3604 
3606 {
3607  return fJetsSCIndex[At];
3608 }
3609 
3611 {
3612  return fIsJetInArray[At];
3613 }
3614 
3616 {
3617  return fJetMaxChargedPt[At];
3618 }
3619 
3621 {
3622  return fNEF;
3623 }
3624 
3625 //
3626 //------------------------------------------------------------------------------------
3627 //------------------------------------------------------------------------------------
3628 //------------------------------------------------------------------------------------
3629 //
3630 
3631 // AlipAJetHistos Class Member Defs
3632 // Constructors
3634 
3635  fOutput(0),
3636 
3637  fh020Rho(0),
3638  fh80100Rho(0),
3639  fhRho(0),
3640  fhRhoCen(0),
3641  fh020BSPt(0),
3642  fh80100BSPt(0),
3643  fhBSPt(0),
3644  fhBSPtCen(0),
3645  fh020BSPtSignal(0),
3646  fh80100BSPtSignal(0),
3647  fhBSPtSignal(0),
3648  fhBSPtCenSignal(0),
3649  fh020DeltaPt(0),
3650  fh80100DeltaPt(0),
3651  fhDeltaPt(0),
3652  fhDeltaPtCen(0),
3653  fh020DeltaPtSignal(0),
3654  fh80100DeltaPtSignal(0),
3655  fhDeltaPtSignal(0),
3656  fhDeltaPtCenSignal(0),
3657  fh020DeltaPtNColl(0),
3658  fh80100DeltaPtNColl(0),
3659  fhDeltaPtNColl(0),
3660  fhDeltaPtCenNColl(0),
3661  fh020BckgFlucPt(0),
3662  fh80100BckgFlucPt(0),
3663  fhBckgFlucPt(0),
3664  fhBckgFlucPtCen(0),
3665 
3666  fpRho(0),
3667  fpLJetRho(0),
3668 
3669  fhJetPtEtaPhi(0),
3670  fhJetPtArea(0),
3671  fhJetConstituentPt(0),
3672  fhJetTracksPt(0),
3673  fhJetClustersPt(0),
3674  fhJetConstituentCounts(0),
3675  fhJetTracksCounts(0),
3676  fhJetClustersCounts(0),
3677  fhJetPtZConstituent(0),
3678  fhJetPtZTrack(0),
3679  fhJetPtZCluster(0),
3680  fhJetPtZLeadingConstituent(0),
3681  fhJetPtZLeadingTrack(0),
3682  fhJetPtZLeadingCluster(0),
3683 
3684  fhEventCentralityVsZNA(0),
3685  fhEventCentralityVsZNAPt(0),
3686 
3687  fNEFOutput(0),
3688  fhJetPtNEF(0),
3689  fhJetNEFInfo(0),
3690  fhJetNEFSignalInfo(0),
3691  fhClusterNEFInfo(0),
3692  fhClusterNEFSignalInfo(0),
3693  fhClusterShapeAll(0),
3694  fhClusterPtCellAll(0),
3695 
3696  fName(0),
3697  fCentralityTag(0),
3698  fCentralityBins(0),
3699  fCentralityLow(0),
3700  fCentralityUp(0),
3701  fPtBins(0),
3702  fPtLow(0),
3703  fPtUp(0),
3704  fRhoPtBins(0),
3705  fRhoPtLow(0),
3706  fRhoPtUp(0),
3707  fDeltaPtBins(0),
3708  fDeltaPtLow(0),
3709  fDeltaPtUp(0),
3710  fBckgFlucPtBins(0),
3711  fBckgFlucPtLow(0),
3712  fBckgFlucPtUp(0),
3713  fLJetPtBins(0),
3714  fLJetPtLow(0),
3715  fLJetPtUp(0),
3716  fRhoValue(0),
3717  fLChargedTrackPtBins(0),
3718  fLChargedTrackPtLow(0),
3719  fLChargedTrackPtUp(0),
3720  fDoNEFQAPlots(0),
3721  fDoNEFSignalOnly(1),
3722  fSignalTrackBias(0),
3723  fDoTHnSparse(0),
3724  fDo3DHistos(0),
3725  fNEFBins(0),
3726  fNEFLow(0),
3727  fNEFUp(0),
3728  fnDimJet(0),
3729  fnDimCluster(0),
3730  fEMCalPhiMin(1.39626),
3731  fEMCalPhiMax(3.26377),
3732  fEMCalEtaMin(-0.7),
3733  fEMCalEtaMax(0.7)
3734 
3735 {
3736  // Dummy constructor ALWAYS needed for I/O.
3737 }
3738 
3740 
3741  fOutput(0),
3742 
3743  fh020Rho(0),
3744  fh80100Rho(0),
3745  fhRho(0),
3746  fhRhoCen(0),
3747  fh020BSPt(0),
3748  fh80100BSPt(0),
3749  fhBSPt(0),
3750  fhBSPtCen(0),
3751  fh020BSPtSignal(0),
3752  fh80100BSPtSignal(0),
3753  fhBSPtSignal(0),
3754  fhBSPtCenSignal(0),
3755  fh020DeltaPt(0),
3756  fh80100DeltaPt(0),
3757  fhDeltaPt(0),
3758  fhDeltaPtCen(0),
3759  fh020DeltaPtSignal(0),
3760  fh80100DeltaPtSignal(0),
3761  fhDeltaPtSignal(0),
3762  fhDeltaPtCenSignal(0),
3763  fh020DeltaPtNColl(0),
3764  fh80100DeltaPtNColl(0),
3765  fhDeltaPtNColl(0),
3766  fhDeltaPtCenNColl(0),
3767  fh020BckgFlucPt(0),
3768  fh80100BckgFlucPt(0),
3769  fhBckgFlucPt(0),
3770  fhBckgFlucPtCen(0),
3771 
3772  fpRho(0),
3773  fpLJetRho(0),
3774 
3775  fhJetPtEtaPhi(0),
3776  fhJetPtArea(0),
3777  fhJetConstituentPt(0),
3778  fhJetTracksPt(0),
3779  fhJetClustersPt(0),
3780  fhJetConstituentCounts(0),
3781  fhJetTracksCounts(0),
3782  fhJetClustersCounts(0),
3783  fhJetPtZConstituent(0),
3784  fhJetPtZTrack(0),
3785  fhJetPtZCluster(0),
3786  fhJetPtZLeadingConstituent(0),
3787  fhJetPtZLeadingTrack(0),
3788  fhJetPtZLeadingCluster(0),
3789 
3790  fhEventCentralityVsZNA(0),
3791  fhEventCentralityVsZNAPt(0),
3792 
3793  fNEFOutput(0),
3794  fhJetPtNEF(0),
3795  fhJetNEFInfo(0),
3796  fhJetNEFSignalInfo(0),
3797  fhClusterNEFInfo(0),
3798  fhClusterNEFSignalInfo(0),
3799  fhClusterShapeAll(0),
3800  fhClusterPtCellAll(0),
3801 
3802  fName(0),
3803  fCentralityTag(0),
3804  fCentralityBins(0),
3805  fCentralityLow(0),
3806  fCentralityUp(0),
3807  fPtBins(0),
3808  fPtLow(0),
3809  fPtUp(0),
3810  fRhoPtBins(0),
3811  fRhoPtLow(0),
3812  fRhoPtUp(0),
3813  fDeltaPtBins(0),
3814  fDeltaPtLow(0),
3815  fDeltaPtUp(0),
3816  fBckgFlucPtBins(0),
3817  fBckgFlucPtLow(0),
3818  fBckgFlucPtUp(0),
3819  fLJetPtBins(0),
3820  fLJetPtLow(0),
3821  fLJetPtUp(0),
3822  fRhoValue(0),
3823  fLChargedTrackPtBins(0),
3824  fLChargedTrackPtLow(0),
3825  fLChargedTrackPtUp(0),
3826  fDoNEFQAPlots(0),
3827  fDoNEFSignalOnly(1),
3828  fSignalTrackBias(0),
3829  fDoTHnSparse(0),
3830  fDo3DHistos(0),
3831  fNEFBins(0),
3832  fNEFLow(0),
3833  fNEFUp(0),
3834  fnDimJet(0),
3835  fnDimCluster(0),
3836  fEMCalPhiMin(1.39626),
3837  fEMCalPhiMax(3.26377),
3838  fEMCalEtaMin(-0.7),
3839  fEMCalEtaMax(0.7)
3840 
3841 {
3842  SetName(name);
3843  SetCentralityTag("ZNA");
3844  SetCentralityRange(100,0,100);
3845  SetPtRange(250,-50,200);
3846  SetRhoPtRange(500,0,50);
3847  SetDeltaPtRange(200,-100,100);
3849  SetLeadingJetPtRange(200,0,200);
3850  SetLeadingChargedTrackPtRange(100,0,100);
3851  SetNEFRange(100,0,1);
3852  DoNEFQAPlots(kFALSE);
3853 
3854  Init();
3855 }
3856 
3857 AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::AlipAJetHistos(const char *name, TString centag, Bool_t doNEF) :
3858 
3859  fOutput(0),
3860 
3861  fh020Rho(0),
3862  fh80100Rho(0),
3863  fhRho(0),
3864  fhRhoCen(0),
3865  fh020BSPt(0),
3866  fh80100BSPt(0),
3867  fhBSPt(0),
3868  fhBSPtCen(0),
3869  fh020BSPtSignal(0),
3870  fh80100BSPtSignal(0),
3871  fhBSPtSignal(0),
3872  fhBSPtCenSignal(0),
3873  fh020DeltaPt(0),
3874  fh80100DeltaPt(0),
3875  fhDeltaPt(0),
3876  fhDeltaPtCen(0),
3877  fh020DeltaPtSignal(0),
3878  fh80100DeltaPtSignal(0),
3879  fhDeltaPtSignal(0),
3880  fhDeltaPtCenSignal(0),
3881  fh020DeltaPtNColl(0),
3882  fh80100DeltaPtNColl(0),
3883  fhDeltaPtNColl(0),
3884  fhDeltaPtCenNColl(0),
3885  fh020BckgFlucPt(0),
3886  fh80100BckgFlucPt(0),
3887  fhBckgFlucPt(0),
3888  fhBckgFlucPtCen(0),
3889 
3890  fpRho(0),
3891  fpLJetRho(0),
3892 
3893  fhJetPtEtaPhi(0),
3894  fhJetPtArea(0),
3895  fhJetConstituentPt(0),
3896  fhJetTracksPt(0),
3897  fhJetClustersPt(0),
3898  fhJetConstituentCounts(0),
3899  fhJetTracksCounts(0),
3900  fhJetClustersCounts(0),
3901  fhJetPtZConstituent(0),
3902  fhJetPtZTrack(0),
3903  fhJetPtZCluster(0),
3904  fhJetPtZLeadingConstituent(0),
3905  fhJetPtZLeadingTrack(0),
3906  fhJetPtZLeadingCluster(0),
3907 
3908  fhEventCentralityVsZNA(0),
3909  fhEventCentralityVsZNAPt(0),
3910 
3911  fNEFOutput(0),
3912  fhJetPtNEF(0),
3913  fhJetNEFInfo(0),
3914  fhJetNEFSignalInfo(0),
3915  fhClusterNEFInfo(0),
3916  fhClusterNEFSignalInfo(0),
3917  fhClusterShapeAll(0),
3918  fhClusterPtCellAll(0),
3919 
3920  fName(0),
3921  fCentralityTag(0),
3922  fCentralityBins(0),
3923  fCentralityLow(0),
3924  fCentralityUp(0),
3925  fPtBins(0),
3926  fPtLow(0),
3927  fPtUp(0),
3928  fRhoPtBins(0),
3929  fRhoPtLow(0),
3930  fRhoPtUp(0),
3931  fDeltaPtBins(0),
3932  fDeltaPtLow(0),
3933  fDeltaPtUp(0),
3934  fBckgFlucPtBins(0),
3935  fBckgFlucPtLow(0),
3936  fBckgFlucPtUp(0),
3937  fLJetPtBins(0),
3938  fLJetPtLow(0),
3939  fLJetPtUp(0),
3940  fRhoValue(0),
3941  fLChargedTrackPtBins(0),
3942  fLChargedTrackPtLow(0),
3943  fLChargedTrackPtUp(0),
3944  fDoNEFQAPlots(0),
3945  fDoNEFSignalOnly(1),
3946  fSignalTrackBias(0),
3947  fDoTHnSparse(0),
3948  fDo3DHistos(0),
3949  fNEFBins(0),
3950  fNEFLow(0),
3951  fNEFUp(0),
3952  fnDimJet(0),
3953  fnDimCluster(0),
3954  fEMCalPhiMin(1.39626),
3955  fEMCalPhiMax(3.26377),
3956  fEMCalEtaMin(-0.7),
3957  fEMCalEtaMax(0.7)
3958 
3959 {
3960  SetName(name);
3961  SetCentralityTag(centag.Data());
3962  SetCentralityRange(100,0,100);
3963  SetPtRange(250,-50,200);
3964  SetRhoPtRange(500,0,50);
3965  SetDeltaPtRange(200,-100,100);
3967  SetLeadingJetPtRange(200,0,200);
3968  SetLeadingChargedTrackPtRange(100,0,100);
3969  SetNEFRange(100,0,1);
3970  DoNEFQAPlots(doNEF);
3971 
3972  Init();
3973 }
3974 
3975 AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::AlipAJetHistos(const char *name, TString centag, Bool_t doNEF, Bool_t doNEFSignalOnly, Bool_t doTHnSparse, Bool_t do3DPlotting) :
3976 
3977  fOutput(0),
3978 
3979  fh020Rho(0),
3980  fh80100Rho(0),
3981  fhRho(0),
3982  fhRhoCen(0),
3983  fh020BSPt(0),
3984  fh80100BSPt(0),
3985  fhBSPt(0),
3986  fhBSPtCen(0),
3987  fh020BSPtSignal(0),
3988  fh80100BSPtSignal(0),
3989  fhBSPtSignal(0),
3990  fhBSPtCenSignal(0),
3991  fh020DeltaPt(0),
3992  fh80100DeltaPt(0),
3993  fhDeltaPt(0),
3994  fhDeltaPtCen(0),
3995  fh020DeltaPtSignal(0),
3996  fh80100DeltaPtSignal(0),
3997  fhDeltaPtSignal(0),
3998  fhDeltaPtCenSignal(0),
3999  fh020DeltaPtNColl(0),
4000  fh80100DeltaPtNColl(0),
4001  fhDeltaPtNColl(0),
4002  fhDeltaPtCenNColl(0),
4003  fh020BckgFlucPt(0),
4004  fh80100BckgFlucPt(0),
4005  fhBckgFlucPt(0),
4006  fhBckgFlucPtCen(0),
4007 
4008  fpRho(0),
4009  fpLJetRho(0),
4010 
4011  fhJetPtEtaPhi(0),
4012  fhJetPtArea(0),
4013  fhJetConstituentPt(0),
4014  fhJetTracksPt(0),
4015  fhJetClustersPt(0),
4016  fhJetConstituentCounts(0),
4017  fhJetTracksCounts(0),
4018  fhJetClustersCounts(0),
4019  fhJetPtZConstituent(0),
4020  fhJetPtZTrack(0),
4021  fhJetPtZCluster(0),
4022  fhJetPtZLeadingConstituent(0),
4023  fhJetPtZLeadingTrack(0),
4024  fhJetPtZLeadingCluster(0),
4025 
4026  fhEventCentralityVsZNA(0),
4027  fhEventCentralityVsZNAPt(0),
4028 
4029  fNEFOutput(0),
4030  fhJetPtNEF(0),
4031  fhJetNEFInfo(0),
4032  fhJetNEFSignalInfo(0),
4033  fhClusterNEFInfo(0),
4034  fhClusterNEFSignalInfo(0),
4035  fhClusterShapeAll(0),
4036  fhClusterPtCellAll(0),
4037 
4038  fName(0),
4039  fCentralityTag(0),
4040  fCentralityBins(0),
4041  fCentralityLow(0),
4042  fCentralityUp(0),
4043  fPtBins(0),
4044  fPtLow(0),
4045  fPtUp(0),
4046  fRhoPtBins(0),
4047  fRhoPtLow(0),
4048  fRhoPtUp(0),
4049  fDeltaPtBins(0),
4050  fDeltaPtLow(0),
4051  fDeltaPtUp(0),
4052  fBckgFlucPtBins(0),
4053  fBckgFlucPtLow(0),
4054  fBckgFlucPtUp(0),
4055  fLJetPtBins(0),
4056  fLJetPtLow(0),
4057  fLJetPtUp(0),
4058  fRhoValue(0),
4059  fLChargedTrackPtBins(0),
4060  fLChargedTrackPtLow(0),
4061  fLChargedTrackPtUp(0),
4062  fDoNEFQAPlots(0),
4063  fDoNEFSignalOnly(1),
4064  fSignalTrackBias(0),
4065  fDoTHnSparse(0),
4066  fDo3DHistos(0),
4067  fNEFBins(0),
4068  fNEFLow(0),
4069  fNEFUp(0),
4070  fnDimJet(0),
4071  fnDimCluster(0),
4072  fEMCalPhiMin(1.39626),
4073  fEMCalPhiMax(3.26377),
4074  fEMCalEtaMin(-0.7),
4075  fEMCalEtaMax(0.7)
4076 
4077 {
4078  SetName(name);
4079  SetCentralityTag(centag.Data());
4080  SetCentralityRange(100,0,100);
4081  SetPtRange(250,-50,200);
4082  SetRhoPtRange(500,0,50);
4083  SetDeltaPtRange(200,-100,100);
4085  SetLeadingJetPtRange(200,0,200);
4086  SetLeadingChargedTrackPtRange(100,0,100);
4087  SetNEFRange(100,0,1);
4088  DoNEFQAPlots(doNEF);
4089  DoNEFSignalOnly(doNEFSignalOnly);
4090  DoTHnSparse(doTHnSparse);
4091  Do3DPlotting(do3DPlotting);
4092 
4093  Init();
4094 }
4095 
4096 // Destructor
4098 {
4099  if (fOutput)
4100  {
4101  delete fOutput;
4102  }
4103 }
4104 
4106 {
4107  Int_t i;
4108 
4109  // Initialize Private Variables
4110  Int_t TCBins = 100;
4111  fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
4112  fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
4113  fEMCalEtaMin=-0.7;
4114  fEMCalEtaMax=0.7;
4115 
4116  fOutput = new TList();
4117  fOutput->SetOwner();
4118  fOutput->SetName(fName);
4119 
4120  TString RhoString="";
4121  TString PtString="";
4122  TString DeltaPtString="";
4123  TString BckgFlucPtString="";
4124  TString CentralityString;
4125  TString EventCentralityString;
4126  CentralityString = Form("Centrality (%s) ",fCentralityTag.Data());
4127  EventCentralityString = Form("%s vs ZNA Centrality ",fCentralityTag.Data());
4128 
4129  // Rho Spectral Plots
4130  RhoString = Form("%d-%d Centrality, Rho Spectrum",0,20);
4131  fh020Rho = new TH1F("fh020Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
4132  fh020Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4133  fh020Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
4134  fh020Rho->Sumw2();
4135 
4136  RhoString = Form("%d-%d Centrality, Rho Spectrum",80,100);
4137  fh80100Rho = new TH1F("fh80100Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
4138  fh80100Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4139  fh80100Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
4140  fh80100Rho->Sumw2();
4141 
4142  RhoString = Form("%d-%d Centrality, Rho Spectrum",0,100);
4143  fhRho = new TH1F("fhRho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
4144  fhRho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4145  fhRho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
4146  fhRho->Sumw2();
4147 
4148  RhoString = "Rho Spectrum vs Centrality";
4149  fhRhoCen = new TH2F("fhRhoCen",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
4150  fhRhoCen->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
4151  fhRhoCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4152  fhRhoCen->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
4153  fhRhoCen->Sumw2();
4154 
4155  // Background Subtracted Plots
4156  PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,20);
4157  fh020BSPt = new TH1F("fh020BSPt",PtString,fPtBins,fPtLow,fPtUp);
4158  fh020BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4159  fh020BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4160  fh020BSPt->Sumw2();
4161 
4162  PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",80,100);
4163  fh80100BSPt = new TH1F("fh80100BSPt",PtString,fPtBins,fPtLow,fPtUp);
4164  fh80100BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4165  fh80100BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4166  fh80100BSPt->Sumw2();
4167 
4168  PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,100);
4169  fhBSPt = new TH1F("fhBSPt",PtString,fPtBins,fPtLow,fPtUp);
4170  fhBSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4171  fhBSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4172  fhBSPt->Sumw2();
4173 
4174  PtString = "Background Subtracted Jet Spectrum vs Centrality";
4175  fhBSPtCen = new TH2F("fhBSPtCen",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
4176  fhBSPtCen->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4177  fhBSPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4178  fhBSPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4179  fhBSPtCen->Sumw2();
4180 
4181  PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,20);
4182  fh020BSPtSignal = new TH1F("fh020BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
4183  fh020BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4184  fh020BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4185  fh020BSPtSignal->Sumw2();
4186 
4187  PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",80,100);
4188  fh80100BSPtSignal = new TH1F("fh80100BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
4189  fh80100BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4190  fh80100BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4191  fh80100BSPtSignal->Sumw2();
4192 
4193  PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,100);
4194  fhBSPtSignal = new TH1F("fhBSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
4195  fhBSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4196  fhBSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4197  fhBSPtSignal->Sumw2();
4198 
4199  PtString = "Background Subtracted Signal Jet Spectrum vs Centrality";
4200  fhBSPtCenSignal = new TH2F("fhBSPtCenSignal",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
4201  fhBSPtCenSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
4202  fhBSPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4203  fhBSPtCenSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4204  fhBSPtCenSignal->Sumw2();
4205 
4206  // Delta Pt Plots with RC at least 2R away from Leading Signal
4207  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
4208  fh020DeltaPt = new TH1F("fh020DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4209  fh020DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4210  fh020DeltaPt->GetYaxis()->SetTitle("Probability Density");
4211  fh020DeltaPt->Sumw2();
4212 
4213  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
4214  fh80100DeltaPt = new TH1F("fh80100DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4215  fh80100DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4216  fh80100DeltaPt->GetYaxis()->SetTitle("Probability Density");
4217  fh80100DeltaPt->Sumw2();
4218 
4219  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
4220  fhDeltaPt = new TH1F("fhDeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4221  fhDeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4222  fhDeltaPt->GetYaxis()->SetTitle("Probability Density");
4223  fhDeltaPt->Sumw2();
4224 
4225  DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
4226  fhDeltaPtCen = new TH2F("fhDeltaPtCen",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
4227  fhDeltaPtCen->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4228  fhDeltaPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4229  fhDeltaPtCen->GetZaxis()->SetTitle("Probability Density");
4230  fhDeltaPtCen->Sumw2();
4231 
4232  // Delta Pt Plots with no spatial restrictions on RC
4233  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
4234  fh020DeltaPtSignal = new TH1F("fh020DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4235  fh020DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4236  fh020DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
4237  fh020DeltaPtSignal->Sumw2();
4238 
4239  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
4240  fh80100DeltaPtSignal = new TH1F("fh80100DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4241  fh80100DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4242  fh80100DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
4243  fh80100DeltaPtSignal->Sumw2();
4244 
4245  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
4246  fhDeltaPtSignal = new TH1F("fhDeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4247  fhDeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4248  fhDeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
4249  fhDeltaPtSignal->Sumw2();
4250 
4251  DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
4252  fhDeltaPtCenSignal = new TH2F("fhDeltaPtCenSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
4253  fhDeltaPtCenSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4254  fhDeltaPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4255  fhDeltaPtCenSignal->GetZaxis()->SetTitle("Probability Density");
4256  fhDeltaPtCenSignal->Sumw2();
4257 
4258  // Delta Pt Plots with NColl restrictions on RC
4259  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
4260  fh020DeltaPtNColl = new TH1F("fh020DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4261  fh020DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4262  fh020DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
4263  fh020DeltaPtNColl->Sumw2();
4264 
4265  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
4266  fh80100DeltaPtNColl = new TH1F("fh80100DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4267  fh80100DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4268  fh80100DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
4269  fh80100DeltaPtNColl->Sumw2();
4270 
4271  DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
4272  fhDeltaPtNColl = new TH1F("fhDeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
4273  fhDeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4274  fhDeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
4275  fhDeltaPtNColl->Sumw2();
4276 
4277  DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
4278  fhDeltaPtCenNColl = new TH2F("fhDeltaPtCenNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
4279  fhDeltaPtCenNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
4280  fhDeltaPtCenNColl->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4281  fhDeltaPtCenNColl->GetZaxis()->SetTitle("Probability Density");
4282  fhDeltaPtCenNColl->Sumw2();
4283 
4284  // Background Fluctuations Pt Plots
4285  BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,20);
4286  fh020BckgFlucPt = new TH1F("fh020BckgFlucPt",PtString,fPtBins,fPtLow,fPtUp);
4287  fh020BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4288  fh020BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4289  fh020BckgFlucPt->Sumw2();
4290 
4291  BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",80,100);
4292  fh80100BckgFlucPt = new TH1F("fh80100BckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
4293  fh80100BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4294  fh80100BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4295  fh80100BckgFlucPt->Sumw2();
4296 
4297  BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,100);
4298  fhBckgFlucPt = new TH1F("fhBckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
4299  fhBckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4300  fhBckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4301  fhBckgFlucPt->Sumw2();
4302 
4303  BckgFlucPtString = "Background Fluctuation p_{T} Spectrum vs Centrality";
4304  fhBckgFlucPtCen = new TH2F("fhBckgFlucPtCen",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
4305  fhBckgFlucPtCen->GetXaxis()->SetTitle("#p_{T} (GeV/c)");
4306  fhBckgFlucPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
4307  fhBckgFlucPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
4308  fhBckgFlucPtCen->Sumw2();
4309 
4310  // Background Density vs Centrality Profile
4311  RhoString = "Background Density vs Centrality";
4312  fpRho = new TProfile("fpRho",RhoString,fCentralityBins,fCentralityLow,fCentralityUp);
4313  fpRho->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
4314  fpRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
4315 
4316  // Background Density vs Leading Jet Profile
4317  fpLJetRho = new TProfile("fpLJetRho","#rho vs Leading Jet p_{T}",fLJetPtBins,fLJetPtLow,fLJetPtUp);
4318  fpLJetRho->GetXaxis()->SetTitle("Leading Jet p_{T}");
4319  fpLJetRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
4320 
4321  // Jet pT vs Area
4322  Int_t JetPtAreaBins=200;
4323  Double_t JetPtAreaLow=0.0;
4324  Double_t JetPtAreaUp=2.0;
4325 
4326  fhJetPtArea = new TH2F("fhJetPtArea","Jet Area Distribution",fPtBins,fPtLow,fPtUp,JetPtAreaBins,JetPtAreaLow,JetPtAreaUp);
4327  fhJetPtArea->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4328  fhJetPtArea->GetYaxis()->SetTitle("A_{jet}");
4329  fhJetPtArea->GetZaxis()->SetTitle("1/N_{Events} dN/dA_{jet}dp_{T}");
4330  fhJetPtArea->Sumw2();
4331 
4332  // Jet pT vs Constituent pT
4333  fhJetConstituentPt = new TH2F("fhJetConstituentPt","Jet constituents p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
4334  fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4335  fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
4336  fhJetConstituentPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,track}");
4337  fhJetConstituentPt->Sumw2();
4338 
4339  fhJetTracksPt = new TH2F("fhJetTracksPt","Jet constituents Tracks p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
4340  fhJetTracksPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4341  fhJetTracksPt->GetYaxis()->SetTitle("Constituent Track p_{T} (GeV/c)");
4342  fhJetTracksPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,track}");
4343  fhJetTracksPt->Sumw2();
4344 
4345  fhJetClustersPt = new TH2F("fhJetClustersPt","Jet constituents Clusters p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
4346  fhJetClustersPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4347  fhJetClustersPt->GetYaxis()->SetTitle("Constituent Cluster p_{T} (GeV/c)");
4348  fhJetClustersPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,cluster}");
4349  fhJetClustersPt->Sumw2();
4350 
4351  // Jet pT vs Constituent Counts
4352  fhJetConstituentCounts = new TH2F("fhJetConstituentCounts","Jet constituents distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
4353  fhJetConstituentCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4354  fhJetConstituentCounts->GetYaxis()->SetTitle("Constituent Count");
4355  fhJetConstituentCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{constituent}");
4356  fhJetConstituentCounts->Sumw2();
4357 
4358  fhJetTracksCounts = new TH2F("fhJetTracksCounts","Jet constituents Tracks distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
4359  fhJetTracksCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4360  fhJetTracksCounts->GetYaxis()->SetTitle("Constituent Track Count");
4361  fhJetTracksCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{track}");
4362  fhJetTracksCounts->Sumw2();
4363 
4364  fhJetClustersCounts = new TH2F("fhJetClustersCounts","Jet constituents Clusters distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
4365  fhJetClustersCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4366  fhJetClustersCounts->GetYaxis()->SetTitle("Constituent Cluster Count");
4367  fhJetClustersCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{cluster}");
4368  fhJetClustersCounts->Sumw2();
4369 
4370  // Jet pT vs z_{constituent/track/cluster}
4371  fhJetPtZConstituent = new TH2F("fhJetPtZConstituent","Jet z_{constituent} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4372  fhJetPtZConstituent->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4373  fhJetPtZConstituent->GetYaxis()->SetTitle("z_{constituent}");
4374  fhJetPtZConstituent->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{constituent}");
4375  fhJetPtZConstituent->Sumw2();
4376 
4377  fhJetPtZTrack = new TH2F("fhJetPtZTrack","Jet z_{track} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4378  fhJetPtZTrack->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4379  fhJetPtZTrack->GetYaxis()->SetTitle("z_{track}");
4380  fhJetPtZTrack->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{track}");
4381  fhJetPtZTrack->Sumw2();
4382 
4383  fhJetPtZCluster = new TH2F("fhJetPtZCluster","Jet z_{cluster} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4384  fhJetPtZCluster->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4385  fhJetPtZCluster->GetYaxis()->SetTitle("z_{cluster}");
4386  fhJetPtZCluster->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{cluster}");
4387  fhJetPtZCluster->Sumw2();
4388 
4389  // Jet pT vs z_Leading{constituent/track/cluster}
4390  fhJetPtZLeadingConstituent = new TH2F("fhJetPtZLeadingConstituent","Jet z_{Leading,constituent} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4391  fhJetPtZLeadingConstituent->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4392  fhJetPtZLeadingConstituent->GetYaxis()->SetTitle("z_{Leading,constituent}");
4393  fhJetPtZLeadingConstituent->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{constituent}");
4394  fhJetPtZLeadingConstituent->Sumw2();
4395 
4396  fhJetPtZLeadingTrack = new TH2F("fhJetPtZLeadingTrack","Jet z_{Leading,track} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4397  fhJetPtZLeadingTrack->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4398  fhJetPtZLeadingTrack->GetYaxis()->SetTitle("z_{Leading,track}");
4399  fhJetPtZLeadingTrack->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{track}");
4400  fhJetPtZLeadingTrack->Sumw2();
4401 
4402  fhJetPtZLeadingCluster = new TH2F("fhJetPtZLeadingCluster","Jet z_{Leading,cluster} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
4403  fhJetPtZLeadingCluster->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4404  fhJetPtZLeadingCluster->GetYaxis()->SetTitle("z_{Leading,cluster}");
4405  fhJetPtZLeadingCluster->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{cluster}");
4406  fhJetPtZLeadingCluster->Sumw2();
4407 
4408  // Event Centralities vs Leading Jet Pt
4409  fhEventCentralityVsZNA = new TH2F("fhEventCentralityVsZNA",EventCentralityString,fCentralityBins,fCentralityLow,fCentralityUp,fCentralityBins,fCentralityLow,fCentralityUp);
4410  fhEventCentralityVsZNA->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
4411  fhEventCentralityVsZNA->GetYaxis()->SetTitle("Centrality (ZNA)");
4412  fhEventCentralityVsZNA->GetZaxis()->SetTitle("Probability Density");
4413  fhEventCentralityVsZNA->Sumw2();
4414 
4415  // 3D Histos
4416  if (fDo3DHistos == kTRUE)
4417  {
4418  // Jet pT, Eta, Phi
4419  fhJetPtEtaPhi = new TH3F("fhJetPtEtaPhi","Jet p_{T} vs #eta-#varphi",fPtBins,fPtLow,fPtUp,TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
4420  fhJetPtEtaPhi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4421  fhJetPtEtaPhi->GetYaxis()->SetTitle("#eta");
4422  fhJetPtEtaPhi->GetZaxis()->SetTitle("#varphi");
4423  fhJetPtEtaPhi->Sumw2();
4424 
4425  EventCentralityString = Form("%s vs ZNA Centrality vs Leading Jet p_{T} ",fCentralityTag.Data());
4426  fhEventCentralityVsZNAPt = new TH3F("fhEventCentralityVsZNAPt",EventCentralityString,fCentralityBins,fCentralityLow,fCentralityUp,fCentralityBins,fCentralityLow,fCentralityUp,fPtBins,fPtLow,fPtUp);
4427  fhEventCentralityVsZNAPt->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
4428  fhEventCentralityVsZNAPt->GetYaxis()->SetTitle("Centrality (ZNA)");
4429  fhEventCentralityVsZNAPt->GetZaxis()->SetTitle("Leading Jet p_{T} (GeV/c)");
4430  fhEventCentralityVsZNAPt->Sumw2();
4431 
4432  fOutput->Add(fhJetPtEtaPhi);
4433  fOutput->Add(fhEventCentralityVsZNAPt);
4434  }
4435 
4436  // Neutral Energy Fraction Histograms & QA
4437  if (fDoNEFQAPlots==kTRUE)
4438  {
4439  fNEFOutput = new TList();
4440  fNEFOutput->SetOwner();
4441  fNEFOutput->SetName("ListNEFQAPlots");
4442 
4443  fhJetPtNEF = new TH2F("fhJetPtNEF","Jet p_{T} vs NEF",fPtBins,fPtLow,fPtUp,fNEFBins,fNEFLow,fNEFUp);
4444  fhJetPtNEF->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
4445  fhJetPtNEF->GetYaxis()->SetTitle("Neutral Energy Fraction");
4446  fhJetPtNEF->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T}dNEF");
4447  fhJetPtNEF->Sumw2();
4448 
4449  SetNEFJetDimensions(10); // Order: nef,Jet Pt,Eta,Phi,Centrality,Constituent mult,Charged mult, Neutral mult, z_leading
4450  SetNEFClusterDimensions(11); // Order: nef,Jet Pt,Eta,Phi,Centrality,F_Cross,z_leading,t_cell,deltat_cell,Cell Count, E_Cluster
4451 
4452  Int_t dimJetBins[fnDimJet];
4453  Double_t dimJetLow[fnDimJet];
4454  Double_t dimJetUp[fnDimJet];
4455 
4456  Int_t dimClusterBins[fnDimCluster];
4457  Double_t dimClusterLow[fnDimCluster];
4458  Double_t dimClusterUp[fnDimCluster];
4459 
4460  // Establish dimensinality bin counts and bounds
4461  // NEF
4462  dimJetBins[0] = dimClusterBins[0] = fNEFBins;
4463  dimJetLow[0] = dimClusterLow[0] = fNEFLow;
4464  dimJetUp[0] = dimClusterUp[0] = fNEFUp;
4465 
4466  // Jet Pt
4467  dimJetBins[1] = dimClusterBins[1] = fPtBins;
4468  dimJetLow[1] = dimClusterLow[1] = fPtLow;
4469  dimJetUp[1] = dimClusterUp[1] = fPtUp;
4470 
4471  // Eta-Phi
4472  dimJetBins[2] = dimJetBins[3] = dimClusterBins[2] = dimClusterBins[3] = TCBins;
4473  dimJetLow[2] = dimClusterLow[2] = fEMCalEtaMin;
4474  dimJetUp[2] = dimClusterUp[2] = fEMCalEtaMax;
4475  dimJetLow[3] = dimClusterLow[3] = fEMCalPhiMin;
4476  dimJetUp[3] = dimClusterUp[3] = fEMCalPhiMax;
4477 
4478  // Centrality
4479  dimJetBins[4] = dimClusterBins[4] = fCentralityBins;
4480  dimJetLow[4] = dimClusterLow[4] = fCentralityLow;
4481  dimJetUp[4] = dimClusterUp[4] = fCentralityUp;
4482 
4483  // z_leading
4484  dimJetBins[5] = dimClusterBins[5] = TCBins;
4485  dimJetLow[5] = dimClusterLow[5] = 0.0;
4486  dimJetUp[5] = dimClusterUp[5] = 1.0;
4487 
4488  // Jets Constituent Multiplicity Info {Total,Charged,Neutral}
4489  for (i=6;i<9;i++)
4490  {
4491  dimJetBins[i] = TCBins;
4492  dimJetLow[i] = 0.0;
4493  dimJetUp[i] = (Double_t)TCBins;
4494  }
4495 
4496  // z_leading^track
4497  dimJetBins[9] = TCBins;
4498  dimJetLow[9] = 0.0;
4499  dimJetUp[9] = 1.0;
4500 
4501  // Cluster E
4502  dimClusterBins[6] = fPtBins;
4503  dimClusterLow[6] = fPtLow;
4504  dimClusterUp[6] = fPtUp;
4505 
4506  // Cluster F_Cross
4507  dimClusterBins[7] = TCBins;
4508  dimClusterLow[7] = 0.0;
4509  dimClusterUp[7] = 1.0;
4510 
4511  // Cluster t_cell
4512  dimClusterBins[8] = 400;
4513  dimClusterLow[8] = -2e-07;
4514  dimClusterUp[8] = 2e-07;
4515 
4516  // Cluster delta t_cell
4517  dimClusterBins[9] = 100;
4518  dimClusterLow[9] = 0.0;
4519  dimClusterUp[9] = 1e-07;
4520 
4521  // Cluster Cell Count
4522  dimClusterBins[10] = TCBins;
4523  dimClusterLow[10] = 0.0;
4524  dimClusterUp[10] = 100.0;
4525 
4526  if (fDoTHnSparse == kTRUE)
4527  {
4528  fhJetNEFSignalInfo = new THnSparseF("fhJetNEFSignalInfo","Signal Jet NEF Information Histogram",fnDimJet,dimJetBins,dimJetLow,dimJetUp);
4529  fhJetNEFSignalInfo->Sumw2();
4530 
4531  fhClusterNEFSignalInfo = new THnSparseF("fhClusterNEFSignalInfo","Signal Jet NEF Cluster Information Histogram",fnDimCluster,dimClusterBins,dimClusterLow,dimClusterUp);
4532  fhClusterNEFSignalInfo->Sumw2();
4533 
4534  // Cluster Shape QA
4535  fhClusterShapeAll = new TH1F("fhClusterShapeAll","Cluster Shape of all CaloClustersCorr",10*TCBins,0,10*TCBins);
4536  fhClusterShapeAll->GetXaxis()->SetTitle("Cells");
4537  fhClusterShapeAll->GetYaxis()->SetTitle("1/N_{Events} dN/dCells");
4538  fhClusterShapeAll->Sumw2();
4539 
4540  fhClusterPtCellAll = new TH2F("fhClusterPtCellAll","Cluster p_{T} vs Cluster Shape of all CaloClustersCorr",fPtBins,fPtLow,fPtUp,10*TCBins,0,10*TCBins);
4541  fhClusterPtCellAll->GetXaxis()->SetTitle("p_{T} (GeV/c)");
4542  fhClusterPtCellAll->GetYaxis()->SetTitle("Cells");
4543  fhClusterPtCellAll->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}dCells");
4544  fhClusterPtCellAll->Sumw2();
4545 
4546  if (fDoNEFSignalOnly == kFALSE)
4547  {
4548  fhJetNEFInfo = new THnSparseF("fhJetNEFInfo","Jet NEF Information Histogram",fnDimJet,dimJetBins,dimJetLow,dimJetUp);
4549  fhJetNEFInfo->Sumw2();
4550 
4551  fhClusterNEFInfo = new THnSparseF("fhClusterNEFInfo","Jet NEF Cluster Information Histogram",fnDimCluster,dimClusterBins,dimClusterLow,dimClusterUp);
4552  fhClusterNEFInfo->Sumw2();
4553 
4554  fNEFOutput->Add(fhJetNEFInfo);
4555  fNEFOutput->Add(fhClusterNEFInfo);
4556  }
4557  fNEFOutput->Add(fhJetNEFSignalInfo);
4558  fNEFOutput->Add(fhClusterNEFSignalInfo);
4559  fNEFOutput->Add(fhClusterShapeAll);
4560  fNEFOutput->Add(fhClusterPtCellAll);
4561  }
4562  fNEFOutput->Add(fhJetPtNEF);
4563  fOutput->Add(fNEFOutput);
4564  }
4565 
4566  // Add Histos & Profiles to List
4567  fOutput->Add(fh020Rho);
4568  fOutput->Add(fh80100Rho);
4569  fOutput->Add(fhRho);
4570  fOutput->Add(fhRhoCen);
4571  fOutput->Add(fh020BSPt);
4572  fOutput->Add(fh80100BSPt);
4573  fOutput->Add(fhBSPt);
4574  fOutput->Add(fhBSPtCen);
4575  fOutput->Add(fh020BSPtSignal);
4576  fOutput->Add(fh80100BSPtSignal);
4577  fOutput->Add(fhBSPtSignal);
4578  fOutput->Add(fhBSPtCenSignal);
4579  fOutput->Add(fh020DeltaPt);
4580  fOutput->Add(fh80100DeltaPt);
4581  fOutput->Add(fhDeltaPt);
4582  fOutput->Add(fhDeltaPtCen);
4583  fOutput->Add(fh020DeltaPtSignal);
4584  fOutput->Add(fh80100DeltaPtSignal);
4585  fOutput->Add(fhDeltaPtSignal);
4586  fOutput->Add(fhDeltaPtCenSignal);
4587  fOutput->Add(fh020DeltaPtNColl);
4588  fOutput->Add(fh80100DeltaPtNColl);
4589  fOutput->Add(fhDeltaPtNColl);
4590  fOutput->Add(fhDeltaPtCenNColl);
4591  fOutput->Add(fh020BckgFlucPt);
4592  fOutput->Add(fh80100BckgFlucPt);
4593  fOutput->Add(fhBckgFlucPt);
4594  fOutput->Add(fhBckgFlucPtCen);
4595  fOutput->Add(fpRho);
4596  fOutput->Add(fpLJetRho);
4597  fOutput->Add(fhJetPtArea);
4598  fOutput->Add(fhJetConstituentPt);
4599  fOutput->Add(fhJetTracksPt);
4600  fOutput->Add(fhJetClustersPt);
4601  fOutput->Add(fhJetConstituentCounts);
4602  fOutput->Add(fhJetTracksCounts);
4603  fOutput->Add(fhJetClustersCounts);
4604  fOutput->Add(fhJetPtZConstituent);
4605  fOutput->Add(fhJetPtZTrack);
4606  fOutput->Add(fhJetPtZCluster);
4607  fOutput->Add(fhJetPtZLeadingConstituent);
4608  fOutput->Add(fhJetPtZLeadingTrack);
4609  fOutput->Add(fhJetPtZLeadingCluster);
4610  fOutput->Add(fhEventCentralityVsZNA);
4611 }
4612 
4614 {
4615  fName = name;
4616 }
4617 
4619 {
4620  fCentralityTag = name.Data();
4621 }
4622 
4624 {
4625  fCentralityBins=bins;
4626  fCentralityLow=low;
4627  fCentralityUp=up;
4628 }
4629 
4630 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::SetPtRange(Int_t bins, Double_t low, Double_t up)
4631 {
4632  fPtBins=bins;
4633  fPtLow=low;
4634  fPtUp=up;
4635 }
4636 
4637 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::SetRhoPtRange(Int_t bins, Double_t low, Double_t up)
4638 {
4639  fRhoPtBins=bins;
4640  fRhoPtLow=low;
4641  fRhoPtUp=up;
4642 }
4643 
4645 {
4646  fDeltaPtBins=bins;
4647  fDeltaPtLow=low;
4648  fDeltaPtUp=up;
4649 }
4650 
4652 {
4653  fBckgFlucPtBins=bins;
4654  fBckgFlucPtLow=low;
4655  fBckgFlucPtUp=up;
4656 }
4657 
4659 {
4660  fLJetPtBins=bins;
4661  fLJetPtLow=low;
4662  fLJetPtUp=up;
4663 }
4664 
4666 {
4667  fLChargedTrackPtBins=bins;
4668  fLChargedTrackPtLow=low;
4669  fLChargedTrackPtUp=up;
4670 }
4671 
4672 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::SetNEFRange(Int_t bins, Double_t low, Double_t up)
4673 {
4674  fNEFBins=bins;
4675  fNEFLow=low;
4676  fNEFUp=up;
4677 }
4678 
4680 {
4681  fSignalTrackBias = chargedBias;
4682 }
4683 
4685 {
4686  fnDimJet = n;
4687 }
4688 
4690 {
4691  fnDimCluster = n;
4692 }
4693 
4695 {
4696  fRhoValue = value;
4697 }
4698 
4700 {
4701  fDoTHnSparse = doTHnSparse;
4702 }
4703 
4705 {
4706  fDo3DHistos = do3DPlotting;
4707 }
4708 
4709 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::FillRho(Double_t eventCentrality, Double_t rho)
4710 {
4711  SetRhoValue(rho);
4712  fhRho->Fill(rho);
4713  fhRhoCen->Fill(rho,eventCentrality);
4714  fpRho->Fill(eventCentrality,rho);
4715 
4716  if (eventCentrality<=20)
4717  {
4718  fh020Rho->Fill(rho);
4719  }
4720  else if (eventCentrality>=80)
4721  {
4722  fh80100Rho->Fill(rho);
4723  }
4724 }
4725 
4726 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::FillBSJS(Double_t eventCentrality, Double_t rho, Double_t signalCut, TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList)
4727 {
4728  Int_t i;
4729  Double_t tempPt=0.0;
4730  Double_t tempChargedHighPt=0.0;
4731 
4732  for (i=0;i<nIndexJetList;i++)
4733  {
4734  //AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
4735  AliEmcalJet *myJet;
4736 
4737  // Should test whether indexJetList[i]==i,
4738  if(indexJetList[i]!=i) cout<<"indexJetList[i]: "<<indexJetList[i]<<", i: "<<i<<endl;
4739  myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
4740  //myJet = (AliEmcalJet*) jetList->At(i);
4741  if(rho!=0){cout<<"Rho in CY: "<<rho<<endl;}
4742  tempPt=myJet->Pt()-rho*myJet->Area();
4743  tempChargedHighPt = myJet->MaxTrackPt();
4744  //ELI here are the background subtracted jets filled!
4745  fhBSPt->Fill(tempPt);
4746  fhBSPtCen->Fill(tempPt,eventCentrality);
4747  if (eventCentrality<=20)
4748  {
4749  fh020BSPt->Fill(tempPt);
4750  }
4751  else if (eventCentrality>=80)
4752  {
4753  fh80100BSPt->Fill(tempPt);
4754  }
4755  if (fSignalTrackBias==kTRUE)
4756  {
4757  if (tempChargedHighPt>=signalCut)
4758  {
4759  fhBSPtSignal->Fill(tempPt);
4760  fhBSPtCenSignal->Fill(tempPt,eventCentrality);
4761  if (eventCentrality<=20)
4762  {
4763  fh020BSPtSignal->Fill(tempPt);
4764  }
4765  else if (eventCentrality>=80)
4766  {
4767  fh80100BSPtSignal->Fill(tempPt);
4768  }
4769  }
4770  }
4771  else
4772  {
4773  if (tempPt>=signalCut)
4774  {
4775  fhBSPtSignal->Fill(tempPt);
4776  fhBSPtCenSignal->Fill(tempPt,eventCentrality);
4777  if (eventCentrality<=20)
4778  {
4779  fh020BSPtSignal->Fill(tempPt);
4780  }
4781  else if (eventCentrality>=80)
4782  {
4783  fh80100BSPtSignal->Fill(tempPt);
4784  }
4785  }
4786  }
4787  tempPt=0.0;
4788  tempChargedHighPt=0.0;
4789  }
4790 }
4791 
4792 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::FillDeltaPt(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
4793 {
4794  Int_t i;
4795  Double_t tempPt=0.0;
4796 
4797  for (i=0;i<nRC;i++)
4798  {
4799  tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
4800  fhDeltaPt->Fill(tempPt);
4801  fhDeltaPtCen->Fill(tempPt,eventCentrality);
4802  if (eventCentrality<=20)
4803  {
4804  fh020DeltaPt->Fill(tempPt);
4805  }
4806  else if (eventCentrality>=80)
4807  {
4808  fh80100DeltaPt->Fill(tempPt);
4809  }
4810  tempPt=0.0;
4811  }
4812 }
4813 
4814 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::FillDeltaPtSignal(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
4815 {
4816  Int_t i;
4817  Double_t tempPt=0.0;
4818 
4819  for (i=0;i<nRC;i++)
4820  {
4821  tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
4822  fhDeltaPtSignal->Fill(tempPt);
4823  fhDeltaPtCenSignal->Fill(tempPt,eventCentrality);
4824  if (eventCentrality<=20)
4825  {
4826  fh020DeltaPtSignal->Fill(tempPt);
4827  }
4828  else if (eventCentrality>=80)
4829  {
4830  fh80100DeltaPtSignal->Fill(tempPt);
4831  }
4832  tempPt=0.0;
4833  }
4834 }
4835 
4836 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::FillDeltaPtNColl(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
4837 {
4838  Int_t i;
4839  Double_t tempPt=0.0;
4840 
4841  for (i=0;i<nRC;i++)
4842  {
4843  tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
4844  fhDeltaPtNColl->Fill(tempPt);
4845  fhDeltaPtCenNColl->Fill(tempPt,eventCentrality);
4846  if (eventCentrality<=20)
4847  {
4848  fh020DeltaPtNColl->Fill(tempPt);
4849  }
4850  else if (eventCentrality>=80)
4851  {
4852  fh80100DeltaPtNColl->Fill(tempPt);
4853  }
4854  tempPt=0.0;
4855  }
4856 }
4857 
4858 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::FillBackgroundFluctuations(Double_t eventCentrality, Double_t rho, Double_t jetRadius)
4859 {
4860  Double_t tempPt=0.0;
4861 
4862  tempPt=rho*TMath::Power(jetRadius,2);
4863  fhBckgFlucPt->Fill(tempPt);
4864  fhBckgFlucPtCen->Fill(tempPt,eventCentrality);
4865  if (eventCentrality<=20)
4866  {
4867  fh020BckgFlucPt->Fill(tempPt);
4868  }
4869  else if (eventCentrality>=80)
4870  {
4871  fh80100BckgFlucPt->Fill(tempPt);
4872  }
4873 }
4874 
4876 {
4877  fpLJetRho->Fill(jetPt,rho);
4878 }
4879 
4881 {
4882  fDoNEFQAPlots = doNEFAna;
4883 }
4884 
4886 {
4887  fDoNEFSignalOnly = doNEFSignalOnly;
4888 }
4889 
4890 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::DoNEFAnalysis(Double_t nefCut, Double_t signalCut, TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList, TObjArray *clusterList, TClonesArray *orgClusterList, AliVEvent *event, AliEMCALGeometry *geometry, AliEMCALRecoUtils *recoUtils, AliVCaloCells *cells)
4891 {
4892  if (fDoNEFQAPlots==kFALSE)
4893  {
4894  return;
4895  }
4896 
4897  if (fDoTHnSparse == kFALSE)
4898  {
4899  Int_t i;
4900 
4901  Double_t nef=0.0;
4902  Double_t jetPt=0.0;
4903  Double_t tempChargedHighPt=0.0;
4904 
4905  for (i=0;i<nIndexJetList;i++)
4906  {
4907  AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
4908  tempChargedHighPt = myJet->MaxTrackPt();
4909  nef=myJet->NEF();
4910  jetPt=myJet->Pt();
4911 
4912  if (fSignalTrackBias==kTRUE)
4913  {
4914  if (tempChargedHighPt>=signalCut && nef<=nefCut)
4915  {
4916  fhJetPtNEF->Fill(jetPt,nef);
4917  }
4918  }
4919  else
4920  {
4921  if (jetPt>=signalCut && nef<=nefCut)
4922  {
4923  fhJetPtNEF->Fill(jetPt,nef);
4924  }
4925  }
4926 
4927  nef=0.0;
4928  jetPt=0.0;
4929  tempChargedHighPt=0.0;
4930  }
4931  return;
4932  }
4933 
4934  Int_t i,j,k;
4935 
4936  Double_t valJet[fnDimJet];
4937  Double_t valCluster[fnDimJet];
4938  for (i=0;i<fnDimJet;i++)
4939  {
4940  valJet[i]=0.0;
4941  valCluster[i]=0.0;
4942  }
4943 
4944  Double_t nef=0.0;
4945  Double_t jetPt=0.0;
4946  Double_t tempChargedHighPt=0.0;
4947  Double_t eta=0.0;
4948  Double_t phi=0.0;
4949  Int_t totalMult=0;
4950  Int_t chargedMult=0;
4951  Int_t neutralMult=0;
4952  Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
4953  Bool_t shared = kFALSE;
4954 
4955  Double_t zLeading=0.0;
4956  Double_t zLeadingTrack=0.0;
4957  Double_t ECluster=0.0;
4958  Double_t eSeed=0.0;
4959  Double_t tCell=0.0;
4960  Double_t eCross=0.0;
4961  Double_t FCross=0.0;
4962 
4963  Double_t lowTime=9.99e99;
4964  Double_t upTime=-9.99e99;
4965  Int_t tempCellID=0;
4966  Double_t tempCellTime=0.0;
4967 
4968  Double_t event_centrality = event->GetCentrality()->GetCentralityPercentile(fCentralityTag);
4969  valJet[4] = valCluster[4] = event_centrality;
4970 
4971  // First, do Jet QA
4972  for (i=0;i<nIndexJetList;i++)
4973  {
4974  AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
4975  tempChargedHighPt = myJet->MaxTrackPt();
4976  nef=myJet->NEF();
4977  jetPt=myJet->Pt();
4978  eta=myJet->Eta();
4979  phi=myJet->Phi();
4980  totalMult=myJet->GetNumberOfConstituents();
4981  chargedMult = myJet->GetNumberOfTracks();
4982  neutralMult=myJet->GetNumberOfClusters();
4983  zLeading=myJet->MaxPartPt()/myJet->Pt();
4984  zLeadingTrack=myJet->MaxTrackPt()/myJet->Pt();
4985 
4986  valJet[0] = valCluster[0] = nef;
4987  valJet[1] = valCluster[1] = jetPt;
4988  valJet[2] = valCluster[2] = eta;
4989  valJet[3] = valCluster[3] = phi;
4990  valJet[5] = valCluster[5] = zLeading;
4991  valJet[9] = zLeadingTrack;
4992 
4993  valJet[6] = totalMult;
4994  valJet[7] = chargedMult;
4995  valJet[8] = neutralMult;
4996 
4997  // Supress filling of this histogram due to memory size of THnSparse when running over large datasets
4998  if (fDoNEFSignalOnly == kFALSE)
4999  {
5000  fhJetNEFInfo->Fill(valJet);
5001  }
5002 
5003  if (fSignalTrackBias==kTRUE)
5004  {
5005  if (tempChargedHighPt>=signalCut && nef<=nefCut)
5006  {
5007  fhJetNEFSignalInfo->Fill(valJet);
5008  fhJetPtNEF->Fill(jetPt,nef);
5009  }
5010  }
5011  else
5012  {
5013  if (jetPt>=signalCut && nef<=nefCut)
5014  {
5015  fhJetNEFSignalInfo->Fill(valJet);
5016  fhJetPtNEF->Fill(jetPt,nef);
5017  }
5018  }
5019 
5020  for (j=0;j<neutralMult;j++)
5021  {
5022  AliVCluster* vcluster = (AliVCluster*) orgClusterList->At(myJet->ClusterAt(j));
5023  ECluster = vcluster->E();
5024  recoUtils->GetMaxEnergyCell(geometry,cells,vcluster,absId,iSupMod,ieta,iphi,shared);
5025  eSeed = cells->GetCellAmplitude(absId);
5026  tCell = cells->GetCellTime(absId);
5027  eCross = recoUtils->GetECross(absId,tCell,cells,event->GetBunchCrossNumber());
5028  FCross = 1 - eCross/eSeed;
5029 
5030  // Obtain Delta t of Cluster
5031  lowTime=9.99e99;
5032  upTime=-9.99e99;
5033  for (k=0;k<vcluster->GetNCells();k++)
5034  {
5035  tempCellID=vcluster->GetCellAbsId(k);
5036  tempCellTime=cells->GetCellTime(tempCellID);
5037  if (tempCellTime>upTime)
5038  {
5039  upTime=tempCellTime;
5040  }
5041  if (tempCellTime<lowTime)
5042  {
5043  lowTime=tempCellTime;
5044  }
5045  }
5046  valCluster[6] = ECluster;
5047  valCluster[7] = FCross;
5048  valCluster[8] = tCell;
5049  valCluster[9] = upTime-lowTime;
5050  valCluster[10] = vcluster->GetNCells();
5051 
5052  if (fDoNEFSignalOnly == kFALSE)
5053  {
5054  fhClusterNEFInfo->Fill(valCluster);
5055  }
5056 
5057  if (fSignalTrackBias==kTRUE)
5058  {
5059  if (tempChargedHighPt>=signalCut && nef<=nefCut)
5060  {
5061  fhClusterNEFSignalInfo->Fill(valCluster);
5062  }
5063  }
5064  else
5065  {
5066  if (myJet->Pt()>=signalCut && nef<=nefCut)
5067  {
5068  fhClusterNEFSignalInfo->Fill(valCluster);
5069  }
5070  }
5071  tempCellID=0;
5072  tempCellTime=0.0;
5073  eSeed=0.0;
5074  tCell=0.0;
5075  eCross=0.0;
5076  FCross=0.0;
5077  iSupMod=-1,absId=-1,ieta=-1,iphi=-1;
5078  }
5079 
5080  nef=0.0;
5081  jetPt=0.0;
5082  tempChargedHighPt=0.0;
5083  eta=0.0;
5084  phi=0.0;
5085  totalMult=0;
5086  chargedMult=0;
5087  neutralMult=0;
5088  zLeading=0.0;
5089  zLeadingTrack=0.0;
5090  ECluster=0.0;
5091  }
5092 
5093  // Now do Cluster QA
5094  for (i=0;i<clusterList->GetEntries();i++)
5095  {
5096  AliVCluster *vcluster = (AliVCluster*) clusterList->At(i);
5097  fhClusterShapeAll->Fill(vcluster->GetNCells());
5098  fhClusterPtCellAll->Fill(vcluster->E(),vcluster->GetNCells());
5099  }
5100 }
5101 
5102 void AliAnalysisTaskFullpAJets_Eli_Mod::AlipAJetHistos::FillMiscJetStats(TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList, TClonesArray *trackList, TClonesArray *clusterList, Double_t *vertex)
5103 {
5104  Int_t i,j;
5105 
5106  TLorentzVector *cluster_vec = new TLorentzVector;
5107  for (i=0;i<nIndexJetList;i++)
5108  {
5109  AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
5110 
5111  if (fDo3DHistos == kTRUE)
5112  {
5113  fhJetPtEtaPhi->Fill(myJet->Pt(),myJet->Eta(),myJet->Phi());
5114  }
5115  fhJetPtArea->Fill(myJet->Pt(),myJet->Area());
5116  fhJetConstituentCounts->Fill(myJet->Pt(),myJet->GetNumberOfConstituents());
5117  fhJetTracksCounts->Fill(myJet->Pt(),myJet->GetNumberOfTracks());
5118  fhJetClustersCounts->Fill(myJet->Pt(),myJet->GetNumberOfClusters());
5119  fhJetPtZLeadingConstituent->Fill(myJet->Pt(),myJet->MaxPartPt()/myJet->Pt());
5120  fhJetPtZLeadingTrack->Fill(myJet->Pt(),myJet->MaxTrackPt()/myJet->Pt());
5121  fhJetPtZLeadingCluster->Fill(myJet->Pt(),myJet->MaxClusterPt()/myJet->Pt());
5122  for (j=0;j<myJet->GetNumberOfTracks();j++)
5123  {
5124  AliVTrack *vtrack = (AliVTrack*) myJet->TrackAt(j,trackList);
5125  fhJetConstituentPt->Fill(myJet->Pt(),vtrack->Pt());
5126  fhJetTracksPt->Fill(myJet->Pt(),vtrack->Pt());
5127  fhJetPtZTrack->Fill(myJet->Pt(),vtrack->Pt()/myJet->Pt());
5128  fhJetPtZConstituent->Fill(myJet->Pt(),vtrack->Pt()/myJet->Pt());
5129  }
5130  for (j=0;j<myJet->GetNumberOfClusters();j++)
5131  {
5132  AliVCluster *vcluster = (AliVCluster*) myJet->ClusterAt(j,clusterList);
5133  vcluster->GetMomentum(*cluster_vec,vertex);
5134  fhJetConstituentPt->Fill(myJet->Pt(),cluster_vec->Pt());
5135  fhJetClustersPt->Fill(myJet->Pt(),vcluster->E());
5136  fhJetPtZCluster->Fill(myJet->Pt(),cluster_vec->Pt()/myJet->Pt());
5137  fhJetPtZConstituent->Fill(myJet->Pt(),cluster_vec->Pt()/myJet->Pt());
5138  }
5139  }
5140  delete cluster_vec;
5141 }
5142 
5144 {
5145  Double_t event_centrality = event->GetCentrality()->GetCentralityPercentile(fCentralityTag);
5146  Double_t event_centrality_ZNA = event->GetCentrality()->GetCentralityPercentile("ZNA");
5147 
5148  fhEventCentralityVsZNA->Fill(event_centrality,event_centrality_ZNA);
5149  if (fDo3DHistos == kTRUE)
5150  {
5151  fhEventCentralityVsZNAPt->Fill(event_centrality,event_centrality_ZNA,leadingJetPt);
5152  }
5153 }
5154 
5156 {
5157  return fOutput;
5158 }
5159 
5161 {
5162  return fRhoValue;
5163 }
5164 
Double_t * fTPCRCBckgFlucNColl
Stores the pT of RC Background clusters in EMCal with no spatial restrictions.
Double_t Area() const
Definition: AliEmcalJet.h:97