#include "AliAnalysisTaskSE.h" #include "AliAnalysisHelperJetTasks.h" #include "THnSparse.h" #include "AliAnalysisManager.h" #include "AliInputEventHandler.h" #include "AliAODHandler.h" #include "AliAODInputHandler.h" #include "AliAODEvent.h" #include "AliESDEvent.h" #include "TAxis.h" #include "AliAODMCParticle.h" #include "AliAnalysisTaskPPJetSpectra.h" #include "AliAnalysisHelperJetTasks.h" ClassImp(AliAnalysisTaskPPJetSpectra) //--------------------------------------------------------------------------------------------------- AliAnalysisTaskPPJetSpectra::AliAnalysisTaskPPJetSpectra() :AliAnalysisTaskSE() ,fOutputList(0) ,fESD(0) ,fAOD(0) ,fAODIn(0) ,fAODOut(0) ,fAODExt(0) ,fNonStdFile("") ,fDebug(0) ,fUseMC(kFALSE) ,fEvtSelectionMask(0) ,fEventClass(-1) ,nVtxContCut(2) ,fVtxZcut(10) ,fVtxRcut(1) ,nTrackFilter(272) ,trackPtMin(0.15) ,trackPtMax(100.) ,trackEtaAbsMax(0.9) ,jetPtMin(0.) ,jetEtaCut(0.) ,jetZmax(0.) ,fhnEvent(0) ,fhnTracks(0) ,fhnMC(0) ,fhnMC2(0) ,fhnRecJetsNoCut(0) ,fhnGenJetsNoCut(0) ,fhnRecJetsCut(0) ,fhnGenJetsCut(0) ,fhnRecBckg(0) ,fhnGenBckg(0) ,fhnRecJetsTrackUEcor(0) ,fhnGenJetsTrackUEcor(0) ,fhnRecJetsBckgUEcor(0) ,fhnGenJetsBckgUEcor(0) ,fhnTrackUE(0) ,fhnParticleUE(0) ,fhnBckgRecUE(0) ,fhnBckgGenUE(0) ,fRecJetBranch("") ,fGenJetBranch("") ,fRecBckgBranch("") ,fGenBckgBranch("") ,fRecJetR(0.) ,fGenJetR(0.) ,fRecBckgR(0.) ,fGenBckgR(0.) ,fTrackType(0) ,fParticleType(0) ,fhnMatching(0) ,fhnTrackCorrMatching(0) ,fhnBckgCorrMatching(0) ,fhnTrackUEanal(0) ,kDoUEanalysis(kFALSE) ,fRejectPileUp(1) { if(fDebug) printf("%s: %d Constructor\n",(char*)__FILE__, __LINE__); //DefineOutput(1, TList::Class()); } //--------------------------------------------------------------------------------------------------- AliAnalysisTaskPPJetSpectra::AliAnalysisTaskPPJetSpectra(const char* name):AliAnalysisTaskSE(name) ,fOutputList(0) ,fESD(0) ,fAOD(0) ,fAODIn(0) ,fAODOut(0) ,fAODExt(0) ,fNonStdFile("") ,fDebug(0) ,fUseMC(kFALSE) ,fEvtSelectionMask(0) ,fEventClass(-1) ,nVtxContCut(2) ,fVtxZcut(10) ,fVtxRcut(1) ,nTrackFilter(272) ,trackPtMin(0.15) ,trackPtMax(100.) ,trackEtaAbsMax(0.9) ,jetPtMin(0.) ,jetEtaCut(0.) ,jetZmax(0.) ,fhnEvent(0) ,fhnTracks(0) ,fhnMC(0) ,fhnMC2(0) ,fhnRecJetsNoCut(0) ,fhnGenJetsNoCut(0) ,fhnRecJetsCut(0) ,fhnGenJetsCut(0) ,fhnRecBckg(0) ,fhnGenBckg(0) ,fhnRecJetsTrackUEcor(0) ,fhnGenJetsTrackUEcor(0) ,fhnRecJetsBckgUEcor(0) ,fhnGenJetsBckgUEcor(0) ,fhnTrackUE(0) ,fhnParticleUE(0) ,fhnBckgRecUE(0) ,fhnBckgGenUE(0) ,fRecJetBranch("") ,fGenJetBranch("") ,fRecBckgBranch("") ,fGenBckgBranch("") ,fRecJetR(0.) ,fGenJetR(0.) ,fRecBckgR(0.) ,fGenBckgR(0.) ,fTrackType(0) ,fParticleType(0) ,fhnMatching(0) ,fhnTrackCorrMatching(0) ,fhnBckgCorrMatching(0) ,fhnTrackUEanal(0) ,kDoUEanalysis(kFALSE) ,fRejectPileUp(1) { if(fDebug) printf("%s: %d Constructor\n",(char*)__FILE__, __LINE__); DefineOutput(1, TList::Class()); } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::Init() { if(fDebug) printf("%s: %d Initialize\n", (char*)__FILE__, __LINE__); return; } //--------------------------------------------------------------------------------------------------- Bool_t AliAnalysisTaskPPJetSpectra::Notify() { if(fDebug) printf("%s: %d Notify\n", (char*)__FILE__, __LINE__); return kTRUE; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::UserCreateOutputObjects() { fOutputList = new TList; PostData(1,fOutputList); fOutputList->SetOwner(1); const Double_t binsPt[] = {0.,2.,4.,6.,8.,10.,12., 14.,16., 18.,20.,24.,28.,32.,38.,44.,50.,58.,66.,76.,86.,100.,114.,128.,150.,250.}; const Int_t nBinsPt = 25; const Double_t trackPtBins[42] = {0,0.2,0.4,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,8.,10.,12.,14.,16.,18.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,95.,100}; const Int_t nTrackPtBins = 41; if(fDebug) printf("%s: %d Creating output objects\n",(char*)__FILE__,__LINE__); // triggered - centrality - z-vtxs - r-vtx - ncontrib - ntracks - IsPileUp?- njetsrec - njetsgen - accepted const Int_t evtsDim = 7; Int_t evtsBins[evtsDim] = { 2, 22, 40, 20, 20, 2, 2}; Double_t evtsMins[evtsDim] = {-0.5, -5., -20., 0., 0., -0.5, -0.5}; Double_t evtsMaxs[evtsDim] = { 1.5, 105., 20., 2., 40., 1.5, 1.5}; // pt - eta - phi - is272 - charge const Int_t trksDim = 7; Int_t trksBins[trksDim] = { nTrackPtBins, 10, 10, 2, 2, 2, 3}; Double_t trksMins[trksDim] = { 0., -1., 0., -0.5,-0.5,-0.5, -1.5}; Double_t trksMaxs[trksDim] = { 0., 1., 2*TMath::Pi(), 1.5, 1.5, 1.5, 1.5}; // pt - eta - phi - AreaCh - AreaN - ptLeadingTrack/ptJet - nTracks const Int_t jetsDim = 7; Int_t jetsBins[jetsDim] = { nBinsPt, 20, 10,3, 5, 20, 30}; Double_t jetsMins[jetsDim] = { 0., -1., 0., 0., 0., 0., 0.}; Double_t jetsMaxs[jetsDim] = { 0., 1., 2*TMath::Pi(), 3., 5., 1.,30.}; // pt - eta - phi - AreaCH - AreaN - UEdensity - ptLeadingTrack/ptJet - ptOrig const Int_t jetsUEcorDim = 8; Int_t jetsUEcorBins[jetsUEcorDim] = { nBinsPt, 10, 10, 3, 5, 100, 20, nBinsPt}; Double_t jetsUEcorMins[jetsUEcorDim] = { 0., -1., 0., 0., 0., 0., 0., 0.}; Double_t jetsUEcorMaxs[jetsUEcorDim] = { 0., 1., 2*TMath::Pi(), 3., 5., 10., 1., 0.}; // leading jet pt - eta - sumAllPt1 - sumAllPt2 - Rpar const Int_t ueDim = 5; Int_t ueBins[ueDim] = { nBinsPt, 20, 100, 100, 100}; Double_t ueMins[ueDim] = { 0., -1, 0., 0., 0.}; Double_t ueMaxs[ueDim] = { 0., 1, 10., 10., 10.}; // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta - dPhi - isSamePtBin const Int_t matchingDim = 10; Int_t matchingBins[matchingDim] = {nBinsPt, nBinsPt, 20, 5, 2, 100, 10, 10, 10, 2 }; Double_t matchingMins[matchingDim] = { 0., 0., 0., 0.05, -0.5, 0, 0, 0, 0,-0.5}; Double_t matchingMaxs[matchingDim] = { 0., 0., 1., 0.55, 1.5, 1,0.5, 0.5,0.5, 1.5}; fhnEvent = new THnSparseF("fhnEvent", "fhnEvent;triggered;centrality;zVtx;rVtx;Ncontrib;IsPileUp;isAccepted", evtsDim, evtsBins, evtsMins, evtsMaxs); fhnTracks = new THnSparseF("fhnTracks", "fhnTraks;pT;eta;phi;isFilterMask272;isMask16;isMask256;charge;", trksDim, trksBins, trksMins, trksMaxs); fhnMC = new THnSparseF("fhnMC", "fhnTraks;pT;eta;phi;isFilterMask272;isMask16;isMask256;charge;", trksDim, trksBins, trksMins, trksMaxs); fhnMC2 = new THnSparseF("fhnMC2", "fhnTraks;pT;eta;phi;isFilterMask272;isMask16;isMask256;charge;", trksDim, trksBins, trksMins, trksMaxs); fhnRecJetsNoCut = new THnSparseF("fhnRecJetsNoCut", "fhnRecJetsNoCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs); fhnGenJetsNoCut = new THnSparseF("fhnGenJetsNoCut", "fhnGenJetsNoCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs); fhnRecJetsCut = new THnSparseF("fhnRecJetsCut", "fhnRecJetsCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs); fhnGenJetsCut = new THnSparseF("fhnGenJetsCut", "fhnGenJetsCut;pT;eta;phi;AreaCH;AreaN;z;nTracks;", jetsDim, jetsBins, jetsMins, jetsMaxs); //UE sparses fhnTrackUE = new THnSparseF("fhnTrackUE", "fhnTrackUE;pT-leading jet;eta;sum1;sum2;avgSum;", ueDim, ueBins, ueMins, ueMaxs); fhnParticleUE = new THnSparseF("fhnParticleUE","fhnParticleUE;pT-leading jet;eta;sum1;sum2;avgSum;", ueDim, ueBins, ueMins, ueMaxs); // pt - eta - phi - AreaCH - AreaN - UEdensity - ptLeadingTrack/ptJet - ptOrig fhnRecJetsTrackUEcor = new THnSparseF("fhnRecJetsTrackUEcor", "fhnRecJetsTrackUEcor;pT;eta;phi;AreaCH;AreaN;UE;z;pTorig;", jetsUEcorDim, jetsUEcorBins, jetsUEcorMins, jetsUEcorMaxs); fhnGenJetsTrackUEcor = new THnSparseF("fhnGenJetsTrackUEcor", "fhnGenJetsTrackUEcor;pT;eta;phi;AreaCH;AreaN;UE;z;pTorig;", jetsUEcorDim, jetsUEcorBins, jetsUEcorMins, jetsUEcorMaxs); // ptGen - ptRec - ptFraction - maxDist - bothWays fhnMatching = new THnSparseF("fhnMatching","fhnMatching;pTgen;pTrec;ptFraction;maxDist;bothWays;1-p_{T,rec}/p_{T,gen}; dR; dEta; dPhi;;isSamePtBin;", matchingDim, matchingBins, matchingMins, matchingMaxs); fhnTrackCorrMatching = new THnSparseF("fhnTrackCorrMatching","fhnTrackCorrMatching;pTgen;pTrec;ptFraction;maxDist;bothWays;1-p_{T,rec}/p_{T,gen}; dR; dEta; dPhi;isSamePtBin;", matchingDim, matchingBins, matchingMins, matchingMaxs); fhnTracks->SetBinEdges(0, trackPtBins); fhnTrackUE->SetBinEdges(0, trackPtBins); fhnParticleUE->SetBinEdges(0, trackPtBins); fhnRecJetsNoCut->SetBinEdges(0, binsPt); fhnGenJetsNoCut->SetBinEdges(0, binsPt); fhnRecJetsCut->SetBinEdges(0, binsPt); fhnGenJetsCut->SetBinEdges(0, binsPt); fhnRecJetsTrackUEcor->SetBinEdges(0,binsPt); fhnRecJetsTrackUEcor->SetBinEdges(7,binsPt); fhnGenJetsTrackUEcor->SetBinEdges(0,binsPt); fhnGenJetsTrackUEcor->SetBinEdges(7,binsPt); fhnMatching->SetBinEdges(0,binsPt); fhnMatching->SetBinEdges(1,binsPt); fhnTrackCorrMatching->SetBinEdges(0,binsPt); fhnTrackCorrMatching->SetBinEdges(1,binsPt); fhnEvent->Sumw2(); fhnTracks->Sumw2(); fhnMC->Sumw2(); fhnMC2->Sumw2(); fhnRecJetsNoCut->Sumw2(); fhnGenJetsNoCut->Sumw2(); fhnRecJetsCut->Sumw2(); fhnGenJetsCut->Sumw2(); fhnRecJetsTrackUEcor->Sumw2(); fhnGenJetsTrackUEcor->Sumw2(); fhnTrackUE->Sumw2(); fhnParticleUE->Sumw2(); fhnMatching->Sumw2(); fhnTrackCorrMatching->Sumw2(); fOutputList->Add(fhnEvent); fOutputList->Add(fhnTracks); fOutputList->Add(fhnMC); fOutputList->Add(fhnMC2); fOutputList->Add(fhnRecJetsNoCut); fOutputList->Add(fhnGenJetsNoCut); fOutputList->Add(fhnRecJetsCut); fOutputList->Add(fhnGenJetsCut); fOutputList->Add(fhnRecJetsTrackUEcor); fOutputList->Add(fhnGenJetsTrackUEcor); fOutputList->Add(fhnTrackUE); fOutputList->Add(fhnParticleUE); fOutputList->Add(fhnMatching); fOutputList->Add(fhnTrackCorrMatching); const Int_t trackUEdim = 11; //max track pt - ch_part_dens - ch_pt_dens - avg_pt_track - zero particles -tMAX_part_dens-tMIN_part_dens-tDIF_part_dens-tMAX_pt_dens-tMIN_pt_dens-tDIF_pt_dens Int_t trackUEbins[trackUEdim] = { 100, 50, 50, 50, 2, 50, 50, 50, 50, 50, 50}; Double_t trackUEmins[trackUEdim] = { 0., 0., 0., 0.,-0.5, 0., 0., 0., 0., 0., 0.}; Double_t trackUEmaxs[trackUEdim] = {100., 5., 5., 5., 1.5, 5., 5., 5., 5., 5., 5.}; fhnTrackUEanal = new THnSparseF("fhnTrackUEanal","fhnTrackUEanal;PTmax;#rho^{ch.part}_{T};#rho^{pT}_{T};avg-pt/track;zero_particles;#rho^{ch.part}_{transMAX};#rho^{ch.part}_{transMIN};#rho^{ch.part}_{transDIF};#rho^{pT}_{transMAX};#rho^{pT}_{transMIN};#rho^{pT}_{transDIF}", trackUEdim, trackUEbins, trackUEmins, trackUEmaxs); fhnTrackUEanal->Sumw2(); fOutputList->Add(fhnTrackUEanal); } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::UserExec(Option_t */*option*/) { Double_t evtContainer[7]; Bool_t isEventSelected = EventSelection(evtContainer); TList trackList; TList particleList; Int_t nTracks = GetListOfTracks(fTrackType, &trackList); Int_t nParticles = GetListOfTracks(fParticleType, &particleList); if(fDebug > 10 ) std::cout<< nTracks<< " " << nParticles< 0) tcaRecJets = dynamic_cast(fAODOut->FindListObject(fRecJetBranch.Data())); if(fAODOut && !tcaGenJets && fGenJetBranch.Length() > 0) tcaGenJets = dynamic_cast(fAODOut->FindListObject(fGenJetBranch.Data())); if(fAODOut && !tcaRecBckg && fRecBckgBranch.Length() > 0) tcaRecBckg = dynamic_cast(fAODOut->FindListObject(fRecBckgBranch.Data())); if(fAODOut && !tcaGenBckg && fGenBckgBranch.Length() > 0) tcaGenBckg = dynamic_cast(fAODOut->FindListObject(fGenBckgBranch.Data())); if(fAODExt) { if(!tcaRecJets && fRecJetBranch.Length() > 0) tcaRecJets = dynamic_cast(fAODExt->GetAOD()->FindListObject(fRecJetBranch.Data())); if(!tcaGenJets && fGenJetBranch.Length() > 0) tcaGenJets = dynamic_cast(fAODExt->GetAOD()->FindListObject(fGenJetBranch.Data())); if(!tcaRecBckg && fRecBckgBranch.Length() > 0) tcaRecBckg = dynamic_cast(fAODExt->GetAOD()->FindListObject(fRecBckgBranch.Data())); if(!tcaGenBckg && fGenBckgBranch.Length() > 0) tcaGenBckg = dynamic_cast(fAODExt->GetAOD()->FindListObject(fGenBckgBranch.Data())); } if(fAODIn) { if(!tcaRecJets && fRecJetBranch.Length() > 0) tcaRecJets = dynamic_cast(fAODIn->FindListObject(fRecJetBranch.Data())); if(!tcaGenJets && fGenJetBranch.Length() > 0) tcaGenJets = dynamic_cast(fAODIn->FindListObject(fGenJetBranch.Data())); if(!tcaRecBckg && fRecBckgBranch.Length() > 0) tcaRecBckg = dynamic_cast(fAODIn->FindListObject(fRecBckgBranch.Data())); if(!tcaGenBckg && fGenBckgBranch.Length() > 0) tcaGenBckg = dynamic_cast(fAODIn->FindListObject(fGenBckgBranch.Data())); } if(!tcaRecJets) rJ = -1; else rJ = tcaRecJets->GetEntries(); if(!tcaGenJets) gJ = -1; else gJ = tcaGenJets->GetEntries(); if(!tcaRecBckg) rB = -1; else rB = tcaRecBckg->GetEntries(); if(!tcaGenBckg) gB = -1; else gB = tcaGenBckg->GetEntries(); TList jetRecListNoCut; TList jetGenListNoCut; Int_t nRecJetsNoCut = (rJ>0?GetListOfJets(tcaRecJets, &jetRecListNoCut, kFALSE):0); Int_t nGenJetsNoCut = (gJ>0?GetListOfJets(tcaGenJets, &jetGenListNoCut, kFALSE):0); TList jetRecListCut; TList jetGenListCut; Int_t nRecJetsCut = (rJ>0?GetListOfJets(tcaRecJets, &jetRecListCut, kTRUE):0); Int_t nGenJetsCut = (gJ>0?GetListOfJets(tcaGenJets, &jetGenListCut, kTRUE):0); TList bckgRecListCut; TList bckgGenListCut; Int_t nRecBckgCut = (rB>0?GetListOfJets(tcaRecBckg, &bckgRecListCut, kTRUE):0); Int_t nGenBckgCut = (gB>0?GetListOfJets(tcaGenBckg, &bckgGenListCut, kTRUE):0); if(fDebug > 10) std::cout<Fill(evtContainer); if(!isEventSelected) { PostData(1,fOutputList); return; } if(kDoUEanalysis) DoUEAnalysis(&trackList, (Double_t)0.15, (Double_t)0.9); Double_t recUEdensity = GetUE(&jetRecListCut, &trackList, fRecJetR, fhnTrackUE); Double_t genUEdensity = GetUE(&jetGenListCut, &particleList, fGenJetR, fhnParticleUE); TList trackUEcorrJetRecListCut; TList trackUEcorrJetGenListCut; CorrectForUE(&jetRecListCut, recUEdensity, &trackUEcorrJetRecListCut, fhnRecJetsTrackUEcor); CorrectForUE(&jetGenListCut, genUEdensity, &trackUEcorrJetGenListCut, fhnGenJetsTrackUEcor); FillJetContainer(&jetRecListNoCut, fhnRecJetsNoCut); FillJetContainer(&jetGenListNoCut, fhnGenJetsNoCut); FillJetContainer(&jetRecListCut, fhnRecJetsCut); FillJetContainer(&jetGenListCut, fhnGenJetsCut); if(trackUEcorrJetRecListCut.GetEntries() != 0 && trackUEcorrJetGenListCut.GetEntries() != 0 ) { MatchJets(kTRUE, &trackUEcorrJetRecListCut, &trackUEcorrJetGenListCut, 0.1, fhnTrackCorrMatching); MatchJets(kTRUE, &trackUEcorrJetRecListCut, &trackUEcorrJetGenListCut, 0.3, fhnTrackCorrMatching); MatchJets(kTRUE, &trackUEcorrJetRecListCut, &trackUEcorrJetGenListCut, 0.5, fhnTrackCorrMatching); } if(jetRecListCut.GetEntries() != 0 && jetGenListCut.GetEntries() != 0 ) { MatchJets(kTRUE, &jetRecListCut, &jetGenListCut, 0.1, fhnMatching); MatchJets(kTRUE, &jetRecListCut, &jetGenListCut, 0.3, fhnMatching); MatchJets(kTRUE, &jetRecListCut, &jetGenListCut, 0.5, fhnMatching); } PostData(1,fOutputList); return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::Terminate(Option_t *) { if(fDebug) printf("%s: %d Terminating\n",(char*)__FILE__, __LINE__); return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::SetRecJetBranch(TString s) { fRecJetBranch = s; if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fRecJetR = 0.2; else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fRecJetR = 0.3; else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fRecJetR = 0.4; else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fRecJetR = 0.5; else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fRecJetR = 0.6; else fRecJetR = 0; return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::SetGenJetBranch(TString s) { fGenJetBranch = s; if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fGenJetR = 0.2; else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fGenJetR = 0.3; else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fGenJetR = 0.4; else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fGenJetR = 0.5; else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fGenJetR = 0.6; else fGenJetR = 0; return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::SetRecBckgBranch(TString s) { fRecBckgBranch = s; if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fRecBckgR = 0.2; else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fRecBckgR = 0.3; else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fRecBckgR = 0.4; else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fRecBckgR = 0.5; else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fRecBckgR = 0.6; else fRecBckgR = 0; return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::SetGenBckgBranch(TString s) { fGenBckgBranch = s; if(s.Contains("KT02") || s.Contains("SISCONE02") || s.Contains("UA102")) fGenBckgR = 0.2; else if(s.Contains("KT03") || s.Contains("SISCONE03") || s.Contains("UA103")) fGenBckgR = 0.3; else if(s.Contains("KT04") || s.Contains("SISCONE04") || s.Contains("UA104")) fGenBckgR = 0.4; else if(s.Contains("KT05") || s.Contains("SISCONE05") || s.Contains("UA105")) fGenBckgR = 0.5; else if(s.Contains("KT06") || s.Contains("SISCONE06") || s.Contains("UA106")) fGenBckgR = 0.6; else fGenBckgR = 0; return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::SetVertexCuts(Int_t nCont, Double_t Zcut, Double_t Rcut) { if(fDebug) printf("%s: %d Setting Vertex Cuts\n",(char*)__FILE__,__LINE__); nVtxContCut = nCont; fVtxZcut = Zcut; fVtxRcut = Rcut; return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::SetTrackCuts(Double_t ptMin, Double_t ptMax, Double_t etaMax) { if(fDebug > 1) printf("%s: %d Track cuts set\n",(char*)__FILE__,__LINE__); trackPtMin = ptMin; trackPtMax = ptMax; trackEtaAbsMax = etaMax; return; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::SetJetCuts(Double_t ptMin, Double_t eta, Double_t z) { jetPtMin = ptMin; jetEtaCut = eta; jetZmax = z; return; } //--------------------------------------------------------------------------------------------------- Bool_t AliAnalysisTaskPPJetSpectra::EventSelection(Double_t evtContainer[6]) { if(fDebug > 1) printf("%s: %d UserExec analysing event %d.\n", (char*)__FILE__, __LINE__, (Int_t)fEntry); // Trigger selection AliInputEventHandler* inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); if(fDebug > 2 ) printf("IsEventSelected: %d eventSelectionMask: %d\n", (Int_t)inputHandler->IsEventSelected(), (Int_t)fEvtSelectionMask); Bool_t isTriggerSelected = inputHandler->IsEventSelected() & fEvtSelectionMask; evtContainer[0] = (Int_t)isTriggerSelected; fESD = dynamic_cast(InputEvent()); if(!fESD && fDebug > 2) printf("%s: %d No ESD event found\n",(char*)__FILE__, __LINE__); if(!fESD) fAODIn = dynamic_cast(InputEvent()); fAODOut = dynamic_cast(AODEvent()); if(!fESD)fAOD = fAODIn; else fAOD = fAODOut; if(fNonStdFile.Length()!=0) { AliAODHandler *aodH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); fAODExt = (aodH?aodH->GetExtension(fNonStdFile.Data()):0); } TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); Float_t centrality = -1; if(fEventClass > 0) { if(handler->InheritsFrom("AliAODHandler")) centrality = fAODIn->GetHeader()->GetCentrality(); else if(fESD) centrality = fESD->GetCentrality()->GetCentralityPercentile("V0M"); else centrality = AliAnalysisHelperJetTasks::EventClass(); } evtContainer[1] = centrality; // Event selection AliAODVertex* primVtx = fAOD->GetPrimaryVertex(); if(!primVtx) { evtContainer[2] = 1e10; evtContainer[3] = 1e10; evtContainer[4] = 1e10; return kFALSE; } else { primVtx->Print(); evtContainer[2] = primVtx->GetZ(); evtContainer[3] = TMath::Sqrt(primVtx->GetX()*primVtx->GetX() + primVtx->GetY()*primVtx->GetY() ); evtContainer[4] = primVtx->GetNContributors(); } if(fRejectPileUp && AliAnalysisHelperJetTasks::IsPileUp()){ if (fDebug > 1) Printf("%s:%d SPD pileup: event REJECTED...", (char*)__FILE__,__LINE__); evtContainer[5] = 0.; return kFALSE; } else evtContainer[5] = 1.; if(evtContainer[0] == 0) { if(fDebug > 2) printf("%s: %d Event rejected: Trigger\n",(char*)__FILE__, __LINE__); return kFALSE; } if(evtContainer[1] >= 0 && evtContainer[1] < 10*((Int_t)fEventClass/10) && evtContainer[1] > 10*((Int_t)fEventClass/10 + 1) ) { if(fDebug > 2) printf("%s: %d Event rejected: Centrality\n",(char*)__FILE__, __LINE__); return kFALSE; } if(TMath::Abs(evtContainer[2]) > fVtxZcut) { if(fDebug > 2) printf("%s: %d Event rejected: Vertex position (z)\n",(char*)__FILE__, __LINE__); return kFALSE; } if(evtContainer[3] > fVtxRcut) { if(fDebug > 2) printf("%s: %d Event rejected: Vertex position (x-y)\n",(char*)__FILE__,__LINE__); return kFALSE; } if(evtContainer[4] < nVtxContCut) { if(fDebug > 2) printf("%s: %d Event rejected: Vertex contributors \n",(char*)__FILE__, __LINE__); return kFALSE; } return kTRUE; } //--------------------------------------------------------------------------------------------------- Int_t AliAnalysisTaskPPJetSpectra::GetListOfTracks(Int_t trackType, TList *trackList) { if(fDebug) printf("%s: %d Getting Tracks\n",(char*)__FILE__,__LINE__); if(trackType == AliAnalysisTaskPPJetSpectra::kNone) { if(fDebug) printf("%s: %d No track type selected", (char*)__FILE__,__LINE__); return 0; } if(1 == trackType) { for(Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) { Bool_t isOK = kTRUE; Bool_t isOK16 = kTRUE; Bool_t isOK256 = kTRUE; AliAODTrack* track = (AliAODTrack*)fAOD->GetTrack(i); if(!track) continue; if(nTrackFilter>0 && !track->TestFilterBit(nTrackFilter) ) isOK = kFALSE; if(nTrackFilter>0 && !track->TestFilterBit(16) ) isOK16 = kFALSE; if(nTrackFilter>0 && !track->TestFilterBit(256) ) isOK256 = kFALSE; if(track->Pt() < trackPtMin || track->Pt() > trackPtMax) isOK = kFALSE; if(trackEtaAbsMax < TMath::Abs(track->Eta())) isOK = kFALSE; Double_t tmpContainer[] = {track->Pt(), track->Eta(), track->Phi(), (Double_t)isOK, (Double_t)isOK16, (Double_t)isOK256, (Double_t)track->Charge()/3}; if(isOK) trackList->Add(track); fhnTracks->Fill(tmpContainer); } } else { TClonesArray *tca=dynamic_cast(fAOD->FindListObject(AliAODMCParticle::StdBranchName())); if(!tca) { if(fDebug > 2) printf("no branch %s\n", AliAODMCParticle::StdBranchName()); } for(Int_t i = 0; i < tca->GetEntries(); i++) { AliAODMCParticle *particle = (AliAODMCParticle*)tca->At(i); if(!particle) continue; if(!particle->IsPhysicalPrimary()) continue; Double_t charge = particle->Charge(); if( charge > 0 ) charge = 1; if( charge < 0 ) charge = -1; if( (trackType == 3) && (charge == 0) ) continue; trackList->Add(particle); Double_t tmpContainer[] = {particle->Pt(), particle->Eta(), particle->Phi(), -999, -999, -999, charge }; if(trackType == 2) fhnMC->Fill(tmpContainer); else fhnMC2->Fill(tmpContainer); } } trackList->Sort(); return trackList->GetEntries(); } //--------------------------------------------------------------------------------------------------- Int_t AliAnalysisTaskPPJetSpectra::GetListOfJets(TClonesArray* jca, TList* jetList, Bool_t useCuts) { if(fDebug) printf("%s: %d Getting list of jets\n",(char*)__FILE__,__LINE__); Double_t ptCut = (useCuts && jetPtMin>0 ?jetPtMin :0.1); Double_t etaCut = (useCuts && jetEtaCut>0?jetEtaCut:0.9); Double_t zCut = (useCuts && jetZmax>0 ?jetZmax :1); if(!jca){ return -2; } for(Int_t i =0 ; i < jca->GetEntries(); i++) { AliAODJet* jet = (AliAODJet*)jca->At(i); if(!jet) continue; Double_t zLeading = -1; // if(jet->GetPtLeading() <= 0) { // TRefArray* jetRef = (TRefArray*)jet->GetRefTracks(); // for(int j = 0; j < jetRef->GetEntries(); j++) { // AliVTrack *track = (AliVTrack*)jet->GetTrack(j); // if(!track) continue; // cout<Pt()<<" "<Pt() > zLeading) zLeading = track->Pt(); // cout<<__LINE__<Pt() > 0) zLeading /= jet->Pt(); // else zLeading = -1; // } //cuty if(ptCut > 0 && ptCut > jet->Pt() ) continue; if(etaCut > 0 && etaCut < TMath::Abs(jet->Eta()) ) continue; if(zCut > 0 && zCut < zLeading && zLeading > 0 ) continue; jetList->Add(jet); } jetList->Sort(); return jetList->GetEntries(); } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::FillJetContainer(TList* jets, THnSparseF* container) { if(!jets || !container) return; for(Int_t i =0 ; i < jets->GetEntries(); i++) { AliAODJet* jet = (AliAODJet*)jets->At(i); if(!jet) continue; // pt - eta - phi - AreaCH - AreaN - ptLeadingTrack/ptJet - ptOrig Double_t tmpContainer[] = {jet->Pt(), jet->Eta(), jet->Phi(), jet->EffectiveAreaCharged(), jet->EffectiveAreaNeutral(), jet->GetPtLeading()/jet->Pt(), (Double_t)((TRefArray*)jet->GetRefTracks())->GetEntries() }; container->Fill(tmpContainer); } return; } //--------------------------------------------------------------------------------------------------- Double_t AliAnalysisTaskPPJetSpectra::GetUE(TList* jets, TList* tracks, Double_t Rpar,THnSparseF* thnCont) { if(!jets || !tracks) { if(fDebug) printf("%s: %d No jets or tracks\n",(char*)__FILE__,__LINE__); return 0; } //find leading jet AliAODJet *jet = 0; for(Int_t iTmp = 0; iTmp < jets->GetEntries(); iTmp++) { AliAODJet* tmp = (AliAODJet*)jets->At(iTmp); if(!tmp) continue; if(!jet) jet = tmp; else if( jet->Pt() < tmp->Pt() ) jet = tmp; } if(!jet) return 0; //get track in perpendicular cone //TVector3 for jet and perpendicular cones with same eta axis TVector3 j(jet->Px(), jet->Py(), jet->Pz()); TVector3 p1(j); TVector3 p2(j); p1.RotateZ(TMath::Pi()/2.); p2.RotateZ(-TMath::Pi()/2.); Double_t sumAllPt1 = 0; Double_t sumAllPt2 = 0; for(int i = 0; i < tracks->GetEntries(); i++) { AliVParticle* tr = (AliVParticle*)tracks->At(i); if(!tr) continue; TVector3 v(tr->Px(), tr->Py(), tr->Pz()); Double_t dR1 = v.DrEtaPhi(p1); Double_t dR2 = v.DrEtaPhi(p2); if(dR1 < Rpar) sumAllPt1+=v.Pt(); if(dR2 < Rpar) sumAllPt2+=v.Pt(); } if(Rpar != 0) { sumAllPt1/=(TMath::Pi()*Rpar*Rpar); sumAllPt2/=(TMath::Pi()*Rpar*Rpar); } Double_t container[] = {j.Pt(), j.Eta(), sumAllPt1, sumAllPt2, (sumAllPt1+sumAllPt2)/2.}; thnCont->Fill(container); return (sumAllPt1+sumAllPt2)/2.; } //--------------------------------------------------------------------------------------------------- Double_t AliAnalysisTaskPPJetSpectra::GetBckgUE(TList* jets, Double_t Rpar, Bool_t isSkipped, THnSparseF* cont) { if(!jets) return 0; Int_t leading01 = -1; Int_t leading02 = -1; Double_t leading01pt = 0; Double_t leading02pt = 0; if(!isSkipped) { for(Int_t i = 0; i < jets->GetEntries(); i++) { AliAODJet* tmp = (AliAODJet*)jets->At(i); if(!tmp) continue; if(tmp->Pt() > leading01pt) { leading02 = leading01; leading01 = i; leading02pt = leading01pt; leading01pt = tmp->Pt(); } else if(tmp->Pt() > leading02) { leading02 = i; leading02pt = tmp->Pt(); } } } Int_t subLeading01 = -1; Int_t subLeading02 = -1; Double_t subLeading01pt = 0; Double_t subLeading02pt = 0; Double_t subLeading01area = 0; Double_t subLeading02area = 0; Double_t avgAll = 0; Double_t areaAll = 0; //search for subleading jets assuming leading 2 (skip02) jets are skipped if(jets->GetEntries() > 0) { for(Int_t i = 0; i < jets->GetEntries(); i++) { if(i == leading01 || i == leading02) continue; AliAODJet* tmp = (AliAODJet*)jets->At(i); if(!tmp) continue; if( subLeading01pt < tmp->Pt() ) { subLeading02 = subLeading01; subLeading01 = i; subLeading02pt = subLeading01pt; subLeading01pt = tmp->Pt(); subLeading02area = subLeading01area; subLeading01area = tmp->EffectiveAreaCharged(); } else if( subLeading02pt < tmp->Pt() ) { subLeading02 = i; subLeading02pt = tmp->Pt(); subLeading02area = tmp->EffectiveAreaCharged(); } } //do avg jet pt Int_t nJets = 0; for(int i = 0; i < jets->GetEntries();i++) { if(i == leading01 || i == leading02) continue; AliAODJet* tmp = (AliAODJet*)jets->At(i); if(!tmp) continue; nJets++; avgAll+= tmp->Pt(); areaAll += tmp->EffectiveAreaCharged(); } } if(avgAll!=0) avgAll/=areaAll; if(subLeading01pt!=0) subLeading01pt/=subLeading01area; if(subLeading02pt!=0)subLeading02pt/=subLeading02area; Double_t tmpContainer[]={subLeading01pt,subLeading02pt,avgAll,Rpar*10}; cont->Fill(tmpContainer); return avgAll; } //--------------------------------------------------------------------------------------------------- Int_t AliAnalysisTaskPPJetSpectra::CorrectForUE(TList* jets, Double_t UE, TList* newList, THnSparseF* container) { if(!jets) return 0; for(Int_t i = 0; i < jets->GetEntries(); i++) { AliAODJet* tmpJet = (AliAODJet*)jets->At(i); if(!tmpJet) continue; AliAODJet* jet = (AliAODJet*)tmpJet->Clone(Form("%s%s",tmpJet->GetName(),"_clone")); Double_t corrUE = UE*(jet->EffectiveAreaCharged() + jet->EffectiveAreaNeutral() ); if(jet->Pt() > corrUE) jet->SetPtEtaPhiM(jet->Pt()- corrUE, jet->Eta(), jet->Phi(), jet->M()); else jet->SetPtEtaPhiM(0, jet->Eta(), jet->Phi(), jet->M()); newList->Add(jet); // pt - eta - phi - AreaCH - AreaN - UEdensity - ptLeadingTrack/ptJet - ptOrig Double_t cont[8]={jet->Pt(), jet->Eta(), jet->Phi(), jet->EffectiveAreaCharged(),jet->EffectiveAreaNeutral(), corrUE, (jet->Pt() > 0?jet->GetPtLeading() / jet->Pt():-1),tmpJet->Pt()}; container->Fill(cont); } newList->Sort(); return newList->GetEntries(); } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::MatchJets(Bool_t doBothWay, TList* recList, TList* genList, Float_t maxDist, THnSparseF* container) { if(!genList || !recList) return; Int_t mode = 1; const Int_t kGenJets = genList->GetEntries(); const Int_t kRecJets = recList->GetEntries(); TArrayI iMatchIndex(kGenJets); TArrayF fPtFraction(kGenJets); TArrayI aGenIndex(kRecJets); TArrayI aRecIndex(kGenJets); Double_t genPt = -1; Double_t recPt = -1; Double_t ptFraction = -1; //-- Closest jets to generated and checked for E-fraction if(!doBothWay) { AliAnalysisHelperJetTasks::GetJetMatching(genList, kGenJets, recList, kRecJets, iMatchIndex, fPtFraction, fDebug, maxDist, mode); for(Int_t ig = 0; ig < kGenJets; ig++) { //get gen-jet AliAODJet* genJet = (AliAODJet*)genList->At(ig); if(!genJet) continue; genPt = genJet->Pt(); //get rec-jet Int_t ir = iMatchIndex[ig]; if(ir<0||ir>=recList->GetEntries()) continue; AliAODJet* recJet = (AliAODJet*)recList->At(ir); if(!recJet) continue; recPt = recJet->Pt(); ptFraction = fPtFraction[ig]; Double_t dR = recJet->DeltaR(genJet); Double_t dEta = recJet->Eta() - genJet->Eta(); Double_t dPhi2 = dR*dR - dEta*dEta; Bool_t isSame = kFALSE; if(CheckPtBin(genPt) == CheckPtBin(recPt)) isSame = kTRUE; // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta - dPhi Double_t cont[] = {genPt,recPt,ptFraction,maxDist,(Float_t)doBothWay, (genPt != 0 ? (1. - recPt/genPt):1), dR, dEta, (dPhi2>=0?TMath::Sqrt(dPhi2):-1), (Double_t)isSame}; container->Fill(cont); } } //-- Closest rec jet to gen and vice-versa & check for E-fraction else { AliAnalysisHelperJetTasks::GetClosestJets(genList, kGenJets, recList, kRecJets, aGenIndex, aRecIndex,fDebug, maxDist); for(Int_t ig = 0; ig < kGenJets; ig++) { //get gen-jet AliAODJet* genJet = (AliAODJet*)genList->At(ig); if(!genJet) continue; genPt = genJet->Pt(); //get rec-jet Int_t ir = aRecIndex[ig]; if(ir>=0&&irGetEntries()) { AliAODJet *recJet = (AliAODJet*)recList->At(ir); if(!recJet) continue; recPt = recJet->Pt(); Double_t dR = recJet->DeltaR(genJet); Double_t dEta = recJet->Eta() - genJet->Eta(); Double_t dPhi2 = dR*dR - dEta*dEta; ptFraction = AliAnalysisHelperJetTasks::GetFractionOfJet(recJet,genJet,mode); Bool_t isSame = kFALSE; if(CheckPtBin(genPt) == CheckPtBin(recPt)) isSame = kTRUE; // ptGen - ptRec - ptFraction - maxDist - bothWays - (1 - rec/gen) - dR - dEta - dPhi Double_t cont[] = {genPt,recPt,ptFraction,maxDist,(Float_t)doBothWay, (genPt != 0 ? (1. - recPt/genPt):1), dR, dEta, (dPhi2>=0?TMath::Sqrt(dPhi2):-1), (Double_t)isSame }; container->Fill(cont); } } } return; } //--------------------------------------------------------------------------------------------------- Int_t AliAnalysisTaskPPJetSpectra::CheckPtBin(Double_t pT) { Double_t bins[] = {0.,2.,4.,6.,8.,10.,12., 14.,16., 18.,20.,24.,28.,32.,38.,44.,50.,58.,66.,76.,86.,100.,114.,128.,150.,250.}; Double_t nBins = 25; for(Int_t i = 0; i <= nBins; i++) { if(pT < bins[i] ) return i; } return nBins+1; } //--------------------------------------------------------------------------------------------------- void AliAnalysisTaskPPJetSpectra::DoUEAnalysis(TList* tracks, Double_t PTcut, Double_t ETAcut) { Double_t PTmax = 0; int iMax = -1; for(int i = 0; i < tracks->GetEntries(); i++) { AliVParticle *part = (AliVParticle*)tracks->At(i); if(!part) continue; if(part->Pt() < PTcut) continue; if(TMath::Abs(part->Eta()) > ETAcut ) continue; if(PTmax < part->Pt()) { PTmax = part->Pt(); iMax = i; } } if(PTmax < trackPtMin) return; Double_t ch_part_density = 0; Double_t ch_part_pt_density = 0; Double_t avg_track_pt = 0; Double_t zero_ch_part = 0; Double_t transMAX_part_density = 0; Double_t transMIN_part_density = 0; Double_t transMAX_pt_density = 0; Double_t transMIN_pt_density = 0; Double_t transDIF_pt_density = 0; Double_t transDIF_part_density = 0; Double_t A1_part_density = 0; Double_t A2_part_density = 0; Double_t A1_pt_density = 0; Double_t A2_pt_density = 0; AliVParticle *maxParticle = (AliVParticle*)tracks->At(iMax); Double_t PHImax = maxParticle->Phi(); for(int i = 0; i < tracks->GetEntries(); i++) { AliVParticle* part = (AliVParticle*)tracks->At(i); if(!part) continue; if(part->Pt() < PTcut) continue; Bool_t isArea1 = kFALSE; Bool_t isArea2 = kFALSE; Double_t dPhi = part->Phi() - PHImax; if(dPhi > TMath::Pi()) dPhi -= TMath::Pi(); if(dPhi < -TMath::Pi()) dPhi +=TMath::Pi(); if(TMath::Abs(dPhi) > TMath::Pi()/3 && TMath::Abs(dPhi) < 2*TMath::Pi()/3 ) { if(dPhi>0) isArea1 = kTRUE; else isArea2 = kTRUE; } if(!isArea1 && !isArea2) continue; ch_part_density += 1; ch_part_pt_density += part->Pt(); avg_track_pt += part->Pt(); if(isArea1) { A1_part_density += 1; A1_pt_density += part->Pt(); } else if(isArea2) { A2_part_density += 1; A2_pt_density += part->Pt(); } } if(ch_part_density > 0) avg_track_pt/=ch_part_density; else zero_ch_part = 1; ch_part_density*=3./(4.*TMath::Pi()*ETAcut); ch_part_pt_density*= 3./(4.*TMath::Pi()*ETAcut); A1_part_density*=3./(4.*TMath::Pi()*ETAcut); A1_pt_density*= 3./(4.*TMath::Pi()*ETAcut); A2_part_density*=3./(4.*TMath::Pi()*ETAcut); A2_pt_density*= 3./(4.*TMath::Pi()*ETAcut); if(A1_part_density < A2_part_density) { transMAX_part_density = A2_part_density; transMIN_part_density = A1_part_density; } else { transMAX_part_density = A1_part_density; transMIN_part_density = A2_part_density; } if(A1_pt_density < A2_pt_density) { transMAX_pt_density = A2_pt_density; transMIN_pt_density = A1_pt_density; } else { transMAX_pt_density = A1_pt_density; transMIN_pt_density = A2_pt_density; } transDIF_pt_density = transMAX_pt_density - transMIN_pt_density; transDIF_part_density = transMAX_part_density - transMIN_part_density; Double_t container[] = {PTmax, ch_part_density, ch_part_pt_density, avg_track_pt, zero_ch_part, transMAX_part_density, transMIN_part_density, transDIF_part_density, transMAX_pt_density, transMIN_pt_density, transDIF_pt_density}; fhnTrackUEanal->Fill(container); }