From 7c9152978e7d6068230a72d7a3ce0e080e7d3bd9 Mon Sep 17 00:00:00 2001 From: kleinb Date: Thu, 28 Feb 2013 15:29:15 +0000 Subject: [PATCH] Updates to pp Full jets, AddTask macro added(R. Ma) --- .../UserTasks/AliAnalysisTaskFullppJet.cxx | 1216 ++++++++++++----- .../UserTasks/AliAnalysisTaskFullppJet.h | 144 +- .../macros/AddTaskAliAnalysisTaskFullppJet.C | 63 + 3 files changed, 1047 insertions(+), 376 deletions(-) create mode 100644 PWGJE/EMCALJetTasks/macros/AddTaskAliAnalysisTaskFullppJet.C diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.cxx index 2a326e7ce7e..4e976d27660 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.cxx +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.cxx @@ -1,6 +1,7 @@ // ************************************** // 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 // ************************************** @@ -35,6 +36,7 @@ #include "AliEMCALRecoUtils.h" #include "TGeoManager.h" #include "AliMCEvent.h" +#include "AliMCParticle.h" #include "AliStack.h" #include "AliGenEventHeader.h" #include "AliGenPythiaEventHeader.h" @@ -51,33 +53,37 @@ using std::vector; 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 @@ -92,15 +98,21 @@ AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet() : 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; r2) 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; i2) 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; iDelete(); delete fJetTCA[s][a][r]; } @@ -312,13 +401,14 @@ AliAnalysisTaskFullppJet::~AliAnalysisTaskFullppJet() 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; rSetOwner(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"); @@ -381,24 +473,38 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() 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); @@ -413,15 +519,16 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() 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) { @@ -447,15 +554,22 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() fOutputList->Add(fEventZ[i]); } - for(Int_t r=0; rAdd(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); @@ -493,8 +607,17 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() 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) { @@ -524,16 +647,42 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() } } } - 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]); + } + } } @@ -541,8 +690,11 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() { if(fStudyMcOverSubE) { - for(Int_t r=0; rAdd(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; @@ -569,14 +760,18 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() { 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; rSetName(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); @@ -590,10 +785,13 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() { 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); @@ -609,7 +807,7 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() } 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++) @@ -638,28 +836,44 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() } } + 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 { @@ -680,12 +894,8 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() 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) { @@ -716,7 +926,6 @@ void AliAnalysisTaskFullppJet::UserCreateOutputObjects() PrintConfig(); BookHistos(); - //PostData(0, AODEvent()); PostData(1, fOutputList); } @@ -738,11 +947,6 @@ void AliAnalysisTaskFullppJet::UserExec(Option_t *) // 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) @@ -750,13 +954,13 @@ void AliAnalysisTaskFullppJet::UserExec(Option_t *) 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(InputEvent()); @@ -772,11 +976,19 @@ void AliAnalysisTaskFullppJet::UserExec(Option_t *) } 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 @@ -799,25 +1011,32 @@ void AliAnalysisTaskFullppJet::UserExec(Option_t *) } 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); @@ -826,52 +1045,68 @@ void AliAnalysisTaskFullppJet::UserExec(Option_t *) 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; rDelete(); - 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 @@ -879,33 +1114,34 @@ void AliAnalysisTaskFullppJet::UserExec(Option_t *) // 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; rGetEntriesFast(); for(Int_t it=0; itAt(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(); @@ -946,12 +1184,13 @@ void AliAnalysisTaskFullppJet::FindDetJets(const Int_t s, const Int_t a, const I Double_t py = pt*TMath::Sin(phi); if(fConstrainChInEMCal && ( TMath::Abs(eta)>0.7 || phi>kPI || phiAddInputVector(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; @@ -970,6 +1209,7 @@ void AliAnalysisTaskFullppJet::FindDetJets(const Int_t s, const Int_t a, const I } if(fVerbosity>5) printf("[i] # of clusters filled: %d\n",countClusters); } + // Run jet finding fDetJetFinder[s][a][r]->Run(); } @@ -980,32 +1220,35 @@ void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper // // 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 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; ij10) 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 constituents = jetFinder->GetJetConstituents(ij); + std::vector 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; ic10 && 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 ++; @@ -1019,7 +1262,7 @@ void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper 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 @@ -1045,7 +1288,10 @@ void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper 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); @@ -1054,21 +1300,40 @@ void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper 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 && ipartGetNumberOfTracks()) + { + 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="<ErrorEffectiveAreaNeutral()<10) cout<<"jet pt="<GetRefTracks()->GetEntries(), (Int_t)constituents.size(), jet->GetNEF()); // For catch good high-E jets @@ -1098,7 +1363,8 @@ void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper } } // End of catching good high-E jets - } + } + if(fVerbosity>5) printf("%s has %d jets\n",fJetArray->GetName(), fJetArray->GetEntries()); } @@ -1139,18 +1405,31 @@ void AliAnalysisTaskFullppJet::ProcessMC(const Int_t r) 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); + } + } } @@ -1180,13 +1459,17 @@ void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t for(UInt_t ij=0; ijFill(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 constituents = jetFinder->GetJetConstituents(ij); + std::vector constituents = jetFinder->GetJetConstituents(ij); Double_t neE=0, leadingZ = 0, maxPt = 0; Int_t constituentType = -1, leadingIndex = 0; for(UInt_t ic=0; ic 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; @@ -1250,20 +1533,24 @@ void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t 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; + } } } } @@ -1371,11 +1658,14 @@ void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t 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; ic0) { for(Int_t it=0; it<3; it++) @@ -1398,54 +1688,290 @@ void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t 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 jetsIncl = jetFinder->GetInclusiveJets(); + for(UInt_t ij=0; ij0.5*etaCut+0.5 || jetPt<1) continue; + std::vector constituents = jetFinder->GetJetConstituents(ij); + Double_t dNKPt = 0, dWPPt = 0, weakPt=0, nkPt=0; + Int_t type = -1; + for(UInt_t ic=0; icGetTrack(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; ipartGetTrack(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 jets_incl_mc = fj.GetInclusiveJets(); + std::vector jets_incl_mc_std = jetFinder->GetInclusiveJets(); + + for(UInt_t ij=0; ij0.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(dRFill(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 jetsIncl = jetFinder->GetInclusiveJets(); Double_t leadpt=0, leadPhi = -999, leadEta = -999; + Int_t leadIndex = -1; for(UInt_t ij=0; ij0.5 ) return; - Double_t ue = 0; - Int_t ntracks = fTrackArray->GetEntriesFast(); - for(Int_t it=0; itAt(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< "<Phi()*TMath::RadToDeg()< 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; itAt(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; icAt(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; iposGetTrack(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) @@ -1469,7 +1995,7 @@ Double_t AliAnalysisTaskFullppJet::GetLeadingZ(const Int_t jetIndex, AliFJWrappe // Double_t z = 0; - std::vector constituents = jetFinder->GetJetConstituents(jetIndex); + std::vector constituents = jetFinder->GetJetConstituents(jetIndex); std::vector jetsIncl = jetFinder->GetInclusiveJets(); Int_t index = -1; Double_t maxPt = 0; @@ -1504,7 +2030,7 @@ Double_t AliAnalysisTaskFullppJet::GetMeasuredJetPtResolution(const Int_t jetInd // Double_t jetSigma2 = 0; - std::vector constituents = jetFinder->GetJetConstituents(jetIndex); + std::vector constituents = jetFinder->GetJetConstituents(jetIndex); for(UInt_t ic=0; ic0) @@ -1575,24 +2101,26 @@ Bool_t AliAnalysisTaskFullppJet::IsGoodMcPartilce(const AliVParticle* vParticle, } //________________________________________________________________________ -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 jetsIncl = jetFinder->GetInclusiveJets(); - for(UInt_t ij=0; ij jets_incl = jetFinder->GetInclusiveJets(); + for(UInt_t ij=0; ij0.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 jetsIncl1 = jetFinder1->GetInclusiveJets(); std::vector jetsIncl2 = jetFinder2->GetInclusiveJets(); - std::vector constituents1 = jetFinder1->GetJetConstituents(index1); + std::vector constituents1 = jetFinder1->GetJetConstituents(index1); Double_t jetPt1 = jetsIncl1[index1].perp(); if(jetPt1<0) return matchedIndex; @@ -1616,7 +2144,7 @@ Int_t AliAnalysisTaskFullppJet::FindEnergyMatchedJet(AliFJWrapper *jetFinder1, c { Double_t jetPt2 = jetsIncl2[ij2].perp(); if(jetPt2<0) return matchedIndex; - std::vector constituents2 = jetFinder2->GetJetConstituents(ij2); + std::vector constituents2 = jetFinder2->GetJetConstituents(ij2); Double_t sharedPt1 = 0., sharedPt2 = 0.; for(UInt_t ic2=0; ic2Pt(); 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; @@ -1795,59 +2323,68 @@ void AliAnalysisTaskFullppJet::GetESDEMCalClusters() 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) @@ -1861,6 +2398,13 @@ void AliAnalysisTaskFullppJet::GetESDEMCalClusters() 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()); @@ -1921,8 +2465,7 @@ Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *clus 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) { @@ -1931,11 +2474,11 @@ Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *clus 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) @@ -1976,10 +2519,8 @@ Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *clus { Double_t neutralE = 0; TArrayI* labels = cluster->GetLabelsArray(); - //cout<GetSize()<GetSize(); il++) { Int_t ipart = labels->At(il); @@ -1996,7 +2537,6 @@ Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *clus mcSubE = cluster->E() - neutralE; if(mcSubE<0) mcSubE=0; - //printf("MIP = %f, subE = %f, clsE = %f, mcSubE = %f\n",MIPE, subE, cluster->E(), mcSubE); } @@ -2011,7 +2551,7 @@ Double_t AliAnalysisTaskFullppJet::GetExoticEnergyFraction(AliESDCaloCluster *cl { // // Exotic fraction: f_cross - // + // Adpated from AliEMCalRecoUtils if(!cluster) return -1; AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells(); @@ -2094,7 +2634,8 @@ void AliAnalysisTaskFullppJet::Terminate(Option_t *) // // Called once at the end of the query // - fhChunkQA->SetBinContent(fEDSFileCounter,fNTracksPerChunk); + Info("Terminate","Terminate"); + AliAnalysisTaskSE::Terminate(); } @@ -2112,6 +2653,7 @@ Int_t AliAnalysisTaskFullppJet::RunOfflineTrigger() Int_t ncl = fESD->GetNumberOfCaloClusters(); for(Int_t icl=0; iclGetCaloCluster(icl); if(!IsGoodCluster(cluster)) continue; Double_t pro = GetOfflineTriggerProbability(cluster); @@ -2119,8 +2661,11 @@ Int_t AliAnalysisTaskFullppJet::RunOfflineTrigger() if(randSetChi2(1); } + else + cluster->SetChi2(0); } return isTrigger; } @@ -2140,7 +2685,7 @@ Double_t AliAnalysisTaskFullppJet::GetOfflineTriggerProbability(AliESDCaloCluste 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 @@ -2149,15 +2694,15 @@ Double_t AliAnalysisTaskFullppJet::GetOfflineTriggerProbability(AliESDCaloCluste // 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; @@ -2226,20 +2771,10 @@ Double_t AliAnalysisTaskFullppJet::GetSmearedTrackPt(AliESDtrack *track) // 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); } @@ -2285,12 +2820,10 @@ void AliAnalysisTaskFullppJet::CheckEventTriggerBit() 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< fTrigCutLow && trigtimes[i] < fTrigCutHigh ) // check if in cut { trigInCut = 1; @@ -2311,8 +2844,10 @@ void AliAnalysisTaskFullppJet::CheckEventTriggerBit() for(Int_t icl=0; iclGetNumberOfCaloClusters(); 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); @@ -2322,7 +2857,7 @@ void AliAnalysisTaskFullppJet::CheckEventTriggerBit() 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 @@ -2363,6 +2898,7 @@ void AliAnalysisTaskFullppJet::PrintConfig() 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"); @@ -2370,8 +2906,11 @@ void AliAnalysisTaskFullppJet::PrintConfig() 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]); @@ -2407,6 +2946,8 @@ void AliAnalysisTaskFullppJet::PrintConfig() 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"); } @@ -2416,7 +2957,7 @@ void AliAnalysisTaskFullppJet::CheckExoticEvent() // // Check if the event containts exotic clusters // - if(!fESD)return; + Double_t leadingE = 0; for(Int_t icl=0; iclGetNumberOfCaloClusters(); icl++) { @@ -2432,23 +2973,25 @@ void AliAnalysisTaskFullppJet::CheckExoticEvent() 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; } @@ -2459,15 +3002,36 @@ Bool_t AliAnalysisTaskFullppJet::IsTPCOnlyVtx() const // 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(fESD->GetPrimaryVertexTracks()); + if(!goodvtx || goodvtx->GetNContributors()<1) + goodvtx = const_cast(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 { diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.h index d09a2622328..9f7c1af5862 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.h +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.h @@ -25,19 +25,19 @@ class AliESDEvent; class AliMCEvent; class AliStack; class AliESDtrack; -class AliESDtrackCuts; +//class AliESDtrackCuts; class AliEMCALGeometry; class AliEMCALRecoUtils; class AliESDCaloCluster; class AliFJWrapper; class AliAODJet; +class TParticle; +#include "AliESDtrackCuts.h" #include "AliAnalysisTaskSE.h" #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" -#include "AliAnalysisTaskSE.h" - class AliAnalysisTaskFullppJet : public AliAnalysisTaskSE { public: AliAnalysisTaskFullppJet(); @@ -53,11 +53,22 @@ class AliAnalysisTaskFullppJet : public AliAnalysisTaskSE { enum {kTPCOnlyVtx = 1<<3, kTrigger = 1<<4, kHighZ = 1<<5, - kSuspicious = 1<<6}; + kSuspicious = 1<<6, + kGluon = 1<<9, + kQuark = 1<<10, + kLeadCh = 1<<11 }; + enum { kHybrid=0, + kTPCOnly=1, + kGlobal=2}; + + void SetAnaType(Int_t ana) { fAnaType=ana; } void SetMCAna(Bool_t mc) { fIsMC=mc; } - void SetMCStandalone(Bool_t mc) { fMCStandalone=mc; } + void SetPhySelForMC(Bool_t phy) { fPhySelForMC=phy; } void SetChargedMC(Bool_t mc) { fChargedMC=mc; } + void SetRejectSPDPileup(Bool_t re) { fRejectPileup=re; } + void SetRejectExoticTrigger(Bool_t re) { fRejectExoticTrigger=re; } + void SetCheckTriggerMask(Bool_t check) { fCheckTriggerMask=check; } void SetRunPeriod(char *p) { fPeriod=p; } void SetOfflineTrigger(Bool_t t) { fOfflineTrigger=t; } void SetXsec(Float_t Xsec) { fXsecScale=Xsec; } @@ -96,6 +107,8 @@ class AliAnalysisTaskFullppJet : public AliAnalysisTaskSE { void SetTrkEffCorrCutZ(Double_t zcut) { fTrkEffCorrCutZ=zcut; } void SetSmearMC(Double_t smear) { fSmearMC=smear; } void SetRunUE(Bool_t run) { fRunUE=run; } + void SetCheckTPCOnlyVtx(Bool_t check) { fCheckTPCOnlyVtx=check; } + void SetRunSecondaries(Bool_t run) { fRunSecondaries=run; } //-------------------------------- // Kinematic cut @@ -129,6 +142,7 @@ class AliAnalysisTaskFullppJet : public AliAnalysisTaskSE { void SetVaryClusterEScale(Double_t vary) { fVaryClusterEScale=vary; } void SetSysClusterERes(Bool_t sys) { fSysClusterERes=sys; } void SetVaryClusterERes(Double_t vary) { fVaryClusterERes=vary; } + void SetSysClusterizer(Bool_t sys) { fSysClusterizer=sys; } protected: @@ -139,9 +153,10 @@ protected: void ProcessMC(const Int_t r=0); void GetMCInfo(); Bool_t IsGoodMcPartilce(const AliVParticle* vParticle, const Int_t ipart); - Int_t FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, const Double_t radius); + Int_t FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR); Int_t FindEnergyMatchedJet(AliFJWrapper *jetFinder1, const Int_t index1, AliFJWrapper *jetFinder2, const Double_t fraction=0.5); - Bool_t IsPrimaryVertexOk(const Double_t trueVtxZ) const; + Bool_t HasPrimaryVertex() const; + Bool_t IsPrimaryVertexOk() const; Bool_t IsTPCOnlyVtx() const; Bool_t IsLEDEvent() const; void CheckExoticEvent(); @@ -156,7 +171,11 @@ protected: void FindDetJets(const Int_t s=0, const Int_t a=0, const Int_t r=0); void FillAODJets(TClonesArray *fJetArray, AliFJWrapper *jetFinder, const Bool_t isTruth = 0); void AnalyzeJets(AliFJWrapper *jetFinder, const Int_t type, const Int_t r); - void RunAnalyzeUE(AliFJWrapper *jetFinder); + void RunAnalyzeUE(AliFJWrapper *jetFinder, const Int_t type, const Bool_t isMCTruth); + void AnalyzeSecondaryContribution(AliFJWrapper *jetFinder, const Int_t r, const Int_t etaCut); + void AnalyzeSecondaryContributionViaMatching(AliFJWrapper *jetFinder, const Int_t r, const Int_t type, const Int_t etaCut); + void GetSecondaryPtFraction(TParticle *particle, Double_t &chPt, Double_t &reGenPt, Double_t &reRecPt); + void CheckTPCOnlyVtx(const UInt_t trigger); Bool_t IsGoodJet(fastjet::PseudoJet jet, Double_t radius); Bool_t IsGoodJet(AliAODJet *jet, Double_t radius); Double_t GetLeadingZ(const Int_t jetIndex, AliFJWrapper *jetFinder); @@ -166,7 +185,6 @@ protected: Double_t GetJetMissPtDueToTrackingEfficiency(const Int_t jetIndex, AliFJWrapper *jetFinder, const Int_t radiusIndex); Double_t GetExoticEnergyFraction(AliESDCaloCluster *cluster); Double_t GetSmearedTrackPt(AliESDtrack *track); - Double_t GetAdditionalSmearing(AliESDtrack *track); enum { kNBins = 20 }; @@ -174,6 +192,9 @@ protected: Int_t fVerbosity; // Control output Int_t fEDSFileCounter; // Keep track of the ESD file inside a chain Int_t fNTracksPerChunk; // Number of tracks per ESD file; used for debugging + Bool_t fRejectPileup; // Flag to reject SPD pileup events + Bool_t fRejectExoticTrigger; // Flag to reject events triggered by exotic clusters + Int_t fAnaType; // 0-local, 1-grid TString fPeriod; // Data period AliESDEvent *fESD; //! ESD object AliAODEvent *fAOD; //! AOD object @@ -182,13 +203,14 @@ protected: TObjArray *fTrackArray; //! Array of input tracks TObjArray *fClusterArray; //! Array of input clusters TArrayI *fMcPartArray; //! Array of MC particles - Bool_t fIsMC; // Flag if analyzing MC data - Bool_t fMCStandalone; // Flag if only analyze Particle-Level MC - Bool_t fChargedMC; // Flag if finding only charged MC jets - Float_t fXsecScale; // Corss section - Double_t fCentrality; //! V0M for current event + Bool_t fIsMC; // Flag of MC data + Bool_t fPhySelForMC; // Flag to run physics selection in case of MC data + Bool_t fChargedMC; // Flag to find only charged MC jets + Float_t fXsecScale; // Corss section for each pT hard bin + Double_t fCentrality; // V0M for current event Double_t fZVtxMax; // Max vertex z cut Int_t fTriggerType; // 0-MB, 1-EMC + Bool_t fCheckTriggerMask; // Flag to check the trigger mask for triggered events Bool_t fIsTPCOnlyVtx; // Flag of events with ONLY TPC vertex Bool_t fIsExoticEvent3GeV; // Flag of events with exotic cluster above 3 GeV Bool_t fIsExoticEvent5GeV; // Flag of events with exotic cluster above 5 GeV @@ -199,9 +221,9 @@ protected: TF1 *fTriggerEfficiency[10]; //! Fit of trigger turn-on curves for EMCal clusters above 4-5 GeV AliEMCALGeometry *fGeom; //! EMCal goemetry utility AliEMCALRecoUtils *fRecoUtil; //! Reco utility - AliESDtrackCuts *fEsdTrackCuts; //! Track cuts for good tracks - AliESDtrackCuts *fHybridTrackCuts1; //! Track cuts for tracks without SPD hit - AliESDtrackCuts *fHybridTrackCuts2; //! Track cuts for tracks witout SPD hit or ITS refit + AliESDtrackCuts *fEsdTrackCuts; // Track cuts for good tracks + AliESDtrackCuts *fHybridTrackCuts1; // Track cuts for tracks without SPD hit + AliESDtrackCuts *fHybridTrackCuts2; // Track cuts for tracks witout SPD hit or ITS refit Int_t fTrackCutsType; // 0-Global track, 1-TPCOnly track Int_t fKinCutType; // 0-cut on track before jet finding, 1-cut on jet with high-pt tracks Double_t fTrkEtaMax; // Max |eta| cut @@ -224,7 +246,6 @@ protected: Float_t fFractionHC; // fraction of hadronic correction Double_t fHCLowerPtCutMIP; // Lower track pt cut for MIP correction TF1 *fClusterEResolution; //! Parameterization of cluster energy resolution from test beam results - TF1 *fMCNonLin; //! Non-linearity of simualtion Double_t fJetNEFMin; // Min jet NEF cut Double_t fJetNEFMax; // Max jet NEF cut Bool_t fSpotGoodJet; // Good jet catching @@ -237,6 +258,9 @@ protected: TH1F *fhCorrTrkEffSample[2][2][kNBins]; //! Tracking efficiency estimated from simulation TRandom3 *fRandomGen; //! Random number generator Bool_t fRunUE; // Run analysis of underlying event + Bool_t fCheckTPCOnlyVtx; // Check events with TPC only vertices + Bool_t fRunSecondaries; // Run analysise for secondary particles + TH2F *fhSecondaryResponse[3]; //! Response matrix for secondary particles Bool_t fSysJetTrigEff; // Flag of systematic uncertainty of jet trigger efficiency Double_t fVaryJetTrigEff; // Variation of cluster E-scale for systematic uncertainty of jet trigger efficiency @@ -248,64 +272,84 @@ protected: Double_t fCutdEta; // Variation of dEta cut Double_t fCutdPhi; // Variation of dPhi cut Bool_t fSysNonLinearity; // Flag of systematic uncertainty of EMCal non-linearity + TF1 *fNonLinear; //! Non-linearity correction functions for data Bool_t fSysClusterEScale; // Flag of systematic uncertainty of EMCal energy scale Double_t fVaryClusterEScale; // Variation of EMCal energy scale Bool_t fSysClusterERes; // Flag of systematic uncertainty of EMCal energy resolution Double_t fVaryClusterERes; // Variation of EMCal energy resolution + Bool_t fSysClusterizer; // Flag of systematic uncertainty on clusterizer - TString fNonStdBranch; //! Non-std branch name for AOD jets - TString fNonStdFile; //! Name of optional file that the non-std branch is written to - TString fAlgorithm; //! name of algorithm - TString fRadius; //! Jet cone radius + TString fNonStdBranch; // Non-std branch name for AOD jets + TString fNonStdFile; // Name of optional file that the non-std branch is written to + TString fAlgorithm; // name of algorithm + TString fRadius; // Jet cone radius Int_t fRecombinationScheme; // Recombination scheme of jet finding - AliFJWrapper *fDetJetFinder[3][2][2]; //! Jet finder - TClonesArray *fJetTCA[3][2][2]; //! TCA of jets: in - akt - 0.4 + AliFJWrapper *fDetJetFinder[3][2][3]; //! Jet finder + TClonesArray *fJetTCA[3][2][3]; //! TCA of jets: in - akt - r Bool_t fConstrainChInEMCal; // Constain charged particle to be in EMCal acceptance Bool_t fRejectNK; // Reject neutron and K_L Bool_t fRejectWD; // Reject primaries, mainly k^0_S, that decay weakly Bool_t fSmearMC; // Flag of smearing tracking resolution in MC to match data. Obselete. TF1 *fTrkPtResData; //! Parameterazation of momentum resolution estimated from data - AliFJWrapper *fTrueJetFinder[2]; //! Jet finder for particle jets - TClonesArray *fMcTruthAntikt[2]; //! TCA of MC truth anti-kt jets + AliFJWrapper *fTrueJetFinder[3]; //! Jet finder for particle jets + TClonesArray *fMcTruthAntikt[3]; //! TCA of MC truth anti-kt jets TList *fOutputList; //! Output list Bool_t fSaveQAHistos; // Flag of saving QA histograms + TH1F *fhJetEventCount; //! Event counts to keep track of rejection criteria TH1F *fhJetEventStat; //! Event counts used for jet analysis - TH1F *fhEventCheck; //! Event statistics for vertex types - TH2F *fhTrkPhiEtaCheck; //! Phi vs Eta of tracks in events with only TPC vertex + TH1F *fhEventStatTPCVtx; //! Event counts for TPC only vertices TH1F *fhChunkQA; //! Check if the chunk is corrupted TH1F *fVertexGenZ[2]; //! Generated event vertex z TH1F *fEventZ[2]; //! reconstructed event vertex z TH1F *fhNTrials[2]; //! # of trials TH1F *fhNMatchedTrack[2]; //! # of matched tracks per cluster TH2F *fhSubEVsTrkPt[2][4]; //! Subtracted energy due to hadronic correction - TH2F *fhNeutralPtInJet[3][2]; //! pt of neutral constituents in jet - TH2F *fhChargedPtInJet[3][2]; //! pt of charged constituents in jet - TH2F *fhLeadNePtInJet[3][2]; //! pt of leading neutral constituents in jet - TH2F *fhLeadChPtInJet[3][2]; //! pt of leading charged constituents in jet - TH2F *fhChLeadZVsJetPt[2][2]; //! Leading charged constituent Z vs jet pt - TH3F *fhJetPtVsZ[3][2]; //! Jet pt vs constituent Z vs constituent type - TH3F *fRelTrkCon[2][2]; //! Jet pt vs track pt contribution vs track class - TH2F *fhJetPtWithTrkThres[2][2]; //! pt of jets containing tracks above certian threshold - TH2F *fhJetPtWithClsThres[2][2]; //! pt of jets containing clusters above certian threshold - TH2F *fhJetPtVsLowPtCons[2][2]; //! Contribution of low pt particles to jet energy - THnSparse *fJetEnergyFraction[3][2]; //! Jet energy fraction - THnSparse *fJetNPartFraction[3][2]; //! Jet NPart fraction - TH1F *fJetCount[3][2]; //! pT distribution of pions detected - TH2F *fhSubClsEVsJetPt[2][2][5]; //! f*subtracted cluster energy vs jet pt - TH2F *fhHCTrkPtClean[2][2][5]; //! Cleanly subtracted charged pt - TH2F *fhHCTrkPtAmbig[2][2][5]; //! Ambiguously subtracted charged pt - TH2F *fHCOverSubE[2][5]; //! Error made by hadronic correction assessed by using particle jet - TH2F *fHCOverSubEFrac[2][5]; //! Error made by hadronic correction assessed by using particle jet - TH3F *fhFcrossVsZleading[2][2]; //! Jet pt vs Fcross vs Zleading - TH1F *fhJetPtInExoticEvent[2][2]; //! Jet pt in exotic events - TH2F *fhJetPtVsUE[2]; //! Underlying event contribution + TH2F *fhNeutralPtInJet[3][3]; //! pt of neutral constituents in jet + TH2F *fhTrigNeuPtInJet[3][3]; //! pt of neutral constituents in jet + TH2F *fhChargedPtInJet[3][3]; //! pt of charged constituents in jet + TH2F *fhLeadNePtInJet[3][3]; //! pt of leading neutral constituents in jet + TH2F *fhLeadChPtInJet[3][3]; //! pt of leading charged constituents in jet + TH2F *fhChLeadZVsJetPt[2][3]; //! Leading charged constituent Z vs jet pt + TH3F *fhJetPtVsZ[3][3]; //! Jet pt vs constituent Z vs constituent type + TH3F *fRelTrkCon[2][3]; //! Jet pt vs track pt contribution vs track class + TH2F *fhJetPtWithTrkThres[2][3]; //! pt of jets containing tracks above certian threshold + TH2F *fhJetPtWithClsThres[2][3]; //! pt of jets containing clusters above certian threshold + TH2F *fhJetPtVsLowPtCons[2][3][2]; //! Contribution of low pt particles to jet energy + THnSparse *fJetEnergyFraction[3][3]; //! Jet energy fraction + THnSparse *fJetNPartFraction[3][3]; //! Jet NPart fraction + TH1F *fJetCount[3][3]; //! pT distribution of pions detected + TH2F *fhSubClsEVsJetPt[2][3][5]; //! f*subtracted cluster energy vs jet pt + TH2F *fhHCTrkPtClean[2][3][5]; //! Cleanly subtracted charged pt + TH2F *fhHCTrkPtAmbig[2][3][5]; //! Ambiguously subtracted charged pt + TH2F *fHCOverSubE[3][5]; //! Error made by hadronic correction assessed by using particle jet + TH2F *fHCOverSubEFrac[3][5]; //! Error made by hadronic correction assessed by using particle jet + TH3F *fhFcrossVsZleading[2][3]; //! Jet pt vs Fcross vs Zleading + TH1F *fhJetPtInExoticEvent[2][3]; //! Jet pt in exotic events + TH2F *fhUEJetPtVsSumPt[3][2][2]; //! Leading jet pt vs underlying event pt + TH2F *fhUEJetPtVsConsPt[3][2][2]; //! Leading jet pt vs constituent pt in underlying event + TH1F *fhUEJetPtNorm[3][2][2]; //! Leading jet normalization TH1F *fhClsE[2]; //! Cluster energy distribution + TH3F *fhJetInTPCOnlyVtx[2][3]; //! Jets in full TPC acceptance in events with TPC only vertex + TH1F *fhSysClusterE[2][2]; //! Cluster energy distribution before and after hadonic correction + TH2F *fhSysNCellVsClsE[2][2]; //! NCell vs cluster energy before and after hadonic correction + + // Secondaries + TH2F *fhNKFracVsJetPt[2][3]; //! Energy fraction lost due to missing neutron and K0L + TH2F *fhWeakFracVsJetPt[2][3]; //! Energy fraction lost due to weakly decaying particles + TH2F *fhJetResponseNK[2][3]; //! Jet response due to missing neutron and K0L using response matrix + TH2F *fhJetResponseWP[2][3]; //! Jet response due to missing weakly decayed particles using response matrix + TH2F *fhJetResolutionNK[2][3]; //! Jet resolution due to missing neutron and K0L using response matrix + TH2F *fhJetResolutionWP[2][3]; //! Jet resolution due to missing weakly decayed particles using response matrix + TH2F *fhJetResponseNKSM[2][3]; //! Jet response due to missing neutron and K0L via matching + TH2F *fhJetResponseWPSM[2][3]; //! Jet response due to missing weakly decayed particles via matching + TH3F *fhJetResolutionNKSM[2][3]; //! Jet resolution due to missing neutron and K0L via matching + TH3F *fhJetResolutionWPSM[2][3]; //! Jet resolution due to missing weakly decayed particles via matching AliAnalysisTaskFullppJet(const AliAnalysisTaskFullppJet&); // not implemented AliAnalysisTaskFullppJet &operator=(const AliAnalysisTaskFullppJet&); // not implemented - ClassDef(AliAnalysisTaskFullppJet, 1); + ClassDef(AliAnalysisTaskFullppJet, 2); }; #endif diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskAliAnalysisTaskFullppJet.C b/PWGJE/EMCALJetTasks/macros/AddTaskAliAnalysisTaskFullppJet.C new file mode 100644 index 00000000000..57dab81515b --- /dev/null +++ b/PWGJE/EMCALJetTasks/macros/AddTaskAliAnalysisTaskFullppJet.C @@ -0,0 +1,63 @@ +AliAnalysisTaskFullppJet *AddTaskAliAnalysisTaskFullppJet(const char *name = "Baseline", + const char *period = "lhc11a", + const Bool_t isMC = kFALSE, + const Bool_t IsPhySelForMC = kFALSE, + const Bool_t offlineTrig = kFALSE, + const Double_t minTrkPt = 0.15, + const Double_t minClsEt = 0.30, + const Bool_t hc = kTRUE, + const Double_t fraction = 1) +{ + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if(!mgr) + { + AliError("No analysise manager is availabe !"); + return NULL; + } + + gROOT->LoadMacro("$ALICE_ROOT/PWGJE/macros/CreateTrackCutsPWGJE.C"); + AliESDtrackCuts *esdTrackCuts = 0x0; + AliESDtrackCuts *hybridTrackCuts1 = 0x0; + AliESDtrackCuts *hybridTrackCuts2 = 0x0; + printf("\n===== Use Hybrid track cuts =====\n"); + esdTrackCuts = CreateTrackCutsPWGJE(10001006); + hybridTrackCuts1 = CreateTrackCutsPWGJE(1006); + hybridTrackCuts2 = CreateTrackCutsPWGJE(10041006); + + AliAnalysisTaskFullppJet *jetTask = new AliAnalysisTaskFullppJet(Form("ppJet_%s_%s",period,name)); + jetTask->SetNonStdBranch(name); + jetTask->SetRunPeriod(period); + jetTask->SetAnaType(1); + jetTask->SetCheckTriggerMask(kFALSE); + jetTask->SetMCAna(isMC); + jetTask->SetPhySelForMC(IsPhySelForMC); + jetTask->SetRejectSPDPileup(kTRUE); + jetTask->SetZvtx(10); + jetTask->SetTrackCutsType(AliAnalysisTaskFullppJet::kHybrid); + jetTask->SetEsdTrackCuts(esdTrackCuts); + jetTask->SetHybridTrackCuts1(hybridTrackCuts1); + jetTask->SetHybridTrackCuts2(hybridTrackCuts2); + jetTask->SetOfflineTrigger(offlineTrig); + jetTask->SetEtaMax(1); + jetTask->SetdEdxRange(75,95); + jetTask->SetEoverPRange(0.8,1.2); + jetTask->SetPtRange(minTrkPt,1e4,minTrkPt,1e4); + jetTask->SetEtRange(minClsEt,1e4,minClsEt,1e4); + jetTask->SetRejectExoticCluster(kTRUE); + jetTask->SetRemoveProblematicSM4(kTRUE); + jetTask->SetStudySubEInHC(kFALSE); + jetTask->SetJetNEFCut(0.02,0.98); + jetTask->SetRejectElectron(hc); + jetTask->SetCorrectHadron(hc); + jetTask->SetHCFraction(fraction); + jetTask->SetRadius("0.4 0.2 0.3"); + jetTask->SetCheckTrkEffCorr(kFALSE); + + mgr->AddTask(jetTask); + TString outfileName = "ppJetOutput.root"; + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(Form("JetOutputList_%s",name), TList::Class(), AliAnalysisManager::kOutputContainer,outfileName.Data()); + mgr->ConnectInput(jetTask,0,mgr->GetCommonInputContainer()); + mgr->ConnectOutput(jetTask, 1, coutput1); + + return jetTask; +} -- 2.43.0