// **************************************
// Full pp jet task - ESD input only
-// Extract the jet spectrum and all the variations
+// Extract the jet spectrum and all the
+// systematic uncertainties
// -R. Ma, Mar 2011
// **************************************
#include "AliEMCALRecoUtils.h"
#include "TGeoManager.h"
#include "AliMCEvent.h"
+#include "AliMCParticle.h"
#include "AliStack.h"
#include "AliGenEventHeader.h"
#include "AliGenPythiaEventHeader.h"
ClassImp(AliAnalysisTaskFullppJet)
-const Float_t kRadius[2] = {0.4,0.2};
+const Float_t kRadius[3] = {0.4,0.2,0.3};
+const Int_t kNRadii = 3;
const Double_t kPI = TMath::Pi();
+const Double_t kdRCut[3] = {0.25,0.1,0.15};
//________________________________________________________________________
AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet() :
AliAnalysisTaskSE("default"),
- fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0),
- fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
- fIsMC(kFALSE), fMCStandalone(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
+ fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE),
+ fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
+ fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
fCentrality(99), fZVtxMax(10),
- fTriggerType(-1), fIsTPCOnlyVtx(kFALSE),
+ fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0),
fGeom(0x0), fRecoUtil(0x0),
fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0),fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9),
fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2),
fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE),
fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE),
- fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0), fMCNonLin(0x0),
+ fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
fJetNEFMin(0.01), fJetNEFMax(0.98),
fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
- fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0),
- fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0), fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
+ fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0),
+ fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0),
+ fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
- fSysNonLinearity(kFALSE), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
+ fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
+ fSysClusterizer(0),
fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
- fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventStat(0x0), fhEventCheck(0x0), fhTrkPhiEtaCheck(0x0), fhChunkQA(0x0)
+ fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
{
// Constructor
fTrkPtMax[i] = 200;
fClsEtMin[i] = 0.15;
fClsEtMax[i] = 200;
+
+ fVertexGenZ[i] = 0x0;
+ fEventZ[i] = 0x0;
+ fhNTrials[i] = 0x0;
+ fhNMatchedTrack[i] = 0x0;
+ fhClsE[i] = 0x0;
+ for(Int_t k=0; k<2; k++)
+ {
+ fhSysClusterE[i][k] = 0x0;
+ fhSysNCellVsClsE[i][k] = 0x0;
+ }
}
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
- fVertexGenZ[r] = 0x0;
- fEventZ[r] = 0x0;
- fhNTrials[r] = 0x0;
- fhNMatchedTrack[r] = 0x0;
- for(Int_t j=0; j<4; j++) fhSubEVsTrkPt[r][j] = 0x0;
for(Int_t j=0; j<5; j++)
{
fHCOverSubE[r][j] = 0x0;
fhHCTrkPtClean[k][r][j] = 0x0;
fhHCTrkPtAmbig[k][r][j] = 0x0;
}
- if(j>2) continue;
- fJetCount[j][r] = 0x0;
- fhNeutralPtInJet[j][r] = 0x0;
- fhChargedPtInJet[j][r] = 0x0;
- fhLeadNePtInJet[j][r] = 0x0;
- fhLeadChPtInJet[j][r] = 0x0;
- fJetEnergyFraction[j][r] = 0x0;
- fJetNPartFraction[j][r] = 0x0;
- if(j==2) continue;
- fRelTrkCon[j][r] = 0x0;
- fhFcrossVsZleading[j][r] = 0x0;
- fhChLeadZVsJetPt[j][r] = 0x0;
- fhJetPtWithTrkThres[j][r]= 0x0;
- fhJetPtWithClsThres[j][r]= 0x0;
- fhJetPtVsLowPtCons[j][r]= 0x0;
- fhJetPtInExoticEvent[r][j] = 0x0;
- }
- for(Int_t j=0; j<2; j++)
- {
- fhCorrTrkEffPtBin[j][r] = 0x0;
- for(Int_t i=0; i<kNBins; i++)
+
+ if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
+
+ if(j<3)
+ {
+ fJetCount[j][r] = 0x0;
+ fhNeutralPtInJet[j][r] = 0x0;
+ fhTrigNeuPtInJet[j][r] = 0x0;
+ fhChargedPtInJet[j][r] = 0x0;
+ fhLeadNePtInJet[j][r] = 0x0;
+ fhLeadChPtInJet[j][r] = 0x0;
+ fJetEnergyFraction[j][r] = 0x0;
+ fJetNPartFraction[j][r] = 0x0;
+ }
+ if(j<2)
{
- fhCorrTrkEffSample[j][r][i] = 0x0;
+ fRelTrkCon[j][r] = 0x0;
+ fhFcrossVsZleading[j][r] = 0x0;
+ fhChLeadZVsJetPt[j][r] = 0x0;
+ fhJetPtWithTrkThres[j][r]= 0x0;
+ fhJetPtWithClsThres[j][r]= 0x0;
+ for(Int_t k=0; k<2; k++)
+ {
+ fhJetPtVsLowPtCons[j][r][k] = 0x0;
+ }
+ fhJetPtInExoticEvent[j][r] = 0x0;
+ fhJetInTPCOnlyVtx[j][r] = 0x0;
+ fhCorrTrkEffPtBin[j][r] = 0x0;
+ for(Int_t i=0; i<kNBins; i++)
+ {
+ fhCorrTrkEffSample[j][r][i] = 0x0;
+ }
+ }
+
+ if(r==0 && j<3)
+ {
+ for(Int_t k=0; k<2; k++)
+ {
+ for(Int_t l=0; l<2; l++)
+ {
+ fhUEJetPtNorm[j][k][l] = 0x0;
+ fhUEJetPtVsSumPt[j][k][l] = 0x0;
+ fhUEJetPtVsConsPt[j][k][l] = 0x0;
+ }
+ }
}
}
- fhJetPtVsUE[r] = 0x0;
- fhClsE[r] = 0x0;
+ for(Int_t i=0; i<2; i++)
+ {
+ fhNKFracVsJetPt[i][r] = 0x0;
+ fhWeakFracVsJetPt[i][r] = 0x0;
+ fhJetResponseNK[i][r] = 0x0;
+ fhJetResponseWP[i][r] = 0x0;
+ fhJetResolutionNK[i][r] = 0x0;
+ fhJetResolutionWP[i][r] = 0x0;
+ fhJetResponseNKSM[i][r] = 0x0;
+ fhJetResponseWPSM[i][r] = 0x0;
+ fhJetResolutionNKSM[i][r] = 0x0;
+ fhJetResolutionWPSM[i][r] = 0x0;
+ }
}
for(Int_t s=0; s<3; s++)
{
for(Int_t a=0; a<2; a++)
{
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
fDetJetFinder[s][a][r] = 0x0;
fJetTCA[s][a][r] = 0x0;
for(Int_t j=0; j<3; j++)
{
fTrkEffFunc[j] = 0x0;
+ fhSecondaryResponse[j] = 0x0;
}
for(Int_t i=0; i<10; i++)
{
//________________________________________________________________________
AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet(const char *name) :
AliAnalysisTaskSE(name),
- fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0),
- fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
- fIsMC(kFALSE), fMCStandalone(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
+ fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE),
+ fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0),
+ fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.),
fCentrality(99), fZVtxMax(10),
- fTriggerType(-1), fIsTPCOnlyVtx(kFALSE),
+ fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0),
fGeom(0x0), fRecoUtil(0x0),
fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0), fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9),
fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2),
fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE),
fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE),
- fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0), fMCNonLin(0x0),
+ fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
fJetNEFMin(0.01), fJetNEFMax(0.98),
fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
- fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0),
- fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0), fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
+ fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0),
+ fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0),
+ fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.),
fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
- fSysNonLinearity(kFALSE), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
+ fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
+ fSysClusterizer(0),
fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
- fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventStat(0x0), fhEventCheck(0x0), fhTrkPhiEtaCheck(0x0), fhChunkQA(0x0)
+ fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
{
// Constructor
fTrkPtMax[i] = 200;
fClsEtMin[i] = 0.15;
fClsEtMax[i] = 200;
+
+ fVertexGenZ[i] = 0x0;
+ fEventZ[i] = 0x0;
+ fhNTrials[i] = 0x0;
+ fhNMatchedTrack[i] = 0x0;
+ fhClsE[i] = 0x0;
+ for(Int_t k=0; k<2; k++)
+ {
+ fhSysClusterE[i][k] = 0x0;
+ fhSysNCellVsClsE[i][k] = 0x0;
+ }
}
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
- fVertexGenZ[r] = 0x0;
- fEventZ[r] = 0x0;
- fhNTrials[r] = 0x0;
- fhNMatchedTrack[r] = 0x0;
- for(Int_t j=0; j<4; j++) fhSubEVsTrkPt[r][j] = 0x0;
for(Int_t j=0; j<5; j++)
{
fHCOverSubE[r][j] = 0x0;
fhHCTrkPtClean[k][r][j] = 0x0;
fhHCTrkPtAmbig[k][r][j] = 0x0;
}
- if(j>2) continue;
- fJetCount[j][r] = 0x0;
- fhNeutralPtInJet[j][r] = 0x0;
- fhChargedPtInJet[j][r] = 0x0;
- fhLeadNePtInJet[j][r] = 0x0;
- fhLeadChPtInJet[j][r] = 0x0;
- fJetEnergyFraction[j][r] = 0x0;
- fJetNPartFraction[j][r] = 0x0;
- if(j==2) continue;
- fRelTrkCon[j][r] = 0x0;
- fhFcrossVsZleading[j][r] = 0x0;
- fhChLeadZVsJetPt[j][r] = 0x0;
- fhJetPtWithTrkThres[j][r]= 0x0;
- fhJetPtWithClsThres[j][r]= 0x0;
- fhJetPtVsLowPtCons[j][r]= 0x0;
- fhJetPtInExoticEvent[r][j] = 0x0;
- }
- for(Int_t j=0; j<2; j++)
- {
- fhCorrTrkEffPtBin[j][r] = 0x0;
- for(Int_t i=0; i<kNBins; i++)
+
+ if(j<4) fhSubEVsTrkPt[r][j] = 0x0;
+
+ if(j<3)
+ {
+ fJetCount[j][r] = 0x0;
+ fhNeutralPtInJet[j][r] = 0x0;
+ fhTrigNeuPtInJet[j][r] = 0x0;
+ fhChargedPtInJet[j][r] = 0x0;
+ fhLeadNePtInJet[j][r] = 0x0;
+ fhLeadChPtInJet[j][r] = 0x0;
+ fJetEnergyFraction[j][r] = 0x0;
+ fJetNPartFraction[j][r] = 0x0;
+ }
+ if(j<2)
{
- fhCorrTrkEffSample[j][r][i] = 0x0;
+ fRelTrkCon[j][r] = 0x0;
+ fhFcrossVsZleading[j][r] = 0x0;
+ fhChLeadZVsJetPt[j][r] = 0x0;
+ fhJetPtWithTrkThres[j][r]= 0x0;
+ fhJetPtWithClsThres[j][r]= 0x0;
+ for(Int_t k=0; k<2; k++)
+ {
+ fhJetPtVsLowPtCons[j][r][k] = 0x0;
+ }
+ fhJetPtInExoticEvent[j][r] = 0x0;
+ fhJetInTPCOnlyVtx[j][r] = 0x0;
+ fhCorrTrkEffPtBin[j][r] = 0x0;
+ for(Int_t i=0; i<kNBins; i++)
+ {
+ fhCorrTrkEffSample[j][r][i] = 0x0;
+ }
+ }
+
+ if(r==0 && j<3)
+ {
+ for(Int_t k=0; k<2; k++)
+ {
+ for(Int_t l=0; l<2; l++)
+ {
+ fhUEJetPtNorm[j][k][l] = 0x0;
+ fhUEJetPtVsSumPt[j][k][l] = 0x0;
+ fhUEJetPtVsConsPt[j][k][l] = 0x0;
+ }
+ }
}
}
- fhJetPtVsUE[r] = 0x0;
- fhClsE[r] = 0x0;
+ for(Int_t i=0; i<2; i++)
+ {
+ fhNKFracVsJetPt[i][r] = 0x0;
+ fhWeakFracVsJetPt[i][r] = 0x0;
+ fhJetResponseNK[i][r] = 0x0;
+ fhJetResponseWP[i][r] = 0x0;
+ fhJetResolutionNK[i][r] = 0x0;
+ fhJetResolutionWP[i][r] = 0x0;
+ fhJetResponseNKSM[i][r] = 0x0;
+ fhJetResponseWPSM[i][r] = 0x0;
+ fhJetResolutionNKSM[i][r] = 0x0;
+ fhJetResolutionWPSM[i][r] = 0x0;
+ }
}
for(Int_t s=0; s<3; s++)
{
for(Int_t a=0; a<2; a++)
{
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
fDetJetFinder[s][a][r] = 0x0;
fJetTCA[s][a][r] = 0x0;
}
}
}
+
for(Int_t j=0; j<3; j++)
{
fTrkEffFunc[j] = 0x0;
+ fhSecondaryResponse[j] = 0x0;
}
for(Int_t i=0; i<10; i++)
{
{
for(Int_t a=0; a<2; a++)
{
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
if(fDetJetFinder[s][a][r]) delete fDetJetFinder[s][a][r];
if(fJetTCA[s][a][r]) { fJetTCA[s][a][r]->Delete(); delete fJetTCA[s][a][r]; }
for(Int_t i=0; i<3; i++)
{
if(fTrkEffFunc[i]) delete fTrkEffFunc[i];
+ if(fhSecondaryResponse[i]) delete fhSecondaryResponse[i];
}
for(Int_t i=0; i<10; i++)
{
if(fTriggerEfficiency[i]) delete fTriggerEfficiency[i];
if(fTriggerCurve[i]) delete fTriggerCurve[i];
}
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
for(Int_t j=0; j<2; j++)
{
if(fRecoUtil) delete fRecoUtil;
if(fClusterEResolution) delete fClusterEResolution;
- if(fMCNonLin) delete fMCNonLin;
+ if(fNonLinear) delete fNonLinear;
if(fTrkPtResData) delete fTrkPtResData;
}
// Create histograms
// Called once
//
+
+ if(fRunUE) fFindChargedOnlyJet = kTRUE;
const Int_t nTrkPtBins = 100;
const Float_t lowTrkPtBin=0, upTrkPtBin=100.;
fOutputList = new TList();
fOutputList->SetOwner(1);
- fhJetEventStat = new TH1F("fhJetEventStat","Event statistics for jet analysis",9,0,9);
+ fhJetEventStat = new TH1F("fhJetEventStat","Event statistics for jet analysis",12,0,12);
fhJetEventStat->GetXaxis()->SetBinLabel(1,"ALL");
fhJetEventStat->GetXaxis()->SetBinLabel(2,"MB");
fhJetEventStat->GetXaxis()->SetBinLabel(3,"MB+vtx+10cm");
fhJetEventStat->GetXaxis()->SetBinLabel(7,"EMC+vtx");
fhJetEventStat->GetXaxis()->SetBinLabel(8,"TriggerBit");
fhJetEventStat->GetXaxis()->SetBinLabel(9,"LED");
+ fhJetEventStat->GetXaxis()->SetBinLabel(10,"MB+TVtx+10cm");
+ fhJetEventStat->GetXaxis()->SetBinLabel(11,"EMC+TVtx+10cm");
+ fhJetEventStat->GetXaxis()->SetBinLabel(12,"ALL-Pileup");
fOutputList->Add(fhJetEventStat);
- fhEventCheck = new TH1F("fhEventCheck","Event statistics for check",11,0,11);
- fhEventCheck->GetXaxis()->SetBinLabel(1,"MB");
- fhEventCheck->GetXaxis()->SetBinLabel(2,"FastOnly");
- fhEventCheck->GetXaxis()->SetBinLabel(3,"FastOnly+TVtx");
- fhEventCheck->GetXaxis()->SetBinLabel(4,"FastOnly+PVtx");
- fhEventCheck->GetXaxis()->SetBinLabel(5,"FastOnly+TVtx+10cm");
- fhEventCheck->GetXaxis()->SetBinLabel(6,"FastOnly+PVtx+10cm");
- fhEventCheck->GetXaxis()->SetBinLabel(7,"!FastOnly");
- fhEventCheck->GetXaxis()->SetBinLabel(8,"!FastOnly+TVtx");
- fhEventCheck->GetXaxis()->SetBinLabel(9,"!FastOnly+PVtx");
- fhEventCheck->GetXaxis()->SetBinLabel(10,"!FastOnly+TVtx+10cm");
- fhEventCheck->GetXaxis()->SetBinLabel(11,"!FastOnly+PVtx+10cm");
- fOutputList->Add(fhEventCheck);
-
- fhTrkPhiEtaCheck = new TH2F("fhTrkPhiEtaCheck","#eta vs #phi of tracks in TPC only vertex events;#eta;#phi",100,-1,1,360,0,360);
- fOutputList->Add(fhTrkPhiEtaCheck);
+ fhJetEventCount = new TH1F("fhJetEventCount","Event statistics for jet analysis",12,0,12);
+ fhJetEventCount->GetXaxis()->SetBinLabel(1,"ALL");
+ fhJetEventCount->GetXaxis()->SetBinLabel(2,"MB");
+ fhJetEventCount->GetXaxis()->SetBinLabel(3,"MB-pileup");
+ fhJetEventCount->GetXaxis()->SetBinLabel(4,"MB+Vtx");
+ fhJetEventCount->GetXaxis()->SetBinLabel(5,"MB+Vtx+10cm");
+ fhJetEventCount->GetXaxis()->SetBinLabel(6,"EMC");
+ fhJetEventCount->GetXaxis()->SetBinLabel(7,"EMC-pileup");
+ fhJetEventCount->GetXaxis()->SetBinLabel(8,"EMC+Vtx");
+ fhJetEventCount->GetXaxis()->SetBinLabel(9,"EMC+Vtx+10cm");
+ fhJetEventCount->GetXaxis()->SetBinLabel(10,"Good EMC");
+ fOutputList->Add(fhJetEventCount);
+
+ if(fCheckTPCOnlyVtx)
+ {
+ fhEventStatTPCVtx = new TH1F("fhEventStatTPCVtx","Event statistics for TPC only vertex",9,0,9);
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(1,"FastOnly");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(2,"FastOnly+PVtx");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(3,"FastOnly+TVtx");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(4,"MB");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(5,"MB+PVtx");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(6,"MB+TVtx");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(7,"EMC");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(8,"EMC+PVtx");
+ fhEventStatTPCVtx->GetXaxis()->SetBinLabel(9,"EMC+TVtx");
+ fOutputList->Add(fhEventStatTPCVtx);
+ }
fhChunkQA = new TH1F("fhChunkQA","# of hybrid tracks per chunk",200,0,200);
fOutputList->Add(fhChunkQA);
Double_t lowBin2[dim2] = {0,0,0,};
Double_t upBin2[dim2] = {200,0.5,50};
- Int_t radiusType = 1;
- if( fRadius.Contains("0.2") ) radiusType = 2;
const char* triggerName[3] = {"MB","EMC","MC"};
const char* triggerTitle[3] = {"MB","EMCal-trigger","MC true"};
const char* fraction[5] = {"MIP","30","50","70","100"};
const char* exotic[2] = {"3GeV","5GeV"};
const char* vertexType[2] = {"All MB","MB with vertex"};
const char *vertexName[2] = {"All","Vertex"};
-
+ const char *clusterizerName[2] = {"before","after"};
+ const char *UEName[2] = {"charged","charged_neutral"};
+ const char *UETitle[2] = {"charged","charged+neutral"};
+ const char *UEEventName[2] = {"LeadingJet","Back-To-Back"};
if(fIsMC)
{
fOutputList->Add(fEventZ[i]);
}
- for(Int_t r=0; r<radiusType; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
+ if(!fRadius.Contains("0.4") && r==0) continue;
+ if(!fRadius.Contains("0.2") && r==1) continue;
+ if(!fRadius.Contains("0.3") && r==2) continue;
+
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);
fOutputList->Add(fJetCount[i][r]);
- 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,lowTrkPtBin,upTrkPtBin);
+ 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);
fOutputList->Add(fhNeutralPtInJet[i][r]);
+
+ 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);
+ fOutputList->Add(fhTrigNeuPtInJet[i][r]);
- 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,lowTrkPtBin,upTrkPtBin);
+ 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);
fOutputList->Add(fhChargedPtInJet[i][r]);
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);
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);
fOutputList->Add(fhJetPtWithClsThres[i][r]);
- fhJetPtVsLowPtCons[i][r] = new TH2F(Form("%s_fhJetPtVsLowPtCons_%1.1f",triggerName[i],kRadius[r]),Form("%s: energy carried by constituents less than 200MeV (EMCal, R=%1.1f,Z<0.98);p_{T}^{low} (GeV/c);p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,100);
- fOutputList->Add(fhJetPtVsLowPtCons[i][r]);
+ 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);
+ fOutputList->Add(fhJetPtVsLowPtCons[i][r][0]);
+
+ 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);
+ fOutputList->Add(fhJetPtVsLowPtCons[i][r][1]);
+
+ if(fCheckTPCOnlyVtx)
+ {
+ 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);
+ fOutputList->Add(fhJetInTPCOnlyVtx[i][r]);
+ }
if(fStudySubEInHC)
{
}
}
}
- if(fRunUE && fFindChargedOnlyJet && i<2)
+
+ if(fRunUE)
{
- fhJetPtVsUE[i] = new TH2F(Form("%s_fhJetPtVsUE",triggerName[i]),Form("%s: jet p_{T} vs UE (R=0.4);p_{T,ch}^{UE} (GeV/c);p_{T,jet} (GeV/c)",triggerTitle[i]),2000,0,20,nbins,xbins);
- fOutputList->Add(fhJetPtVsUE[i]);
+ for(Int_t k=0; k<2; k++)
+ {
+ for(Int_t l=0; l<2; l++)
+ {
+ 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);
+ fOutputList->Add(fhUEJetPtNorm[i][k][l]);
+
+ 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);
+ fOutputList->Add(fhUEJetPtVsSumPt[i][k][l]);
+
+ 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);
+ fOutputList->Add(fhUEJetPtVsConsPt[i][k][l]);
+ }
+ }
}
+
if(fSysJetTrigEff)
{
fhClsE[i] = new TH1F(Form("%s_fhClsE",triggerName[i]),Form("%s: cluster E;E (GeV)",triggerTitle[i]),1000,0,100);
fOutputList->Add(fhClsE[i]);
}
+
+ if(fSysClusterizer)
+ {
+ for(Int_t k=0; k<2; k++)
+ {
+ 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);
+ fOutputList->Add(fhSysClusterE[i][k]);
+
+ 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);
+ fOutputList->Add(fhSysNCellVsClsE[i][k]);
+ }
+ }
}
{
if(fStudyMcOverSubE)
{
- for(Int_t r=0; r<radiusType; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
+ if(!fRadius.Contains("0.4") && r==0) continue;
+ if(!fRadius.Contains("0.2") && r==1) continue;
+ if(!fRadius.Contains("0.3") && r==2) continue;
for(Int_t i=0; i<5; i++)
{
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);
}
}
}
+ if(fRunSecondaries && fAnaType==0)
+ {
+ for(Int_t i=0; i<2; i++)
+ {
+ for(Int_t r=0; r<3; r++)
+ {
+ 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);
+ fOutputList->Add(fhNKFracVsJetPt[i][r]);
+
+ 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);
+ fOutputList->Add(fhWeakFracVsJetPt[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResponseNK[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResponseWP[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResolutionNK[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResolutionWP[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResponseNKSM[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResponseWPSM[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResolutionNKSM[i][r]);
+
+ 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);
+ fOutputList->Add(fhJetResolutionWPSM[i][r]);
+ }
+ }
+ }
}
printf("\n=======================================\n");
printf("===== Jet task configuration ==========\n");
+
if(fNonStdBranch.Length()!=0)
{
const char* species[3] = {"in","ch","ne"};
const char* algorithm[2] = {"akt","kt"};
- const char* radii[2] = {"04","02"};
+ const char* radii[kNRadii] = {"04","02","03"};
for(Int_t s=0; s<3; s++)
{
if(!fFindChargedOnlyJet && s==1) continue;
{
if(!fAlgorithm.Contains("aKt") && a==0) continue;
if(!fAlgorithm.Contains("kt") && a==1) continue;
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
if(!fRadius.Contains("0.4") && r==0) continue;
if(!fRadius.Contains("0.2") && r==1) continue;
- fJetTCA[s][a][r] = new TClonesArray("AliAODJet",0);
- fJetTCA[s][a][r]->SetName(Form("Jet_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
- AddAODBranch("TClonesArray",&fJetTCA[s][a][r],fNonStdFile.Data());
- printf("Add branch: Jet_%s_%s_%s_%s\n",species[s],algorithm[a],radii[r],fNonStdBranch.Data());
+ if(!fRadius.Contains("0.3") && r==2) continue;
+ if(fAnaType==0)
+ {
+ fJetTCA[s][a][r] = new TClonesArray("AliAODJet",0);
+ fJetTCA[s][a][r]->SetName(Form("Jet_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()));
+ AddAODBranch("TClonesArray",&fJetTCA[s][a][r],fNonStdFile.Data());
+ printf("Add branch: Jet_%s_%s_%s_%s\n",species[s],algorithm[a],radii[r],fNonStdBranch.Data());
+ }
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()));
if(a==0) fDetJetFinder[s][a][r]->SetAlgorithm(fastjet::antikt_algorithm);
{
if(fIsMC && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase))
{
- fMcTruthAntikt[r] = new TClonesArray("AliAODJet",0);
- fMcTruthAntikt[r]->SetName(Form("Jet_mc_truth_in_akt_%s_%s",radii[r],fNonStdBranch.Data()));
- AddAODBranch("TClonesArray",&fMcTruthAntikt[r],fNonStdFile.Data());
- printf("Add branch: Jet_mc_truth_in_akt_%s_%s\n",radii[r],fNonStdBranch.Data());
+ if(fAnaType==0)
+ {
+ fMcTruthAntikt[r] = new TClonesArray("AliAODJet",0);
+ fMcTruthAntikt[r]->SetName(Form("Jet_mc_truth_in_akt_%s_%s",radii[r],fNonStdBranch.Data()));
+ AddAODBranch("TClonesArray",&fMcTruthAntikt[r],fNonStdFile.Data());
+ printf("Add branch: Jet_mc_truth_in_akt_%s_%s\n",radii[r],fNonStdBranch.Data());
+ }
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()));
fTrueJetFinder[r]->SetAlgorithm(fastjet::antikt_algorithm);
}
fRandomGen = new TRandom3(0);
- if(fCheckTrkEffCorr)
+ if(fCheckTrkEffCorr && fAnaType==0)
{
TFile f ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffFit.root","read");
for(Int_t i=0; i<3; i++)
}
}
+ if(fRunSecondaries && fAnaType==0)
+ {
+ const char *secondaryName[3] = {"k0S","lamda","hyperon"};
+ TFile f2 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/SecondaryResponse.root","read");
+ for(Int_t j=0; j<3; j++)
+ {
+ fhSecondaryResponse[j] = new TH2F(*((TH2F*)f2.Get(Form("DetectorReponse_%s",secondaryName[j]))));
+ }
+ }
- TFile f ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TriggerCurve.root","read");
- fTriggerMask = new TH2F(*((TH2F*)f.Get("lhc11a_TriggerMask")));
- if(fOfflineTrigger)
+ if(fCheckTriggerMask)
{
- for(Int_t i=0; i<10; i++)
+ char *name = "TriggerCurve.root";
+ if(fAnaType==0) name = "/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TriggerCurve.root";
+ TFile f (name,"read");
+ fTriggerMask = new TH2F(*((TH2F*)f.Get("lhc11a_TriggerMask")));
+ if(fOfflineTrigger)
{
- fTriggerEfficiency[i] = new TF1(*((TF1*)f.Get(Form("lhc11a_TriggerEfficiency_SM%d_fit",i))));
- fTriggerCurve[i] = new TH1D(*((TH1D*)f.Get(Form("lhc11a_TriggerCurve_SM%d",i))));
- }
+ for(Int_t i=0; i<10; i++)
+ {
+ fTriggerEfficiency[i] = new TF1(*((TF1*)f.Get(Form("lhc11a_TriggerEfficiency_SM%d_fit",i))));
+ fTriggerCurve[i] = new TH1D(*((TH1D*)f.Get(Form("lhc11a_TriggerCurve_SM%d",i))));
+ }
+ }
+ f.Close();
}
- f.Close();
fClusterEResolution = new TF1("fClusterEResolution","sqrt([0]^2+[1]^2*x+([2]*x)^2)*0.01");
fClusterEResolution->SetParameters(4.35,9.07,1.63);
- fMCNonLin = new TF1("f1","([0]/((x+[1])^[2]))+1",0.1,200.0);
- fMCNonLin->SetParameters(3.11111e-02, -5.71666e-02, 5.67995e-01);
-
fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
+ if(!fGeom)
+ {
+ AliError("No EMCal geometry is available!");
+ return;
+ }
fRecoUtil = new AliEMCALRecoUtils();
- if(fRejectExoticCluster)
+ if(fRejectExoticCluster || fRejectExoticTrigger)
fRecoUtil->SwitchOnRejectExoticCluster();
else
{
fRecoUtil->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrected);
if(fSysNonLinearity)
{
- fRecoUtil->SetNonLinearityParam(0,9.88165e-01);
- fRecoUtil->SetNonLinearityParam(1,3.07553e-01);
- fRecoUtil->SetNonLinearityParam(2,4.91690e-01);
- fRecoUtil->SetNonLinearityParam(3,1.07910e-01);
- fRecoUtil->SetNonLinearityParam(4,1.54119e+02);
- fRecoUtil->SetNonLinearityParam(5,2.91142e+01);
+ fNonLinear = new TF1("TB_oldBest","([0])*(1./(1.+[1]*exp(-x/[2]))*1./(1.+[3]*exp((x-[4])/[5])))*1/([6])",0.1,110.);
+ fNonLinear->SetParameters(0.99078, 0.161499, 0.655166, 0.134101, 163.282, 23.6904, 0.978);
}
if(fSysTrkClsMth)
{
PrintConfig();
BookHistos();
- //PostData(0, AODEvent());
PostData(1, fOutputList);
}
// Main loop, called for each event.
//
- if(fPeriod.CompareTo("lhc11b10a",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11b10b",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc12a15a",TString::kIgnoreCase)==0)
- {
- fIsMC = kTRUE;
- }
-
// Get event pointers, check for signs of life
Double_t vtxTrueZ = -100;
if(fIsMC)
fMC = MCEvent();
if (!fMC)
{
- Printf("ERROR: Could not retrieve MC event");
+ printf("ERROR: Could not retrieve MC event");
return;
}
fStack = fMC->Stack();
TParticle *particle = fStack->Particle(0);
- if(particle) vtxTrueZ = particle->Vz();
- //printf("Generated vertex coordinate: (x,y,z) = (%4.2f, %4.2f, %4.2f)\n", particle->Vx(), particle->Vy(), particle->Vz());
+ if(particle) vtxTrueZ = particle->Vz(); // True vertex z in MC events
+ if(fVerbosity>10) printf("Generated vertex coordinate: (x,y,z) = (%4.2f, %4.2f, %4.2f)\n", particle->Vx(), particle->Vy(), particle->Vz());
}
fESD = dynamic_cast<AliESDEvent*>(InputEvent());
}
fhJetEventStat->Fill(0.5);
+ fhJetEventCount->Fill(0.5);
+ if(fIsMC)
+ {
+ GetMCInfo();
+ if(fVertexGenZ[0]) fVertexGenZ[0]->Fill(vtxTrueZ);
+ }
+
// Centrality, vertex, other event variables...
if(fESD)
{
UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
if (trigger==0) return;
+ if(fCheckTPCOnlyVtx) CheckTPCOnlyVtx(trigger);
if(!fIsMC)
{
if (trigger & AliVEvent::kFastOnly) return; // Reject fast trigger cluster
}
else
{
- if (trigger & AliVEvent::kMB)
+ if(!fPhySelForMC) fTriggerType = 0;
+ else if (trigger & AliVEvent::kAnyINT) fTriggerType = 0;
+ else fTriggerType = -1;
+
+ if(fOfflineTrigger)
{
- fTriggerType = 0;
- if(fOfflineTrigger) fTriggerType = RunOfflineTrigger();
+ RunOfflineTrigger();
+ if(fIsEventTriggerBit) fTriggerType = 1;
}
- else fTriggerType = -1;
}
if(fTriggerType==-1)
{
- //printf("Error: worng trigger type %s\n",(fESD->GetFiredTriggerClasses()).Data());
+ if(fVerbosity>10) printf("Error: worng trigger type %s\n",(fESD->GetFiredTriggerClasses()).Data());
return;
}
}
- if(fIsMC)
- {
- GetMCInfo();
- if(fVertexGenZ[0]) fVertexGenZ[0]->Fill(vtxTrueZ);
- }
+
+ fhJetEventCount->Fill(1.5+fTriggerType*4);
+ if(fRejectPileup && fESD->IsPileupFromSPD() ) return; // reject pileup
+ fhJetEventStat->Fill(11.5);
+ fhJetEventCount->Fill(2.5+fTriggerType*4);
+
+
+ fIsTPCOnlyVtx = 0;
+ // Reject LED events
if (IsLEDEvent())
{
fhJetEventStat->Fill(8.5);
fhJetEventStat->Fill(1.5+fTriggerType*2);
- if(!IsPrimaryVertexOk(vtxTrueZ)) return;
+ // Check if primary vertex exists
+ if (!HasPrimaryVertex()) return;
+ fhJetEventStat->Fill(5.5+fTriggerType);
+ fhJetEventCount->Fill(3.5+fTriggerType*4);
+
+ const AliESDVertex* vtx = fESD->GetPrimaryVertex();
+ Double_t zVertex = vtx->GetZ();
+ if(fEventZ[fTriggerType]) fEventZ[fTriggerType]->Fill(zVertex);
+ if(fVertexGenZ[1]) fVertexGenZ[1]->Fill(vtxTrueZ);
+
+ // Check if |Z_vtx|<10cm
+ if( TMath::Abs(zVertex) > fZVtxMax ) return;
+ fhJetEventCount->Fill(4.5+fTriggerType*4);
+
+ // Check if only TPC vertex exists primitive
fIsTPCOnlyVtx = IsTPCOnlyVtx();
if(!fIsMC && fTriggerType==1)
{
+ // Check if event has valid trigger bit
CheckEventTriggerBit();
- if(fIsEventTriggerBit) fhJetEventStat->Fill(7.5);
+ if(fIsEventTriggerBit)
+ {
+ fhJetEventStat->Fill(7.5);
+ fhJetEventCount->Fill(9.5);
+ }
else return;
}
fhJetEventStat->Fill(2.5+fTriggerType*2);
+ if(fIsTPCOnlyVtx) fhJetEventStat->Fill(9.5+fTriggerType);
+ // Check if event contains exotic clusters
CheckExoticEvent();
- AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
+ // Write jet tree into AOD output in local analysis mode
+ if(fAnaType==0) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
- // Fast Jet calls BEG --------------------------------------------------------
-
+ // Clean up arrays
for(Int_t s=0; s<3; s++)
{
for(Int_t a=0; a<2; a++)
{
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
- if(fJetTCA[s][a][r])
- {
- fJetTCA[s][a][r]->Delete();
- fDetJetFinder[s][a][r]->Clear();
- }
+ if(fJetTCA[s][a][r]) fJetTCA[s][a][r]->Delete();
+ if(fDetJetFinder[s][a][r]) fDetJetFinder[s][a][r]->Clear();
if(s==0 && a==0)
{
- if(fMcTruthAntikt[r])
- {
- fMcTruthAntikt[r]->Delete();
- fTrueJetFinder[r]->Clear();
- }
+ if(fMcTruthAntikt[r]) fMcTruthAntikt[r]->Delete();
+ if(fTrueJetFinder[r]) fTrueJetFinder[r]->Clear();
}
}
}
}
- if(fVerbosity>10) printf("# of jets after clear: %d\n",fJetTCA[0][0][0]->GetEntries());
+ if(fVerbosity>5) printf("# of jets after clear: %d\n",fJetTCA[0][0][0]->GetEntries());
if(fTrackArray) fTrackArray->Delete();
if(fClusterArray) fClusterArray->Delete();
fMcPartArray->Reset();
+
if (fESD)
{
// get the tracks and fill the input vector for the jet finders
// get EMCal clusters and fill the input vector for the jet finders
GetESDEMCalClusters();
-
+
for(Int_t s=0; s<3; s++)
{
for(Int_t a=0; a<2; a++)
{
- for(Int_t r=0; r<2; r++)
+ for(Int_t r=0; r<kNRadii; r++)
{
//Detector jets
- if(fJetTCA[s][a][r])
+ if(fDetJetFinder[s][a][r])
{
FindDetJets(s,a,r);
- FillAODJets(fJetTCA[s][a][r], fDetJetFinder[s][a][r], 0);
+ if(fJetTCA[s][a][r]) FillAODJets(fJetTCA[s][a][r], fDetJetFinder[s][a][r], 0);
if(s==0 && a==0 && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase))
{
- if(fSaveQAHistos) AnalyzeJets(fDetJetFinder[s][a][r],fTriggerType, r);
- // s == 1 cannot be true here if(fRunUE && fFindChargedOnlyJet && s==1 && a==0 && r==0) RunAnalyzeUE(fDetJetFinder[s][a][r]);
+ if(fSaveQAHistos) AnalyzeJets(fDetJetFinder[s][a][r],fTriggerType, r); // run analysis on jets
}
}
}
}
}
+
+ if(fRunUE && fDetJetFinder[1][0][0]) RunAnalyzeUE(fDetJetFinder[1][0][0],fTriggerType,0); // Run analysis on underlying event
}
if(fIsMC)
{
- for(Int_t r=0; r<2; r++)
- if(fMcTruthAntikt[r]) ProcessMC(r);
+ for(Int_t r=0; r<kNRadii; r++)
+ if(fTrueJetFinder[r]) ProcessMC(r); // find particle level jets
}
// Fast Jet calls END --------------------------------------------------------
if(isCh)
{
+ // Feed in charged tracks
Int_t countTracks = 0;
Int_t ntracks = fTrackArray->GetEntriesFast();
for(Int_t it=0; it<ntracks; it++)
{
AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
- if(!t || t->Pt() < fTrkPtMin[fTriggerType]) continue;
+ if(!t) continue;
countTracks++;
Double_t e = t->P();
Double_t pt = t->Pt();
+ if(s==1 && fRunUE && pt<1) continue;
if(fRecombinationScheme==0) e = TMath::Sqrt(t->P()*t->P()+0.139*0.139);
if(fSysTrkPtRes) pt += GetSmearedTrackPt(t);
Double_t phi = t->Phi();
Double_t py = pt*TMath::Sin(phi);
if(fConstrainChInEMCal && ( TMath::Abs(eta)>0.7 || phi>kPI || phi<TMath::DegToRad()*80) ) continue;
fDetJetFinder[s][a][r]->AddInputVector(px,py,t->Pz(), e, it+1);
- //printf("%d: m_{ch}=%f\n",it+1,t->P()*t->P()-t->Px()*t->Px()-t->Py()*t->Py()-t->Pz()*t->Pz());
+ 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());
}
if(fVerbosity>5) printf("[i] # of tracks filled: %d\n",countTracks);
}
if(isNe)
{
+ // Feed in EMCal clusters
Double_t vertex[3] = {0, 0, 0};
fESD->GetVertex()->GetXYZ(vertex) ;
TLorentzVector gamma;
}
if(fVerbosity>5) printf("[i] # of clusters filled: %d\n",countClusters);
}
+ // Run jet finding
fDetJetFinder[s][a][r]->Run();
}
//
// Fill the found jets into AOD output
// Only consider jets pointing to EMCal and with pt above 1GeV/c
+ //
Int_t radiusIndex = 0;
TString arrayName = fJetArray->GetName();
if(arrayName.Contains("_02_")) radiusIndex = 1;
+ if(arrayName.Contains("_03_")) radiusIndex = 2;
std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
- if(fVerbosity>0) printf("[i] # of jets in %s : %d\n",fJetArray->GetName(),(Int_t)jetsIncl.size());
+ if(fVerbosity>5 && radiusIndex==0) printf("[i] # of jets in %s : %d\n",fJetArray->GetName(),(Int_t)jetsIncl.size());
AliAODJet *jet = 0;
Int_t jetCount = 0;
for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
{
- //printf("fastjet: eta=%f, phi=%f\n",jetsIncl[ij].eta(),jetsIncl[ij].phi()*TMath::RadToDeg());
+ if(fVerbosity>10) printf("fastjet: eta=%f, phi=%f\n",jetsIncl[ij].eta(),jetsIncl[ij].phi()*TMath::RadToDeg());
if(!IsGoodJet(jetsIncl[ij],0)) continue;
AliAODJet tmpJet (jetsIncl[ij].px(), jetsIncl[ij].py(), jetsIncl[ij].pz(), jetsIncl[ij].E());
jet = new ((*fJetArray)[jetCount]) AliAODJet(tmpJet);
jetCount++;
- //printf("AOD jet: ij=%d, pt=%f, eta=%f, phi=%f\n",ij, jet->Pt(), jet->Eta(),jet->Phi()*TMath::RadToDeg());
+ 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());
jet->GetRefTracks()->Clear();
- std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
+ std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
Double_t totE=0, totPt=0, totChPt=0, leadChPt=0, neE=0, totNePt=0, leadNePt=0, leadPt=0;
Double_t leadTrkType=0, nChPart = 0, nNePart = 0;
Bool_t isHighPtTrigger = kFALSE, isTriggering = kFALSE, isHyperTrack = kFALSE;
+ Int_t leadIndex = -1;
for(UInt_t ic=0; ic<constituents.size(); ic++)
{
- //printf("ic=%d: user_index=%d, E=%f\n",ic,constituents[ic].user_index(),constituents[ic].E());
+ if(fVerbosity>10 && radiusIndex==0) printf("ic=%d: user_index=%d, E=%f\n",ic,constituents[ic].user_index(),constituents[ic].E());
if(constituents[ic].user_index()<0)
{
nNePart ++;
if(!isTruth)
{
AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
- if(cluster->Chi2()==1) isTriggering = kTRUE;
+ if(cluster->Chi2()>0.5) isTriggering = kTRUE;
}
}
else
totE += constituents[ic].E();
totPt += constituents[ic].perp();
if(constituents[ic].perp()>leadPt)
- leadPt = constituents[ic].perp();
+ {
+ leadPt = constituents[ic].perp();
+ leadIndex = ic;
+ }
}
if(fIsEventTriggerBit) jet->SetTrigger(AliAODJet::kTRDTriggered);
if(isHighPtTrigger) jet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
if(isTriggering) jet->SetTrigger(AliAnalysisTaskFullppJet::kTrigger);
if(isHyperTrack) jet->SetTrigger(AliAnalysisTaskFullppJet::kSuspicious);
if(fIsTPCOnlyVtx) jet->SetTrigger(AliAnalysisTaskFullppJet::kTPCOnlyVtx);
+ if(constituents[leadIndex].user_index()>0) jet->SetTrigger(AliAnalysisTaskFullppJet::kLeadCh);
+
if(jetsIncl[ij].E()>0) jet->SetNEF(neE/jetsIncl[ij].E());
+ jet->SetPtLeading(leadPt);
+ jet->SetPtSubtracted(leadTrkType,0);
Double_t effAErrCh = 0, effAErrNe = leadTrkType;
- Double_t chBgEnergy = 0, neBgEnergy = nNePart;
+ Double_t chBgEnergy = 10, neBgEnergy = nNePart;
if(!isTruth)
{
effAErrCh = GetMeasuredJetPtResolution(ij,jetFinder);
- if(fCheckTrkEffCorr)
- chBgEnergy = GetJetMissPtDueToTrackingEfficiency(ij,jetFinder,radiusIndex);
+ if(fIsMC && constituents[leadIndex].user_index()>0)
+ {
+ AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[leadIndex].user_index()-1);
+ Double_t pt = track->Pt();
+ Int_t ipart = track->GetLabel();
+ if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
+ {
+ AliVParticle* vParticle = fMC->GetTrack(ipart);
+ if(vParticle)
+ {
+ Double_t truePt = vParticle->Pt();
+ chBgEnergy = (truePt-pt)/truePt;
+ }
+ }
+ }
}
+
jet->SetEffArea(leadPt, jetFinder->GetJetArea(ij),effAErrCh,effAErrNe);
jet->SetBgEnergy(chBgEnergy,neBgEnergy);
- //cout<<"jet pt="<<jetsIncl[ij].perp()<<", nef="<<jet->GetNEF()<<", trk eff corr="<<chBgEnergy<<endl;
+ if(fVerbosity>10 && isTruth) cout<<jet->ErrorEffectiveAreaCharged()<<" "<<jet->ErrorEffectiveAreaNeutral()<<endl;
+ if(fVerbosity>10) cout<<"jet pt="<<jetsIncl[ij].perp()<<", nef="<<jet->GetNEF()<<", trk eff corr="<<chBgEnergy<<endl;
- if(fVerbosity>0)
+ if(fVerbosity>5)
printf("# of ref tracks: %d = %d, and nef=%f\n",jet->GetRefTracks()->GetEntries(), (Int_t)constituents.size(), jet->GetNEF());
// For catch good high-E jets
}
}
// End of catching good high-E jets
- }
+ }
+ if(fVerbosity>5) printf("%s has %d jets\n",fJetArray->GetName(), fJetArray->GetEntries());
}
if( fChargedMC && vParticle->Charge()==0 ) continue;
Int_t index = 1;
- if(vParticle->Charge()==0) { index=-1; }
- fTrueJetFinder[r]->AddInputVector(vParticle->Px(), vParticle->Py(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
- //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());
-
+ if(vParticle->Charge()==0) { index=-1; }
fMcPartArray->AddAt(ipart, countPart);
countPart++;
+
+ fTrueJetFinder[r]->AddInputVector(vParticle->Px(), vParticle->Py(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
+ if(fVerbosity>10)
+ 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());
}
fMcPartArray->Set(countPart);
fTrueJetFinder[r]->Run();
- FillAODJets(fMcTruthAntikt[r], fTrueJetFinder[r], 1);
- if(fSaveQAHistos) AnalyzeJets(fTrueJetFinder[r], 2, r);
+ if(fMcTruthAntikt[r]) FillAODJets(fMcTruthAntikt[r], fTrueJetFinder[r], 1);
+ if(fSaveQAHistos) AnalyzeJets(fTrueJetFinder[r], 2, r);
+ if(fRunUE && r==0) RunAnalyzeUE(fTrueJetFinder[r], 2, 1);
+
+ // Run analysis on secondary particles
+ if(fRunSecondaries && fAnaType==0)
+ {
+ for(Int_t i=0; i<2; i++)
+ {
+ AnalyzeSecondaryContribution(fTrueJetFinder[r],r,i);
+ AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,0,i);
+ AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,1,i);
+ }
+ }
}
for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
{
+ if(type<2 && fCheckTPCOnlyVtx && fIsTPCOnlyVtx)
+ {
+ fhJetInTPCOnlyVtx[type][r]->Fill(jetsIncl[ij].perp(),jetsIncl[ij].phi()*TMath::RadToDeg(),jetsIncl[ij].eta());
+ }
if(!IsGoodJet(jetsIncl[ij],kRadius[r])) continue; // Fidiual cut
Float_t jetEta = jetsIncl[ij].eta();
Float_t jetPhi = jetsIncl[ij].phi();
Float_t jetE = jetsIncl[ij].E();
Float_t jetPt = jetsIncl[ij].perp();
- std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
+ std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
Double_t neE=0, leadingZ = 0, maxPt = 0;
Int_t constituentType = -1, leadingIndex = 0;
for(UInt_t ic=0; ic<constituents.size(); ic++)
if(leadingZ > 0.98) continue; // Z cut
fJetCount[type][r]->Fill(jetPt);
- if(type==1 && fIsExoticEvent3GeV) fhJetPtInExoticEvent[r][0]->Fill(jetPt);
- if(type==1 && fIsExoticEvent5GeV) fhJetPtInExoticEvent[r][1]->Fill(jetPt);
+ if(type==1 && fIsExoticEvent3GeV) fhJetPtInExoticEvent[0][r]->Fill(jetPt);
+ if(type==1 && fIsExoticEvent5GeV) fhJetPtInExoticEvent[1][r]->Fill(jetPt);
Double_t totTrkPt[3] = {0.,0.,0.};
Double_t chPt = 0;
fhNeutralPtInJet[type][r]->Fill(jetPt, partPt);
if(partPt>leadNePt)
leadNePt = partPt;
- if((fStudySubEInHC || fStudyMcOverSubE) && type<2)
+ if(type<2)
{
AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
Double_t clsE = cluster->E();
- cluster->GetMomentum(gamma, vertex);
- Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
-
- Double_t subEtmp = cluster->GetDispersion();
- Double_t mipETmp = cluster->GetEmcCpvDistance();
- mcSubE += cluster->GetDistanceToBadChannel()*sinTheta;
- subClsE[0] += ((mipETmp>clsE)?clsE:mipETmp)*sinTheta;
- for(Int_t j=1; j<5; j++)
+ if(cluster->Chi2()>0.5) fhTrigNeuPtInJet[type][r]->Fill(jetPt, partPt);
+ if(fStudySubEInHC || fStudyMcOverSubE)
{
- subClsE[j] += ((fraction[j]*subEtmp>clsE)?clsE:fraction[j]*subEtmp)*sinTheta;
+ cluster->GetMomentum(gamma, vertex);
+ Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
+
+ Double_t subEtmp = cluster->GetDispersion();
+ Double_t mipETmp = cluster->GetEmcCpvDistance();
+ mcSubE += cluster->GetDistanceToBadChannel()*sinTheta;
+ subClsE[0] += ((mipETmp>clsE)?clsE:mipETmp)*sinTheta;
+ for(Int_t j=1; j<5; j++)
+ {
+ subClsE[j] += ((fraction[j]*subEtmp>clsE)?clsE:fraction[j]*subEtmp)*sinTheta;
+ }
}
}
}
Double_t thres[3] = {15,25,40};
Int_t okCh[3] = {0,0,0};
Int_t okNe[3] = {0,0,0};
- Double_t lowPt = 0;
+ Double_t lowPt[2] = {0.,0.};
for(UInt_t ic=0; ic<constituents.size(); ic++)
{
Float_t partPt = constituents[ic].perp();
- if(partPt<0.2) lowPt += partPt;
+
+ if(partPt<0.3) lowPt[0] += partPt;
+ else if(partPt<0.5) lowPt[1] += partPt;
+
if(constituents[ic].user_index()>0)
{
for(Int_t it=0; it<3; it++)
if(okNe[i]==1)
fhJetPtWithClsThres[type][r]->Fill(i,jetPt);
}
- fhJetPtVsLowPtCons[type][r]->Fill(jetPt,lowPt);
+ for(Int_t k=0; k<2; k++) fhJetPtVsLowPtCons[type][r][k]->Fill(jetPt,lowPt[k]/jetPt);
}
}
}
//________________________________________________________________________
-void AliAnalysisTaskFullppJet::RunAnalyzeUE(AliFJWrapper *jetFinder)
+void AliAnalysisTaskFullppJet::AnalyzeSecondaryContribution(AliFJWrapper *jetFinder, const Int_t r, const Int_t etaCut)
+{
+ //
+ // Analyze secondaries
+ //
+
+ std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
+ for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
+ {
+ Float_t jetEta = jetsIncl[ij].eta();
+ Float_t jetPt = jetsIncl[ij].perp();
+ if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<1) continue;
+ std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
+ Double_t dNKPt = 0, dWPPt = 0, weakPt=0, nkPt=0;
+ Int_t type = -1;
+ for(UInt_t ic=0; ic<constituents.size(); ic++)
+ {
+ Int_t ipart = TMath::Abs(constituents[ic].user_index())-1;
+ AliMCParticle* mcParticle = (AliMCParticle*)fMC->GetTrack(ipart);
+ if(!mcParticle) continue;
+ Int_t pdg = mcParticle->PdgCode();
+ if(pdg==310 || (pdg<=3400 && pdg>=3100)) weakPt += mcParticle->Pt();
+ if(pdg==130 || pdg==2112) nkPt += mcParticle->Pt();
+
+ // Weak decaying particles
+ if(pdg==310) type=0;
+ else if (pdg==3122) type=1;
+ else if (pdg<=3400 && pdg>=3100) type=2;
+ else type=-1;
+ if(type>-1)
+ {
+ Int_t binx = fhSecondaryResponse[type]->GetXaxis()->FindFixBin(mcParticle->Pt());
+ TH1F *htmp = (TH1F*)fhSecondaryResponse[type]->ProjectionY(Form("hpro_%d_%d_%d",(Int_t)Entry(),ij,ic),binx,binx);
+ Double_t pro = htmp->GetRandom();
+ dWPPt += pro*mcParticle->Pt() - mcParticle->Pt();
+ delete htmp;
+ }
+
+ // Missing neutron and K0L
+ if(pdg==130 || pdg==2112) dNKPt -= mcParticle->Pt();
+ }
+
+ fhNKFracVsJetPt[etaCut][r]->Fill(jetPt,nkPt/jetPt);
+ fhWeakFracVsJetPt[etaCut][r]->Fill(jetPt,weakPt/jetPt);
+
+ Double_t jetNewWPPt = jetPt + dWPPt;
+ fhJetResponseWP[etaCut][r]->Fill(jetPt,jetNewWPPt);
+ fhJetResolutionWP[etaCut][r]->Fill(jetPt,(jetPt-jetNewWPPt)/jetPt);
+
+ Double_t jetNewNKPt = jetPt + dNKPt;
+ fhJetResponseNK[etaCut][r]->Fill(jetPt,jetNewNKPt);
+ fhJetResolutionNK[etaCut][r]->Fill(jetPt,(jetPt-jetNewNKPt)/jetPt);
+ }
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskFullppJet::AnalyzeSecondaryContributionViaMatching(AliFJWrapper *jetFinder, const Int_t r, const Int_t type, const Int_t etaCut)
+{
+ //
+ // Estimate the contribution of missing energies via matching
+ // type = 0 -> NK
+ // type = 1 -> WP
+ //
+
+ if(type!=0 && type!=1) return;
+
+ AliFJWrapper fj("fj","fj");
+ fj.CopySettingsFrom(*jetFinder);
+
+ Int_t npart = fMC->GetNumberOfTracks();
+ Int_t pType = -1;
+ for(Int_t ipart=0; ipart<npart; ipart++)
+ {
+ AliVParticle* vParticle = fMC->GetTrack(ipart);
+ if(!IsGoodMcPartilce(vParticle,ipart)) continue;
+ Int_t pdgCode = vParticle->PdgCode();
+ Double_t pt = vParticle->Pt();
+ if( type==0 && (pdgCode==130 || pdgCode==2112) ) continue; // Reject NK
+ if(type==1)
+ {
+ if(pdgCode==310) pType=0;
+ else if (pdgCode==3122) pType=1;
+ else if (pdgCode<=3400 && pdgCode>=3100) pType=2;
+ else pType=-1;
+ if(pType>-1)
+ {
+ Int_t binx = fhSecondaryResponse[pType]->GetXaxis()->FindFixBin(vParticle->Pt());
+ TH1F *htmp = (TH1F*)fhSecondaryResponse[pType]->ProjectionY(Form("hpro_%d_%d",(Int_t)Entry(),ipart),binx,binx);
+ Double_t pro = htmp->GetRandom();
+ pt = pro * vParticle->Pt();
+ delete htmp;
+ }
+ }
+ Int_t index = 1;
+ if(vParticle->Charge()==0) { index=-1; }
+ fj.AddInputVector(vParticle->Px()*pt/vParticle->Pt(), vParticle->Py()*pt/vParticle->Pt(), vParticle->Pz(), vParticle->P(), (ipart+1)*index);
+ }
+
+ fj.Run();
+ std::vector<fastjet::PseudoJet> jets_incl_mc = fj.GetInclusiveJets();
+ std::vector<fastjet::PseudoJet> jets_incl_mc_std = jetFinder->GetInclusiveJets();
+
+ for(UInt_t ij=0; ij<jets_incl_mc.size(); ij++)
+ {
+ Float_t jetEta = jets_incl_mc[ij].eta();
+ Float_t jetPt = jets_incl_mc[ij].perp();
+ if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<5) continue;
+ Double_t dEta, dPhi;
+ Int_t index = FindSpatialMatchedJet(jets_incl_mc[ij], jetFinder, dEta, dPhi, 1);
+ if(index>-1)
+ {
+ Double_t jetPtStd = jets_incl_mc_std[index].perp();
+ Double_t dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
+ if(type==0) fhJetResolutionNKSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR);
+ if(type==1) fhJetResolutionWPSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR);
+
+ if(dR<kdRCut[r])
+ {
+ if(type==0) fhJetResponseNKSM[etaCut][r]->Fill(jetPtStd,jetPt);
+ if(type==1) fhJetResponseWPSM[etaCut][r]->Fill(jetPtStd,jetPt);
+ }
+ }
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskFullppJet::RunAnalyzeUE(AliFJWrapper *jetFinder, const Int_t type, const Bool_t isMCTruth)
{
//
// Run analysis to estimate the underlying event
- // Adapt from Oliver
+ // Adapted from Oliver
std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
Double_t leadpt=0, leadPhi = -999, leadEta = -999;
+ Int_t leadIndex = -1;
for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
{
- Double_t jetPt = jetsIncl[ij].perp();
Double_t jetEta = jetsIncl[ij].eta();
+ Double_t jetPt = jetsIncl[ij].perp();
Double_t jetPhi = jetsIncl[ij].phi();
- if(TMath::Abs(jetEta)<0.5) continue;
if(leadpt<jetPt)
{
leadpt = jetPt;
leadPhi = jetPhi;
leadEta = jetEta;
+ leadIndex = ij;
}
}
+ if(leadpt<1 || TMath::Abs(leadEta)>0.5 ) return;
- Double_t ue = 0;
- Int_t ntracks = fTrackArray->GetEntriesFast();
- for(Int_t it=0; it<ntracks; it++)
+ Int_t backIndex = -1;
+ Bool_t isBackToBack = kFALSE;
+ for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
{
- AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
- if(!t) continue;
- Double_t axis = leadPhi + 0.5 * kPI;
- if(axis > 2*kPI) axis -= 2*kPI;
- Double_t dPhi = TMath::Abs(axis-t->Phi());
- Double_t dEta = TMath::Abs(leadEta-t->Eta());
- cout<<leadPhi*TMath::RadToDeg()<<" -> "<<t->Phi()*TMath::RadToDeg()<<endl;
- if(dPhi > 2*kPI) dPhi -= 2*kPI;
- if(TMath::Abs(dPhi*dPhi+dEta*dEta)<0.4)
+ Double_t dPhi = TMath::Abs(jetsIncl[ij].phi()-leadPhi);
+ if(dPhi > kPI) dPhi = 2*kPI - dPhi;
+ if(dPhi>150*TMath::DegToRad() && jetsIncl[ij].perp()/leadpt > 0.8)
{
- ue += t->Pt();
+ backIndex = ij;
+ isBackToBack = kTRUE;
}
}
- fhJetPtVsUE[fTriggerType]->Fill(ue,leadpt);
-}
+ for(Int_t ij=0; ij<(Int_t)jetsIncl.size(); ij++)
+ {
+ if(ij==leadIndex || ij==backIndex) continue;
+ if(TMath::Abs(jetsIncl[ij].eta())>0.5) continue;
+ if(jetsIncl[ij].perp()>15) isBackToBack = kFALSE;
+ }
+
+
+ Double_t axis[2] = {leadPhi+0.5 * kPI, leadPhi-0.5 * kPI};
+ if(axis[0]>2*kPI) axis[0] -= 2*kPI;
+ if(axis[1]<0) axis[1] += 2*kPI;
+
+ fhUEJetPtNorm[type][0][0]->Fill(leadpt);
+ if(isBackToBack) fhUEJetPtNorm[type][0][1]->Fill(leadpt);
+
+ if(!isMCTruth)
+ {
+ Int_t ntracks = fTrackArray->GetEntriesFast();
+ Double_t vertex[3] = {0, 0, 0};
+ fESD->GetVertex()->GetXYZ(vertex) ;
+ TLorentzVector gamma;
+ Int_t nclusters = fClusterArray->GetEntriesFast();
+
+ for(Int_t j=0; j<2; j++)
+ {
+ Double_t ueCh = 0, ueChNe = 0.;
+ for(Int_t it=0; it<ntracks; it++)
+ {
+ AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it);
+ if(!t) continue;
+ Double_t dPhi = TMath::Abs(axis[j]-t->Phi());
+ Double_t dEta = TMath::Abs(leadEta-t->Eta());
+ if(dPhi > kPI) dPhi = 2*kPI - dPhi;
+ if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
+ {
+ fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,t->Pt());
+ if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,t->Pt());
+ ueCh += t->Pt();
+ if(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4)
+ {
+ fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,t->Pt());
+ if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,t->Pt());
+ ueChNe += t->Pt();
+ }
+ }
+ }
+ fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh);
+ if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh);
+
+ if(!(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4)) continue;
+ fhUEJetPtNorm[type][1][0]->Fill(leadpt);
+ if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt);
+ for(Int_t ic=0; ic<nclusters; ic++)
+ {
+ AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic);
+ if(!cl) continue;
+ cl->GetMomentum(gamma, vertex);
+ Double_t clsPhi = gamma.Phi();
+ if(clsPhi<0) clsPhi += 2*kPI;
+ Double_t dPhi = TMath::Abs(axis[j]-clsPhi);
+ Double_t dEta = TMath::Abs(leadEta-gamma.Eta());
+ if(dPhi > kPI) dPhi = 2*kPI - dPhi;
+ if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
+ {
+ ueChNe += gamma.Pt();
+ fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,gamma.Pt());
+ if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,gamma.Pt());
+ }
+ }
+ fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe);
+ if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe);
+ }
+ }
+ else
+ {
+ fhUEJetPtNorm[type][1][0]->Fill(leadpt);
+ if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt);
+ Int_t npart = fMcPartArray->GetSize();
+ for(Int_t j=0; j<2; j++)
+ {
+ Double_t ueCh = 0, ueChNe = 0., ueChNe2 = 0.;
+ for(Int_t ipos=0; ipos<npart; ipos++)
+ {
+ AliVParticle* vParticle = fMC->GetTrack(fMcPartArray->At(ipos));
+ if(!vParticle) continue;
+ Double_t dPhi = TMath::Abs(axis[j]-vParticle->Phi());
+ Double_t dEta = TMath::Abs(leadEta-vParticle->Eta());
+ if(dPhi > kPI) dPhi = 2*kPI - dPhi;
+ if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4)
+ {
+ Double_t pt = vParticle->Pt();
+ ueChNe += pt;
+ fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,pt);
+ if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,pt);
+ if( vParticle->Charge()!=0 )
+ {
+ fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,pt);
+ if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,pt);
+ ueCh += pt;
+ ueChNe2 += pt;
+ }
+ else
+ {
+ if(pt>0.5) ueChNe2 += pt;
+ }
+ }
+ }
+ fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh);
+ fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe);
+ if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh);
+ if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe);
+ }
+ }
+}
//________________________________________________________________________
Bool_t AliAnalysisTaskFullppJet::IsGoodJet(AliAODJet *jet, Double_t rad)
//
Double_t z = 0;
- std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
+ std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
Int_t index = -1;
Double_t maxPt = 0;
//
Double_t jetSigma2 = 0;
- std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
+ std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
for(UInt_t ic=0; ic<constituents.size(); ic++)
{
if(constituents[ic].user_index()>0)
}
//________________________________________________________________________
-Int_t AliAnalysisTaskFullppJet::FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, const Double_t rad)
+Int_t AliAnalysisTaskFullppJet::FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR)
{
//
// Find spatially matched detector and particle jets
//
- Double_t maxR=0.25;
Int_t index=-1;
- std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
- for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
+ dEta=-999, dPhi=-999;
+ std::vector<fastjet::PseudoJet> jets_incl = jetFinder->GetInclusiveJets();
+ for(UInt_t ij=0; ij<jets_incl.size(); ij++)
{
- if(!IsGoodJet(jetsIncl[ij],rad)) continue;
- if(GetLeadingZ(ij, jetFinder)>0.98) continue;
- Double_t tmpR = TMath::Sqrt( TMath::Power(jet.eta()-jetsIncl[ij].eta(), 2) + TMath::Power(jet.phi()-jetsIncl[ij].phi(), 2) );
+ if(jets_incl[ij].perp()<5) continue;
+ if(TMath::Abs(jets_incl[ij].eta())>1) continue;
+ Double_t tmpR = TMath::Sqrt( TMath::Power(jet.eta()-jets_incl[ij].eta(), 2) + TMath::Power(jet.phi()-jets_incl[ij].phi(), 2) );
if(tmpR<maxR)
{
maxR=tmpR;
index=ij;
+ dEta = jet.eta()-jets_incl[ij].eta();
+ dPhi = jet.phi()-jets_incl[ij].phi();
}
}
return index;
Int_t matchedIndex = -1;
std::vector<fastjet::PseudoJet> jetsIncl1 = jetFinder1->GetInclusiveJets();
std::vector<fastjet::PseudoJet> jetsIncl2 = jetFinder2->GetInclusiveJets();
- std::vector<fastjet::PseudoJet> constituents1 = jetFinder1->GetJetConstituents(index1);
+ std::vector<fastjet::PseudoJet> constituents1 = jetFinder1->GetJetConstituents(index1);
Double_t jetPt1 = jetsIncl1[index1].perp();
if(jetPt1<0) return matchedIndex;
{
Double_t jetPt2 = jetsIncl2[ij2].perp();
if(jetPt2<0) return matchedIndex;
- std::vector<fastjet::PseudoJet> constituents2 = jetFinder2->GetJetConstituents(ij2);
+ std::vector<fastjet::PseudoJet> constituents2 = jetFinder2->GetJetConstituents(ij2);
Double_t sharedPt1 = 0., sharedPt2 = 0.;
for(UInt_t ic2=0; ic2<constituents2.size(); ic2++)
{
Double_t pt = newtrack->Pt();
Double_t eta = newtrack->Eta();
- if (pt < 0.15 || pt > fTrkPtMax[fTriggerType] || TMath::Abs(eta) > fTrkEtaMax)
+ if (pt < fTrkPtMin[fTriggerType] || pt > fTrkPtMax[fTriggerType] || TMath::Abs(eta) > fTrkEtaMax)
{
delete newtrack;
continue;
fRecoUtil->SetTracksMatchedToCluster(fESD);
}
- Double_t vertex[3] = {0, 0, 0};
- fESD->GetVertex()->GetXYZ(vertex) ;
- TLorentzVector gamma;
-
const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
for(Int_t i = 0 ; i < nCaloClusters; i++)
{
AliESDCaloCluster * cl = (AliESDCaloCluster *) fESD->GetCaloCluster(i);
if(!IsGoodCluster(cl)) continue;
AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cl);
+ if (newCluster->E() < fClsEtMin[fTriggerType] || newCluster->E() > fClsEtMax[fTriggerType])
+ {delete newCluster; continue;}
- // Trigger efficiency
- if(fSysJetTrigEff)
+ // Absolute scale
+ if(fPeriod.Contains("lhc11a",TString::kIgnoreCase))
{
- Double_t newE = newCluster->E() * (1+fVaryJetTrigEff);
+ Double_t newE = newCluster->E() * 1.02;
newCluster->SetE(newE);
- if(fTriggerType==0) fhClsE[fTriggerType]->Fill(newCluster->E());
- if(fTriggerType==1 && newCluster->Chi2()==1) fhClsE[fTriggerType]->Fill(newCluster->E());
}
-
- // non-linearity correction
- if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0)
- {
- Double_t correctedEnergy = fRecoUtil->CorrectClusterEnergyLinearity(newCluster);
- newCluster->SetE(correctedEnergy);
- }
- if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase))
+ // Trigger efficiency systematic uncertainty
+ if(fSysJetTrigEff)
{
- Double_t correctedEnergy = newCluster->E() * fMCNonLin->Eval(newCluster->E());
- newCluster->SetE(correctedEnergy);
+ Double_t newE = newCluster->E() * (1+fVaryJetTrigEff);
+ newCluster->SetE(newE);
+ if(fTriggerType==0) fhClsE[fTriggerType]->Fill(newCluster->E());
+ if(fTriggerType==1)
+ {
+ if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase)) fhClsE[0]->Fill(newCluster->E());
+ if(newCluster->Chi2()>0.5) fhClsE[fTriggerType]->Fill(newCluster->E());
+ }
}
- // Cluster energy scale
+ // Cluster energy scale systematic uncertainty
if(fSysClusterEScale)
{
Double_t newE = newCluster->E() * (1+fVaryClusterEScale);
newCluster->SetE(newE);
}
- // Cluster energy resolution
+ // Cluster energy resolution systematic uncertainty
if(fSysClusterERes)
+ {
+ Double_t oldE = newCluster->E();
+ Double_t resolution = fClusterEResolution->Eval(oldE);
+ Double_t smear = resolution * TMath::Sqrt((1+fVaryClusterERes)*(1+fVaryClusterERes)-1);
+ Double_t newE = oldE + fRandomGen->Gaus(0, smear);
+ newCluster->SetE(newE);
+ }
+
+ // non-linearity systematic uncertainty
+ if(fSysNonLinearity)
{
Double_t oldE = newCluster->E();
- Double_t resolution = fClusterEResolution->Eval(oldE);
- Double_t smear = resolution * fVaryClusterERes;
- Double_t newE = oldE + fRandomGen->Gaus(0, smear);
+ Double_t newE = oldE/fNonLinear->Eval(oldE) * 1.012;
newCluster->SetE(newE);
}
- newCluster->GetMomentum(gamma, vertex);
- if (gamma.Pt() < fClsEtMin[fTriggerType] || gamma.Pt() > fClsEtMax[fTriggerType]) {delete newCluster; continue;}
+ // clusterizer systematic uncertainty
+ if(fSysClusterizer)
+ {
+ fhSysClusterE[fTriggerType][0]->Fill(newCluster->E());
+ fhSysNCellVsClsE[fTriggerType][0]->Fill(newCluster->E(),newCluster->GetNCells());
+ }
+ // hadronic correction
Double_t clsE = newCluster->E();
Double_t subE = 0., eRes = 0., mcSubE = 0, MIPE = 0.;
if(fElectronRejection || fHadronicCorrection)
newCluster->SetEmcCpvDistance(MIPE);
newCluster->SetDistanceToBadChannel(mcSubE);
fClusterArray->Add(newCluster);
+
+ // clusterizer systematic uncertainty
+ if(fSysClusterizer)
+ {
+ fhSysClusterE[fTriggerType][1]->Fill(newCluster->E());
+ fhSysNCellVsClsE[fTriggerType][1]->Fill(newCluster->E(),newCluster->GetNCells());
+ }
}
if(fVerbosity>5) printf("[i] # of EMCal clusters in event: %d\n", fClusterArray->GetEntries());
Double_t trkP = track->P();
sumTrkPt += track->Pt(); sumTrkP += track->P();
if(fSysTrkPtRes) trkP = TMath::Sqrt(TMath::Power(track->Pt()+GetSmearedTrackPt(track),2)+TMath::Power(track->Pz(),2));
-
- if(IsElectron(track,cluster->E()))
+ if(IsElectron(track,cluster->E())) //electrons
{
if(fElectronRejection)
{
eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2);
}
}
- else
+ else //hadrons
{
if(fHadronicCorrection)
{
- MIPE += (trkP>0.27)?0.27:trkP;
+ MIPE += (trkP>0.27)?0.27:trkP; //MIP correction
if(fFractionHC>2)
{
if(trkP>0.27)
{
Double_t neutralE = 0;
TArrayI* labels = cluster->GetLabelsArray();
- //cout<<labels<<endl;
if(labels)
{
- //cout<<"# of MC particles: "<<labels->GetSize()<<endl;
for(Int_t il=0; il<labels->GetSize(); il++)
{
Int_t ipart = labels->At(il);
mcSubE = cluster->E() - neutralE;
if(mcSubE<0)
mcSubE=0;
- //printf("MIP = %f, subE = %f, clsE = %f, mcSubE = %f\n",MIPE, subE, cluster->E(), mcSubE);
}
{
//
// Exotic fraction: f_cross
- //
+ // Adpated from AliEMCalRecoUtils
if(!cluster) return -1;
AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
//
// Called once at the end of the query
//
- fhChunkQA->SetBinContent(fEDSFileCounter,fNTracksPerChunk);
+ Info("Terminate","Terminate");
+ AliAnalysisTaskSE::Terminate();
}
Int_t ncl = fESD->GetNumberOfCaloClusters();
for(Int_t icl=0; icl<ncl; icl++)
{
+ // Check every cluster
AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
if(!IsGoodCluster(cluster)) continue;
Double_t pro = GetOfflineTriggerProbability(cluster);
if(rand<pro)
{
isTrigger = 1;
+ fIsEventTriggerBit = 1;
cluster->SetChi2(1);
}
+ else
+ cluster->SetChi2(0);
}
return isTrigger;
}
AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells();
Int_t iSupMod = -1, absID = -1, ieta = -1, iphi = -1,iTower = -1, iIphi = -1, iIeta = -1;
Bool_t share = kFALSE;
- fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share);
+ fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share); // Get the position of the most energetic cell in the cluster
fGeom->GetCellIndex(absID,iSupMod,iTower,iIphi,iIeta);
fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
// convert co global phi eta
// get corresponding FALTRO
Int_t fphi = gphi / 2;
Int_t feta = geta / 2;
- if(fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5 && iSupMod>-1)
+ if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5 && iSupMod>-1) // check the trigger mask
{
Double_t clsE = cluster->E();
- if(fSysJetTrigEff) clsE = clsE * (1+fVaryJetTrigEff);
- if(clsE>5) pro = 1;
+ if(fSysClusterEScale) clsE = clsE * (1+fVaryClusterEScale); // Used for systematic uncertainty. Not needed for regular analysis
+ if(clsE>10) pro = 1; // Probability is 1 at high E
else
{
Int_t bin = fTriggerCurve[iSupMod]->FindFixBin(clsE);
- pro = fTriggerCurve[iSupMod]->GetBinContent(bin)/fTriggerEfficiency[iSupMod]->Eval(10);
+ pro = fTriggerCurve[iSupMod]->GetBinContent(bin)/fTriggerEfficiency[iSupMod]->Eval(10); // Read the probability from trigger turn-on curves
}
}
return pro;
// Smear track momentum
//
- Double_t resolution = track->Pt()*TMath::Sqrt(track->GetSigma1Pt2());
- Double_t smear = (resolution*fVaryTrkPtRes)*track->Pt();
+ Double_t resolution = track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2());
+ Double_t smear = resolution*TMath::Sqrt((1+fVaryTrkPtRes)*(1+fVaryTrkPtRes)-1);
return fRandomGen->Gaus(0, smear);
-}
-//________________________________________________________________________
-Double_t AliAnalysisTaskFullppJet::GetAdditionalSmearing(AliESDtrack *track)
-{
- //
- // Smear track momentum
- //
-
- Double_t resolution = fTrkPtResData->Eval(track->Pt())*track->Pt();
- return fRandomGen->Gaus(0, resolution);
}
fCaloTrigger->GetNL0Times( ntimes ); // get dimension of time arrays
if( ntimes < 1 ) continue; // no L0s in this channel. Presence of the channel in the iterator still does not guarantee that L0 was produced!!
fCaloTrigger->GetL0Times( trigtimes ); // get timing array
- if(fTriggerMask->GetBinContent(globCol+1,globRow+1)<0.5) continue;
+ if(fCheckTriggerMask && fTriggerMask->GetBinContent(globCol+1,globRow+1)<0.5) continue;
trigInCut = 0;
- //cout<<globCol<<" "<<globRow<<" "<<ntimes<<endl;
for( Int_t i = 0; i < ntimes; i++ )
{
- //printf("trigger times: %d\n",trigtimes[i]);
if( trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh ) // check if in cut
{
trigInCut = 1;
for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
{
AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
- if(!IsGoodCluster(cluster)) continue;
+ if(!cluster || !cluster->IsEMCAL()) continue;
cluster->SetChi2(0);
+ if(!IsGoodCluster(cluster)) continue;
+
//Clusters with most energetic cell in dead region can't be triggered
fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, nSupMod, ieta, iphi, share);
fGeom->GetCellIndex(absID,nSupMod,nModule, nIphi, nIeta);
fphi = gphi / 2;
feta = geta / 2;
- if(fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5)
+ if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5)
{
nCell = cluster->GetNCells(); // get cluster cells
cellAddrs = cluster->GetCellsAbsId(); // get the cell addresses
const char *triggerType[2] = {"MB","EMC"};
const char *decision[2] = {"no","yes"};
const char *recombination[] = {"E_scheme","pt_scheme","pt2_scheme","Et_scheme","Et2_scheme","BIpt_scheme","BIpt2_scheme"};
+ const char *type[2] = {"local","grid"};
if(fStudySubEInHC || fStudyMcOverSubE)
{
printf("\n\n=================================\n");
printf("====== NOT for PHYSICS! ======\n\n");
}
printf("Run period: %s\n",fPeriod.Data());
+ printf("Reject SPD pileup: %s\n",decision[fRejectPileup]);
+ printf("Reject exotic triggered events: %s\n",decision[fRejectExoticTrigger]);
printf("Is this MC data: %s\n",decision[fIsMC]);
printf("Only find charged jets in MC: %s\n",decision[fChargedMC]);
+ printf("Analyze on local or grid? %s\n",type[fAnaType]);
printf("Run offline trigger on MC: %s\n",decision[fOfflineTrigger]);
if(fIsMC)
printf("Is K0 and n included: %s\n",decision[1-fRejectNK]);
printf("Systematics: EMCal energy scale: %s with uncertainty of %1.0f%%\n",decision[fSysClusterEScale],fVaryClusterEScale*100);
printf("Systematics: EMCal energy resolution: %s with uncertainty of %1.0f%%\n",decision[fSysClusterERes],fVaryClusterERes*100);
printf("Smear lhc12a15a: %s\n",decision[fSmearMC]);
+ printf("Run UE analysis: %s\n",decision[fRunUE]);
+ printf("Run secondaries: %s\n",decision[fRunSecondaries]);
printf("=======================================\n\n");
}
//
// Check if the event containts exotic clusters
//
- if(!fESD)return;
+
Double_t leadingE = 0;
for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
{
if(leadingE>5) fIsExoticEvent5GeV = kTRUE;
}
+
//________________________________________________________________________
-Bool_t AliAnalysisTaskFullppJet::IsPrimaryVertexOk(const Double_t trueVtxZ) const
+Bool_t AliAnalysisTaskFullppJet::HasPrimaryVertex() const
{
//
- // Check if the event vertex is good
+ // Check if the primary vertex exists
//
-
const AliESDVertex* vtx = fESD->GetPrimaryVertex();
if (!vtx || vtx->GetNContributors()<1) return kFALSE;
- fhJetEventStat->Fill(5.5+fTriggerType);
-
- Double_t zVertex = vtx->GetZ();
- if(fEventZ[fTriggerType]) fEventZ[fTriggerType]->Fill(zVertex);
- if(fVertexGenZ[1]) fVertexGenZ[1]->Fill(trueVtxZ);
+ else return kTRUE;
+}
- if( TMath::Abs(zVertex) > fZVtxMax ) return kFALSE;
+//________________________________________________________________________
+Bool_t AliAnalysisTaskFullppJet::IsPrimaryVertexOk() const
+{
+ //
+ // Check if the event vertex is good
+ //
return kTRUE;
}
// Check if the event only has valid TPC vertex
//
- const AliESDVertex* vertex1 = fESD->GetPrimaryVertexTracks();
- const AliESDVertex* vertex2 = fESD->GetPrimaryVertexSPD();
- const AliESDVertex* vertex3 = fESD->GetPrimaryVertexTPC();
- if(vertex1 && vertex1->GetNContributors()==0 && vertex2 && vertex2->GetStatus()==0 && vertex3 && vertex3->GetNContributors()>0 )
+ const Bool_t isPrimaryVtx = HasPrimaryVertex();
+
+ Bool_t isGoodVtx = kTRUE;
+ AliESDVertex* goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexTracks());
+ if(!goodvtx || goodvtx->GetNContributors()<1)
+ goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexSPD()); // SPD vertex
+ if(!goodvtx || !goodvtx->GetStatus()) isGoodVtx = kFALSE; // Event rejected
+
+ if( isPrimaryVtx && !isGoodVtx )
return kTRUE;
else
return kFALSE;
}
+//________________________________________________________________________
+void AliAnalysisTaskFullppJet::CheckTPCOnlyVtx(const UInt_t trigger)
+{
+ //
+ // Check the fraction of accepted events that have only TPC vertex
+ //
+ Int_t lTriggerType = -1;
+ if (trigger & AliVEvent::kMB) lTriggerType = 1;
+ if (trigger & AliVEvent::kFastOnly) lTriggerType = 0;
+ if (trigger & AliVEvent::kEMC1) lTriggerType = 2;
+ if(lTriggerType==-1) return;
+ fhEventStatTPCVtx->Fill(0.5+lTriggerType*3);
+ if(HasPrimaryVertex()) fhEventStatTPCVtx->Fill(1.5+lTriggerType*3);
+ if(IsTPCOnlyVtx()) fhEventStatTPCVtx->Fill(2.5+lTriggerType*3);
+}
+
//________________________________________________________________________
Bool_t AliAnalysisTaskFullppJet::IsLEDEvent() const
{
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff