AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskFullppJet.cxx
Go to the documentation of this file.
1 // **************************************
2 // Full pp jet task - ESD input only
3 // Extract the jet spectrum and all the
4 // systematic uncertainties
5 // -R. Ma, Mar 2011
6 // **************************************
7 
8 #include <TCanvas.h>
9 #include <TChain.h>
10 #include <TFormula.h>
11 #include <TH1.h>
12 #include <TH2.h>
13 #include <TH3.h>
14 #include <TProfile2D.h>
15 #include <THnSparse.h>
16 #include <TROOT.h>
17 #include <TTree.h>
18 #include <TArrayI.h>
19 #include <TClonesArray.h>
20 #include <TRandom3.h>
21 #include <TFile.h>
22 #include <TF1.h>
23 
24 #include "AliAODEvent.h"
25 #include "AliAODInputHandler.h"
26 #include "AliAnalysisManager.h"
27 #include "AliAnalysisTask.h"
28 #include "AliCentrality.h"
30 #include "AliESDEvent.h"
31 #include "AliESDInputHandler.h"
32 #include "AliESDtrackCuts.h"
33 #include "AliVParticle.h"
34 #include "AliInputEventHandler.h"
35 #include "AliEMCALGeometry.h"
36 #include "AliEMCALRecoUtils.h"
37 #include "TGeoManager.h"
38 #include "AliMCEvent.h"
39 #include "AliMCParticle.h"
40 #include "AliStack.h"
41 #include "AliGenEventHeader.h"
42 #include "AliGenPythiaEventHeader.h"
43 #include "TGeoGlobalMagField.h"
44 
45 #include "AliAODJet.h"
46 #include "AliFJWrapper.h"
47 
48 #include <iostream>
49 using std::cout;
50 using std::cerr;
51 using std::endl;
52 using std::vector;
53 
55 
56 const Float_t kRadius[3] = {0.4,0.2,0.3};
57 const Double_t kPI = TMath::Pi();
58 const Double_t kdRCut[3] = {0.25,0.1,0.15};
59 
60 //________________________________________________________________________
62  AliAnalysisTaskSE("default"),
63  fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE),
64  fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
65  fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
66  fCentrality(99), fZVtxMax(10),
67  fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
68  fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0),
69  fGeom(0x0), fRecoUtil(0x0),
70  fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0),fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9),
71  fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2),
72  fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE),
73  fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE),
74  fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
75  fJetNEFMin(0.01), fJetNEFMax(0.98),
76  fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
77  fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0),
78  fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0),
79  fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
80  fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
81  fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
82  fSysClusterizer(0),
83  fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
84  fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
85  fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
86 {
87  // Constructor
88 
89  // Define input and output slots here
90  // Input slot #0 works with a TChain
91  DefineInput(0, TChain::Class());
92  DefineOutput(1, TList::Class());
93  // Output slot #0 id reserved by the base class for AOD
94  for(Int_t i=0; i<2; i++)
95  {
96  fTrkPtMin[i] = 0.15;
97  fTrkPtMax[i] = 200;
98  fClsEtMin[i] = 0.15;
99  fClsEtMax[i] = 200;
100 
101  fVertexGenZ[i] = 0x0;
102  fEventZ[i] = 0x0;
103  fhNTrials[i] = 0x0;
104  fhNMatchedTrack[i] = 0x0;
105  fhClsE[i] = 0x0;
106  for(Int_t k=0; k<2; k++)
107  {
108  fhSysClusterE[i][k] = 0x0;
109  fhSysNCellVsClsE[i][k] = 0x0;
110  }
111  }
112 
113  for(Int_t r=0; r<kNRadii; r++)
114  {
115  for(Int_t j=0; j<5; j++)
116  {
117  fHCOverSubE[r][j] = 0x0;
118  fHCOverSubEFrac[r][j] = 0x0;
119  for(Int_t k=0; k<2; k++)
120  {
121  fhSubClsEVsJetPt[k][r][j] = 0x0;
122  fhHCTrkPtClean[k][r][j] = 0x0;
123  fhHCTrkPtAmbig[k][r][j] = 0x0;
124  }
125 
126  if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
127 
128  if(j<3)
129  {
130  fJetCount[j][r] = 0x0;
131  fhNeutralPtInJet[j][r] = 0x0;
132  fhTrigNeuPtInJet[j][r] = 0x0;
133  fhChargedPtInJet[j][r] = 0x0;
134  fhLeadNePtInJet[j][r] = 0x0;
135  fhLeadChPtInJet[j][r] = 0x0;
136  fJetEnergyFraction[j][r] = 0x0;
137  fJetNPartFraction[j][r] = 0x0;
138  }
139  if(j<2)
140  {
141  fRelTrkCon[j][r] = 0x0;
142  fhFcrossVsZleading[j][r] = 0x0;
143  fhChLeadZVsJetPt[j][r] = 0x0;
144  fhJetPtWithTrkThres[j][r]= 0x0;
145  fhJetPtWithClsThres[j][r]= 0x0;
146  for(Int_t k=0; k<2; k++)
147  {
148  fhJetPtVsLowPtCons[j][r][k] = 0x0;
149  }
150  fhJetPtInExoticEvent[j][r] = 0x0;
151  fhJetInTPCOnlyVtx[j][r] = 0x0;
152  fhCorrTrkEffPtBin[j][r] = 0x0;
153  for(Int_t i=0; i<kNBins; i++)
154  {
155  fhCorrTrkEffSample[j][r][i] = 0x0;
156  }
157  }
158 
159  if(r==0 && j<3)
160  {
161  for(Int_t k=0; k<2; k++)
162  {
163  for(Int_t l=0; l<2; l++)
164  {
165  fhUEJetPtNorm[j][k][l] = 0x0;
166  fhUEJetPtVsSumPt[j][k][l] = 0x0;
167  fhUEJetPtVsConsPt[j][k][l] = 0x0;
168  }
169  }
170  }
171  }
172  for(Int_t i=0; i<2; i++)
173  {
174  fhNKFracVsJetPt[i][r] = 0x0;
175  fhWeakFracVsJetPt[i][r] = 0x0;
176  fhJetResponseNK[i][r] = 0x0;
177  fhJetResponseWP[i][r] = 0x0;
178  fhJetResolutionNK[i][r] = 0x0;
179  fhJetResolutionWP[i][r] = 0x0;
180  fhJetResponseNKSM[i][r] = 0x0;
181  fhJetResponseWPSM[i][r] = 0x0;
182  fhJetResolutionNKSM[i][r] = 0x0;
183  fhJetResolutionWPSM[i][r] = 0x0;
184  }
185  }
186 
187  for(Int_t s=0; s<3; s++)
188  {
189  for(Int_t a=0; a<2; a++)
190  {
191  for(Int_t r=0; r<kNRadii; r++)
192  {
193  fDetJetFinder[s][a][r] = 0x0;
194  fJetTCA[s][a][r] = 0x0;
195  if(s==0 && a==0)
196  {
197  fTrueJetFinder[r] = 0x0;
198  fMcTruthAntikt[r] = 0x0;
199  }
200  }
201  }
202  }
203 
204  for(Int_t j=0; j<3; j++)
205  {
206  fTrkEffFunc[j] = 0x0;
207  fhSecondaryResponse[j] = 0x0;
208  }
209  for(Int_t i=0; i<10; i++)
210  {
211  fTriggerCurve[i] = 0x0;
212  fTriggerEfficiency[i] = 0x0;
213  }
214 
215 }
216 
217 //________________________________________________________________________
219  AliAnalysisTaskSE(name),
220  fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE),
221  fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
222  fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
223  fCentrality(99), fZVtxMax(10),
224  fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
225  fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0),
226  fGeom(0x0), fRecoUtil(0x0),
227  fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0), fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9),
228  fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2),
229  fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE),
230  fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE),
231  fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
232  fJetNEFMin(0.01), fJetNEFMax(0.98),
233  fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
234  fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0),
235  fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0),
236  fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
237  fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
238  fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
239  fSysClusterizer(0),
240  fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
241  fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
242  fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
243 {
244  // Constructor
245 
246  // Define input and output slots here
247  // Input slot #0 works with a TChain
248  DefineInput(0, TChain::Class());
249  DefineOutput(1, TList::Class());
250  // Output slot #0 id reserved by the base class for AOD
251  for(Int_t i=0; i<2; i++)
252  {
253  fTrkPtMin[i] = 0.15;
254  fTrkPtMax[i] = 200;
255  fClsEtMin[i] = 0.15;
256  fClsEtMax[i] = 200;
257 
258  fVertexGenZ[i] = 0x0;
259  fEventZ[i] = 0x0;
260  fhNTrials[i] = 0x0;
261  fhNMatchedTrack[i] = 0x0;
262  fhClsE[i] = 0x0;
263  for(Int_t k=0; k<2; k++)
264  {
265  fhSysClusterE[i][k] = 0x0;
266  fhSysNCellVsClsE[i][k] = 0x0;
267  }
268  }
269 
270  for(Int_t r=0; r<kNRadii; r++)
271  {
272  for(Int_t j=0; j<5; j++)
273  {
274  fHCOverSubE[r][j] = 0x0;
275  fHCOverSubEFrac[r][j] = 0x0;
276  for(Int_t k=0; k<2; k++)
277  {
278  fhSubClsEVsJetPt[k][r][j] = 0x0;
279  fhHCTrkPtClean[k][r][j] = 0x0;
280  fhHCTrkPtAmbig[k][r][j] = 0x0;
281  }
282 
283  if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
284 
285  if(j<3)
286  {
287  fJetCount[j][r] = 0x0;
288  fhNeutralPtInJet[j][r] = 0x0;
289  fhTrigNeuPtInJet[j][r] = 0x0;
290  fhChargedPtInJet[j][r] = 0x0;
291  fhLeadNePtInJet[j][r] = 0x0;
292  fhLeadChPtInJet[j][r] = 0x0;
293  fJetEnergyFraction[j][r] = 0x0;
294  fJetNPartFraction[j][r] = 0x0;
295  }
296  if(j<2)
297  {
298  fRelTrkCon[j][r] = 0x0;
299  fhFcrossVsZleading[j][r] = 0x0;
300  fhChLeadZVsJetPt[j][r] = 0x0;
301  fhJetPtWithTrkThres[j][r]= 0x0;
302  fhJetPtWithClsThres[j][r]= 0x0;
303  for(Int_t k=0; k<2; k++)
304  {
305  fhJetPtVsLowPtCons[j][r][k] = 0x0;
306  }
307  fhJetPtInExoticEvent[j][r] = 0x0;
308  fhJetInTPCOnlyVtx[j][r] = 0x0;
309  fhCorrTrkEffPtBin[j][r] = 0x0;
310  for(Int_t i=0; i<kNBins; i++)
311  {
312  fhCorrTrkEffSample[j][r][i] = 0x0;
313  }
314  }
315 
316  if(r==0 && j<3)
317  {
318  for(Int_t k=0; k<2; k++)
319  {
320  for(Int_t l=0; l<2; l++)
321  {
322  fhUEJetPtNorm[j][k][l] = 0x0;
323  fhUEJetPtVsSumPt[j][k][l] = 0x0;
324  fhUEJetPtVsConsPt[j][k][l] = 0x0;
325  }
326  }
327  }
328  }
329  for(Int_t i=0; i<2; i++)
330  {
331  fhNKFracVsJetPt[i][r] = 0x0;
332  fhWeakFracVsJetPt[i][r] = 0x0;
333  fhJetResponseNK[i][r] = 0x0;
334  fhJetResponseWP[i][r] = 0x0;
335  fhJetResolutionNK[i][r] = 0x0;
336  fhJetResolutionWP[i][r] = 0x0;
337  fhJetResponseNKSM[i][r] = 0x0;
338  fhJetResponseWPSM[i][r] = 0x0;
339  fhJetResolutionNKSM[i][r] = 0x0;
340  fhJetResolutionWPSM[i][r] = 0x0;
341  }
342  }
343 
344  for(Int_t s=0; s<3; s++)
345  {
346  for(Int_t a=0; a<2; a++)
347  {
348  for(Int_t r=0; r<kNRadii; r++)
349  {
350  fDetJetFinder[s][a][r] = 0x0;
351  fJetTCA[s][a][r] = 0x0;
352  if(s==0 && a==0)
353  {
354  fTrueJetFinder[r] = 0x0;
355  fMcTruthAntikt[r] = 0x0;
356  }
357  }
358  }
359  }
360 
361  for(Int_t j=0; j<3; j++)
362  {
363  fTrkEffFunc[j] = 0x0;
364  fhSecondaryResponse[j] = 0x0;
365  }
366  for(Int_t i=0; i<10; i++)
367  {
368  fTriggerCurve[i] = 0x0;
369  fTriggerEfficiency[i] = 0x0;
370  }
371 }
372 
373 //________________________________________________________________________
375 {
376  //Destructor
377 
378  if(fEsdTrackCuts) delete fEsdTrackCuts;
381  if(fOutputList)
382  { fOutputList->Delete(); delete fOutputList;}
383  for(Int_t s=0; s<3; s++)
384  {
385  for(Int_t a=0; a<2; a++)
386  {
387  for(Int_t r=0; r<kNRadii; r++)
388  {
389  if(fDetJetFinder[s][a][r]) delete fDetJetFinder[s][a][r];
390  if(fJetTCA[s][a][r]) { fJetTCA[s][a][r]->Delete(); delete fJetTCA[s][a][r]; }
391  if(s==0 && a==0)
392  {
393  if(fTrueJetFinder[r]) delete fTrueJetFinder[r];
394  if(fMcTruthAntikt[r]) { fMcTruthAntikt[r]->Delete(); delete fMcTruthAntikt[r]; }
395  }
396  }
397  }
398  }
399  if(fRandomGen) delete fRandomGen;
400  for(Int_t i=0; i<3; i++)
401  {
402  if(fTrkEffFunc[i]) delete fTrkEffFunc[i];
403  if(fhSecondaryResponse[i]) delete fhSecondaryResponse[i];
404  }
405  for(Int_t i=0; i<10; i++)
406  {
407  if(fTriggerEfficiency[i]) delete fTriggerEfficiency[i];
408  if(fTriggerCurve[i]) delete fTriggerCurve[i];
409  }
410  for(Int_t r=0; r<kNRadii; r++)
411  {
412  for(Int_t j=0; j<2; j++)
413  {
414  if(fhCorrTrkEffPtBin[j][r]) delete fhCorrTrkEffPtBin[j][r];
415  for(Int_t i=0; i<kNBins; i++)
416  {
417  if(fhCorrTrkEffSample[j][r][i]) delete fhCorrTrkEffSample[j][r][i];
418  }
419  }
420  }
421  if(fTrackArray) { fTrackArray->Delete(); delete fTrackArray; }
422  if(fClusterArray) { fClusterArray->Delete(); delete fClusterArray; }
423  if(fMcPartArray) { fMcPartArray->Reset(); delete fMcPartArray; }
424 
425  if(fRecoUtil) delete fRecoUtil;
427  if(fNonLinear) delete fNonLinear;
428  if(fTrkPtResData) delete fTrkPtResData;
429 }
430 
431 //________________________________________________________________________
433 {
434  //
435  // Fill the number of tracks per chunk
436  //
437 
439  fNTracksPerChunk = 0;
440  fEDSFileCounter++;
441  return kTRUE;
442 }
443 
444 //________________________________________________________________________
446 {
447  // Create histograms
448  // Called once
449  //
450 
451  if(fRunUE) fFindChargedOnlyJet = kTRUE;
452 
453  const Int_t nTrkPtBins = 100;
454  const Float_t lowTrkPtBin=0, upTrkPtBin=100.;
455 
456  const Int_t nbins = 220;
457  Double_t xbins[221];
458  for(Int_t i=0; i<nbins+1; i++)
459  xbins[i] = i;
460 
461  OpenFile(1);
462  fOutputList = new TList();
463  fOutputList->SetOwner(1);
464 
465  fhJetEventStat = new TH1F("fhJetEventStat","Event statistics for jet analysis",12,0,12);
466  fhJetEventStat->GetXaxis()->SetBinLabel(1,"ALL");
467  fhJetEventStat->GetXaxis()->SetBinLabel(2,"MB");
468  fhJetEventStat->GetXaxis()->SetBinLabel(3,"MB+vtx+10cm");
469  fhJetEventStat->GetXaxis()->SetBinLabel(4,"EMC");
470  fhJetEventStat->GetXaxis()->SetBinLabel(5,"EMC+vtx+10cm");
471  fhJetEventStat->GetXaxis()->SetBinLabel(6,"MB+vtx");
472  fhJetEventStat->GetXaxis()->SetBinLabel(7,"EMC+vtx");
473  fhJetEventStat->GetXaxis()->SetBinLabel(8,"TriggerBit");
474  fhJetEventStat->GetXaxis()->SetBinLabel(9,"LED");
475  fhJetEventStat->GetXaxis()->SetBinLabel(10,"MB+TVtx+10cm");
476  fhJetEventStat->GetXaxis()->SetBinLabel(11,"EMC+TVtx+10cm");
477  fhJetEventStat->GetXaxis()->SetBinLabel(12,"ALL-Pileup");
479 
480  fhJetEventCount = new TH1F("fhJetEventCount","Event statistics for jet analysis",12,0,12);
481  fhJetEventCount->GetXaxis()->SetBinLabel(1,"ALL");
482  fhJetEventCount->GetXaxis()->SetBinLabel(2,"MB");
483  fhJetEventCount->GetXaxis()->SetBinLabel(3,"MB-pileup");
484  fhJetEventCount->GetXaxis()->SetBinLabel(4,"MB+Vtx");
485  fhJetEventCount->GetXaxis()->SetBinLabel(5,"MB+Vtx+10cm");
486  fhJetEventCount->GetXaxis()->SetBinLabel(6,"EMC");
487  fhJetEventCount->GetXaxis()->SetBinLabel(7,"EMC-pileup");
488  fhJetEventCount->GetXaxis()->SetBinLabel(8,"EMC+Vtx");
489  fhJetEventCount->GetXaxis()->SetBinLabel(9,"EMC+Vtx+10cm");
490  fhJetEventCount->GetXaxis()->SetBinLabel(10,"Good EMC");
492 
493  if(fCheckTPCOnlyVtx)
494  {
495  fhEventStatTPCVtx = new TH1F("fhEventStatTPCVtx","Event statistics for TPC only vertex",9,0,9);
496  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(1,"FastOnly");
497  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(2,"FastOnly+PVtx");
498  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(3,"FastOnly+TVtx");
499  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(4,"MB");
500  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(5,"MB+PVtx");
501  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(6,"MB+TVtx");
502  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(7,"EMC");
503  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(8,"EMC+PVtx");
504  fhEventStatTPCVtx->GetXaxis()->SetBinLabel(9,"EMC+TVtx");
506  }
507 
508  fhChunkQA = new TH1F("fhChunkQA","# of hybrid tracks per chunk",200,0,200);
509  fOutputList->Add(fhChunkQA);
510 
511  const Int_t dim1 = 3;
512  Int_t nBins1[dim1] = {200,50,110};
513  Double_t lowBin1[dim1] = {0,0,0,};
514  Double_t upBin1[dim1] = {200,0.5,1.1};
515 
516  const Int_t dim2 = 3;
517  Int_t nBins2[dim2] = {200,50,50};
518  Double_t lowBin2[dim2] = {0,0,0,};
519  Double_t upBin2[dim2] = {200,0.5,50};
520 
521  const char* triggerName[3] = {"MB","EMC","MC"};
522  const char* triggerTitle[3] = {"MB","EMCal-trigger","MC true"};
523  const char* fraction[5] = {"MIP","30","50","70","100"};
524  const char* exotic[2] = {"3GeV","5GeV"};
525  const char* vertexType[2] = {"All MB","MB with vertex"};
526  const char *vertexName[2] = {"All","Vertex"};
527  const char *clusterizerName[2] = {"before","after"};
528  const char *UEName[2] = {"charged","charged_neutral"};
529  const char *UETitle[2] = {"charged","charged+neutral"};
530  const char *UEEventName[2] = {"LeadingJet","Back-To-Back"};
531 
532  if(fIsMC)
533  {
534  for(Int_t i=0; i<2; i++)
535  {
536  fhNTrials[i] = new TH1F(Form("MC_%s_fhNTrials",triggerName[i]),Form("MC-%s: # of trials",triggerName[i]),1,0,1);
537  fOutputList->Add(fhNTrials[i]);
538 
539  fVertexGenZ[i] = new TH1F(Form("%s_fVertexGenZ",vertexName[i]),Form("Distribution of vertex z (%s); z (cm)",vertexType[i]),60,-30,30);
540  fOutputList->Add(fVertexGenZ[i]);
541  }
542  }
543 
544  for(Int_t i=0; i<3; i++)
545  {
546  if(!fIsMC && i==2) continue;
547 
548  if(fSaveQAHistos)
549  {
550  if(i<2)
551  {
552  fEventZ[i] = new TH1F(Form("%s_fEventZ",triggerName[i]),Form("%s: Distribution of vertex z; z (cm)",triggerTitle[i]),60,-30,30);
553  fOutputList->Add(fEventZ[i]);
554  }
555 
556  for(Int_t r=0; r<kNRadii; r++)
557  {
558  if(!fRadius.Contains("0.4") && r==0) continue;
559  if(!fRadius.Contains("0.2") && r==1) continue;
560  if(!fRadius.Contains("0.3") && r==2) continue;
561 
562  fJetCount[i][r] = new TH1F(Form("%s_fJetCount_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins);
563  fOutputList->Add(fJetCount[i][r]);
564 
565  fhNeutralPtInJet[i][r] = new TH2F(Form("%s_fhNeutralPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of neutral constituents vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin);
566  fOutputList->Add(fhNeutralPtInJet[i][r]);
567 
568  fhTrigNeuPtInJet[i][r] = new TH2F(Form("%s_fhTrigNeuPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of triggered neutral constituents vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin);
569  fOutputList->Add(fhTrigNeuPtInJet[i][r]);
570 
571  fhChargedPtInJet[i][r] = new TH2F(Form("%s_fhChargedPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of charged constituents vs jet p_{T} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin);
572  fOutputList->Add(fhChargedPtInJet[i][r]);
573 
574  fhLeadNePtInJet[i][r] = new TH2F(Form("%s_fhLeadNePtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of leading neutral constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T}^{ne} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,100);
575  fOutputList->Add(fhLeadNePtInJet[i][r]);
576 
577  fhLeadChPtInJet[i][r] = new TH2F(Form("%s_fhLeadChPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of leading charged constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T}^{ch} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,100);
578  fOutputList->Add(fhLeadChPtInJet[i][r]);
579 
580  fhJetPtVsZ[i][r] = new TH3F(Form("%s_fhJetPtVsZ_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs Z_{h} vs constituent type in EMCal (R=%1.1f);p_{T}^{jet} (GeV/c); Z_{h};constituent type",triggerTitle[i],kRadius[r]),200,0,200,110,-0.05,1.05,4,0,4);
581  fOutputList->Add(fhJetPtVsZ[i][r]);
582 
583  fJetEnergyFraction[i][r] = new THnSparseF(Form("%s_fJetEnergyFraction_%1.1f",triggerName[i],kRadius[r]),Form("%s: Jet p_{T} vs radius vs energy fraction in EMCal (R=%1.1f,Z<0.98); p_{T};R;Fraction",triggerName[i],kRadius[r]),dim1,nBins1,lowBin1,upBin1);
584  fOutputList->Add(fJetEnergyFraction[i][r]);
585 
586  fJetNPartFraction[i][r] = new THnSparseF(Form("%s_fJetNPartFraction_%1.1f",triggerName[i],kRadius[r]),Form("%s: Jet p_{T} vs radius vs NPart in EMCal (R=%1.1f,Z<0.98); p_{T};R;NPart",triggerName[i],kRadius[r]),dim2,nBins2,lowBin2,upBin2);
587  fOutputList->Add(fJetNPartFraction[i][r]);
588 
589  if(i<2)
590  {
591  fhJetPtInExoticEvent[i][r] = new TH1F(Form("EMC_fhJetPtInExoticEvent_%1.1f_%s",kRadius[r],exotic[i]),Form("EMC: jet p_{T} in events with exotic cluster > %s (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c)",exotic[i],kRadius[r]),nbins,xbins);
592  fOutputList->Add(fhJetPtInExoticEvent[i][r]);
593 
594  fRelTrkCon[i][r] = new TH3F(Form("%s_fRelTrkCon_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs (sum p_{T}^{ch})/p_{T}^{jet} vs track class in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); (sum p_{T}^{ch})/p_{T}^{jet}; track class",triggerTitle[i],kRadius[r]),200,0,200,110,-0.05,1.05,3,0,3);
595  fOutputList->Add(fRelTrkCon[i][r]);
596 
597  fhFcrossVsZleading[i][r] = new TH3F(Form("%s_fhFcrossVsZleading_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs F_{cross} vs Z_{leading}^{ne} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);F_{cross};Z_{leading}^{ne}",triggerTitle[i],kRadius[r]),200,0,200,55,-0.05,1.05,55,-0.05,1.05);
598  fOutputList->Add(fhFcrossVsZleading[i][r]);
599 
600  fhChLeadZVsJetPt[i][r] = new TH2F(Form("%s_fhChLeadZVsJetPt_%1.1f",triggerName[i],kRadius[r]),Form("%s: Z of leading charged constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); Z_{leading}^{ch}",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1);
601  fOutputList->Add(fhChLeadZVsJetPt[i][r]);
602 
603  fhJetPtWithTrkThres[i][r] = new TH2F(Form("%s_fhJetPtWithTrkThres_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of jets containing tracks above certain threshold (15,25,40 GeV/c) (EMCal, R=%1.1f,Z<0.98);Threshold type;p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),3,0,3,nbins,xbins);
604  fOutputList->Add(fhJetPtWithTrkThres[i][r]);
605 
606  fhJetPtWithClsThres[i][r] = new TH2F(Form("%s_fhJetPtWithClsThres_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of jets containing clusters above certain threshold (15,25,40 GeV) (EMCal, R=%1.1f,Z<0.98);Threshold type;p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),3,0,3,nbins,xbins);
607  fOutputList->Add(fhJetPtWithClsThres[i][r]);
608 
609  fhJetPtVsLowPtCons[i][r][0] = new TH2F(Form("%s_fhJetPtVsLowPtCons_150-300MeV_%1.1f",triggerName[i],kRadius[r]),Form("%s: energy carried by constituents in 150-300MeV (EMCal, R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);p_{T,em}^{low} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1);
610  fOutputList->Add(fhJetPtVsLowPtCons[i][r][0]);
611 
612  fhJetPtVsLowPtCons[i][r][1] = new TH2F(Form("%s_fhJetPtVsLowPtCons_300-500MeV_%1.1f",triggerName[i],kRadius[r]),Form("%s: energy carried by constituents in 300-500MeV (EMCal, R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);p_{T,em}^{low} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1);
613  fOutputList->Add(fhJetPtVsLowPtCons[i][r][1]);
614 
615  if(fCheckTPCOnlyVtx)
616  {
617  fhJetInTPCOnlyVtx[i][r] = new TH3F(Form("%s_fhJetInTPCOnlyVtx_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet pt in events with only TPC vertex (Full, R=%1.1f);p_{T}^{jet} (GeV/c);#phi;#eta",triggerTitle[i],kRadius[r]),20,0,100,36,0,360,20,-1,1);
618  fOutputList->Add(fhJetInTPCOnlyVtx[i][r]);
619  }
620 
621  if(fStudySubEInHC)
622  {
623  for(Int_t k=0; k<5; k++)
624  {
625  fhSubClsEVsJetPt[i][r][k] = new TH2F(Form("%s_fhSubClsEVsJetPt_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: relative %s%% subtracted cluster Et vs jet pt (R=%1.1f);jet p_{T} (GeV/c);E_{t}^{sub}/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,50,0,0.5);
626  fOutputList->Add(fhSubClsEVsJetPt[i][r][k]);
627 
628  fhHCTrkPtClean[i][r][k] = new TH2F(Form("%s_fhHCTrkPtClean_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: sum of track p_{T} that are cleanly subtracted vs jet pt (%s%%, R=%1.1f);jet p_{T} (GeV/c);#sum(p_{T,trk}^{clean})/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,100,-0.005,0.995);
629  fOutputList->Add(fhHCTrkPtClean[i][r][k]);
630 
631  fhHCTrkPtAmbig[i][r][k] = new TH2F(Form("%s_fhHCTrkPtAmbig_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: sum of track p_{T} that are ambiguously subtracted vs jet pt (%s%%, R=%1.1f);jet p_{T} (GeV/c);#sum(p_{T,trk}^{ambig})/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,100,-0.005,0.995);
632  fOutputList->Add(fhHCTrkPtAmbig[i][r][k]);
633  }
634  }
635  }
636  }
637  if(i<2)
638  {
639  fhNMatchedTrack[i] = new TH1F(Form("%s_fhNMatchedTrack",triggerName[i]),Form("%s: # of matched tracks per cluster; N_{mth}",triggerTitle[i]),5,0,5);
640  fOutputList->Add(fhNMatchedTrack[i]);
641 
642  for(Int_t j=0; j<4; j++)
643  {
644  fhSubEVsTrkPt[i][j] = new TH2F(Form("%s_fhSubEVsTrkPt_%s",triggerName[i],fraction[j+1]),Form("%s: fractional subtracted energy (%s%% HC);#sum(p_{ch,T}^{mth}) (GeV/c);E_{sub}/#sum(P_{ch}^{mth})",triggerTitle[i],fraction[j+1]),50,0,50,110,0,1.1);
645  fOutputList->Add(fhSubEVsTrkPt[i][j]);
646  }
647  }
648  }
649 
650  if(fRunUE)
651  {
652  for(Int_t k=0; k<2; k++)
653  {
654  for(Int_t l=0; l<2; l++)
655  {
656  fhUEJetPtNorm[i][k][l] = new TH1F(Form("%s_fhUEJetPtNorm_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} in TPC (%s in %s event);p_{T,jet}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins);
657  fOutputList->Add(fhUEJetPtNorm[i][k][l]);
658 
659  fhUEJetPtVsSumPt[i][k][l] = new TH2F(Form("%s_fhUEJetPtVsSumPt_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} vs underlying event contribution (R=0.4,%s in %s event);p_{T,jet}^{ch} (GeV/c);p_{T,UE}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins,500,0,50);
660  fOutputList->Add(fhUEJetPtVsSumPt[i][k][l]);
661 
662  fhUEJetPtVsConsPt[i][k][l] = new TH2F(Form("%s_fhUEJetPtVsConsPt_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} vs constituent pt in UE (R=0.4,%s in %s event);p_{T,jet}^{ch} (GeV/c);p_{T,cons}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins,500,0,50);
663  fOutputList->Add(fhUEJetPtVsConsPt[i][k][l]);
664  }
665  }
666  }
667 
668  if(fSysJetTrigEff)
669  {
670  fhClsE[i] = new TH1F(Form("%s_fhClsE",triggerName[i]),Form("%s: cluster E;E (GeV)",triggerTitle[i]),1000,0,100);
671  fOutputList->Add(fhClsE[i]);
672  }
673 
674  if(fSysClusterizer)
675  {
676  for(Int_t k=0; k<2; k++)
677  {
678  fhSysClusterE[i][k] = new TH1F(Form("%s_fhSysClusterE_%sHC",triggerName[i],clusterizerName[k]),Form("%s: cluster E %s hadronic correction;E (GeV)",triggerTitle[i],clusterizerName[k]),100,0,100);
679  fOutputList->Add(fhSysClusterE[i][k]);
680 
681  fhSysNCellVsClsE[i][k] = new TH2F(Form("%s_fhSysNCellVsClsE_%sHC",triggerName[i],clusterizerName[k]),Form("%s: NCell vs cluster E %s hadronic correction;E (GeV);NCell",triggerTitle[i],clusterizerName[k]),100,0,100,50,0,50);
682  fOutputList->Add(fhSysNCellVsClsE[i][k]);
683  }
684  }
685  }
686 
687 
688  if(fIsMC)
689  {
690  if(fStudyMcOverSubE)
691  {
692  for(Int_t r=0; r<kNRadii; r++)
693  {
694  if(!fRadius.Contains("0.4") && r==0) continue;
695  if(!fRadius.Contains("0.2") && r==1) continue;
696  if(!fRadius.Contains("0.3") && r==2) continue;
697  for(Int_t i=0; i<5; i++)
698  {
699  fHCOverSubE[r][i] = new TH2F(Form("%s_HC_over_sub_e_%s_%1.1f",triggerName[2],fraction[i],kRadius[r]),Form("%s: oversubtracted neutral Et by %s%% HC (R=%1.1f);particle jet p_{T} (GeV/c);#DeltaE_{t} (GeV)",triggerName[2],fraction[i],kRadius[r]),nbins,xbins,200,-49.75,50.25);
700  fOutputList->Add(fHCOverSubE[r][i]);
701  fHCOverSubEFrac[r][i] = new TH2F(Form("%s_HC_over_sub_e_frac_%s_%1.1f",triggerName[2],fraction[i],kRadius[r]),Form("%s: relative oversubtracted neutral Et fraction by %s%% HC (R=%1.1f);jet p_{T} (GeV/c);#DeltaE_{t}/p_{T}^{jet}",triggerName[2],fraction[i],kRadius[r]),nbins,xbins,200,-0.995,1.005);
702  fOutputList->Add(fHCOverSubEFrac[r][i]);
703  }
704  }
705  }
706  if(fRunSecondaries && fAnaType==0)
707  {
708  for(Int_t i=0; i<2; i++)
709  {
710  for(Int_t r=0; r<3; r++)
711  {
712  fhNKFracVsJetPt[i][r] = new TH2F(Form("%s_fhNKFracVsJetPt_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: energy fraction carried by n,k^{0}_{L} vs jet p_{T} (R=%1.1f,|#eta|<%1.1f);p_{T,jet} (GeV/c);fraction",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,0,1);
713  fOutputList->Add(fhNKFracVsJetPt[i][r]);
714 
715  fhWeakFracVsJetPt[i][r] = new TH2F(Form("%s_fhWeakFracVsJetPt_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: energy fraction carried by k^{0}_{S},hyperon vs jet p_{T} (R=%1.1f,|#eta|<%1.1f);p_{T,jet} (GeV/c);fraction",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,0,1);
716  fOutputList->Add(fhWeakFracVsJetPt[i][r]);
717 
718  fhJetResponseNK[i][r] = new TH2F(Form("%s_fhJetResponseNK_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
719  fOutputList->Add(fhJetResponseNK[i][r]);
720 
721  fhJetResponseWP[i][r] = new TH2F(Form("%s_fhJetResponseWP_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
722  fOutputList->Add(fhJetResponseWP[i][r]);
723 
724  fhJetResolutionNK[i][r] = new TH2F(Form("%s_fhJetResolutionNK_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L}: (p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{rec} vs p_{T,jet}^{rec} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T}/p_{T}",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,-0.995,1.005);
725  fOutputList->Add(fhJetResolutionNK[i][r]);
726 
727  fhJetResolutionWP[i][r] = new TH2F(Form("%s_fhJetResolutionWP_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon: (p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{rec} vs p_{T,jet}^{rec} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T}/p_{T}",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,-0.995,1.005);
728  fOutputList->Add(fhJetResolutionWP[i][r]);
729 
730  fhJetResponseNKSM[i][r] = new TH2F(Form("%s_fhJetResponseNKSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L} via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
731  fOutputList->Add(fhJetResponseNKSM[i][r]);
732 
733  fhJetResponseWPSM[i][r] = new TH2F(Form("%s_fhJetResponseWPSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon via matching(R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins);
734  fOutputList->Add(fhJetResponseWPSM[i][r]);
735 
736  fhJetResolutionNKSM[i][r] = new TH3F(Form("%s_fhJetResolutionNKSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet resolution due to missing n and k^{0}_{L} via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T,jet}/p_{T,jet}^{rec};dR",triggerName[2],kRadius[r],i*0.5+0.5),220,0,220,200,-0.995,1.005,100,0,1);
737  fOutputList->Add(fhJetResolutionNKSM[i][r]);
738 
739  fhJetResolutionWPSM[i][r] = new TH3F(Form("%s_fhJetResolutionWPSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet resolution due tok^{0}_{S}, #Lambda and hyperon via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T,jet}/p_{T,jet}^{rec};dR",triggerName[2],kRadius[r],i*0.5+0.5),220,0,220,200,-0.995,1.005,100,0,1);
740  fOutputList->Add(fhJetResolutionWPSM[i][r]);
741  }
742  }
743  }
744  }
745 
746  printf("\n=======================================\n");
747  printf("===== Jet task configuration ==========\n");
748 
749  if(fNonStdBranch.Length()!=0)
750  {
751  const char* species[3] = {"in","ch","ne"};
752  const char* algorithm[2] = {"akt","kt"};
753  const char* radii[kNRadii] = {"04","02","03"};
754  for(Int_t s=0; s<3; s++)
755  {
756  if(!fFindChargedOnlyJet && s==1) continue;
757  if(!fFindNeutralOnlyJet && s==2) continue;
758  for(Int_t a=0; a<2; a++)
759  {
760  if(!fAlgorithm.Contains("aKt") && a==0) continue;
761  if(!fAlgorithm.Contains("kt") && a==1) continue;
762  for(Int_t r=0; r<kNRadii; r++)
763  {
764  if(!fRadius.Contains("0.4") && r==0) continue;
765  if(!fRadius.Contains("0.2") && r==1) continue;
766  if(!fRadius.Contains("0.3") && r==2) continue;
767  if(fAnaType==0)
768  {
769  fJetTCA[s][a][r] = new TClonesArray("AliAODJet",0);
770  fJetTCA[s][a][r]->SetName(Form("Jet_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
771  AddAODBranch("TClonesArray",&fJetTCA[s][a][r],fNonStdFile.Data());
772  printf("Add branch: Jet_%s_%s_%s_%s\n",species[s],algorithm[a],radii[r],fNonStdBranch.Data());
773  }
774 
775  fDetJetFinder[s][a][r] = new AliFJWrapper(Form("DetJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()),Form("DetJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
776  if(a==0) fDetJetFinder[s][a][r]->SetAlgorithm(fastjet::antikt_algorithm);
777  if(a==1) fDetJetFinder[s][a][r]->SetAlgorithm(fastjet::kt_algorithm);
778  if(fRecombinationScheme==0) fDetJetFinder[s][a][r]->SetRecombScheme(fastjet::E_scheme);
779  fDetJetFinder[s][a][r]->SetR(kRadius[r]);
780  fDetJetFinder[s][a][r]->SetMaxRap(0.9);
781  fDetJetFinder[s][a][r]->Clear();
782 
783  if(s==0 && a==0)
784  {
785  if(fIsMC && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase))
786  {
787  if(fAnaType==0)
788  {
789  fMcTruthAntikt[r] = new TClonesArray("AliAODJet",0);
790  fMcTruthAntikt[r]->SetName(Form("Jet_mc_truth_in_akt_%s_%s",radii[r],fNonStdBranch.Data()));
791  AddAODBranch("TClonesArray",&fMcTruthAntikt[r],fNonStdFile.Data());
792  printf("Add branch: Jet_mc_truth_in_akt_%s_%s\n",radii[r],fNonStdBranch.Data());
793  }
794 
795  fTrueJetFinder[r] = new AliFJWrapper(Form("TrueJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()),Form("TrueJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
796  fTrueJetFinder[r]->SetAlgorithm(fastjet::antikt_algorithm);
797  fTrueJetFinder[r]->SetR(kRadius[r]);
798  fTrueJetFinder[r]->SetMaxRap(0.9);
799  if(fRecombinationScheme==0) fTrueJetFinder[r]->SetRecombScheme(fastjet::E_scheme);
800  fTrueJetFinder[r]->Clear();
801  }
802  }
803  }
804  }
805  }
806  }
807 
808  fRandomGen = new TRandom3(0);
809  if(fCheckTrkEffCorr && fAnaType==0)
810  {
811  TFile f ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffFit.root","read");
812  for(Int_t i=0; i<3; i++)
813  {
814  fTrkEffFunc[i] = new TF1(*((TF1*)f.Get(Form("Trk_eff_fit_%d",i+1))));
815  }
816  f.Close();
817 
818  if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0)
819  {
820  TFile f1 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffSampling.root","read");
821  for(Int_t j=0; j<2; j++)
822  {
823  Int_t tmp = j;
824  if(fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0) tmp = 2;
825  for(Int_t r=0; r<kNRadii; r++)
826  {
827  fhCorrTrkEffPtBin[j][r] = new TH1F(*((TH1F*)f1.Get(Form("%s_%s_NTrackPerPtBin_%1.1f",fPeriod.Data(),triggerName[tmp],kRadius[r]))));
828  for(Int_t i=0; i<kNBins; i++)
829  {
830  fhCorrTrkEffSample[j][r][i] = new TH1F(*((TH1F*)f1.Get(Form("%s_%s_ChTrkPt_Bin%d_%1.1f",fPeriod.Data(),triggerName[tmp],i+1,kRadius[r]))));
831  }
832  }
833  }
834  f1.Close();
835  }
836  }
837 
838  if(fRunSecondaries && fAnaType==0)
839  {
840  const char *secondaryName[3] = {"k0S","lamda","hyperon"};
841  TFile f2 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/SecondaryResponse.root","read");
842  for(Int_t j=0; j<3; j++)
843  {
844  fhSecondaryResponse[j] = new TH2F(*((TH2F*)f2.Get(Form("DetectorReponse_%s",secondaryName[j]))));
845  }
846  }
847 
849  {
850  TString name = "TriggerCurve.root";
851  if(fAnaType==0) name = "/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TriggerCurve.root";
852  TFile f (name.Data(),"read");
853  fTriggerMask = new TH2F(*((TH2F*)f.Get("lhc11a_TriggerMask")));
854  if(fOfflineTrigger)
855  {
856  for(Int_t i=0; i<10; i++)
857  {
858  fTriggerEfficiency[i] = new TF1(*((TF1*)f.Get(Form("lhc11a_TriggerEfficiency_SM%d_fit",i))));
859  fTriggerCurve[i] = new TH1D(*((TH1D*)f.Get(Form("lhc11a_TriggerCurve_SM%d",i))));
860  }
861  }
862  f.Close();
863  }
864 
865  fClusterEResolution = new TF1("fClusterEResolution","sqrt([0]^2+[1]^2*x+([2]*x)^2)*0.01");
866  fClusterEResolution->SetParameters(4.35,9.07,1.63);
867 
868  fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
869  if(!fGeom)
870  {
871  AliError("No EMCal geometry is available!");
872  return;
873  }
874  fRecoUtil = new AliEMCALRecoUtils();
876  fRecoUtil->SwitchOnRejectExoticCluster();
877  else
878  {
879  fRecoUtil->SwitchOffRejectExoticCluster();
880  fRecoUtil->SwitchOffRejectExoticCell();
881  }
883  {
884  fRecoUtil->SwitchOnBadChannelsRemoval();
885  // Remove the problematic SM4 region due to wrong event sequence
886  for(Int_t ieta=0; ieta<36; ieta++)
887  for(Int_t iphi=0; iphi<8; iphi++)
888  fRecoUtil->SetEMCALChannelStatus(4,ieta,iphi);
889  }
890  else
891  fRecoUtil->SwitchOffBadChannelsRemoval();
892 
893  fRecoUtil->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrected);
894  if(fSysNonLinearity)
895  {
896  fNonLinear = new TF1("TB_oldBest","([0])*(1./(1.+[1]*exp(-x/[2]))*1./(1.+[3]*exp((x-[4])/[5])))*1/([6])",0.1,110.);
897  fNonLinear->SetParameters(0.99078, 0.161499, 0.655166, 0.134101, 163.282, 23.6904, 0.978);
898  }
899  if(fSysTrkClsMth)
900  {
901  fRecoUtil->SwitchOnCutEtaPhiSeparate();
902  fRecoUtil->SetCutEta(fCutdEta);
903  fRecoUtil->SetCutPhi(fCutdPhi);
904  }
905 
906  fTrackArray = new TObjArray();
907  fTrackArray->SetOwner(1);
908  fClusterArray = new TObjArray();
909  fClusterArray->SetOwner(1);
910  fMcPartArray = new TArrayI();
911 
912 
913  //error calculation in THnSparse
914  Int_t nObj = fOutputList->GetEntries();
915  for(Int_t i=0; i<nObj; i++)
916  {
917  TObject *obj = (TObject*) fOutputList->At(i);
918  if (obj->IsA()->InheritsFrom( "THnSparse" ))
919  {
920  THnSparseF *hn = (THnSparseF*)obj;
921  hn->Sumw2();
922  }
923  }
924 
925 
926  PrintConfig();
927  BookHistos();
928  PostData(1, fOutputList);
929 }
930 
931 
932 
933 //________________________________________________________________________
935 {
936  // Book histograms.
937 
938  if(fVerbosity>10) printf("[i] Booking histograms \n");
939  return;
940 }
941 
942 //________________________________________________________________________
944 {
945  //
946  // Main loop, called for each event.
947  //
948 
949  // Get event pointers, check for signs of life
950  Double_t vtxTrueZ = -100;
951  if(fIsMC)
952  {
953  fMC = MCEvent();
954  if (!fMC)
955  {
956  printf("ERROR: Could not retrieve MC event");
957  return;
958  }
959  fStack = fMC->Stack();
960  TParticle *particle = fStack->Particle(0);
961  if(particle)
962  {
963  vtxTrueZ = particle->Vz(); // True vertex z in MC events
964  if(fVerbosity>10) printf("Generated vertex coordinate: (x,y,z) = (%4.2f, %4.2f, %4.2f)\n", particle->Vx(), particle->Vy(), particle->Vz());
965  }
966  }
967 
968  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
969  if (!fESD)
970  {
971  AliError("fESD is not available");
972  return;
973  }
974 
975  fhJetEventStat->Fill(0.5);
976  fhJetEventCount->Fill(0.5);
977  if(fIsMC)
978  {
979  GetMCInfo();
980  if(fVertexGenZ[0]) fVertexGenZ[0]->Fill(vtxTrueZ);
981  }
982 
983  // Centrality, vertex, other event variables...
984  UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
985  if (trigger==0) return;
986  if(fCheckTPCOnlyVtx) CheckTPCOnlyVtx(trigger);
987  if(!fIsMC)
988  {
989  if (trigger & AliVEvent::kFastOnly) return; // Reject fast trigger cluster
990  if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0)
991  {
992  if (trigger & AliVEvent::kMB) fTriggerType = 0;
993  else if(trigger & AliVEvent::kEMC1) fTriggerType = 1;
994  else fTriggerType = -1;
995  }
996  else if (fPeriod.CompareTo("lhc11c",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11d",TString::kIgnoreCase)==0)
997  {
998  if (trigger & AliVEvent::kINT7) fTriggerType = 0;
999  else if(trigger & AliVEvent::kEMC7) fTriggerType = 1;
1000  else fTriggerType = -1;
1001  }
1002  else
1003  {
1004  return;
1005  }
1006  }
1007  else
1008  {
1009  if(!fPhySelForMC) fTriggerType = 0;
1010  else if (trigger & AliVEvent::kAnyINT) fTriggerType = 0;
1011  else fTriggerType = -1;
1012 
1013  if(fOfflineTrigger)
1014  {
1017  }
1018  }
1019 
1020  if(fTriggerType==-1)
1021  {
1022  if(fVerbosity>10) printf("Error: worng trigger type %s\n",(fESD->GetFiredTriggerClasses()).Data());
1023  return;
1024  }
1025 
1026  fhJetEventCount->Fill(1.5+fTriggerType*4);
1027  if(fRejectPileup && fESD->IsPileupFromSPD() ) return; // reject pileup
1028  fhJetEventStat->Fill(11.5);
1029  fhJetEventCount->Fill(2.5+fTriggerType*4);
1030 
1031 
1032  fIsTPCOnlyVtx = 0;
1033  // Reject LED events
1034  if (IsLEDEvent())
1035  {
1036  fhJetEventStat->Fill(8.5);
1037  return;
1038  }
1039 
1040  fhJetEventStat->Fill(1.5+fTriggerType*2);
1041 
1042  // Check if primary vertex exists
1043  if (!HasPrimaryVertex()) return;
1044  fhJetEventStat->Fill(5.5+fTriggerType);
1045  fhJetEventCount->Fill(3.5+fTriggerType*4);
1046 
1047  const AliESDVertex* vtx = fESD->GetPrimaryVertex();
1048  Double_t zVertex = vtx->GetZ();
1049  if(fEventZ[fTriggerType]) fEventZ[fTriggerType]->Fill(zVertex);
1050  if(fVertexGenZ[1]) fVertexGenZ[1]->Fill(vtxTrueZ);
1051 
1052  // Check if |Z_vtx|<10cm
1053  if( TMath::Abs(zVertex) > fZVtxMax ) return;
1054  fhJetEventCount->Fill(4.5+fTriggerType*4);
1055 
1056  // Check if only TPC vertex exists primitive
1058 
1059  if(!fIsMC && fTriggerType==1)
1060  {
1061  // Check if event has valid trigger bit
1063  if(fIsEventTriggerBit)
1064  {
1065  fhJetEventStat->Fill(7.5);
1066  fhJetEventCount->Fill(9.5);
1067  }
1068  else return;
1069  }
1070  fhJetEventStat->Fill(2.5+fTriggerType*2);
1071  if(fIsTPCOnlyVtx) fhJetEventStat->Fill(9.5+fTriggerType);
1072 
1073  // Check if event contains exotic clusters
1074  CheckExoticEvent();
1075 
1076  // Write jet tree into AOD output in local analysis mode
1077  if(fAnaType==0) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
1078 
1079  // Clean up arrays
1080  for(Int_t s=0; s<3; s++)
1081  {
1082  for(Int_t a=0; a<2; a++)
1083  {
1084  for(Int_t r=0; r<kNRadii; r++)
1085  {
1086  if(fJetTCA[s][a][r]) fJetTCA[s][a][r]->Delete();
1087  if(fDetJetFinder[s][a][r]) fDetJetFinder[s][a][r]->Clear();
1088 
1089  if(s==0 && a==0)
1090  {
1091  if(fMcTruthAntikt[r]) fMcTruthAntikt[r]->Delete();
1092  if(fTrueJetFinder[r]) fTrueJetFinder[r]->Clear();
1093  }
1094  }
1095  }
1096  }
1097 
1098  if(fVerbosity>5) printf("# of jets after clear: %d\n",fJetTCA[0][0][0]->GetEntries());
1099 
1100  if(fTrackArray) fTrackArray->Delete();
1101  if(fClusterArray) fClusterArray->Delete();
1102  fMcPartArray->Reset();
1103 
1104 
1105  // get the tracks and fill the input vector for the jet finders
1106  GetESDTrax();
1107 
1108  // get EMCal clusters and fill the input vector for the jet finders
1110 
1111  for(Int_t s=0; s<3; s++)
1112  {
1113  for(Int_t a=0; a<2; a++)
1114  {
1115  for(Int_t r=0; r<kNRadii; r++)
1116  {
1117  //Detector jets
1118  if(fDetJetFinder[s][a][r])
1119  {
1120  FindDetJets(s,a,r);
1121  if(fJetTCA[s][a][r]) FillAODJets(fJetTCA[s][a][r], fDetJetFinder[s][a][r], 0);
1122  if(s==0 && a==0 && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase))
1123  {
1124  if(fSaveQAHistos) AnalyzeJets(fDetJetFinder[s][a][r],fTriggerType, r); // run analysis on jets
1125  }
1126  }
1127  }
1128  }
1129  }
1130  if(fRunUE && fDetJetFinder[1][0][0]) RunAnalyzeUE(fDetJetFinder[1][0][0],fTriggerType,0); // Run analysis on underlying event
1131 
1132  if(fIsMC)
1133  {
1134  for(Int_t r=0; r<kNRadii; r++)
1135  if(fTrueJetFinder[r]) ProcessMC(r); // find particle level jets
1136  }
1137 
1138  // Fast Jet calls END --------------------------------------------------------
1139 
1140  PostData(1, fOutputList);
1141  return;
1142 }
1143 
1144 
1145 //________________________________________________________________________
1147 {
1148  //
1149  // Jet finding is happening here
1150  //
1151 
1152  Bool_t isCh = kTRUE;
1153  Bool_t isNe = kTRUE;
1154  if(s==1) isNe = kFALSE;
1155  if(s==2) isCh = kFALSE;
1156 
1157  if(isCh)
1158  {
1159  // Feed in charged tracks
1160  Int_t countTracks = 0;
1161  Int_t ntracks = fTrackArray->GetEntriesFast();
1162  for(Int_t it=0; it<ntracks; it++)
1163  {
1164  AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
1165  if(!t) continue;
1166  countTracks++;
1167  Double_t e = t->P();
1168  Double_t pt = t->Pt();
1169  if(s==1 && fRunUE && pt<1) continue;
1170  if(fRecombinationScheme==0) e = TMath::Sqrt(t->P()*t->P()+0.139*0.139);
1171  if(fSysTrkPtRes) pt += GetSmearedTrackPt(t);
1172  Double_t phi = t->Phi();
1173  Double_t eta = t->Eta();
1174  Double_t px = pt*TMath::Cos(phi);
1175  Double_t py = pt*TMath::Sin(phi);
1176  if(fConstrainChInEMCal && ( TMath::Abs(eta)>0.7 || phi>kPI || phi<TMath::DegToRad()*80) ) continue;
1177  fDetJetFinder[s][a][r]->AddInputVector(px,py,t->Pz(), e, it+1);
1178  if(fVerbosity>10) printf("%d: m_{ch}=%f\n",it+1,t->P()*t->P()-t->Px()*t->Px()-t->Py()*t->Py()-t->Pz()*t->Pz());
1179  }
1180  if(fVerbosity>5) printf("[i] # of tracks filled: %d\n",countTracks);
1181  }
1182  if(isNe)
1183  {
1184  // Feed in EMCal clusters
1185  Double_t vertex[3] = {0, 0, 0};
1186  fESD->GetVertex()->GetXYZ(vertex) ;
1187  TLorentzVector gamma;
1188 
1189  Int_t countClusters = 0;
1190  Int_t nclusters = fClusterArray->GetEntriesFast();
1191  for (Int_t ic = 0; ic < nclusters; ic++)
1192  {
1193  AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic);
1194  if(!cl) continue;
1195  cl->GetMomentum(gamma, vertex);
1196  countClusters++;
1197  Double_t e = gamma.P();
1198  if(fRecombinationScheme==0) e = TMath::Sqrt(gamma.P()*gamma.P()+0.139*0.139);
1199  fDetJetFinder[s][a][r]->AddInputVector(gamma.Px(), gamma.Py(), gamma.Pz(), e, (ic+1)*(-1));
1200  }
1201  if(fVerbosity>5) printf("[i] # of clusters filled: %d\n",countClusters);
1202  }
1203  // Run jet finding
1204  fDetJetFinder[s][a][r]->Run();
1205 }
1206 
1207 
1208 //________________________________________________________________________
1209 void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper *jetFinder, const Bool_t isTruth)
1210 {
1211  //
1212  // Fill the found jets into AOD output
1213  // Only consider jets pointing to EMCal and with pt above 1GeV/c
1214  //
1215 
1216  Int_t radiusIndex = 0;
1217  TString arrayName = fJetArray->GetName();
1218  if(arrayName.Contains("_02_")) radiusIndex = 1;
1219  if(arrayName.Contains("_03_")) radiusIndex = 2;
1220  std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1221  if(fVerbosity>5 && radiusIndex==0) printf("[i] # of jets in %s : %d\n",fJetArray->GetName(),(Int_t)jetsIncl.size());
1222  AliAODJet *jet = 0;
1223  Int_t jetCount = 0;
1224  for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1225  {
1226  if(fVerbosity>10) printf("fastjet: eta=%f, phi=%f\n",jetsIncl[ij].eta(),jetsIncl[ij].phi()*TMath::RadToDeg());
1227  if(!IsGoodJet(jetsIncl[ij],0)) continue;
1228 
1229  AliAODJet tmpJet (jetsIncl[ij].px(), jetsIncl[ij].py(), jetsIncl[ij].pz(), jetsIncl[ij].E());
1230  jet = new ((*fJetArray)[jetCount]) AliAODJet(tmpJet);
1231  jetCount++;
1232  if(fVerbosity>10 && radiusIndex==0) printf("AOD jet: ij=%d, pt=%f, eta=%f, phi=%f\n",ij, jet->Pt(), jet->Eta(),jet->Phi()*TMath::RadToDeg());
1233 
1234  jet->GetRefTracks()->Clear();
1235  std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
1236  Double_t totE=0, totPt=0, totChPt=0, leadChPt=0, neE=0, totNePt=0, leadNePt=0, leadPt=0;
1237  Double_t leadTrkType=0, nChPart = 0, nNePart = 0;
1238  Bool_t isHighPtTrigger = kFALSE, isTriggering = kFALSE, isHyperTrack = kFALSE;
1239  Int_t leadIndex = -1;
1240  for(UInt_t ic=0; ic<constituents.size(); ic++)
1241  {
1242  if(fVerbosity>10 && radiusIndex==0) printf("ic=%d: user_index=%d, E=%f\n",ic,constituents[ic].user_index(),constituents[ic].E());
1243  if(constituents[ic].user_index()<0)
1244  {
1245  nNePart ++;
1246  totNePt += constituents[ic].perp();
1247  if(constituents[ic].perp()>leadNePt)
1248  leadNePt = constituents[ic].perp();
1249 
1250  neE += constituents[ic].E();
1251  if(constituents[ic].perp()>fClsEtMax[fTriggerType])
1252  isHighPtTrigger = kTRUE;
1253  if(!isTruth)
1254  {
1255  AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1256  if(cluster->Chi2()>0.5) isTriggering = kTRUE;
1257  }
1258  }
1259  else
1260  {
1261  totChPt += constituents[ic].perp();
1262  nChPart ++;
1263  if(constituents[ic].perp()>leadChPt)
1264  {
1265  leadChPt = constituents[ic].perp();
1266  if(!isTruth)
1267  {
1268  AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1269  leadTrkType = track->GetTRDQuality();
1270  if(track->GetIntegratedLength()>500) isHyperTrack = kTRUE;
1271  }
1272  }
1273  if(constituents[ic].perp()>fTrkPtMax[fTriggerType])
1274  isHighPtTrigger = kTRUE;
1275  }
1276  TParticle part;
1277  part.SetMomentum(constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),constituents[ic].E());
1278  jet->AddTrack(&part); //The references are not usable, this line is aimed to count the number of contituents
1279  totE += constituents[ic].E();
1280  totPt += constituents[ic].perp();
1281  if(constituents[ic].perp()>leadPt)
1282  {
1283  leadPt = constituents[ic].perp();
1284  leadIndex = ic;
1285  }
1286  }
1287  if(fIsEventTriggerBit) jet->SetTrigger(AliAODJet::kTRDTriggered);
1288  if(isHighPtTrigger) jet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1289  if(fTriggerType==1) jet->SetTrigger(AliAODJet::kEMCALTriggered);
1290  if(GetLeadingZ(ij,jetFinder) > 0.98 ) jet->SetTrigger(AliAnalysisTaskFullppJet::kHighZ);
1291  if(isTriggering) jet->SetTrigger(AliAnalysisTaskFullppJet::kTrigger);
1292  if(isHyperTrack) jet->SetTrigger(AliAnalysisTaskFullppJet::kSuspicious);
1294  if(constituents[leadIndex].user_index()>0) jet->SetTrigger(AliAnalysisTaskFullppJet::kLeadCh);
1295 
1296  if(jetsIncl[ij].E()>0) jet->SetNEF(neE/jetsIncl[ij].E());
1297  jet->SetPtLeading(leadPt);
1298  jet->SetPtSubtracted(leadTrkType,0);
1299 
1300  Double_t effAErrCh = 0, effAErrNe = leadTrkType;
1301  Double_t chBgEnergy = 10, neBgEnergy = nNePart;
1302  if(!isTruth)
1303  {
1304  effAErrCh = GetMeasuredJetPtResolution(ij,jetFinder);
1305  if(fIsMC && constituents[leadIndex].user_index()>0)
1306  {
1307  AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[leadIndex].user_index()-1);
1308  Double_t pt = track->Pt();
1309  Int_t ipart = track->GetLabel();
1310  if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
1311  {
1312  AliVParticle* vParticle = fMC->GetTrack(ipart);
1313  if(vParticle)
1314  {
1315  Double_t truePt = vParticle->Pt();
1316  chBgEnergy = (truePt-pt)/truePt;
1317  }
1318  }
1319  }
1320  }
1321 
1322  jet->SetEffArea(leadPt, jetFinder->GetJetArea(ij),effAErrCh,effAErrNe);
1323  jet->SetBgEnergy(chBgEnergy,neBgEnergy);
1324  if(fVerbosity>10 && isTruth) cout<<jet->ErrorEffectiveAreaCharged()<<" "<<jet->ErrorEffectiveAreaNeutral()<<endl;
1325  if(fVerbosity>10) cout<<"jet pt="<<jetsIncl[ij].perp()<<", nef="<<jet->GetNEF()<<", trk eff corr="<<chBgEnergy<<endl;
1326 
1327  if(fVerbosity>5)
1328  printf("# of ref tracks: %d = %d, and nef=%f\n",jet->GetRefTracks()->GetEntries(), (Int_t)constituents.size(), jet->GetNEF());
1329 
1330  // For catch good high-E jets
1331  if(fSpotGoodJet && !isTruth)
1332  {
1333  if(jetsIncl[ij].perp()>100)
1334  {
1335  printf("\n\n--- HAHAHA: High pt jets ---\n");
1336  printf("File: %s, event = %d, pt=%f\n",CurrentFileName(),(Int_t)Entry(),jetsIncl[ij].perp());
1337  printf("%s , pt < %f\n", fJetArray->GetName(), fTrkPtMax[1]);
1338  printf("Jet: E=%f, eta=%f, phi=%f, # of constituents=%d, nef=%f\n",jetsIncl[ij].E(),jetsIncl[ij].eta(), jetsIncl[ij].phi()*TMath::RadToDeg(), (Int_t)constituents.size(),jet->GetNEF());
1339  for(UInt_t ic=0; ic<constituents.size(); ic++)
1340  {
1341  if(constituents[ic].user_index()<0)
1342  {
1343  AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1344  printf("id = %d, cluster with pt = %f, ncell = %d\n",ic,constituents[ic].perp(), cluster->GetNCells());
1345  }
1346  else
1347  {
1348  AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1349  printf("id = %d, track with pt = %f, track class = %d, originalPt = %f\n",ic,constituents[ic].perp(),(Int_t)track->GetTRDQuality(), track->GetIntegratedLength());
1350 
1351  }
1352  }
1353  printf("==============================\n\n");
1354  }
1355  }
1356  // End of catching good high-E jets
1357  }
1358  if(fVerbosity>5) printf("%s has %d jets\n",fJetArray->GetName(), fJetArray->GetEntries());
1359 }
1360 
1361 //________________________________________________________________________
1363 {
1364  //
1365  // Find the parent parton for a given MC particle by tracing back the parent
1366  // DUMMY; NEED improvement
1367 
1368  return ipart;
1369 }
1370 
1371 
1372 //________________________________________________________________________
1374 {
1375  //
1376  // Check if the jet pt and direction fulfill the requirement
1377  //
1378 
1379  if(jet.perp()<1) return kFALSE;
1380  if(TMath::Abs(jet.eta())>(0.7-rad)) return kFALSE;
1381  if(jet.phi() < (80*TMath::DegToRad()+rad) || jet.phi() > (180*TMath::DegToRad()-rad) ) return kFALSE;
1382  return kTRUE;
1383 }
1384 
1385 
1386 //________________________________________________________________________
1388 {
1389  //
1390  // Main function for jet finding in MC
1391  //
1392 
1393  Int_t npart = fMC->GetNumberOfTracks();
1394  fMcPartArray->Set(npart);
1395  Int_t countPart = 0;
1396  for(Int_t ipart=0; ipart<npart; ipart++)
1397  {
1398  AliVParticle* vParticle = fMC->GetTrack(ipart);
1399  if(!IsGoodMcPartilce(vParticle,ipart)) continue;
1400 
1401  Int_t pdgCode = vParticle->PdgCode();
1402  if( fRejectNK && (pdgCode==130 || pdgCode==2112) ) continue;
1403 
1404  if( fRejectWD && (pdgCode==310 || pdgCode==3112 || pdgCode==3122 || pdgCode==3222 || pdgCode==3312 || pdgCode==3322 || pdgCode==3334)) continue;
1405 
1406  if( fChargedMC && vParticle->Charge()==0 ) continue;
1407 
1408  Int_t index = 1;
1409  if(vParticle->Charge()==0) { index=-1; }
1410  fMcPartArray->AddAt(ipart, countPart);
1411  countPart++;
1412 
1413  fTrueJetFinder[r]->AddInputVector(vParticle->Px(), vParticle->Py(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
1414  if(fVerbosity>10)
1415  printf("Input particle: ipart=%d, pdg=%d, species=%s, charge=%d, E=%4.3f\n",(ipart+1)*index,pdgCode,vParticle->GetName(), vParticle->Charge(),vParticle->P());
1416  }
1417  fMcPartArray->Set(countPart);
1418  fTrueJetFinder[r]->Run();
1419 
1422  if(fRunUE && r==0) RunAnalyzeUE(fTrueJetFinder[r], 2, 1);
1423 
1424  // Run analysis on secondary particles
1425  if(fRunSecondaries && fAnaType==0)
1426  {
1427  for(Int_t i=0; i<2; i++)
1428  {
1432  }
1433  }
1434 }
1435 
1436 
1437 //________________________________________________________________________
1438 void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t type, const Int_t r)
1439 {
1440  //
1441  // Fill all the QA plots for jets
1442  // Especailly all the constituents plot must be filled here since the information
1443  // is not available in the output AOD
1444 
1445  Double_t vertex[3] = {0, 0, 0};
1446  fESD->GetVertex()->GetXYZ(vertex) ;
1447  TLorentzVector gamma;
1448 
1449  const Int_t nBins = fJetEnergyFraction[type][r]->GetAxis(1)->GetNbins();
1450  Float_t radiusCut[nBins];
1451  Float_t eFraction[nBins];
1452  Int_t nPart[nBins];
1453  for(Int_t i=0; i<nBins; i++)
1454  {
1455  radiusCut[i] = (fJetEnergyFraction[type][r]->GetAxis(1)->GetBinWidth(1)) * (i+1);
1456  eFraction[i] = 0.;
1457  nPart[nBins] = 0;
1458  }
1459  std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1460 
1461  for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1462  {
1463  if(type<2 && fCheckTPCOnlyVtx && fIsTPCOnlyVtx)
1464  {
1465  fhJetInTPCOnlyVtx[type][r]->Fill(jetsIncl[ij].perp(),jetsIncl[ij].phi()*TMath::RadToDeg(),jetsIncl[ij].eta());
1466  }
1467  if(!IsGoodJet(jetsIncl[ij],kRadius[r])) continue; // Fidiual cut
1468  Float_t jetEta = jetsIncl[ij].eta();
1469  Float_t jetPhi = jetsIncl[ij].phi();
1470  Float_t jetE = jetsIncl[ij].E();
1471  Float_t jetPt = jetsIncl[ij].perp();
1472 
1473  std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
1474  Double_t neE=0, leadingZ = 0, maxPt = 0;
1475  Int_t constituentType = -1, leadingIndex = 0;
1476  for(UInt_t ic=0; ic<constituents.size(); ic++)
1477  {
1478  if(constituents[ic].perp()>maxPt)
1479  {
1480  maxPt = constituents[ic].perp();
1481  leadingIndex = constituents[ic].user_index();
1482  }
1483 
1484  if(constituents[ic].user_index()<0)
1485  {
1486  neE += constituents[ic].E();
1487  constituentType = 3;
1488  }
1489  else
1490  {
1491  if(type==2) constituentType = 0;
1492  else
1493  {
1494  AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1495  constituentType = (Int_t)track->GetTRDQuality();
1496  }
1497  }
1498  Double_t cz = GetZ(constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),jetsIncl[ij].px(),jetsIncl[ij].py(),jetsIncl[ij].pz());
1499  fhJetPtVsZ[type][r]->Fill(jetPt,cz, constituentType);
1500  if(cz>leadingZ) leadingZ = cz;
1501  }
1502 
1503  if(type<2 && leadingIndex<0)
1504  {
1505  AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(leadingIndex*(-1)-1);
1506  fhFcrossVsZleading[type][r]->Fill(jetPt,GetExoticEnergyFraction(cluster),leadingZ);
1507  }
1508 
1509  if(leadingZ > 0.98) continue; // Z cut
1510 
1511  fJetCount[type][r]->Fill(jetPt);
1512  if(type==1 && fIsExoticEvent3GeV) fhJetPtInExoticEvent[0][r]->Fill(jetPt);
1513  if(type==1 && fIsExoticEvent5GeV) fhJetPtInExoticEvent[1][r]->Fill(jetPt);
1514 
1515  Double_t totTrkPt[3] = {0.,0.,0.};
1516  Double_t chPt = 0;
1517  for(Int_t i=0; i<nBins; i++) { eFraction[i] = 0.; nPart[i] = 0;}
1518  Double_t mcSubE=0.;
1519  Double_t subClsE[5] = {0,0,0,0,0};
1520  Double_t subTrkPtClean[5] = {0,0,0,0,0};
1521  Double_t subTrkPtAmbig[5] = {0,0,0,0,0};
1522  Double_t fraction[5] = {270,0.3,0.5,0.7,1};
1523  Double_t leadChPt=0., leadNePt=0.;
1524  Int_t leadChIndex = -1;
1525  for(UInt_t ic=0; ic<constituents.size(); ic++)
1526  {
1527  Float_t partEta = constituents[ic].eta();
1528  Float_t partPhi = constituents[ic].phi();
1529  Float_t partE = constituents[ic].E();
1530  Float_t partPt = constituents[ic].perp();
1531 
1532  if(constituents[ic].user_index()<0)
1533  {
1534  fhNeutralPtInJet[type][r]->Fill(jetPt, partPt);
1535  if(partPt>leadNePt)
1536  leadNePt = partPt;
1537  if(type<2)
1538  {
1539  AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1540  Double_t clsE = cluster->E();
1541  if(cluster->Chi2()>0.5) fhTrigNeuPtInJet[type][r]->Fill(jetPt, partPt);
1543  {
1544  cluster->GetMomentum(gamma, vertex);
1545  Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
1546 
1547  Double_t subEtmp = cluster->GetDispersion();
1548  Double_t mipETmp = cluster->GetEmcCpvDistance();
1549  mcSubE += cluster->GetDistanceToBadChannel()*sinTheta;
1550  subClsE[0] += ((mipETmp>clsE)?clsE:mipETmp)*sinTheta;
1551  for(Int_t j=1; j<5; j++)
1552  {
1553  subClsE[j] += ((fraction[j]*subEtmp>clsE)?clsE:fraction[j]*subEtmp)*sinTheta;
1554  }
1555  }
1556  }
1557  }
1558  else
1559  {
1560  if(partPt>leadChPt) {leadChPt = partPt;leadChIndex=ic;}
1561  fhChargedPtInJet[type][r]->Fill(jetPt, partPt);
1562  chPt += constituents[ic].perp();
1563  if(type<2)
1564  {
1565  AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1566  Int_t trkType = (Int_t)track->GetTRDQuality();
1567  totTrkPt[trkType] += partPt;
1568  Int_t clusterIndex = track->GetEMCALcluster();
1569  Int_t clusterPos = -1;
1570  Bool_t isExist = kFALSE;
1571  for(Int_t j=0; j<fClusterArray->GetEntriesFast(); j++)
1572  {
1573  AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(j);
1574  if( clusterIndex == cluster->GetID() )
1575  {
1576  isExist = kTRUE;
1577  clusterPos = j;
1578  break;
1579  }
1580  }
1581  if(isExist)
1582  {
1583  AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(clusterPos);
1584  Double_t subEtmp = cluster->GetDispersion();
1585  Double_t mipETmp = cluster->GetEmcCpvDistance();
1586  for(Int_t k=0; k<5; k++)
1587  {
1588  if(k==0)
1589  {
1590  if(mipETmp>cluster->E()) subTrkPtClean[k] += partPt;
1591  else subTrkPtAmbig[k] += partPt;
1592  }
1593  else
1594  {
1595  if(fraction[k]*subEtmp>cluster->E()) subTrkPtClean[k] += partPt;
1596  else subTrkPtAmbig[k] += partPt;
1597  }
1598  }
1599  }
1600  }
1601  }
1602 
1603  Float_t dR = TMath::Sqrt( (partEta-jetEta)*(partEta-jetEta) + (partPhi-jetPhi)*(partPhi-jetPhi) );
1604  for(Int_t i=0; i<nBins; i++)
1605  {
1606  if( dR < radiusCut[i] )
1607  {
1608  eFraction[i] += partE;
1609  nPart[i]++;
1610  }
1611  }
1612  }
1613 
1614  fhLeadNePtInJet[type][r]->Fill(jetPt, leadNePt);
1615  fhLeadChPtInJet[type][r]->Fill(jetPt, leadChPt);
1616  if(type<2 && leadChIndex>-1)
1617  {
1618  Double_t cz = GetZ(constituents[leadChIndex].px(),constituents[leadChIndex].py(),constituents[leadChIndex].pz(),jetsIncl[ij].px(),jetsIncl[ij].py(),jetsIncl[ij].pz());
1619  fhChLeadZVsJetPt[type][r]->Fill(jetPt, cz);
1620  }
1621 
1622  if((fStudySubEInHC||fStudyMcOverSubE) && type<2)
1623  {
1624  for(Int_t i=0; i<5; i++)
1625  {
1626  if(fStudySubEInHC)
1627  {
1628  fhSubClsEVsJetPt[type][r][i]->Fill(jetPt-subClsE[4],subClsE[i]/(jetPt-subClsE[4]));
1629  fhHCTrkPtClean[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtClean[i]/(jetPt-subClsE[4]));
1630  fhHCTrkPtAmbig[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtAmbig[i]/(jetPt-subClsE[4]));
1631  }
1632  if(type==0 && fIsMC && fStudyMcOverSubE)
1633  {
1634  fHCOverSubE[r][i]->Fill(jetPt-subClsE[4],subClsE[i]-mcSubE);
1635  fHCOverSubEFrac[r][i]->Fill(jetPt-subClsE[4],(subClsE[i]-mcSubE)/(jetPt-subClsE[4]));
1636  }
1637  }
1638  }
1639  if(type<2 && chPt>0)
1640  {
1641  for(Int_t i=0; i<3; i++)
1642  {
1643  fRelTrkCon[type][r]->Fill(jetPt,totTrkPt[i]/chPt,i);
1644  }
1645  }
1646 
1647 
1648  for(Int_t ibin=0; ibin<nBins; ibin++)
1649  {
1650  Double_t fill1[3] = {jetPt,radiusCut[ibin]-0.005,eFraction[ibin]/jetE};
1651  fJetEnergyFraction[type][r]->Fill(fill1);
1652  Double_t fill2[3] = {jetPt,radiusCut[ibin]-0.005,static_cast<Double_t>(nPart[ibin])};
1653  fJetNPartFraction[type][r]->Fill(fill2);
1654  }
1655 
1656  // Get the jet pt containing tracks or clusters above some threshold
1657  if(type<2)
1658  {
1659  Double_t thres[3] = {15,25,40};
1660  Int_t okCh[3] = {0,0,0};
1661  Int_t okNe[3] = {0,0,0};
1662  Double_t lowPt[2] = {0.,0.};
1663  for(UInt_t ic=0; ic<constituents.size(); ic++)
1664  {
1665  Float_t partPt = constituents[ic].perp();
1666 
1667  if(partPt<0.3) lowPt[0] += partPt;
1668  else if(partPt<0.5) lowPt[1] += partPt;
1669 
1670  if(constituents[ic].user_index()>0)
1671  {
1672  for(Int_t it=0; it<3; it++)
1673  {
1674  if(partPt>thres[it]) okCh[it] = 1;
1675  }
1676  }
1677  else
1678  {
1679  for(Int_t icl=0; icl<3; icl++)
1680  {
1681  if(partPt>thres[icl]) okNe[icl] = 1;
1682  }
1683  }
1684  }
1685  for(Int_t i=0; i<3; i++)
1686  {
1687  if(okCh[i]==1)
1688  fhJetPtWithTrkThres[type][r]->Fill(i,jetPt);
1689  if(okNe[i]==1)
1690  fhJetPtWithClsThres[type][r]->Fill(i,jetPt);
1691  }
1692  for(Int_t k=0; k<2; k++) fhJetPtVsLowPtCons[type][r][k]->Fill(jetPt,lowPt[k]/jetPt);
1693  }
1694  }
1695 }
1696 
1697 //________________________________________________________________________
1699 {
1700  //
1701  // Analyze secondaries
1702  //
1703 
1704  std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1705  for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1706  {
1707  Float_t jetEta = jetsIncl[ij].eta();
1708  Float_t jetPt = jetsIncl[ij].perp();
1709  if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<1) continue;
1710  std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
1711  Double_t dNKPt = 0, dWPPt = 0, weakPt=0, nkPt=0;
1712  Int_t type = -1;
1713  for(UInt_t ic=0; ic<constituents.size(); ic++)
1714  {
1715  Int_t ipart = TMath::Abs(constituents[ic].user_index())-1;
1716  AliMCParticle* mcParticle = (AliMCParticle*)fMC->GetTrack(ipart);
1717  if(!mcParticle) continue;
1718  Int_t pdg = mcParticle->PdgCode();
1719  if(pdg==310 || (pdg<=3400 && pdg>=3100)) weakPt += mcParticle->Pt();
1720  if(pdg==130 || pdg==2112) nkPt += mcParticle->Pt();
1721 
1722  // Weak decaying particles
1723  if(pdg==310) type=0;
1724  else if (pdg==3122) type=1;
1725  else if (pdg<=3400 && pdg>=3100) type=2;
1726  else type=-1;
1727  if(type>-1)
1728  {
1729  Int_t binx = fhSecondaryResponse[type]->GetXaxis()->FindFixBin(mcParticle->Pt());
1730  TH1F *htmp = (TH1F*)fhSecondaryResponse[type]->ProjectionY(Form("hpro_%d_%d_%d",(Int_t)Entry(),ij,ic),binx,binx);
1731  Double_t pro = htmp->GetRandom();
1732  dWPPt += pro*mcParticle->Pt() - mcParticle->Pt();
1733  delete htmp;
1734  }
1735 
1736  // Missing neutron and K0L
1737  if(pdg==130 || pdg==2112) dNKPt -= mcParticle->Pt();
1738  }
1739 
1740  fhNKFracVsJetPt[etaCut][r]->Fill(jetPt,nkPt/jetPt);
1741  fhWeakFracVsJetPt[etaCut][r]->Fill(jetPt,weakPt/jetPt);
1742 
1743  Double_t jetNewWPPt = jetPt + dWPPt;
1744  fhJetResponseWP[etaCut][r]->Fill(jetPt,jetNewWPPt);
1745  fhJetResolutionWP[etaCut][r]->Fill(jetPt,(jetPt-jetNewWPPt)/jetPt);
1746 
1747  Double_t jetNewNKPt = jetPt + dNKPt;
1748  fhJetResponseNK[etaCut][r]->Fill(jetPt,jetNewNKPt);
1749  fhJetResolutionNK[etaCut][r]->Fill(jetPt,(jetPt-jetNewNKPt)/jetPt);
1750  }
1751 }
1752 
1753 
1754 //________________________________________________________________________
1756 {
1757  //
1758  // Estimate the contribution of missing energies via matching
1759  // type = 0 -> NK
1760  // type = 1 -> WP
1761  //
1762 
1763  if(type!=0 && type!=1) return;
1764 
1765  AliFJWrapper fj("fj","fj");
1766  fj.CopySettingsFrom(*jetFinder);
1767 
1768  Int_t npart = fMC->GetNumberOfTracks();
1769  Int_t pType = -1;
1770  for(Int_t ipart=0; ipart<npart; ipart++)
1771  {
1772  AliVParticle* vParticle = fMC->GetTrack(ipart);
1773  if(!IsGoodMcPartilce(vParticle,ipart)) continue;
1774  Int_t pdgCode = vParticle->PdgCode();
1775  Double_t pt = vParticle->Pt();
1776  if( type==0 && (pdgCode==130 || pdgCode==2112) ) continue; // Reject NK
1777  if(type==1)
1778  {
1779  if(pdgCode==310) pType=0;
1780  else if (pdgCode==3122) pType=1;
1781  else if (pdgCode<=3400 && pdgCode>=3100) pType=2;
1782  else pType=-1;
1783  if(pType>-1)
1784  {
1785  Int_t binx = fhSecondaryResponse[pType]->GetXaxis()->FindFixBin(vParticle->Pt());
1786  TH1F *htmp = (TH1F*)fhSecondaryResponse[pType]->ProjectionY(Form("hpro_%d_%d",(Int_t)Entry(),ipart),binx,binx);
1787  Double_t pro = htmp->GetRandom();
1788  pt = pro * vParticle->Pt();
1789  delete htmp;
1790  }
1791  }
1792  Int_t index = 1;
1793  if(vParticle->Charge()==0) { index=-1; }
1794  fj.AddInputVector(vParticle->Px()*pt/vParticle->Pt(), vParticle->Py()*pt/vParticle->Pt(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
1795  }
1796 
1797  fj.Run();
1798  std::vector<fastjet::PseudoJet> jets_incl_mc = fj.GetInclusiveJets();
1799  std::vector<fastjet::PseudoJet> jets_incl_mc_std = jetFinder->GetInclusiveJets();
1800 
1801  for(UInt_t ij=0; ij<jets_incl_mc.size(); ij++)
1802  {
1803  Float_t jetEta = jets_incl_mc[ij].eta();
1804  Float_t jetPt = jets_incl_mc[ij].perp();
1805  if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<5) continue;
1806  Double_t dEta, dPhi;
1807  Int_t index = FindSpatialMatchedJet(jets_incl_mc[ij], jetFinder, dEta, dPhi, 1);
1808  if(index>-1)
1809  {
1810  Double_t jetPtStd = jets_incl_mc_std[index].perp();
1811  Double_t dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
1812  if(type==0) fhJetResolutionNKSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR);
1813  if(type==1) fhJetResolutionWPSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR);
1814 
1815  if(dR<kdRCut[r])
1816  {
1817  if(type==0) fhJetResponseNKSM[etaCut][r]->Fill(jetPtStd,jetPt);
1818  if(type==1) fhJetResponseWPSM[etaCut][r]->Fill(jetPtStd,jetPt);
1819  }
1820  }
1821  }
1822 }
1823 
1824 //________________________________________________________________________
1825 void AliAnalysisTaskFullppJet::RunAnalyzeUE(AliFJWrapper *jetFinder, const Int_t type, const Bool_t isMCTruth)
1826 {
1827  //
1828  // Run analysis to estimate the underlying event
1829  // Adapted from Oliver
1830 
1831  std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1832  Double_t leadpt=0, leadPhi = -999, leadEta = -999;
1833  Int_t leadIndex = -1;
1834  for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1835  {
1836  Double_t jetEta = jetsIncl[ij].eta();
1837  Double_t jetPt = jetsIncl[ij].perp();
1838  Double_t jetPhi = jetsIncl[ij].phi();
1839  if(leadpt<jetPt)
1840  {
1841  leadpt = jetPt;
1842  leadPhi = jetPhi;
1843  leadEta = jetEta;
1844  leadIndex = ij;
1845  }
1846  }
1847  if(leadpt<1 || TMath::Abs(leadEta)>0.5 ) return;
1848 
1849  Int_t backIndex = -1;
1850  Bool_t isBackToBack = kFALSE;
1851  for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1852  {
1853  Double_t dPhi = TMath::Abs(jetsIncl[ij].phi()-leadPhi);
1854  if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1855  if(dPhi>150*TMath::DegToRad() && jetsIncl[ij].perp()/leadpt > 0.8)
1856  {
1857  backIndex = ij;
1858  isBackToBack = kTRUE;
1859  }
1860  }
1861 
1862  for(Int_t ij=0; ij<(Int_t)jetsIncl.size(); ij++)
1863  {
1864  if(ij==leadIndex || ij==backIndex) continue;
1865  if(TMath::Abs(jetsIncl[ij].eta())>0.5) continue;
1866  if(jetsIncl[ij].perp()>15) isBackToBack = kFALSE;
1867  }
1868 
1869 
1870  Double_t axis[2] = {leadPhi+0.5 * kPI, leadPhi-0.5 * kPI};
1871  if(axis[0]>2*kPI) axis[0] -= 2*kPI;
1872  if(axis[1]<0) axis[1] += 2*kPI;
1873 
1874  fhUEJetPtNorm[type][0][0]->Fill(leadpt);
1875  if(isBackToBack) fhUEJetPtNorm[type][0][1]->Fill(leadpt);
1876 
1877  if(!isMCTruth)
1878  {
1879  Int_t ntracks = fTrackArray->GetEntriesFast();
1880  Double_t vertex[3] = {0, 0, 0};
1881  fESD->GetVertex()->GetXYZ(vertex) ;
1882  TLorentzVector gamma;
1883  Int_t nclusters = fClusterArray->GetEntriesFast();
1884 
1885  for(Int_t j=0; j<2; j++)
1886  {
1887  Double_t ueCh = 0, ueChNe = 0.;
1888  for(Int_t it=0; it<ntracks; it++)
1889  {
1890  AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
1891  if(!t) continue;
1892  Double_t dPhi = TMath::Abs(axis[j]-t->Phi());
1893  Double_t dEta = TMath::Abs(leadEta-t->Eta());
1894  if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1895  if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
1896  {
1897  fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,t->Pt());
1898  if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,t->Pt());
1899  ueCh += t->Pt();
1900  if(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4)
1901  {
1902  fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,t->Pt());
1903  if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,t->Pt());
1904  ueChNe += t->Pt();
1905  }
1906  }
1907  }
1908  fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh);
1909  if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh);
1910 
1911  if(!(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4)) continue;
1912  fhUEJetPtNorm[type][1][0]->Fill(leadpt);
1913  if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt);
1914  for(Int_t ic=0; ic<nclusters; ic++)
1915  {
1916  AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic);
1917  if(!cl) continue;
1918  cl->GetMomentum(gamma, vertex);
1919  Double_t clsPhi = gamma.Phi();
1920  if(clsPhi<0) clsPhi += 2*kPI;
1921  Double_t dPhi = TMath::Abs(axis[j]-clsPhi);
1922  Double_t dEta = TMath::Abs(leadEta-gamma.Eta());
1923  if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1924  if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
1925  {
1926  ueChNe += gamma.Pt();
1927  fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,gamma.Pt());
1928  if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,gamma.Pt());
1929  }
1930  }
1931  fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe);
1932  if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe);
1933  }
1934  }
1935  else
1936  {
1937  fhUEJetPtNorm[type][1][0]->Fill(leadpt);
1938  if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt);
1939  Int_t npart = fMcPartArray->GetSize();
1940  for(Int_t j=0; j<2; j++)
1941  {
1942  Double_t ueCh = 0, ueChNe = 0., ueChNe2 = 0.;
1943  for(Int_t ipos=0; ipos<npart; ipos++)
1944  {
1945  AliVParticle* vParticle = fMC->GetTrack(fMcPartArray->At(ipos));
1946  if(!vParticle) continue;
1947  Double_t dPhi = TMath::Abs(axis[j]-vParticle->Phi());
1948  Double_t dEta = TMath::Abs(leadEta-vParticle->Eta());
1949  if(dPhi > kPI) dPhi = 2*kPI - dPhi;
1950  if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
1951  {
1952  Double_t pt = vParticle->Pt();
1953  ueChNe += pt;
1954  fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,pt);
1955  if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,pt);
1956  if( vParticle->Charge()!=0 )
1957  {
1958  fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,pt);
1959  if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,pt);
1960  ueCh += pt;
1961  ueChNe2 += pt;
1962  }
1963  else
1964  {
1965  if(pt>0.5) ueChNe2 += pt;
1966  }
1967  }
1968  }
1969  fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh);
1970  fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe);
1971  if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh);
1972  if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe);
1973  }
1974  }
1975 }
1976 
1977 //________________________________________________________________________
1979 {
1980  //
1981  // Check if it is a good jet
1982  //
1983 
1984  if(jet->Pt()<1) return kFALSE;
1985  if(TMath::Abs(jet->Eta())>(0.7-rad)) return kFALSE;
1986  if(jet->Phi() < (80*TMath::DegToRad()+rad) || jet->Phi() > (180*TMath::DegToRad()-rad) ) return kFALSE;
1987  return kTRUE;
1988 }
1989 
1990 
1991 //________________________________________________________________________
1993 {
1994  //
1995  // Get the leading z of the input jet
1996  //
1997 
1998  Double_t z = 0;
1999  std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
2000  std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
2001  Int_t index = -1;
2002  Double_t maxPt = 0;
2003  for(UInt_t ic=0; ic<constituents.size(); ic++)
2004  {
2005  if(constituents[ic].perp()>maxPt)
2006  {
2007  maxPt = constituents[ic].perp();
2008  index = ic;
2009  }
2010  }
2011  if(index>-1)
2012  z = GetZ(constituents[index].px(),constituents[index].py(),constituents[index].pz(),jetsIncl[jetIndex].px(),jetsIncl[jetIndex].py(),jetsIncl[jetIndex].pz());
2013  return z;
2014 }
2015 
2016 //________________________________________________________________________
2017 Double_t AliAnalysisTaskFullppJet::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
2018 {
2019  //
2020  // Get the z of a constituent inside of a jet
2021  //
2022 
2023  return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
2024 }
2025 
2026 //________________________________________________________________________
2028 {
2029  //
2030  // Get jet energy resoultion due to intrinsic detector effects
2031  //
2032 
2033  Double_t jetSigma2 = 0;
2034  std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
2035  for(UInt_t ic=0; ic<constituents.size(); ic++)
2036  {
2037  if(constituents[ic].user_index()>0)
2038  {
2039  AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
2040  jetSigma2 += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2041  }
2042  else
2043  {
2044  AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
2045  jetSigma2 += TMath::Power(cluster->GetTOF(),2);
2046  }
2047  }
2048  return TMath::Sqrt(jetSigma2);
2049 }
2050 
2051 
2052 
2053 //________________________________________________________________________
2055 {
2056  //
2057  // Correct the tracking inefficiency explicitly
2058  //
2059 
2060  Double_t misspt = 0;
2061  std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
2062  Double_t jetPt = jetsIncl[jetIndex].perp();
2063  if(!fhCorrTrkEffPtBin[fTriggerType][radiusIndex])
2064  {
2065  printf("Warning: can't get the mean # of tracks per jet with pt=%f in: trigger=%d, radiusIndex=%d\n",jetPt,fTriggerType,radiusIndex);
2066  return 0;
2067  }
2068  Int_t ibin = fhCorrTrkEffPtBin[fTriggerType][radiusIndex]->FindFixBin(jetPt);
2069  if(!fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1] || fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1]->Integral()<0.001)
2070  {
2071  printf("Warning: no sampling distrubtion for jet with pt=%f\n",jetPt);
2072  return 0;
2073  }
2074  Double_t ntrack = fhCorrTrkEffPtBin[fTriggerType][radiusIndex]->GetBinContent(ibin);
2075  Int_t nTrk = (Int_t) ntrack;
2076  Double_t res = ntrack-nTrk;
2077  Double_t pro1 = fRandomGen->Uniform();
2078  if(pro1<res) nTrk++;
2079  for(Int_t itry=0; itry<nTrk; itry++)
2080  {
2081  Double_t trkPt = fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1]->GetRandom();
2082  if(trkPt/jetPt>fTrkEffCorrCutZ) continue;
2083  Double_t eff = GetTrkEff(trkPt);
2084  Double_t pro = fRandomGen->Uniform();
2085  if(pro>eff) misspt += trkPt;
2086  }
2087  return misspt;
2088 }
2089 
2090 
2091 //________________________________________________________________________
2092 Bool_t AliAnalysisTaskFullppJet::IsGoodMcPartilce(const AliVParticle* vParticle, const Int_t ipart)
2093 {
2094  //
2095  // Select good primary particles to feed into jet finder
2096  //
2097 
2098  if(!vParticle) return kFALSE;
2099  if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;
2100  if (TMath::Abs(vParticle->Eta())>2) return kFALSE;
2101  return kTRUE;
2102 }
2103 
2104 //________________________________________________________________________
2105 Int_t AliAnalysisTaskFullppJet::FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR)
2106 {
2107  //
2108  // Find spatially matched detector and particle jets
2109  //
2110 
2111  Int_t index=-1;
2112  dEta=-999, dPhi=-999;
2113  std::vector<fastjet::PseudoJet> jets_incl = jetFinder->GetInclusiveJets();
2114  for(UInt_t ij=0; ij<jets_incl.size(); ij++)
2115  {
2116  if(jets_incl[ij].perp()<5) continue;
2117  if(TMath::Abs(jets_incl[ij].eta())>1) continue;
2118  Double_t tmpR = TMath::Sqrt( TMath::Power(jet.eta()-jets_incl[ij].eta(), 2) + TMath::Power(jet.phi()-jets_incl[ij].phi(), 2) );
2119  if(tmpR<maxR)
2120  {
2121  maxR=tmpR;
2122  index=ij;
2123  dEta = jet.eta()-jets_incl[ij].eta();
2124  dPhi = jet.phi()-jets_incl[ij].phi();
2125  }
2126  }
2127  return index;
2128 }
2129 
2130 //________________________________________________________________________
2131 Int_t AliAnalysisTaskFullppJet::FindEnergyMatchedJet(AliFJWrapper *jetFinder1, const Int_t index1, AliFJWrapper *jetFinder2, const Double_t fraction)
2132 {
2133  //
2134  // Find matched detector and particle jets based on shared constituents
2135  //
2136 
2137  Int_t matchedIndex = -1;
2138  std::vector<fastjet::PseudoJet> jetsIncl1 = jetFinder1->GetInclusiveJets();
2139  std::vector<fastjet::PseudoJet> jetsIncl2 = jetFinder2->GetInclusiveJets();
2140  std::vector<fastjet::PseudoJet> constituents1 = jetFinder1->GetJetConstituents(index1);
2141  Double_t jetPt1 = jetsIncl1[index1].perp();
2142  if(jetPt1<0) return matchedIndex;
2143 
2144  for(UInt_t ij2=0; ij2<jetsIncl2.size(); ij2++)
2145  {
2146  Double_t jetPt2 = jetsIncl2[ij2].perp();
2147  if(jetPt2<0) return matchedIndex;
2148  std::vector<fastjet::PseudoJet> constituents2 = jetFinder2->GetJetConstituents(ij2);
2149  Double_t sharedPt1 = 0., sharedPt2 = 0.;
2150  for(UInt_t ic2=0; ic2<constituents2.size(); ic2++)
2151  {
2152  Int_t mcLabel = constituents2[ic2].user_index()-1;
2153  Double_t consPt2 = constituents2[ic2].perp();
2154  for(UInt_t ic1=0; ic1<constituents1.size(); ic1++)
2155  {
2156  Double_t consPt1 = constituents1[ic1].perp();
2157  if(constituents1[ic1].user_index()>0)
2158  {
2159  AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents1[ic1].user_index()-1);
2160  if(track->GetLabel()==mcLabel) {sharedPt2 += consPt2; sharedPt1 += consPt1; cout<<"Found a matched track"<<endl;break;}
2161  }
2162  else
2163  {
2164  AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents1[ic1].user_index()*(-1)-1);
2165  if(cluster->GetLabel()==mcLabel) {sharedPt2 += consPt2; sharedPt1 += consPt1; cout<<"Found a matched cluster"<<endl;break;}
2166  }
2167  }
2168  }
2169  cout<<sharedPt1/jetPt1<<" "<<sharedPt2/jetPt2<<endl;
2170  if(sharedPt1/jetPt1 > fraction && sharedPt2/jetPt2 > fraction)
2171  {
2172  matchedIndex = ij2;
2173  break;
2174  }
2175  }
2176  return matchedIndex;
2177 }
2178 
2179 //________________________________________________________________________
2181 {
2182  //
2183  // Get esd tracks.
2184  //
2185 
2186  Int_t nTrax = fESD->GetNumberOfTracks();
2187  if(fVerbosity>5) printf("[i] # of tracks in event: %d\n",nTrax);
2188 
2189  for (Int_t i = 0; i < nTrax; ++i)
2190  {
2191  AliESDtrack* esdtrack = fESD->GetTrack(i);
2192  if (!esdtrack)
2193  {
2194  AliError(Form("Couldn't get ESD track %d\n", i));
2195  continue;
2196  }
2197 
2198  AliESDtrack *newtrack = GetAcceptTrack(esdtrack);
2199  if(!newtrack) continue;
2200 
2201  Double_t pt = newtrack->Pt();
2202  Double_t eta = newtrack->Eta();
2203  if (pt < fTrkPtMin[fTriggerType] || pt > fTrkPtMax[fTriggerType] || TMath::Abs(eta) > fTrkEtaMax)
2204  {
2205  delete newtrack;
2206  continue;
2207  }
2208 
2209  if(fSysTrkEff)
2210  {
2211  Double_t rand = fRandomGen->Uniform();
2212  if(rand<fVaryTrkEff)
2213  {
2214  delete newtrack;
2215  continue;
2216  }
2217  }
2218  newtrack->SetIntegratedLength(esdtrack->Pt());
2219  fTrackArray->Add(newtrack);
2220  }
2221  if(fVerbosity>5) printf("[i] # of tracks in event: %d\n", fTrackArray->GetEntries());
2222  fNTracksPerChunk += fTrackArray->GetEntries();
2223 }
2224 
2225 //
2226 //________________________________________________________________________
2227 //
2228 AliESDtrack *AliAnalysisTaskFullppJet::GetAcceptTrack(AliESDtrack *esdtrack)
2229 {
2230  //
2231  // Get the hybrid tracks
2232  //
2233 
2234  AliESDtrack *newTrack = 0x0;
2235 
2236  if(fTrackCutsType==0 || fTrackCutsType==3)
2237  {
2238  if(fEsdTrackCuts->AcceptTrack(esdtrack))
2239  {
2240  newTrack = new AliESDtrack(*esdtrack);
2241  newTrack->SetTRDQuality(0);
2242  }
2243  else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
2244  {
2245  if(esdtrack->GetConstrainedParam())
2246  {
2247  newTrack = new AliESDtrack(*esdtrack);
2248  const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
2249  newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
2250  newTrack->SetTRDQuality(1);
2251  }
2252  else
2253  return 0x0;
2254  }
2255  else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
2256  {
2257  if(esdtrack->GetConstrainedParam())
2258  {
2259  newTrack = new AliESDtrack(*esdtrack);
2260  const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
2261  newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
2262  newTrack->SetTRDQuality(2);
2263  }
2264  else
2265  return 0x0;
2266  }
2267  else
2268  {
2269  return 0x0;
2270  }
2271  }
2272  else if(fTrackCutsType==1)
2273  {
2274  if(fEsdTrackCuts->AcceptTrack(esdtrack))
2275  {
2276  newTrack = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());// use TPC only tracks with non default SPD vertex
2277  if(!newTrack) return 0x0;
2278  AliExternalTrackParam exParam;
2279  Bool_t relate = newTrack->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam); //constrain to SPD vertex
2280  if( !relate )
2281  {
2282  delete newTrack;
2283  return 0x0;
2284  }
2285  newTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
2286  newTrack->SetTRDQuality(1);
2287  }
2288  else
2289  {
2290  return 0x0;
2291  }
2292  }
2293  else if (fTrackCutsType==2)
2294  {
2295  if(fEsdTrackCuts->AcceptTrack(esdtrack))
2296  {
2297  newTrack = new AliESDtrack(*esdtrack);
2298  newTrack->SetTRDQuality(0);
2299  }
2300  else
2301  return 0x0;
2302  }
2303  else
2304  {
2305  printf("Unknown track cuts type: %d\n",fTrackCutsType);
2306  return 0x0;
2307  }
2308 
2309  return newTrack;
2310 }
2311 
2312 //________________________________________________________________________
2314 {
2315  //
2316  // Get emcal clusters - selected
2317  //
2318 
2319  if(fSysTrkClsMth)
2320  {
2321  if (!TGeoGlobalMagField::Instance()->GetField()) fESD->InitMagneticField();
2322  fRecoUtil->FindMatches(fESD,0x0,fGeom);
2323  fRecoUtil->SetClusterMatchedToTrack(fESD);
2324  fRecoUtil->SetTracksMatchedToCluster(fESD);
2325  }
2326 
2327  const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
2328  for(Int_t i = 0 ; i < nCaloClusters; i++)
2329  {
2330  AliESDCaloCluster * cl = (AliESDCaloCluster *) fESD->GetCaloCluster(i);
2331  if(!IsGoodCluster(cl)) continue;
2332  AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cl);
2333  if (newCluster->E() < fClsEtMin[fTriggerType] || newCluster->E() > fClsEtMax[fTriggerType])
2334  {delete newCluster; continue;}
2335 
2336  // Absolute scale
2337  if(fPeriod.Contains("lhc11a",TString::kIgnoreCase))
2338  {
2339  Double_t newE = newCluster->E() * 1.02;
2340  newCluster->SetE(newE);
2341  }
2342 
2343  // Trigger efficiency systematic uncertainty
2344  if(fSysJetTrigEff)
2345  {
2346  Double_t newE = newCluster->E() * (1+fVaryJetTrigEff);
2347  newCluster->SetE(newE);
2348  if(fTriggerType==0) fhClsE[fTriggerType]->Fill(newCluster->E());
2349  if(fTriggerType==1)
2350  {
2351  if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase)) fhClsE[0]->Fill(newCluster->E());
2352  if(newCluster->Chi2()>0.5) fhClsE[fTriggerType]->Fill(newCluster->E());
2353  }
2354  }
2355 
2356  // Cluster energy scale systematic uncertainty
2357  if(fSysClusterEScale)
2358  {
2359  Double_t newE = newCluster->E() * (1+fVaryClusterEScale);
2360  newCluster->SetE(newE);
2361  }
2362 
2363  // Cluster energy resolution systematic uncertainty
2364  if(fSysClusterERes)
2365  {
2366  Double_t oldE = newCluster->E();
2367  Double_t resolution = fClusterEResolution->Eval(oldE);
2368  Double_t smear = resolution * TMath::Sqrt((1+fVaryClusterERes)*(1+fVaryClusterERes)-1);
2369  Double_t newE = oldE + fRandomGen->Gaus(0, smear);
2370  newCluster->SetE(newE);
2371  }
2372 
2373  // non-linearity systematic uncertainty
2374  if(fSysNonLinearity)
2375  {
2376  Double_t oldE = newCluster->E();
2377  Double_t newE = oldE/fNonLinear->Eval(oldE) * 1.012;
2378  newCluster->SetE(newE);
2379  }
2380 
2381  // clusterizer systematic uncertainty
2382  if(fSysClusterizer)
2383  {
2384  fhSysClusterE[fTriggerType][0]->Fill(newCluster->E());
2385  fhSysNCellVsClsE[fTriggerType][0]->Fill(newCluster->E(),newCluster->GetNCells());
2386  }
2387 
2388  // hadronic correction
2389  Double_t clsE = newCluster->E();
2390  Double_t subE = 0., eRes = 0., mcSubE = 0, MIPE = 0.;
2392  subE= SubtractClusterEnergy(newCluster, eRes, MIPE, mcSubE);
2393 
2394  if(!fStudySubEInHC && !fStudyMcOverSubE) clsE -= subE;
2395  if(clsE<0) {delete newCluster; continue;}
2396  newCluster->SetE(clsE);
2397  newCluster->SetTOF(eRes);
2398  newCluster->SetDispersion(subE);
2399  newCluster->SetEmcCpvDistance(MIPE);
2400  newCluster->SetDistanceToBadChannel(mcSubE);
2401  fClusterArray->Add(newCluster);
2402 
2403  // clusterizer systematic uncertainty
2404  if(fSysClusterizer)
2405  {
2406  fhSysClusterE[fTriggerType][1]->Fill(newCluster->E());
2407  fhSysNCellVsClsE[fTriggerType][1]->Fill(newCluster->E(),newCluster->GetNCells());
2408  }
2409  }
2410 
2411  if(fVerbosity>5) printf("[i] # of EMCal clusters in event: %d\n", fClusterArray->GetEntries());
2412 
2413 }
2414 
2415 
2416 //________________________________________________________________________
2418 {
2419  //
2420  // Select good clusters
2421  //
2422 
2423  if(!cluster) return kFALSE;
2424  if (!cluster->IsEMCAL()) return kFALSE;
2425  if(fRejectExoticCluster && fRecoUtil->IsExoticCluster(cluster, (AliVCaloCells*)fESD->GetEMCALCells())) return kFALSE;
2426  if(fRemoveBadChannel && fRecoUtil->ClusterContainsBadChannel(fGeom, cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE;
2427  return kTRUE;
2428 }
2429 
2430 
2431 //________________________________________________________________________
2432 Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *cluster, Double_t &eRes, Double_t &MIPE, Double_t &mcSubE)
2433 {
2434  //
2435  // Hadronic correction
2436  //
2437 
2438  mcSubE = 0;
2439  eRes = 0;
2440  MIPE = 0;
2441 
2442  eRes += TMath::Power(fClusterEResolution->Eval(cluster->E()),2);
2443  Double_t subE = 0., sumTrkPt = 0., sumTrkP = 0.;
2444  Int_t nTrack = 0;
2445  TArrayI *matched = cluster->GetTracksMatched();
2446  if(!matched) return 0;
2447  for(Int_t im=0; im<matched->GetSize(); im++)
2448  {
2449  Int_t trkIndex = matched->At(im);
2450  if(trkIndex<0 || trkIndex>=fESD->GetNumberOfTracks()) continue;
2451  Bool_t isSelected = kFALSE;
2452  Int_t index = -1;
2453  for(Int_t j=0; j<fTrackArray->GetEntriesFast(); j++)
2454  {
2455  AliESDtrack *tr = (AliESDtrack*)fTrackArray->At(j);
2456  if( trkIndex == tr->GetID() )
2457  {
2458  isSelected = kTRUE;
2459  index = j;
2460  break;
2461  }
2462  }
2463  if(!isSelected) continue;
2464  nTrack++;
2465  AliESDtrack *track = (AliESDtrack*)fTrackArray->At(index);
2466  Double_t trkP = track->P();
2467  sumTrkPt += track->Pt(); sumTrkP += track->P();
2468  if(fSysTrkPtRes) trkP = TMath::Sqrt(TMath::Power(track->Pt()+GetSmearedTrackPt(track),2)+TMath::Power(track->Pz(),2));
2469  if(IsElectron(track,cluster->E())) //electrons
2470  {
2471  if(fElectronRejection)
2472  {
2473  subE+= trkP;
2474  MIPE+= trkP;
2475  eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2476  }
2477  }
2478  else //hadrons
2479  {
2481  {
2482  MIPE += (trkP>0.27)?0.27:trkP; //MIP correction
2483  if(fFractionHC>2)
2484  {
2485  if(trkP>0.27)
2486  {
2487  subE+=0.27;
2488  }
2489  else
2490  {
2491  subE+=trkP;
2492  eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2493  }
2494  }
2495  else
2496  {
2497  if(trkP>fHCLowerPtCutMIP) subE += 0.27;
2498  else subE+=fFractionHC*trkP;
2499  eRes += TMath::Power(fFractionHC*track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2500  }
2501  }
2502  }
2503  }
2504 
2505  if(fSaveQAHistos) fhNMatchedTrack[fTriggerType]->Fill(nTrack);
2506  if(nTrack>0)
2507  {
2508  Double_t fraction[4] = {0.3,0.5,0.7,1.0};
2509  for(Int_t j=0; j<4; j++)
2510  {
2511  Double_t subETmp = sumTrkP*fraction[j];
2512  if(subETmp>cluster->E()) subETmp = cluster->E();
2513  if(fSaveQAHistos) fhSubEVsTrkPt[fTriggerType][j]->Fill(sumTrkPt,subETmp/sumTrkP);
2514  }
2515  }
2516 
2517  eRes = TMath::Sqrt(eRes);
2518 
2519  if(fIsMC && nTrack>0)
2520  {
2521  Double_t neutralE = 0;
2522  TArrayI* labels = cluster->GetLabelsArray();
2523  if(labels)
2524  {
2525  for(Int_t il=0; il<labels->GetSize(); il++)
2526  {
2527  Int_t ipart = labels->At(il);
2528  if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
2529  {
2530  AliVParticle* vParticle = fMC->GetTrack(ipart);
2531  if(vParticle->Charge()==0)
2532  {
2533  neutralE += vParticle->E();
2534  }
2535  }
2536  }
2537  }
2538  mcSubE = cluster->E() - neutralE;
2539  if(mcSubE<0)
2540  mcSubE=0;
2541  }
2542 
2543 
2544  return subE;
2545 }
2546 
2547 
2548 //
2549 //________________________________________________________________________
2550 //
2552 {
2553  //
2554  // Exotic fraction: f_cross
2555  // Adpated from AliEMCalRecoUtils
2556 
2557  if(!cluster) return -1;
2558  AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
2559  if(!cells) return -1;
2560 
2561  // Get highest energy tower
2562  Int_t iSupMod = -1, absID = -1, ieta = -1, iphi = -1,iTower = -1, iIphi = -1, iIeta = -1;
2563  Bool_t share = kFALSE;
2564  fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share);
2565  fGeom->GetCellIndex(absID,iSupMod,iTower,iIphi,iIeta);
2566  fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2567 
2568  Int_t absID1 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi+1, ieta);
2569  Int_t absID2 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi-1, ieta);
2570  Int_t absID3 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi, ieta+1);
2571  Int_t absID4 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi, ieta-1);
2572 
2573  Float_t ecell = 0, ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
2574  Double_t tcell = 0, tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0;
2575  Bool_t accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0;
2576  const Int_t bc = 0;
2577 
2578  accept = fRecoUtil->AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
2579 
2580  if(!accept) return -1;
2581 
2582  if(ecell < 0.5) return -1;
2583 
2584  accept1 = fRecoUtil->AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells);
2585  accept2 = fRecoUtil->AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells);
2586  accept3 = fRecoUtil->AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells);
2587  accept4 = fRecoUtil->AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells);
2588 
2589  const Double_t exoticCellDiffTime = 1e6;
2590  if(TMath::Abs(tcell-tcell1)*1.e9 > exoticCellDiffTime) ecell1 = 0 ;
2591  if(TMath::Abs(tcell-tcell2)*1.e9 > exoticCellDiffTime) ecell2 = 0 ;
2592  if(TMath::Abs(tcell-tcell3)*1.e9 > exoticCellDiffTime) ecell3 = 0 ;
2593  if(TMath::Abs(tcell-tcell4)*1.e9 > exoticCellDiffTime) ecell4 = 0 ;
2594 
2595  Float_t eCross = ecell1+ecell2+ecell3+ecell4;
2596 
2597  return 1-eCross/ecell;
2598 }
2599 
2600 
2601 //________________________________________________________________________
2603 {
2604  //
2605  // Get # of trials per ESD event
2606  //
2607 
2608  AliStack *stack = fMC->Stack();
2609  if(stack)
2610  {
2611  AliGenEventHeader *head = dynamic_cast<AliGenEventHeader*>(fMC->GenEventHeader());
2612  if (head == 0x0)
2613  {
2614  AliError("Could not get the event header");
2615  return;
2616  }
2617 
2618  AliGenPythiaEventHeader *headPy = dynamic_cast<AliGenPythiaEventHeader*>(head);
2619  if (headPy != 0x0)
2620  {
2621  if (headPy->Trials() > 0)
2622  {
2623  fhNTrials[0]->Fill(0.5,headPy->Trials());
2624 
2625  }
2626  }
2627  }
2628 }
2629 
2630 
2631 
2632 //________________________________________________________________________
2634 {
2635  //
2636  // Called once at the end of the query
2637  //
2638  Info("Terminate","Terminate");
2639  AliAnalysisTaskSE::Terminate();
2640 
2641 }
2642 
2643 //
2644 //________________________________________________________________________
2645 //
2647 {
2648  //
2649  // Run trigger offline
2650  //
2651 
2652  fIsEventTriggerBit = 0;
2653  Int_t isTrigger = 0;
2654  Int_t ncl = fESD->GetNumberOfCaloClusters();
2655  for(Int_t icl=0; icl<ncl; icl++)
2656  {
2657  // Check every cluster
2658  AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2659  if(!IsGoodCluster(cluster)) continue;
2660  Double_t pro = GetOfflineTriggerProbability(cluster);
2661  Double_t rand = fRandomGen->Uniform();
2662  if(rand<pro)
2663  {
2664  isTrigger = 1;
2665  fIsEventTriggerBit = 1;
2666  cluster->SetChi2(1);
2667  }
2668  else
2669  cluster->SetChi2(0);
2670  }
2671  return isTrigger;
2672 }
2673 
2674 
2675 //
2676 //________________________________________________________________________
2677 //
2679 {
2680  //
2681  // Get the probablity of the given cluster to trigger the event
2682  //
2683 
2684  Double_t pro = 0;
2685  // Check the trigger mask
2686  AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
2687  Int_t iSupMod = -1, absID = -1, ieta = -1, iphi = -1,iTower = -1, iIphi = -1, iIeta = -1;
2688  Bool_t share = kFALSE;
2689  fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share); // Get the position of the most energetic cell in the cluster
2690  fGeom->GetCellIndex(absID,iSupMod,iTower,iIphi,iIeta);
2691  fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
2692  // convert co global phi eta
2693  Int_t gphi = iphi + 24*(iSupMod/2);
2694  Int_t geta = ieta + 48*(iSupMod%2);
2695  // get corresponding FALTRO
2696  Int_t fphi = gphi / 2;
2697  Int_t feta = geta / 2;
2698  if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5 && iSupMod>-1) // check the trigger mask
2699  {
2700  Double_t clsE = cluster->E();
2701  if(fSysClusterEScale) clsE = clsE * (1+fVaryClusterEScale); // Used for systematic uncertainty. Not needed for regular analysis
2702  if(clsE>10) pro = 1; // Probability is 1 at high E
2703  else
2704  {
2705  Int_t bin = fTriggerCurve[iSupMod]->FindFixBin(clsE);
2706  pro = fTriggerCurve[iSupMod]->GetBinContent(bin)/fTriggerEfficiency[iSupMod]->Eval(10); // Read the probability from trigger turn-on curves
2707  }
2708  }
2709  return pro;
2710 }
2711 
2712 
2713 
2714 //
2715 //________________________________________________________________________
2716 //
2718 {
2719  //
2720  // Return the given cluster supermodule
2721  //
2722 
2723  Float_t pos[3];
2724  cluster->GetPosition(pos);
2725  TVector3 clsVec(pos[0],pos[1],pos[2]);
2726 
2727  Int_t sMod=-1;
2728  fGeom->SuperModuleNumberFromEtaPhi(clsVec.Eta(), clsVec.Phi(), sMod);
2729  return sMod;
2730 }
2731 
2732 
2733 //________________________________________________________________________
2735 {
2736  //
2737  // Check if the given track is an electron candidate based on de/dx
2738  //
2739 
2740  if(track->GetTPCsignal()<=fdEdxMax && track->GetTPCsignal()>=fdEdxMin && (clsE/track->P())<=fEoverPMax && (clsE/track->P())>=fEoverPMin )
2741  return kTRUE;
2742  else
2743  return kFALSE;
2744 }
2745 
2746 
2747 //________________________________________________________________________
2749 {
2750  //
2751  // Get tracking efficiency estimated from simulation
2752  //
2753 
2754  Double_t eff = 1;
2755  Double_t ptBound[4] = {0, 0.5, 3.8, 300};
2756 
2757  for(Int_t i=0; i<3; i++)
2758  {
2759  if( inPt < ptBound[i+1] && inPt >= ptBound[i])
2760  {
2761  eff = fTrkEffFunc[i]->Eval(inPt);
2762  break;
2763  }
2764  }
2765  return eff;
2766 }
2767 
2768 //________________________________________________________________________
2770 {
2771  //
2772  // Smear track momentum
2773  //
2774 
2775  Double_t resolution = track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2());
2776  Double_t smear = resolution*TMath::Sqrt((1+fVaryTrkPtRes)*(1+fVaryTrkPtRes)-1);
2777  return fRandomGen->Gaus(0, smear);
2778 
2779 }
2780 
2781 
2782 //________________________________________________________________________
2784 {
2785  //
2786  // Check if the triggered events have correct trigger bit
2787  //
2788 
2789  fIsEventTriggerBit = 0;
2790  // constants
2791  const Int_t nColsModule = 2;
2792  const Int_t nRowsModule = 5;
2793  const Int_t nRowsFeeModule = 24;
2794  const Int_t nColsFeeModule = 48;
2795  const Int_t nColsFaltroModule = 24;
2796  const Int_t nRowsFaltroModule = 12;
2797 
2798  // part 1, trigger extraction -------------------------------------
2799  Int_t trigtimes[30], globCol, globRow, ntimes;
2800  Int_t trigger[nColsFaltroModule*nColsModule][nRowsFaltroModule*nRowsModule];
2801 
2802  // erase trigger maps
2803  for( Int_t i = 0; i < nColsFaltroModule*nColsModule; i++ )
2804  {
2805  for( Int_t j = 0; j < nRowsFaltroModule*nRowsModule; j++ )
2806  {
2807  trigger[i][j] = 0;
2808  }
2809  }
2810 
2811  Int_t fTrigCutLow = 7, fTrigCutHigh = 10, trigInCut = 0;
2812  AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
2813  // go through triggers
2814  if( fCaloTrigger->GetEntries() > 0 )
2815  {
2816  // needs reset
2817  fCaloTrigger->Reset();
2818  while( fCaloTrigger->Next() )
2819  {
2820  fCaloTrigger->GetPosition( globCol, globRow ); // get position in global 2x2 tower coordinates
2821  fCaloTrigger->GetNL0Times( ntimes ); // get dimension of time arrays
2822  if( ntimes < 1 ) continue; // no L0s in this channel. Presence of the channel in the iterator still does not guarantee that L0 was produced!!
2823  fCaloTrigger->GetL0Times( trigtimes ); // get timing array
2824  if(fCheckTriggerMask && fTriggerMask->GetBinContent(globCol+1,globRow+1)<0.5) continue;
2825  trigInCut = 0;
2826  for( Int_t i = 0; i < ntimes; i++ )
2827  {
2828  if( trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh ) // check if in cut
2829  {
2830  trigInCut = 1;
2831  }
2832  }
2833  if(trigInCut==1) trigger[globCol][globRow] = 1;
2834  } // calo trigger entries
2835  } // has calo trigger entries
2836 
2837 
2838 
2839  // part 2 go through the clusters here -----------------------------------
2840  Int_t nCell, iCell;
2841  UShort_t *cellAddrs;
2842  Int_t absID = -1, nSupMod=-1, nModule=-1, nIphi=-1, nIeta=-1, iphi=-1, ieta=-1, gphi=-1, geta=-1, feta=-1, fphi=-1;
2843  Bool_t share = kFALSE;
2844  AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
2845  for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2846  {
2847  AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2848  if(!cluster || !cluster->IsEMCAL()) continue;
2849  cluster->SetChi2(0);
2850  if(!IsGoodCluster(cluster)) continue;
2851 
2852  //Clusters with most energetic cell in dead region can't be triggered
2853  fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, nSupMod, ieta, iphi, share);
2854  fGeom->GetCellIndex(absID,nSupMod,nModule, nIphi, nIeta);
2855  fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2856  gphi = iphi + nRowsFeeModule*(nSupMod/2);
2857  geta = ieta + nColsFeeModule*(nSupMod%2);
2858  fphi = gphi / 2;
2859  feta = geta / 2;
2860 
2861  if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5)
2862  {
2863  nCell = cluster->GetNCells(); // get cluster cells
2864  cellAddrs = cluster->GetCellsAbsId(); // get the cell addresses
2865  for( iCell = 0; iCell < nCell; iCell++ )
2866  {
2867  // get cell position
2868  fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2869  fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2870 
2871  // convert co global phi eta
2872  gphi = iphi + nRowsFeeModule*(nSupMod/2);
2873  geta = ieta + nColsFeeModule*(nSupMod%2);
2874 
2875  // get corresponding FALTRO
2876  fphi = gphi / 2;
2877  feta = geta / 2;
2878  // try to match with a triggered
2879  if( trigger[feta][fphi] == 1)
2880  {
2881  cluster->SetChi2(1);
2882  fIsEventTriggerBit = 1;
2883  //break;
2884  }
2885  } // cells
2886  }
2887  } // clusters
2888 }
2889 
2890 
2891 //________________________________________________________________________
2893 {
2894  //
2895  // Print configuration
2896  //
2897 
2898  const char *trackCutName[3] = {"Hybrid tracks","TPCOnly tracks","Golden tracks"};
2899  const char *triggerType[2] = {"MB","EMC"};
2900  const char *decision[2] = {"no","yes"};
2901  const char *recombination[] = {"E_scheme","pt_scheme","pt2_scheme","Et_scheme","Et2_scheme","BIpt_scheme","BIpt2_scheme"};
2902  const char *type[2] = {"local","grid"};
2904  {
2905  printf("\n\n=================================\n");
2906  printf("======WARNING: HC is ingored!======\n");
2907  printf("====== NOT for PHYSICS! ======\n\n");
2908  }
2909  printf("Run period: %s\n",fPeriod.Data());
2910  printf("Reject SPD pileup: %s\n",decision[fRejectPileup]);
2911  printf("Reject exotic triggered events: %s\n",decision[fRejectExoticTrigger]);
2912  printf("Is this MC data: %s\n",decision[fIsMC]);
2913  printf("Only find charged jets in MC: %s\n",decision[fChargedMC]);
2914  printf("Analyze on local or grid? %s\n",type[fAnaType]);
2915  printf("Run offline trigger on MC: %s\n",decision[fOfflineTrigger]);
2916  if(fIsMC)
2917  printf("Is K0 and n included: %s\n",decision[1-fRejectNK]);
2918  printf("Constrain tracks in EMCal acceptance: %s\n",decision[fConstrainChInEMCal]);
2919  printf("Track selection: %s, |eta| < %2.1f\n",trackCutName[fTrackCutsType], fTrkEtaMax);
2920  for(Int_t i=0; i<2; i++)
2921  {
2922  printf("Track pt cut: %s -> %2.2f < pT < %2.1f\n",triggerType[i], fTrkPtMin[i], fTrkPtMax[i]);
2923  }
2924  for(Int_t i=0; i<2; i++)
2925  {
2926  printf("Cluster selection: %s -> %2.2f < Et < %2.1f\n",triggerType[i],fClsEtMin[i], fClsEtMax[i]);
2927  }
2928  printf("Electron selectoin: %2.0f < dE/dx < %2.0f, %1.1f < E/P < %1.1f\n",fdEdxMin, fdEdxMax, fEoverPMin, fEoverPMax);
2929  printf("Reject exotic cluster: %s\n",decision[fRejectExoticCluster]);
2930  printf("Remove problematic region in SM4: %s\n",decision[fRemoveBadChannel]);
2931  printf("Use only good SM (1,2,6,7,8,9) for trigger: %s\n",decision[fUseGoodSM]);
2932  printf("Reject electron: %s\n", decision[fElectronRejection]);
2933  printf("Correct hadron: %s\n",decision[fHadronicCorrection]);
2934  printf("HC fraction: %2.1f up to %2.0f GeV/c\n",fFractionHC,fHCLowerPtCutMIP);
2935  printf("Find charged jets: %s\n", decision[fFindChargedOnlyJet]);
2936  printf("Find netural jets: %s\n", decision[fFindNeutralOnlyJet]);
2937  printf("Find good jets: %s\n",decision[fSpotGoodJet]);
2938  printf("Jet radius: %s\n",fRadius.Data());
2939  printf("Jet recombination scheme: %s\n",recombination[fRecombinationScheme]);
2940  printf("Correct tracking efficiency: %s\n",decision[fCheckTrkEffCorr]);
2941  printf("Save jet QA histos: %s\n",decision[fSaveQAHistos]);
2942  printf("Systematics: jet efficiency: %s with variation %1.0f%%\n",decision[fSysJetTrigEff],fVaryJetTrigEff*100);
2943  printf("Systematics: tracking efficiency: %s with variation %1.0f%%\n",decision[fSysTrkEff],fVaryTrkEff*100);
2944  printf("Systematics: track pt resolution: %s with variation %1.0f%%\n",decision[fSysTrkPtRes],fVaryTrkPtRes*100);
2945  printf("Systematics: track-cluster matching: %s with |dEta|<%2.3f, |dPhi|<%2.3f\n",decision[fSysTrkClsMth],fCutdEta,fCutdPhi);
2946  printf("Systematics: EMCal non-linearity: %s\n",decision[fSysNonLinearity]);
2947  printf("Systematics: EMCal energy scale: %s with uncertainty of %1.0f%%\n",decision[fSysClusterEScale],fVaryClusterEScale*100);
2948  printf("Systematics: EMCal energy resolution: %s with uncertainty of %1.0f%%\n",decision[fSysClusterERes],fVaryClusterERes*100);
2949  printf("Smear lhc12a15a: %s\n",decision[fSmearMC]);
2950  printf("Run UE analysis: %s\n",decision[fRunUE]);
2951  printf("Run secondaries: %s\n",decision[fRunSecondaries]);
2952  printf("=======================================\n\n");
2953 }
2954 
2955 //________________________________________________________________________
2957 {
2958  //
2959  // Check if the event containts exotic clusters
2960  //
2961 
2962  Double_t leadingE = 0;
2963  Int_t nCls = fESD->GetNumberOfCaloClusters();
2964  for(Int_t icl=0; icl<nCls; icl++)
2965  {
2966  AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2967  if(!cluster) continue;
2968  if (!cluster->IsEMCAL()) continue;
2969  if(fRecoUtil->IsExoticCluster(cluster, (AliVCaloCells*)fESD->GetEMCALCells()) && cluster->E()>leadingE)
2970  {
2971  leadingE = cluster->E();
2972  }
2973  }
2974  if(leadingE>3) fIsExoticEvent3GeV = kTRUE;
2975  if(leadingE>5) fIsExoticEvent5GeV = kTRUE;
2976 }
2977 
2978 
2979 //________________________________________________________________________
2981 {
2982  //
2983  // Check if the primary vertex exists
2984  //
2985  const AliESDVertex* vtx = fESD->GetPrimaryVertex();
2986  if (!vtx || vtx->GetNContributors()<1) return kFALSE;
2987  else return kTRUE;
2988 }
2989 
2990 
2991 //________________________________________________________________________
2993 {
2994  //
2995  // Check if the event vertex is good
2996  //
2997  return kTRUE;
2998 }
2999 
3000 //________________________________________________________________________
3002 {
3003  //
3004  // Check if the event only has valid TPC vertex
3005  //
3006 
3007  const Bool_t isPrimaryVtx = HasPrimaryVertex();
3008 
3009  Bool_t isGoodVtx = kTRUE;
3010  AliESDVertex* goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexTracks());
3011  if(!goodvtx || goodvtx->GetNContributors()<1)
3012  goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexSPD()); // SPD vertex
3013  if(!goodvtx || !goodvtx->GetStatus()) isGoodVtx = kFALSE; // Event rejected
3014 
3015  if( isPrimaryVtx && !isGoodVtx )
3016  return kTRUE;
3017  else
3018  return kFALSE;
3019 }
3020 
3021 //________________________________________________________________________
3023 {
3024  //
3025  // Check the fraction of accepted events that have only TPC vertex
3026  //
3027  Int_t lTriggerType = -1;
3028  if (trigger & AliVEvent::kMB) lTriggerType = 1;
3029  if (trigger & AliVEvent::kFastOnly) lTriggerType = 0;
3030  if (trigger & AliVEvent::kEMC1) lTriggerType = 2;
3031  if(lTriggerType==-1) return;
3032  fhEventStatTPCVtx->Fill(0.5+lTriggerType*3);
3033  if(HasPrimaryVertex()) fhEventStatTPCVtx->Fill(1.5+lTriggerType*3);
3034  if(IsTPCOnlyVtx()) fhEventStatTPCVtx->Fill(2.5+lTriggerType*3);
3035 }
3036 
3037 //________________________________________________________________________
3039 {
3040  //
3041  // Check if the event is contaminated by LED signal
3042  //
3043 
3044  AliESDCaloCells *cells = fESD->GetEMCALCells();
3045  Short_t nCells = cells->GetNumberOfCells();
3046  Int_t nCellCount[2] = {0,0};
3047  for(Int_t iCell=0; iCell<nCells; iCell++)
3048  {
3049  Int_t cellId = cells->GetCellNumber(iCell);
3050  Double_t cellE = cells->GetCellAmplitude(cellId);
3051  Int_t sMod = fGeom->GetSuperModuleNumber(cellId);
3052 
3053  if(sMod==3 || sMod==4)
3054  {
3055  if(cellE>0.1)
3056  nCellCount[sMod-3]++;
3057  }
3058  }
3059  Bool_t isLED=kFALSE;
3060 
3061  if(fPeriod.CompareTo("lhc11a")==0)
3062  {
3063  if (nCellCount[1] > 100)
3064  isLED = kTRUE;
3065  Int_t runN = fESD->GetRunNumber();
3066  if (runN>=146858 && runN<=146860)
3067  {
3068  if(fTriggerType==0 && nCellCount[0]>=21) isLED=kTRUE;
3069  if(fTriggerType==1 && nCellCount[0]>=35) isLED=kTRUE;
3070  }
3071  }
3072 
3073  return isLED;
3074 }
Int_t pdg
Bool_t fRunUE
Random number generator.
TH2F * fhJetPtWithClsThres[2][kNRadii]
pt of jets containing tracks above certian threshold
virtual Int_t Run()
double Double_t
Definition: External.C:58
Definition: External.C:260
Bool_t IsGoodCluster(AliESDCaloCluster *cluster)
TH2F * fhJetResolutionNK[2][kNRadii]
Jet response due to missing weakly decayed particles using response matrix.
TH2F * fhHCTrkPtAmbig[2][kNRadii][5]
Cleanly subtracted charged pt.
TH1F * fhUEJetPtNorm[3][2][2]
Leading jet pt vs constituent pt in underlying event.
Definition: External.C:236
TH3F * fhJetPtVsZ[3][kNRadii]
Leading charged constituent Z vs jet pt.
TH1F * fhJetPtInExoticEvent[2][kNRadii]
Jet pt vs Fcross vs Zleading.
TH2F * fhJetResponseNKSM[2][kNRadii]
Jet resolution due to missing weakly decayed particles using response matrix.
TH2F * fhHCTrkPtClean[2][kNRadii][5]
f*subtracted cluster energy vs jet pt
AliESDtrackCuts * fEsdTrackCuts
Reco utility.
void AnalyzeSecondaryContributionViaMatching(AliFJWrapper *jetFinder, const Int_t r, const Int_t type, const Int_t etaCut)
TH1F * fEventZ[2]
Generated event vertex z.
TH1F * fhEventStatTPCVtx
Event counts used for jet analysis.
TClonesArray * fMcTruthAntikt[kNRadii]
Jet finder for particle jets.
virtual void CopySettingsFrom(const AliFJWrapper &wrapper)
TH1F * fhSysClusterE[2][2]
Jets in full TPC acceptance in events with TPC only vertex.
const Float_t kRadius[3]
TH3F * fhFcrossVsZleading[2][kNRadii]
Error made by hadronic correction assessed by using particle jet.
AliESDtrack * GetAcceptTrack(AliESDtrack *esdtrack)
TH2F * fhNeutralPtInJet[3][kNRadii]
Subtracted energy due to hadronic correction.
TH2F * fhJetResponseNK[2][kNRadii]
Energy fraction lost due to weakly decaying particles.
TH2F * fhChargedPtInJet[3][kNRadii]
pt of neutral constituents in jet
TH2F * fhTrigNeuPtInJet[3][kNRadii]
pt of neutral constituents in jet
void SetMaxRap(Double_t maxrap)
Definition: AliFJWrapper.h:126
TH2F * fhUEJetPtVsSumPt[3][2][2]
Jet pt in exotic events.
Double_t GetLeadingZ(const Int_t jetIndex, AliFJWrapper *jetFinder)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
const Double_t kPI
AliStack * stack
TF1 * fTriggerEfficiency[10]
Trigger turn-on curves of EMCal clusters.
void FillAODJets(TClonesArray *fJetArray, AliFJWrapper *jetFinder, const Bool_t isTruth=0)
TH2F * fhSysNCellVsClsE[2][2]
Cluster energy distribution before and after hadonic correction.
Bool_t IsElectron(AliESDtrack *track, Double_t clsE) const
AliFJWrapper * fDetJetFinder[3][2][kNRadii]
void FindDetJets(const Int_t s=0, const Int_t a=0, const Int_t r=0)
TH2F * fhJetResponseWP[2][kNRadii]
Jet response due to missing neutron and K0L using response matrix.
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
TH2F * fhSubClsEVsJetPt[2][kNRadii][5]
pT distribution of pions detected
TH1F * fJetCount[3][kNRadii]
Jet NPart fraction.
Double_t GetSmearedTrackPt(AliESDtrack *track)
Int_t GetClusterSuperModule(AliESDCaloCluster *cluster)
THnSparse * fJetEnergyFraction[3][kNRadii]
Contribution of low pt particles to jet energy.
TClonesArray * fJetTCA[3][2][kNRadii]
Jet finder.
Double_t GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
int Int_t
Definition: External.C:63
Bool_t fIsMC
Array of MC particles.
TH1F * fhClsE[2]
Leading jet normalization.
Double_t GetExoticEnergyFraction(AliESDCaloCluster *cluster)
unsigned int UInt_t
Definition: External.C:33
TH1F * fhNTrials[2]
reconstructed event vertex z
float Float_t
Definition: External.C:68
Double_t GetJetArea(UInt_t idx) const
TH1D * fTriggerCurve[10]
Offline trigger mask.
TObjArray * fClusterArray
Array of input tracks.
void CheckTPCOnlyVtx(const UInt_t trigger)
TH2F * fhJetResponseWPSM[2][kNRadii]
Jet response due to missing neutron and K0L via matching.
void RunAnalyzeUE(AliFJWrapper *jetFinder, const Int_t type, const Bool_t isMCTruth)
Double_t SubtractClusterEnergy(AliESDCaloCluster *cluster, Double_t &eRes, Double_t &MIPE, Double_t &MCsubE)
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
Double_t GetJetMissPtDueToTrackingEfficiency(const Int_t jetIndex, AliFJWrapper *jetFinder, const Int_t radiusIndex)
TH2F * fhJetResolutionWP[2][kNRadii]
Jet resolution due to missing neutron and K0L using response matrix.
Definition: External.C:212
TH1F * fhJetEventStat
Event counts to keep track of rejection criteria.
Bool_t fSysJetTrigEff
Response matrix for secondary particles.
AliEMCALGeometry * fGeom
Fit of trigger turn-on curves for EMCal clusters above 4-5 GeV.
TList * fOutputList
TCA of MC truth anti-kt jets.
AliEMCALRecoUtils * fRecoUtil
EMCal goemetry utility.
TH2F * fhJetPtWithTrkThres[2][kNRadii]
Jet pt vs track pt contribution vs track class.
virtual void Clear(const Option_t *="")
AliMCEvent * fMC
AOD object.
AliFJWrapper * fTrueJetFinder[kNRadii]
Parameterazation of momentum resolution estimated from data.
Double_t GetOfflineTriggerProbability(AliESDCaloCluster *cluster)
TH2F * fhChLeadZVsJetPt[2][kNRadii]
pt of leading charged constituents in jet
TH2F * fhUEJetPtVsConsPt[3][2][2]
Leading jet pt vs underlying event pt.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TH3F * fRelTrkCon[2][kNRadii]
Jet pt vs constituent Z vs constituent type.
short Short_t
Definition: External.C:23
TH2F * fhNKFracVsJetPt[2][kNRadii]
NCell vs cluster energy before and after hadonic correction.
TH3F * fhJetResolutionWPSM[2][kNRadii]
Jet resolution due to missing neutron and K0L via matching.
TArrayI * fMcPartArray
Array of input clusters.
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
TH2F * fhSubEVsTrkPt[2][4]
of matched tracks per cluster
TObjArray * fTrackArray
MC stack.
TH1F * fVertexGenZ[2]
Check if the chunk is corrupted.
Int_t GetParentParton(const Int_t ipart)
TH3F * fhJetResolutionNKSM[2][kNRadii]
Jet response due to missing weakly decayed particles via matching.
Bool_t fSysClusterEScale
Non-linearity correction functions for data.
Int_t FindEnergyMatchedJet(AliFJWrapper *jetFinder1, const Int_t index1, AliFJWrapper *jetFinder2, const Double_t fraction=0.5)
TH2F * fhLeadNePtInJet[3][kNRadii]
pt of charged constituents in jet
TH2F * fhJetPtVsLowPtCons[2][kNRadii][2]
pt of jets containing clusters above certian threshold
Bool_t fConstrainChInEMCal
TCA of jets: in - akt - r.
unsigned short UShort_t
Definition: External.C:28
TH1F * fhChunkQA
Event counts for TPC only vertices.
Bool_t IsGoodJet(fastjet::PseudoJet jet, Double_t radius)
Int_t FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR)
const char Option_t
Definition: External.C:48
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
TH1F * fhCorrTrkEffPtBin[2][kNRadii]
Fit function of tracking efficiency.
const Int_t nbins
bool Bool_t
Definition: External.C:53
TH3F * fhJetInTPCOnlyVtx[2][kNRadii]
Cluster energy distribution.
Bool_t IsGoodMcPartilce(const AliVParticle *vParticle, const Int_t ipart)
TH2F * fHCOverSubE[kNRadii][5]
Ambiguously subtracted charged pt.
Double_t GetMeasuredJetPtResolution(const Int_t jetIndex, AliFJWrapper *jetFinder)
TH1F * fhCorrTrkEffSample[2][kNRadii][kNBins]
Number of tracks per jet pt bin, used to correct the tracking efficiency explicitly.
const Double_t kdRCut[3]
TH2F * fHCOverSubEFrac[kNRadii][5]
Error made by hadronic correction assessed by using particle jet.
void AnalyzeSecondaryContribution(AliFJWrapper *jetFinder, const Int_t r, const Int_t etaCut)
TH2F * fhWeakFracVsJetPt[2][kNRadii]
Energy fraction lost due to missing neutron and K0L.
TRandom3 * fRandomGen
Tracking efficiency estimated from simulation.
THnSparse * fJetNPartFraction[3][kNRadii]
Jet energy fraction.
TH2F * fhLeadChPtInJet[3][kNRadii]
pt of leading neutral constituents in jet
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
void AnalyzeJets(AliFJWrapper *jetFinder, const Int_t type, const Int_t r)