1 // **************************************
2 // Full pp jet task - ESD input only
3 // Extract the jet spectrum and all the
4 // systematic uncertainties
6 // **************************************
14 #include <TProfile2D.h>
15 #include <THnSparse.h>
19 #include <TClonesArray.h>
24 #include "AliAODEvent.h"
25 #include "AliAODInputHandler.h"
26 #include "AliAnalysisManager.h"
27 #include "AliAnalysisTask.h"
28 #include "AliCentrality.h"
29 #include "AliAnalysisTaskFullppJet.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"
41 #include "AliGenEventHeader.h"
42 #include "AliGenPythiaEventHeader.h"
43 #include "TGeoGlobalMagField.h"
45 #include "AliAODJet.h"
46 #include "AliFJWrapper.h"
54 ClassImp(AliAnalysisTaskFullppJet)
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};
60 //________________________________________________________________________
61 AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet() :
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),
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)
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++)
101 fVertexGenZ[i] = 0x0;
104 fhNMatchedTrack[i] = 0x0;
106 for(Int_t k=0; k<2; k++)
108 fhSysClusterE[i][k] = 0x0;
109 fhSysNCellVsClsE[i][k] = 0x0;
113 for(Int_t r=0; r<kNRadii; r++)
115 for(Int_t j=0; j<5; j++)
117 fHCOverSubE[r][j] = 0x0;
118 fHCOverSubEFrac[r][j] = 0x0;
119 for(Int_t k=0; k<2; k++)
121 fhSubClsEVsJetPt[k][r][j] = 0x0;
122 fhHCTrkPtClean[k][r][j] = 0x0;
123 fhHCTrkPtAmbig[k][r][j] = 0x0;
126 if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
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;
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++)
148 fhJetPtVsLowPtCons[j][r][k] = 0x0;
150 fhJetPtInExoticEvent[j][r] = 0x0;
151 fhJetInTPCOnlyVtx[j][r] = 0x0;
152 fhCorrTrkEffPtBin[j][r] = 0x0;
153 for(Int_t i=0; i<kNBins; i++)
155 fhCorrTrkEffSample[j][r][i] = 0x0;
161 for(Int_t k=0; k<2; k++)
163 for(Int_t l=0; l<2; l++)
165 fhUEJetPtNorm[j][k][l] = 0x0;
166 fhUEJetPtVsSumPt[j][k][l] = 0x0;
167 fhUEJetPtVsConsPt[j][k][l] = 0x0;
172 for(Int_t i=0; i<2; i++)
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;
187 for(Int_t s=0; s<3; s++)
189 for(Int_t a=0; a<2; a++)
191 for(Int_t r=0; r<kNRadii; r++)
193 fDetJetFinder[s][a][r] = 0x0;
194 fJetTCA[s][a][r] = 0x0;
197 fTrueJetFinder[r] = 0x0;
198 fMcTruthAntikt[r] = 0x0;
204 for(Int_t j=0; j<3; j++)
206 fTrkEffFunc[j] = 0x0;
207 fhSecondaryResponse[j] = 0x0;
209 for(Int_t i=0; i<10; i++)
211 fTriggerCurve[i] = 0x0;
212 fTriggerEfficiency[i] = 0x0;
217 //________________________________________________________________________
218 AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet(const char *name) :
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),
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)
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++)
258 fVertexGenZ[i] = 0x0;
261 fhNMatchedTrack[i] = 0x0;
263 for(Int_t k=0; k<2; k++)
265 fhSysClusterE[i][k] = 0x0;
266 fhSysNCellVsClsE[i][k] = 0x0;
270 for(Int_t r=0; r<kNRadii; r++)
272 for(Int_t j=0; j<5; j++)
274 fHCOverSubE[r][j] = 0x0;
275 fHCOverSubEFrac[r][j] = 0x0;
276 for(Int_t k=0; k<2; k++)
278 fhSubClsEVsJetPt[k][r][j] = 0x0;
279 fhHCTrkPtClean[k][r][j] = 0x0;
280 fhHCTrkPtAmbig[k][r][j] = 0x0;
283 if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
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;
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++)
305 fhJetPtVsLowPtCons[j][r][k] = 0x0;
307 fhJetPtInExoticEvent[j][r] = 0x0;
308 fhJetInTPCOnlyVtx[j][r] = 0x0;
309 fhCorrTrkEffPtBin[j][r] = 0x0;
310 for(Int_t i=0; i<kNBins; i++)
312 fhCorrTrkEffSample[j][r][i] = 0x0;
318 for(Int_t k=0; k<2; k++)
320 for(Int_t l=0; l<2; l++)
322 fhUEJetPtNorm[j][k][l] = 0x0;
323 fhUEJetPtVsSumPt[j][k][l] = 0x0;
324 fhUEJetPtVsConsPt[j][k][l] = 0x0;
329 for(Int_t i=0; i<2; i++)
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;
344 for(Int_t s=0; s<3; s++)
346 for(Int_t a=0; a<2; a++)
348 for(Int_t r=0; r<kNRadii; r++)
350 fDetJetFinder[s][a][r] = 0x0;
351 fJetTCA[s][a][r] = 0x0;
354 fTrueJetFinder[r] = 0x0;
355 fMcTruthAntikt[r] = 0x0;
361 for(Int_t j=0; j<3; j++)
363 fTrkEffFunc[j] = 0x0;
364 fhSecondaryResponse[j] = 0x0;
366 for(Int_t i=0; i<10; i++)
368 fTriggerCurve[i] = 0x0;
369 fTriggerEfficiency[i] = 0x0;
373 //________________________________________________________________________
374 AliAnalysisTaskFullppJet::~AliAnalysisTaskFullppJet()
378 if(fEsdTrackCuts) delete fEsdTrackCuts;
379 if(fHybridTrackCuts1) delete fHybridTrackCuts1;
380 if(fHybridTrackCuts2) delete fHybridTrackCuts2;
382 { fOutputList->Delete(); delete fOutputList;}
383 for(Int_t s=0; s<3; s++)
385 for(Int_t a=0; a<2; a++)
387 for(Int_t r=0; r<kNRadii; r++)
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]; }
393 if(fTrueJetFinder[r]) delete fTrueJetFinder[r];
394 if(fMcTruthAntikt[r]) { fMcTruthAntikt[r]->Delete(); delete fMcTruthAntikt[r]; }
399 if(fRandomGen) delete fRandomGen;
400 for(Int_t i=0; i<3; i++)
402 if(fTrkEffFunc[i]) delete fTrkEffFunc[i];
403 if(fhSecondaryResponse[i]) delete fhSecondaryResponse[i];
405 for(Int_t i=0; i<10; i++)
407 if(fTriggerEfficiency[i]) delete fTriggerEfficiency[i];
408 if(fTriggerCurve[i]) delete fTriggerCurve[i];
410 for(Int_t r=0; r<kNRadii; r++)
412 for(Int_t j=0; j<2; j++)
414 if(fhCorrTrkEffPtBin[j][r]) delete fhCorrTrkEffPtBin[j][r];
415 for(Int_t i=0; i<kNBins; i++)
417 if(fhCorrTrkEffSample[j][r][i]) delete fhCorrTrkEffSample[j][r][i];
421 if(fTrackArray) { fTrackArray->Delete(); delete fTrackArray; }
422 if(fClusterArray) { fClusterArray->Delete(); delete fClusterArray; }
423 if(fMcPartArray) { fMcPartArray->Reset(); delete fMcPartArray; }
425 if(fRecoUtil) delete fRecoUtil;
426 if(fClusterEResolution) delete fClusterEResolution;
427 if(fNonLinear) delete fNonLinear;
428 if(fTrkPtResData) delete fTrkPtResData;
431 //________________________________________________________________________
432 Bool_t AliAnalysisTaskFullppJet::Notify()
435 // Fill the number of tracks per chunk
438 fhChunkQA->SetBinContent(fEDSFileCounter,fNTracksPerChunk);
439 fNTracksPerChunk = 0;
444 //________________________________________________________________________
445 void AliAnalysisTaskFullppJet::UserCreateOutputObjects()
451 if(fRunUE) fFindChargedOnlyJet = kTRUE;
453 const Int_t nTrkPtBins = 100;
454 const Float_t lowTrkPtBin=0, upTrkPtBin=100.;
456 const Int_t nbins = 220;
458 for(Int_t i=0; i<nbins+1; i++)
462 fOutputList = new TList();
463 fOutputList->SetOwner(1);
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");
478 fOutputList->Add(fhJetEventStat);
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");
491 fOutputList->Add(fhJetEventCount);
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");
505 fOutputList->Add(fhEventStatTPCVtx);
508 fhChunkQA = new TH1F("fhChunkQA","# of hybrid tracks per chunk",200,0,200);
509 fOutputList->Add(fhChunkQA);
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};
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};
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"};
534 for(Int_t i=0; i<2; i++)
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]);
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]);
544 for(Int_t i=0; i<3; i++)
546 if(!fIsMC && i==2) continue;
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]);
556 for(Int_t r=0; r<kNRadii; r++)
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;
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
623 for(Int_t k=0; k<5; k++)
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]);
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]);
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]);
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]);
642 for(Int_t j=0; j<4; j++)
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]);
652 for(Int_t k=0; k<2; k++)
654 for(Int_t l=0; l<2; l++)
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]);
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]);
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]);
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]);
676 for(Int_t k=0; k<2; k++)
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]);
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]);
692 for(Int_t r=0; r<kNRadii; r++)
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++)
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]);
706 if(fRunSecondaries && fAnaType==0)
708 for(Int_t i=0; i<2; i++)
710 for(Int_t r=0; r<3; r++)
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
746 printf("\n=======================================\n");
747 printf("===== Jet task configuration ==========\n");
749 if(fNonStdBranch.Length()!=0)
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++)
756 if(!fFindChargedOnlyJet && s==1) continue;
757 if(!fFindNeutralOnlyJet && s==2) continue;
758 for(Int_t a=0; a<2; a++)
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++)
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;
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());
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();
785 if(fIsMC && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase))
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());
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();
808 fRandomGen = new TRandom3(0);
809 if(fCheckTrkEffCorr && fAnaType==0)
811 TFile f ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffFit.root","read");
812 for(Int_t i=0; i<3; i++)
814 fTrkEffFunc[i] = new TF1(*((TF1*)f.Get(Form("Trk_eff_fit_%d",i+1))));
818 if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0)
820 TFile f1 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffSampling.root","read");
821 for(Int_t j=0; j<2; j++)
824 if(fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0) tmp = 2;
825 for(Int_t r=0; r<kNRadii; r++)
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++)
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]))));
838 if(fRunSecondaries && fAnaType==0)
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++)
844 fhSecondaryResponse[j] = new TH2F(*((TH2F*)f2.Get(Form("DetectorReponse_%s",secondaryName[j]))));
848 if(fCheckTriggerMask)
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")));
856 for(Int_t i=0; i<10; i++)
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))));
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);
868 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
871 AliError("No EMCal geometry is available!");
874 fRecoUtil = new AliEMCALRecoUtils();
875 if(fRejectExoticCluster || fRejectExoticTrigger)
876 fRecoUtil->SwitchOnRejectExoticCluster();
879 fRecoUtil->SwitchOffRejectExoticCluster();
880 fRecoUtil->SwitchOffRejectExoticCell();
882 if(fRemoveBadChannel)
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);
891 fRecoUtil->SwitchOffBadChannelsRemoval();
893 fRecoUtil->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrected);
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);
901 fRecoUtil->SwitchOnCutEtaPhiSeparate();
902 fRecoUtil->SetCutEta(fCutdEta);
903 fRecoUtil->SetCutPhi(fCutdPhi);
906 fTrackArray = new TObjArray();
907 fTrackArray->SetOwner(1);
908 fClusterArray = new TObjArray();
909 fClusterArray->SetOwner(1);
910 fMcPartArray = new TArrayI();
913 //error calculation in THnSparse
914 Int_t nObj = fOutputList->GetEntries();
915 for(Int_t i=0; i<nObj; i++)
917 TObject *obj = (TObject*) fOutputList->At(i);
918 if (obj->IsA()->InheritsFrom( "THnSparse" ))
920 THnSparseF *hn = (THnSparseF*)obj;
928 PostData(1, fOutputList);
933 //________________________________________________________________________
934 void AliAnalysisTaskFullppJet::BookHistos()
938 if(fVerbosity>10) printf("[i] Booking histograms \n");
942 //________________________________________________________________________
943 void AliAnalysisTaskFullppJet::UserExec(Option_t *)
946 // Main loop, called for each event.
949 // Get event pointers, check for signs of life
950 Double_t vtxTrueZ = -100;
956 printf("ERROR: Could not retrieve MC event");
959 fStack = fMC->Stack();
960 TParticle *particle = fStack->Particle(0);
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());
968 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
971 AliError("fESD is not available");
975 fhJetEventStat->Fill(0.5);
976 fhJetEventCount->Fill(0.5);
980 if(fVertexGenZ[0]) fVertexGenZ[0]->Fill(vtxTrueZ);
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);
989 if (trigger & AliVEvent::kFastOnly) return; // Reject fast trigger cluster
990 if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0)
992 if (trigger & AliVEvent::kMB) fTriggerType = 0;
993 else if(trigger & AliVEvent::kEMC1) fTriggerType = 1;
994 else fTriggerType = -1;
996 else if (fPeriod.CompareTo("lhc11c",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11d",TString::kIgnoreCase)==0)
998 if (trigger & AliVEvent::kINT7) fTriggerType = 0;
999 else if(trigger & AliVEvent::kEMC7) fTriggerType = 1;
1000 else fTriggerType = -1;
1009 if(!fPhySelForMC) fTriggerType = 0;
1010 else if (trigger & AliVEvent::kAnyINT) fTriggerType = 0;
1011 else fTriggerType = -1;
1015 RunOfflineTrigger();
1016 if(fIsEventTriggerBit) fTriggerType = 1;
1020 if(fTriggerType==-1)
1022 if(fVerbosity>10) printf("Error: worng trigger type %s\n",(fESD->GetFiredTriggerClasses()).Data());
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);
1033 // Reject LED events
1036 fhJetEventStat->Fill(8.5);
1040 fhJetEventStat->Fill(1.5+fTriggerType*2);
1042 // Check if primary vertex exists
1043 if (!HasPrimaryVertex()) return;
1044 fhJetEventStat->Fill(5.5+fTriggerType);
1045 fhJetEventCount->Fill(3.5+fTriggerType*4);
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);
1052 // Check if |Z_vtx|<10cm
1053 if( TMath::Abs(zVertex) > fZVtxMax ) return;
1054 fhJetEventCount->Fill(4.5+fTriggerType*4);
1056 // Check if only TPC vertex exists primitive
1057 fIsTPCOnlyVtx = IsTPCOnlyVtx();
1059 if(!fIsMC && fTriggerType==1)
1061 // Check if event has valid trigger bit
1062 CheckEventTriggerBit();
1063 if(fIsEventTriggerBit)
1065 fhJetEventStat->Fill(7.5);
1066 fhJetEventCount->Fill(9.5);
1070 fhJetEventStat->Fill(2.5+fTriggerType*2);
1071 if(fIsTPCOnlyVtx) fhJetEventStat->Fill(9.5+fTriggerType);
1073 // Check if event contains exotic clusters
1076 // Write jet tree into AOD output in local analysis mode
1077 if(fAnaType==0) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
1080 for(Int_t s=0; s<3; s++)
1082 for(Int_t a=0; a<2; a++)
1084 for(Int_t r=0; r<kNRadii; r++)
1086 if(fJetTCA[s][a][r]) fJetTCA[s][a][r]->Delete();
1087 if(fDetJetFinder[s][a][r]) fDetJetFinder[s][a][r]->Clear();
1091 if(fMcTruthAntikt[r]) fMcTruthAntikt[r]->Delete();
1092 if(fTrueJetFinder[r]) fTrueJetFinder[r]->Clear();
1098 if(fVerbosity>5) printf("# of jets after clear: %d\n",fJetTCA[0][0][0]->GetEntries());
1100 if(fTrackArray) fTrackArray->Delete();
1101 if(fClusterArray) fClusterArray->Delete();
1102 fMcPartArray->Reset();
1105 // get the tracks and fill the input vector for the jet finders
1108 // get EMCal clusters and fill the input vector for the jet finders
1109 GetESDEMCalClusters();
1111 for(Int_t s=0; s<3; s++)
1113 for(Int_t a=0; a<2; a++)
1115 for(Int_t r=0; r<kNRadii; r++)
1118 if(fDetJetFinder[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))
1124 if(fSaveQAHistos) AnalyzeJets(fDetJetFinder[s][a][r],fTriggerType, r); // run analysis on jets
1130 if(fRunUE && fDetJetFinder[1][0][0]) RunAnalyzeUE(fDetJetFinder[1][0][0],fTriggerType,0); // Run analysis on underlying event
1134 for(Int_t r=0; r<kNRadii; r++)
1135 if(fTrueJetFinder[r]) ProcessMC(r); // find particle level jets
1138 // Fast Jet calls END --------------------------------------------------------
1140 PostData(1, fOutputList);
1145 //________________________________________________________________________
1146 void AliAnalysisTaskFullppJet::FindDetJets(const Int_t s, const Int_t a, const Int_t r)
1149 // Jet finding is happening here
1152 Bool_t isCh = kTRUE;
1153 Bool_t isNe = kTRUE;
1154 if(s==1) isNe = kFALSE;
1155 if(s==2) isCh = kFALSE;
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++)
1164 AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
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());
1180 if(fVerbosity>5) printf("[i] # of tracks filled: %d\n",countTracks);
1184 // Feed in EMCal clusters
1185 Double_t vertex[3] = {0, 0, 0};
1186 fESD->GetVertex()->GetXYZ(vertex) ;
1187 TLorentzVector gamma;
1189 Int_t countClusters = 0;
1190 Int_t nclusters = fClusterArray->GetEntriesFast();
1191 for (Int_t ic = 0; ic < nclusters; ic++)
1193 AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic);
1195 cl->GetMomentum(gamma, vertex);
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));
1201 if(fVerbosity>5) printf("[i] # of clusters filled: %d\n",countClusters);
1204 fDetJetFinder[s][a][r]->Run();
1208 //________________________________________________________________________
1209 void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper *jetFinder, const Bool_t isTruth)
1212 // Fill the found jets into AOD output
1213 // Only consider jets pointing to EMCal and with pt above 1GeV/c
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());
1224 for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
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;
1229 AliAODJet tmpJet (jetsIncl[ij].px(), jetsIncl[ij].py(), jetsIncl[ij].pz(), jetsIncl[ij].E());
1230 jet = new ((*fJetArray)[jetCount]) AliAODJet(tmpJet);
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());
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++)
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)
1246 totNePt += constituents[ic].perp();
1247 if(constituents[ic].perp()>leadNePt)
1248 leadNePt = constituents[ic].perp();
1250 neE += constituents[ic].E();
1251 if(constituents[ic].perp()>fClsEtMax[fTriggerType])
1252 isHighPtTrigger = kTRUE;
1255 AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1256 if(cluster->Chi2()>0.5) isTriggering = kTRUE;
1261 totChPt += constituents[ic].perp();
1263 if(constituents[ic].perp()>leadChPt)
1265 leadChPt = constituents[ic].perp();
1268 AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1269 leadTrkType = track->GetTRDQuality();
1270 if(track->GetIntegratedLength()>500) isHyperTrack = kTRUE;
1273 if(constituents[ic].perp()>fTrkPtMax[fTriggerType])
1274 isHighPtTrigger = kTRUE;
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)
1283 leadPt = constituents[ic].perp();
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);
1293 if(fIsTPCOnlyVtx) jet->SetTrigger(AliAnalysisTaskFullppJet::kTPCOnlyVtx);
1294 if(constituents[leadIndex].user_index()>0) jet->SetTrigger(AliAnalysisTaskFullppJet::kLeadCh);
1296 if(jetsIncl[ij].E()>0) jet->SetNEF(neE/jetsIncl[ij].E());
1297 jet->SetPtLeading(leadPt);
1298 jet->SetPtSubtracted(leadTrkType,0);
1300 Double_t effAErrCh = 0, effAErrNe = leadTrkType;
1301 Double_t chBgEnergy = 10, neBgEnergy = nNePart;
1304 effAErrCh = GetMeasuredJetPtResolution(ij,jetFinder);
1305 if(fIsMC && constituents[leadIndex].user_index()>0)
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())
1312 AliVParticle* vParticle = fMC->GetTrack(ipart);
1315 Double_t truePt = vParticle->Pt();
1316 chBgEnergy = (truePt-pt)/truePt;
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;
1328 printf("# of ref tracks: %d = %d, and nef=%f\n",jet->GetRefTracks()->GetEntries(), (Int_t)constituents.size(), jet->GetNEF());
1330 // For catch good high-E jets
1331 if(fSpotGoodJet && !isTruth)
1333 if(jetsIncl[ij].perp()>100)
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++)
1341 if(constituents[ic].user_index()<0)
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());
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());
1353 printf("==============================\n\n");
1356 // End of catching good high-E jets
1358 if(fVerbosity>5) printf("%s has %d jets\n",fJetArray->GetName(), fJetArray->GetEntries());
1361 //________________________________________________________________________
1362 Int_t AliAnalysisTaskFullppJet::GetParentParton(const Int_t ipart)
1365 // Find the parent parton for a given MC particle by tracing back the parent
1366 // DUMMY; NEED improvement
1372 //________________________________________________________________________
1373 Bool_t AliAnalysisTaskFullppJet::IsGoodJet(fastjet::PseudoJet jet, Double_t rad)
1376 // Check if the jet pt and direction fulfill the requirement
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;
1386 //________________________________________________________________________
1387 void AliAnalysisTaskFullppJet::ProcessMC(const Int_t r)
1390 // Main function for jet finding in MC
1393 Int_t npart = fMC->GetNumberOfTracks();
1394 fMcPartArray->Set(npart);
1395 Int_t countPart = 0;
1396 for(Int_t ipart=0; ipart<npart; ipart++)
1398 AliVParticle* vParticle = fMC->GetTrack(ipart);
1399 if(!IsGoodMcPartilce(vParticle,ipart)) continue;
1401 Int_t pdgCode = vParticle->PdgCode();
1402 if( fRejectNK && (pdgCode==130 || pdgCode==2112) ) continue;
1404 if( fRejectWD && (pdgCode==310 || pdgCode==3112 || pdgCode==3122 || pdgCode==3222 || pdgCode==3312 || pdgCode==3322 || pdgCode==3334)) continue;
1406 if( fChargedMC && vParticle->Charge()==0 ) continue;
1409 if(vParticle->Charge()==0) { index=-1; }
1410 fMcPartArray->AddAt(ipart, countPart);
1413 fTrueJetFinder[r]->AddInputVector(vParticle->Px(), vParticle->Py(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
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());
1417 fMcPartArray->Set(countPart);
1418 fTrueJetFinder[r]->Run();
1420 if(fMcTruthAntikt[r]) FillAODJets(fMcTruthAntikt[r], fTrueJetFinder[r], 1);
1421 if(fSaveQAHistos) AnalyzeJets(fTrueJetFinder[r], 2, r);
1422 if(fRunUE && r==0) RunAnalyzeUE(fTrueJetFinder[r], 2, 1);
1424 // Run analysis on secondary particles
1425 if(fRunSecondaries && fAnaType==0)
1427 for(Int_t i=0; i<2; i++)
1429 AnalyzeSecondaryContribution(fTrueJetFinder[r],r,i);
1430 AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,0,i);
1431 AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,1,i);
1437 //________________________________________________________________________
1438 void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t type, const Int_t r)
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
1445 Double_t vertex[3] = {0, 0, 0};
1446 fESD->GetVertex()->GetXYZ(vertex) ;
1447 TLorentzVector gamma;
1449 const Int_t nBins = fJetEnergyFraction[type][r]->GetAxis(1)->GetNbins();
1450 Float_t radiusCut[nBins];
1451 Float_t eFraction[nBins];
1453 for(Int_t i=0; i<nBins; i++)
1455 radiusCut[i] = (fJetEnergyFraction[type][r]->GetAxis(1)->GetBinWidth(1)) * (i+1);
1459 std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1461 for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1463 if(type<2 && fCheckTPCOnlyVtx && fIsTPCOnlyVtx)
1465 fhJetInTPCOnlyVtx[type][r]->Fill(jetsIncl[ij].perp(),jetsIncl[ij].phi()*TMath::RadToDeg(),jetsIncl[ij].eta());
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();
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++)
1478 if(constituents[ic].perp()>maxPt)
1480 maxPt = constituents[ic].perp();
1481 leadingIndex = constituents[ic].user_index();
1484 if(constituents[ic].user_index()<0)
1486 neE += constituents[ic].E();
1487 constituentType = 3;
1491 if(type==2) constituentType = 0;
1494 AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1495 constituentType = (Int_t)track->GetTRDQuality();
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;
1503 if(type<2 && leadingIndex<0)
1505 AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(leadingIndex*(-1)-1);
1506 fhFcrossVsZleading[type][r]->Fill(jetPt,GetExoticEnergyFraction(cluster),leadingZ);
1509 if(leadingZ > 0.98) continue; // Z cut
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);
1515 Double_t totTrkPt[3] = {0.,0.,0.};
1517 for(Int_t i=0; i<nBins; i++) { eFraction[i] = 0.; nPart[i] = 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++)
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();
1532 if(constituents[ic].user_index()<0)
1534 fhNeutralPtInJet[type][r]->Fill(jetPt, partPt);
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);
1542 if(fStudySubEInHC || fStudyMcOverSubE)
1544 cluster->GetMomentum(gamma, vertex);
1545 Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
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++)
1553 subClsE[j] += ((fraction[j]*subEtmp>clsE)?clsE:fraction[j]*subEtmp)*sinTheta;
1560 if(partPt>leadChPt) {leadChPt = partPt;leadChIndex=ic;}
1561 fhChargedPtInJet[type][r]->Fill(jetPt, partPt);
1562 chPt += constituents[ic].perp();
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++)
1573 AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(j);
1574 if( clusterIndex == cluster->GetID() )
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++)
1590 if(mipETmp>cluster->E()) subTrkPtClean[k] += partPt;
1591 else subTrkPtAmbig[k] += partPt;
1595 if(fraction[k]*subEtmp>cluster->E()) subTrkPtClean[k] += partPt;
1596 else subTrkPtAmbig[k] += partPt;
1603 Float_t dR = TMath::Sqrt( (partEta-jetEta)*(partEta-jetEta) + (partPhi-jetPhi)*(partPhi-jetPhi) );
1604 for(Int_t i=0; i<nBins; i++)
1606 if( dR < radiusCut[i] )
1608 eFraction[i] += partE;
1614 fhLeadNePtInJet[type][r]->Fill(jetPt, leadNePt);
1615 fhLeadChPtInJet[type][r]->Fill(jetPt, leadChPt);
1616 if(type<2 && leadChIndex>-1)
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);
1622 if((fStudySubEInHC||fStudyMcOverSubE) && type<2)
1624 for(Int_t i=0; i<5; i++)
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]));
1632 if(type==0 && fIsMC && fStudyMcOverSubE)
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]));
1639 if(type<2 && chPt>0)
1641 for(Int_t i=0; i<3; i++)
1643 fRelTrkCon[type][r]->Fill(jetPt,totTrkPt[i]/chPt,i);
1648 for(Int_t ibin=0; ibin<nBins; ibin++)
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);
1656 // Get the jet pt containing tracks or clusters above some threshold
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++)
1665 Float_t partPt = constituents[ic].perp();
1667 if(partPt<0.3) lowPt[0] += partPt;
1668 else if(partPt<0.5) lowPt[1] += partPt;
1670 if(constituents[ic].user_index()>0)
1672 for(Int_t it=0; it<3; it++)
1674 if(partPt>thres[it]) okCh[it] = 1;
1679 for(Int_t icl=0; icl<3; icl++)
1681 if(partPt>thres[icl]) okNe[icl] = 1;
1685 for(Int_t i=0; i<3; i++)
1688 fhJetPtWithTrkThres[type][r]->Fill(i,jetPt);
1690 fhJetPtWithClsThres[type][r]->Fill(i,jetPt);
1692 for(Int_t k=0; k<2; k++) fhJetPtVsLowPtCons[type][r][k]->Fill(jetPt,lowPt[k]/jetPt);
1697 //________________________________________________________________________
1698 void AliAnalysisTaskFullppJet::AnalyzeSecondaryContribution(AliFJWrapper *jetFinder, const Int_t r, const Int_t etaCut)
1701 // Analyze secondaries
1704 std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1705 for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
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;
1713 for(UInt_t ic=0; ic<constituents.size(); ic++)
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();
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;
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();
1736 // Missing neutron and K0L
1737 if(pdg==130 || pdg==2112) dNKPt -= mcParticle->Pt();
1740 fhNKFracVsJetPt[etaCut][r]->Fill(jetPt,nkPt/jetPt);
1741 fhWeakFracVsJetPt[etaCut][r]->Fill(jetPt,weakPt/jetPt);
1743 Double_t jetNewWPPt = jetPt + dWPPt;
1744 fhJetResponseWP[etaCut][r]->Fill(jetPt,jetNewWPPt);
1745 fhJetResolutionWP[etaCut][r]->Fill(jetPt,(jetPt-jetNewWPPt)/jetPt);
1747 Double_t jetNewNKPt = jetPt + dNKPt;
1748 fhJetResponseNK[etaCut][r]->Fill(jetPt,jetNewNKPt);
1749 fhJetResolutionNK[etaCut][r]->Fill(jetPt,(jetPt-jetNewNKPt)/jetPt);
1754 //________________________________________________________________________
1755 void AliAnalysisTaskFullppJet::AnalyzeSecondaryContributionViaMatching(AliFJWrapper *jetFinder, const Int_t r, const Int_t type, const Int_t etaCut)
1758 // Estimate the contribution of missing energies via matching
1763 if(type!=0 && type!=1) return;
1765 AliFJWrapper fj("fj","fj");
1766 fj.CopySettingsFrom(*jetFinder);
1768 Int_t npart = fMC->GetNumberOfTracks();
1770 for(Int_t ipart=0; ipart<npart; ipart++)
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
1779 if(pdgCode==310) pType=0;
1780 else if (pdgCode==3122) pType=1;
1781 else if (pdgCode<=3400 && pdgCode>=3100) pType=2;
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();
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);
1798 std::vector<fastjet::PseudoJet> jets_incl_mc = fj.GetInclusiveJets();
1799 std::vector<fastjet::PseudoJet> jets_incl_mc_std = jetFinder->GetInclusiveJets();
1801 for(UInt_t ij=0; ij<jets_incl_mc.size(); ij++)
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);
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);
1817 if(type==0) fhJetResponseNKSM[etaCut][r]->Fill(jetPtStd,jetPt);
1818 if(type==1) fhJetResponseWPSM[etaCut][r]->Fill(jetPtStd,jetPt);
1824 //________________________________________________________________________
1825 void AliAnalysisTaskFullppJet::RunAnalyzeUE(AliFJWrapper *jetFinder, const Int_t type, const Bool_t isMCTruth)
1828 // Run analysis to estimate the underlying event
1829 // Adapted from Oliver
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++)
1836 Double_t jetEta = jetsIncl[ij].eta();
1837 Double_t jetPt = jetsIncl[ij].perp();
1838 Double_t jetPhi = jetsIncl[ij].phi();
1847 if(leadpt<1 || TMath::Abs(leadEta)>0.5 ) return;
1849 Int_t backIndex = -1;
1850 Bool_t isBackToBack = kFALSE;
1851 for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
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)
1858 isBackToBack = kTRUE;
1862 for(Int_t ij=0; ij<(Int_t)jetsIncl.size(); ij++)
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;
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;
1874 fhUEJetPtNorm[type][0][0]->Fill(leadpt);
1875 if(isBackToBack) fhUEJetPtNorm[type][0][1]->Fill(leadpt);
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();
1885 for(Int_t j=0; j<2; j++)
1887 Double_t ueCh = 0, ueChNe = 0.;
1888 for(Int_t it=0; it<ntracks; it++)
1890 AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
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)
1897 fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,t->Pt());
1898 if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,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)
1902 fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,t->Pt());
1903 if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,t->Pt());
1908 fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh);
1909 if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh);
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++)
1916 AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic);
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)
1926 ueChNe += gamma.Pt();
1927 fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,gamma.Pt());
1928 if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,gamma.Pt());
1931 fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe);
1932 if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe);
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++)
1942 Double_t ueCh = 0, ueChNe = 0., ueChNe2 = 0.;
1943 for(Int_t ipos=0; ipos<npart; ipos++)
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)
1952 Double_t pt = vParticle->Pt();
1954 fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,pt);
1955 if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,pt);
1956 if( vParticle->Charge()!=0 )
1958 fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,pt);
1959 if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,pt);
1965 if(pt>0.5) ueChNe2 += pt;
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);
1977 //________________________________________________________________________
1978 Bool_t AliAnalysisTaskFullppJet::IsGoodJet(AliAODJet *jet, Double_t rad)
1981 // Check if it is a good jet
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;
1991 //________________________________________________________________________
1992 Double_t AliAnalysisTaskFullppJet::GetLeadingZ(const Int_t jetIndex, AliFJWrapper *jetFinder)
1995 // Get the leading z of the input jet
1999 std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
2000 std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
2003 for(UInt_t ic=0; ic<constituents.size(); ic++)
2005 if(constituents[ic].perp()>maxPt)
2007 maxPt = constituents[ic].perp();
2012 z = GetZ(constituents[index].px(),constituents[index].py(),constituents[index].pz(),jetsIncl[jetIndex].px(),jetsIncl[jetIndex].py(),jetsIncl[jetIndex].pz());
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
2020 // Get the z of a constituent inside of a jet
2023 return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
2026 //________________________________________________________________________
2027 Double_t AliAnalysisTaskFullppJet::GetMeasuredJetPtResolution(const Int_t jetIndex, AliFJWrapper *jetFinder)
2030 // Get jet energy resoultion due to intrinsic detector effects
2033 Double_t jetSigma2 = 0;
2034 std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
2035 for(UInt_t ic=0; ic<constituents.size(); ic++)
2037 if(constituents[ic].user_index()>0)
2039 AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
2040 jetSigma2 += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2044 AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
2045 jetSigma2 += TMath::Power(cluster->GetTOF(),2);
2048 return TMath::Sqrt(jetSigma2);
2053 //________________________________________________________________________
2054 Double_t AliAnalysisTaskFullppJet::GetJetMissPtDueToTrackingEfficiency(const Int_t jetIndex, AliFJWrapper *jetFinder, const Int_t radiusIndex)
2057 // Correct the tracking inefficiency explicitly
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])
2065 printf("Warning: can't get the mean # of tracks per jet with pt=%f in: trigger=%d, radiusIndex=%d\n",jetPt,fTriggerType,radiusIndex);
2068 Int_t ibin = fhCorrTrkEffPtBin[fTriggerType][radiusIndex]->FindFixBin(jetPt);
2069 if(!fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1] || fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1]->Integral()<0.001)
2071 printf("Warning: no sampling distrubtion for jet with pt=%f\n",jetPt);
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++)
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;
2091 //________________________________________________________________________
2092 Bool_t AliAnalysisTaskFullppJet::IsGoodMcPartilce(const AliVParticle* vParticle, const Int_t ipart)
2095 // Select good primary particles to feed into jet finder
2098 if(!vParticle) return kFALSE;
2099 if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;
2100 if (TMath::Abs(vParticle->Eta())>2) return kFALSE;
2104 //________________________________________________________________________
2105 Int_t AliAnalysisTaskFullppJet::FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR)
2108 // Find spatially matched detector and particle jets
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++)
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) );
2123 dEta = jet.eta()-jets_incl[ij].eta();
2124 dPhi = jet.phi()-jets_incl[ij].phi();
2130 //________________________________________________________________________
2131 Int_t AliAnalysisTaskFullppJet::FindEnergyMatchedJet(AliFJWrapper *jetFinder1, const Int_t index1, AliFJWrapper *jetFinder2, const Double_t fraction)
2134 // Find matched detector and particle jets based on shared constituents
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;
2144 for(UInt_t ij2=0; ij2<jetsIncl2.size(); ij2++)
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++)
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++)
2156 Double_t consPt1 = constituents1[ic1].perp();
2157 if(constituents1[ic1].user_index()>0)
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;}
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;}
2169 cout<<sharedPt1/jetPt1<<" "<<sharedPt2/jetPt2<<endl;
2170 if(sharedPt1/jetPt1 > fraction && sharedPt2/jetPt2 > fraction)
2176 return matchedIndex;
2179 //________________________________________________________________________
2180 void AliAnalysisTaskFullppJet::GetESDTrax()
2186 Int_t nTrax = fESD->GetNumberOfTracks();
2187 if(fVerbosity>5) printf("[i] # of tracks in event: %d\n",nTrax);
2189 for (Int_t i = 0; i < nTrax; ++i)
2191 AliESDtrack* esdtrack = fESD->GetTrack(i);
2194 AliError(Form("Couldn't get ESD track %d\n", i));
2198 AliESDtrack *newtrack = GetAcceptTrack(esdtrack);
2199 if(!newtrack) continue;
2201 Double_t pt = newtrack->Pt();
2202 Double_t eta = newtrack->Eta();
2203 if (pt < fTrkPtMin[fTriggerType] || pt > fTrkPtMax[fTriggerType] || TMath::Abs(eta) > fTrkEtaMax)
2211 Double_t rand = fRandomGen->Uniform();
2212 if(rand<fVaryTrkEff)
2218 newtrack->SetIntegratedLength(esdtrack->Pt());
2219 fTrackArray->Add(newtrack);
2221 if(fVerbosity>5) printf("[i] # of tracks in event: %d\n", fTrackArray->GetEntries());
2222 fNTracksPerChunk += fTrackArray->GetEntries();
2226 //________________________________________________________________________
2228 AliESDtrack *AliAnalysisTaskFullppJet::GetAcceptTrack(AliESDtrack *esdtrack)
2231 // Get the hybrid tracks
2234 AliESDtrack *newTrack = 0x0;
2236 if(fTrackCutsType==0 || fTrackCutsType==3)
2238 if(fEsdTrackCuts->AcceptTrack(esdtrack))
2240 newTrack = new AliESDtrack(*esdtrack);
2241 newTrack->SetTRDQuality(0);
2243 else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
2245 if(esdtrack->GetConstrainedParam())
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);
2255 else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
2257 if(esdtrack->GetConstrainedParam())
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);
2272 else if(fTrackCutsType==1)
2274 if(fEsdTrackCuts->AcceptTrack(esdtrack))
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
2285 newTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
2286 newTrack->SetTRDQuality(1);
2293 else if (fTrackCutsType==2)
2295 if(fEsdTrackCuts->AcceptTrack(esdtrack))
2297 newTrack = new AliESDtrack(*esdtrack);
2298 newTrack->SetTRDQuality(0);
2305 printf("Unknown track cuts type: %d\n",fTrackCutsType);
2312 //________________________________________________________________________
2313 void AliAnalysisTaskFullppJet::GetESDEMCalClusters()
2316 // Get emcal clusters - selected
2321 if (!TGeoGlobalMagField::Instance()->GetField()) fESD->InitMagneticField();
2322 fRecoUtil->FindMatches(fESD,0x0,fGeom);
2323 fRecoUtil->SetClusterMatchedToTrack(fESD);
2324 fRecoUtil->SetTracksMatchedToCluster(fESD);
2327 const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
2328 for(Int_t i = 0 ; i < nCaloClusters; i++)
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;}
2337 if(fPeriod.Contains("lhc11a",TString::kIgnoreCase))
2339 Double_t newE = newCluster->E() * 1.02;
2340 newCluster->SetE(newE);
2343 // Trigger efficiency systematic uncertainty
2346 Double_t newE = newCluster->E() * (1+fVaryJetTrigEff);
2347 newCluster->SetE(newE);
2348 if(fTriggerType==0) fhClsE[fTriggerType]->Fill(newCluster->E());
2351 if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase)) fhClsE[0]->Fill(newCluster->E());
2352 if(newCluster->Chi2()>0.5) fhClsE[fTriggerType]->Fill(newCluster->E());
2356 // Cluster energy scale systematic uncertainty
2357 if(fSysClusterEScale)
2359 Double_t newE = newCluster->E() * (1+fVaryClusterEScale);
2360 newCluster->SetE(newE);
2363 // Cluster energy resolution systematic uncertainty
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);
2373 // non-linearity systematic uncertainty
2374 if(fSysNonLinearity)
2376 Double_t oldE = newCluster->E();
2377 Double_t newE = oldE/fNonLinear->Eval(oldE) * 1.012;
2378 newCluster->SetE(newE);
2381 // clusterizer systematic uncertainty
2384 fhSysClusterE[fTriggerType][0]->Fill(newCluster->E());
2385 fhSysNCellVsClsE[fTriggerType][0]->Fill(newCluster->E(),newCluster->GetNCells());
2388 // hadronic correction
2389 Double_t clsE = newCluster->E();
2390 Double_t subE = 0., eRes = 0., mcSubE = 0, MIPE = 0.;
2391 if(fElectronRejection || fHadronicCorrection)
2392 subE= SubtractClusterEnergy(newCluster, eRes, MIPE, mcSubE);
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);
2403 // clusterizer systematic uncertainty
2406 fhSysClusterE[fTriggerType][1]->Fill(newCluster->E());
2407 fhSysNCellVsClsE[fTriggerType][1]->Fill(newCluster->E(),newCluster->GetNCells());
2411 if(fVerbosity>5) printf("[i] # of EMCal clusters in event: %d\n", fClusterArray->GetEntries());
2416 //________________________________________________________________________
2417 Bool_t AliAnalysisTaskFullppJet::IsGoodCluster(AliESDCaloCluster *cluster)
2420 // Select good clusters
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;
2431 //________________________________________________________________________
2432 Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *cluster, Double_t &eRes, Double_t &MIPE, Double_t &mcSubE)
2435 // Hadronic correction
2442 eRes += TMath::Power(fClusterEResolution->Eval(cluster->E()),2);
2443 Double_t subE = 0., sumTrkPt = 0., sumTrkP = 0.;
2445 TArrayI *matched = cluster->GetTracksMatched();
2446 if(!matched) return 0;
2447 for(Int_t im=0; im<matched->GetSize(); im++)
2449 Int_t trkIndex = matched->At(im);
2450 if(trkIndex<0 || trkIndex>=fESD->GetNumberOfTracks()) continue;
2451 Bool_t isSelected = kFALSE;
2453 for(Int_t j=0; j<fTrackArray->GetEntriesFast(); j++)
2455 AliESDtrack *tr = (AliESDtrack*)fTrackArray->At(j);
2456 if( trkIndex == tr->GetID() )
2463 if(!isSelected) continue;
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
2471 if(fElectronRejection)
2475 eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
2480 if(fHadronicCorrection)
2482 MIPE += (trkP>0.27)?0.27:trkP; //MIP correction
2492 eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
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);
2505 if(fSaveQAHistos) fhNMatchedTrack[fTriggerType]->Fill(nTrack);
2508 Double_t fraction[4] = {0.3,0.5,0.7,1.0};
2509 for(Int_t j=0; j<4; j++)
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);
2517 eRes = TMath::Sqrt(eRes);
2519 if(fIsMC && nTrack>0)
2521 Double_t neutralE = 0;
2522 TArrayI* labels = cluster->GetLabelsArray();
2525 for(Int_t il=0; il<labels->GetSize(); il++)
2527 Int_t ipart = labels->At(il);
2528 if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
2530 AliVParticle* vParticle = fMC->GetTrack(ipart);
2531 if(vParticle->Charge()==0)
2533 neutralE += vParticle->E();
2538 mcSubE = cluster->E() - neutralE;
2549 //________________________________________________________________________
2551 Double_t AliAnalysisTaskFullppJet::GetExoticEnergyFraction(AliESDCaloCluster *cluster)
2554 // Exotic fraction: f_cross
2555 // Adpated from AliEMCalRecoUtils
2557 if(!cluster) return -1;
2558 AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
2559 if(!cells) return -1;
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);
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);
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;
2578 accept = fRecoUtil->AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells);
2580 if(!accept) return -1;
2582 if(ecell < 0.5) return -1;
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);
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 ;
2595 Float_t eCross = ecell1+ecell2+ecell3+ecell4;
2597 return 1-eCross/ecell;
2601 //________________________________________________________________________
2602 void AliAnalysisTaskFullppJet::GetMCInfo()
2605 // Get # of trials per ESD event
2608 AliStack *stack = fMC->Stack();
2611 AliGenEventHeader *head = dynamic_cast<AliGenEventHeader*>(fMC->GenEventHeader());
2614 AliError("Could not get the event header");
2618 AliGenPythiaEventHeader *headPy = dynamic_cast<AliGenPythiaEventHeader*>(head);
2621 if (headPy->Trials() > 0)
2623 fhNTrials[0]->Fill(0.5,headPy->Trials());
2632 //________________________________________________________________________
2633 void AliAnalysisTaskFullppJet::Terminate(Option_t *)
2636 // Called once at the end of the query
2638 Info("Terminate","Terminate");
2639 AliAnalysisTaskSE::Terminate();
2644 //________________________________________________________________________
2646 Int_t AliAnalysisTaskFullppJet::RunOfflineTrigger()
2649 // Run trigger offline
2652 fIsEventTriggerBit = 0;
2653 Int_t isTrigger = 0;
2654 Int_t ncl = fESD->GetNumberOfCaloClusters();
2655 for(Int_t icl=0; icl<ncl; icl++)
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();
2665 fIsEventTriggerBit = 1;
2666 cluster->SetChi2(1);
2669 cluster->SetChi2(0);
2676 //________________________________________________________________________
2678 Double_t AliAnalysisTaskFullppJet::GetOfflineTriggerProbability(AliESDCaloCluster *cluster)
2681 // Get the probablity of the given cluster to trigger the event
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
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
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
2715 //________________________________________________________________________
2717 Int_t AliAnalysisTaskFullppJet::GetClusterSuperModule(AliESDCaloCluster *cluster)
2720 // Return the given cluster supermodule
2724 cluster->GetPosition(pos);
2725 TVector3 clsVec(pos[0],pos[1],pos[2]);
2728 fGeom->SuperModuleNumberFromEtaPhi(clsVec.Eta(), clsVec.Phi(), sMod);
2733 //________________________________________________________________________
2734 Bool_t AliAnalysisTaskFullppJet::IsElectron(AliESDtrack *track, Double_t clsE) const
2737 // Check if the given track is an electron candidate based on de/dx
2740 if(track->GetTPCsignal()<=fdEdxMax && track->GetTPCsignal()>=fdEdxMin && (clsE/track->P())<=fEoverPMax && (clsE/track->P())>=fEoverPMin )
2747 //________________________________________________________________________
2748 Double_t AliAnalysisTaskFullppJet::GetTrkEff(Double_t inPt)
2751 // Get tracking efficiency estimated from simulation
2755 Double_t ptBound[4] = {0, 0.5, 3.8, 300};
2757 for(Int_t i=0; i<3; i++)
2759 if( inPt < ptBound[i+1] && inPt >= ptBound[i])
2761 eff = fTrkEffFunc[i]->Eval(inPt);
2768 //________________________________________________________________________
2769 Double_t AliAnalysisTaskFullppJet::GetSmearedTrackPt(AliESDtrack *track)
2772 // Smear track momentum
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);
2782 //________________________________________________________________________
2783 void AliAnalysisTaskFullppJet::CheckEventTriggerBit()
2786 // Check if the triggered events have correct trigger bit
2789 fIsEventTriggerBit = 0;
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;
2798 // part 1, trigger extraction -------------------------------------
2799 Int_t trigtimes[30], globCol, globRow, ntimes;
2800 Int_t trigger[nColsFaltroModule*nColsModule][nRowsFaltroModule*nRowsModule];
2802 // erase trigger maps
2803 for( Int_t i = 0; i < nColsFaltroModule*nColsModule; i++ )
2805 for( Int_t j = 0; j < nRowsFaltroModule*nRowsModule; j++ )
2811 Int_t fTrigCutLow = 7, fTrigCutHigh = 10, trigInCut = 0;
2812 AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
2813 // go through triggers
2814 if( fCaloTrigger->GetEntries() > 0 )
2817 fCaloTrigger->Reset();
2818 while( fCaloTrigger->Next() )
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;
2826 for( Int_t i = 0; i < ntimes; i++ )
2828 if( trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh ) // check if in cut
2833 if(trigInCut==1) trigger[globCol][globRow] = 1;
2834 } // calo trigger entries
2835 } // has calo trigger entries
2839 // part 2 go through the clusters here -----------------------------------
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++)
2847 AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2848 if(!cluster || !cluster->IsEMCAL()) continue;
2849 cluster->SetChi2(0);
2850 if(!IsGoodCluster(cluster)) continue;
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);
2861 if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5)
2863 nCell = cluster->GetNCells(); // get cluster cells
2864 cellAddrs = cluster->GetCellsAbsId(); // get the cell addresses
2865 for( iCell = 0; iCell < nCell; iCell++ )
2867 // get cell position
2868 fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2869 fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2871 // convert co global phi eta
2872 gphi = iphi + nRowsFeeModule*(nSupMod/2);
2873 geta = ieta + nColsFeeModule*(nSupMod%2);
2875 // get corresponding FALTRO
2878 // try to match with a triggered
2879 if( trigger[feta][fphi] == 1)
2881 cluster->SetChi2(1);
2882 fIsEventTriggerBit = 1;
2891 //________________________________________________________________________
2892 void AliAnalysisTaskFullppJet::PrintConfig()
2895 // Print configuration
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"};
2903 if(fStudySubEInHC || fStudyMcOverSubE)
2905 printf("\n\n=================================\n");
2906 printf("======WARNING: HC is ingored!======\n");
2907 printf("====== NOT for PHYSICS! ======\n\n");
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]);
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++)
2922 printf("Track pt cut: %s -> %2.2f < pT < %2.1f\n",triggerType[i], fTrkPtMin[i], fTrkPtMax[i]);
2924 for(Int_t i=0; i<2; i++)
2926 printf("Cluster selection: %s -> %2.2f < Et < %2.1f\n",triggerType[i],fClsEtMin[i], fClsEtMax[i]);
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");
2955 //________________________________________________________________________
2956 void AliAnalysisTaskFullppJet::CheckExoticEvent()
2959 // Check if the event containts exotic clusters
2962 Double_t leadingE = 0;
2963 Int_t nCls = fESD->GetNumberOfCaloClusters();
2964 for(Int_t icl=0; icl<nCls; icl++)
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)
2971 leadingE = cluster->E();
2974 if(leadingE>3) fIsExoticEvent3GeV = kTRUE;
2975 if(leadingE>5) fIsExoticEvent5GeV = kTRUE;
2979 //________________________________________________________________________
2980 Bool_t AliAnalysisTaskFullppJet::HasPrimaryVertex() const
2983 // Check if the primary vertex exists
2985 const AliESDVertex* vtx = fESD->GetPrimaryVertex();
2986 if (!vtx || vtx->GetNContributors()<1) return kFALSE;
2991 //________________________________________________________________________
2992 Bool_t AliAnalysisTaskFullppJet::IsPrimaryVertexOk() const
2995 // Check if the event vertex is good
3000 //________________________________________________________________________
3001 Bool_t AliAnalysisTaskFullppJet::IsTPCOnlyVtx() const
3004 // Check if the event only has valid TPC vertex
3007 const Bool_t isPrimaryVtx = HasPrimaryVertex();
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
3015 if( isPrimaryVtx && !isGoodVtx )
3021 //________________________________________________________________________
3022 void AliAnalysisTaskFullppJet::CheckTPCOnlyVtx(const UInt_t trigger)
3025 // Check the fraction of accepted events that have only TPC vertex
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);
3037 //________________________________________________________________________
3038 Bool_t AliAnalysisTaskFullppJet::IsLEDEvent() const
3041 // Check if the event is contaminated by LED signal
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++)
3049 Int_t cellId = cells->GetCellNumber(iCell);
3050 Double_t cellE = cells->GetCellAmplitude(cellId);
3051 Int_t sMod = fGeom->GetSuperModuleNumber(cellId);
3053 if(sMod==3 || sMod==4)
3056 nCellCount[sMod-3]++;
3059 Bool_t isLED=kFALSE;
3061 if(fPeriod.CompareTo("lhc11a")==0)
3063 if (nCellCount[1] > 100)
3065 Int_t runN = fESD->GetRunNumber();
3066 if (runN>=146858 && runN<=146860)
3068 if(fTriggerType==0 && nCellCount[0]>=21) isLED=kTRUE;
3069 if(fTriggerType==1 && nCellCount[0]>=35) isLED=kTRUE;