X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliAnalysisTaskJetCluster.cxx;h=8173f51041f1c3a62c4f2234669ba11073332811;hb=9188451d40f248875b09fa4c77f6dbb5438a9473;hp=583128727d58f601f6712a74290c5b8e2e8a35b0;hpb=ea950747c684d5786929bf2db5b1ff50fca2ae2d;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx index 583128727d5..8173f51041f 100644 --- a/JETAN/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/AliAnalysisTaskJetCluster.cxx @@ -32,10 +32,12 @@ #include #include #include +#include #include #include #include #include "TDatabasePDG.h" +#include #include "AliAnalysisTaskJetCluster.h" #include "AliAnalysisManager.h" @@ -65,34 +67,50 @@ // get info on how fastjet was configured #include "fastjet/config.h" +using std::vector; ClassImp(AliAnalysisTaskJetCluster) AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){ + // + // Destructor + // + delete fRef; delete fRandom; if(fTCAJetsOut)fTCAJetsOut->Delete(); delete fTCAJetsOut; + if(fTCAJetsOutRan)fTCAJetsOutRan->Delete(); delete fTCAJetsOutRan; + if(fTCARandomConesOut)fTCARandomConesOut->Delete(); delete fTCARandomConesOut; + if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete(); delete fTCARandomConesOutRan; + if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset(); delete fAODJetBackgroundOut; } -AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), +AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): + AliAnalysisTaskSE(), fAOD(0x0), fAODExtension(0x0), fRef(new TRefArray), fUseAODTrackInput(kFALSE), fUseAODMCInput(kFALSE), - fUseGlobalSelection(kFALSE), fUseBackgroundCalc(kFALSE), + fEventSelection(kFALSE), + fRequireVZEROAC(kFALSE), + fRequireTZEROvtx(kFALSE), + fUseHFcuts(kFALSE), fFilterMask(0), + fFilterMaskBestPt(0), + fFilterType(0), + fJetTypes(kJet), fTrackTypeRec(kTrackUndef), fTrackTypeGen(kTrackUndef), fNSkipLeadingRan(0), @@ -101,15 +119,39 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fAvgTrials(1), fExternalWeight(1), fTrackEtaWindow(0.9), + fRequireITSRefit(0), + fApplySharedClusterCut(0), fRecEtaWindow(0.5), fTrackPtCut(0.), fJetOutputMinPt(0.150), + fMaxTrackPtInJet(100.), fJetTriggerPtCut(0), + fVtxZCut(8), + fVtxR2Cut(1), fCentCutUp(0), fCentCutLo(0), + fStoreRhoLeadingTrackCorr(kFALSE), fNonStdBranch(""), fBackgroundBranch(""), fNonStdFile(""), + fMomResH1(0x0), + fMomResH2(0x0), + fMomResH3(0x0), + fMomResH1Fit(0x0), + fMomResH2Fit(0x0), + fMomResH3Fit(0x0), + fhEffH1(0x0), + fhEffH2(0x0), + fhEffH3(0x0), + fUseTrPtResolutionSmearing(kFALSE), + fUseDiceEfficiency(kFALSE), + fDiceEfficiencyMinPt(-1.), + fUseTrPtResolutionFromOADB(kFALSE), + fUseTrEfficiencyFromOADB(kFALSE), + fPathTrPtResolution(""), + fPathTrEfficiency(""), + fChangeEfficiencyFraction(0.), + fEfficiencyFixed(1.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -178,8 +220,39 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fh2PtNchNRan(0x0), fh2TracksLeadingJetPhiPtRan(0x0), fh2TracksLeadingJetPhiPtWRan(0x0), + fh3CentvsRhoLeadingTrackPt(0x0), + fh3CentvsSigmaLeadingTrackPt(0x0), + fh3MultvsRhoLeadingTrackPt(0x0), + fh3MultvsSigmaLeadingTrackPt(0x0), + fh3CentvsRhoLeadingTrackPtQ1(0x0), + fh3CentvsRhoLeadingTrackPtQ2(0x0), + fh3CentvsRhoLeadingTrackPtQ3(0x0), + fh3CentvsRhoLeadingTrackPtQ4(0x0), + fh3CentvsSigmaLeadingTrackPtQ1(0x0), + fh3CentvsSigmaLeadingTrackPtQ2(0x0), + fh3CentvsSigmaLeadingTrackPtQ3(0x0), + fh3CentvsSigmaLeadingTrackPtQ4(0x0), + fh3MultvsRhoLeadingTrackPtQ1(0x0), + fh3MultvsRhoLeadingTrackPtQ2(0x0), + fh3MultvsRhoLeadingTrackPtQ3(0x0), + fh3MultvsRhoLeadingTrackPtQ4(0x0), + fh3MultvsSigmaLeadingTrackPtQ1(0x0), + fh3MultvsSigmaLeadingTrackPtQ2(0x0), + fh3MultvsSigmaLeadingTrackPtQ3(0x0), + fh3MultvsSigmaLeadingTrackPtQ4(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0), + fh2PtGenPtSmeared(0x0), + fp1Efficiency(0x0), + fp1PtResolution(0x0), fHistList(0x0) { + // + // Constructor + // + for(int i = 0;i<3;i++){ fh1BiARandomCones[i] = 0; fh1BiARandomConesRan[i] = 0; @@ -199,9 +272,15 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fRef(new TRefArray), fUseAODTrackInput(kFALSE), fUseAODMCInput(kFALSE), - fUseGlobalSelection(kFALSE), fUseBackgroundCalc(kFALSE), + fEventSelection(kFALSE), + fRequireVZEROAC(kFALSE), + fRequireTZEROvtx(kFALSE), + fUseHFcuts(kFALSE), fFilterMask(0), + fFilterMaskBestPt(0), + fFilterType(0), + fJetTypes(kJet), fTrackTypeRec(kTrackUndef), fTrackTypeGen(kTrackUndef), fNSkipLeadingRan(0), @@ -210,15 +289,39 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fAvgTrials(1), fExternalWeight(1), fTrackEtaWindow(0.9), + fRequireITSRefit(0), + fApplySharedClusterCut(0), fRecEtaWindow(0.5), fTrackPtCut(0.), fJetOutputMinPt(0.150), + fMaxTrackPtInJet(100.), fJetTriggerPtCut(0), + fVtxZCut(8), + fVtxR2Cut(1), fCentCutUp(0), fCentCutLo(0), + fStoreRhoLeadingTrackCorr(kFALSE), fNonStdBranch(""), fBackgroundBranch(""), fNonStdFile(""), + fMomResH1(0x0), + fMomResH2(0x0), + fMomResH3(0x0), + fMomResH1Fit(0x0), + fMomResH2Fit(0x0), + fMomResH3Fit(0x0), + fhEffH1(0x0), + fhEffH2(0x0), + fhEffH3(0x0), + fUseTrPtResolutionSmearing(kFALSE), + fUseDiceEfficiency(kFALSE), + fDiceEfficiencyMinPt(-1.), + fUseTrPtResolutionFromOADB(kFALSE), + fUseTrEfficiencyFromOADB(kFALSE), + fPathTrPtResolution(""), + fPathTrEfficiency(""), + fChangeEfficiencyFraction(0.), + fEfficiencyFixed(1.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -287,8 +390,39 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fh2PtNchNRan(0x0), fh2TracksLeadingJetPhiPtRan(0x0), fh2TracksLeadingJetPhiPtWRan(0x0), + fh3CentvsRhoLeadingTrackPt(0x0), + fh3CentvsSigmaLeadingTrackPt(0x0), + fh3MultvsRhoLeadingTrackPt(0x0), + fh3MultvsSigmaLeadingTrackPt(0x0), + fh3CentvsRhoLeadingTrackPtQ1(0x0), + fh3CentvsRhoLeadingTrackPtQ2(0x0), + fh3CentvsRhoLeadingTrackPtQ3(0x0), + fh3CentvsRhoLeadingTrackPtQ4(0x0), + fh3CentvsSigmaLeadingTrackPtQ1(0x0), + fh3CentvsSigmaLeadingTrackPtQ2(0x0), + fh3CentvsSigmaLeadingTrackPtQ3(0x0), + fh3CentvsSigmaLeadingTrackPtQ4(0x0), + fh3MultvsRhoLeadingTrackPtQ1(0x0), + fh3MultvsRhoLeadingTrackPtQ2(0x0), + fh3MultvsRhoLeadingTrackPtQ3(0x0), + fh3MultvsRhoLeadingTrackPtQ4(0x0), + fh3MultvsSigmaLeadingTrackPtQ1(0x0), + fh3MultvsSigmaLeadingTrackPtQ2(0x0), + fh3MultvsSigmaLeadingTrackPtQ3(0x0), + fh3MultvsSigmaLeadingTrackPtQ4(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0), + fh2PtGenPtSmeared(0x0), + fp1Efficiency(0x0), + fp1PtResolution(0x0), fHistList(0x0) { + // + // named ctor + // + for(int i = 0;i<3;i++){ fh1BiARandomCones[i] = 0; fh1BiARandomConesRan[i] = 0; @@ -299,7 +433,7 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fh2TracksLeadingJetPhiPtC[i] = 0; fh2TracksLeadingJetPhiPtWC[i] = 0; } - DefineOutput(1, TList::Class()); + DefineOutput(1, TList::Class()); } @@ -337,38 +471,61 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() // -> cleared in the UserExec.... // here we can also have the case that the brnaches are written to a separate file - fTCAJetsOut = new TClonesArray("AliAODJet", 0); - fTCAJetsOut->SetName(fNonStdBranch.Data()); - AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data()); - - - - fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); - fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); - AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data()); - + if(fJetTypes&kJet){ + fTCAJetsOut = new TClonesArray("AliAODJet", 0); + fTCAJetsOut->SetName(fNonStdBranch.Data()); + AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data()); + } + + if(fJetTypes&kJetRan){ + fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); + fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) { + if( fEfficiencyFixed < 1.) + fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.))); + else + fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); + } + AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data()); + } + if(fUseBackgroundCalc){ if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){ fAODJetBackgroundOut = new AliAODJetEventBackground(); fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) { + if( fEfficiencyFixed < 1.) + fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.))); + else + fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); + } AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data()); } } // create the branch for the random cones with the same R TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone); - + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) { + if( fEfficiencyFixed < 1.) + cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone); + else + cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone); + } if(fNRandomCones>0){ - if(!AODEvent()->FindListObject(cName.Data())){ - fTCARandomConesOut = new TClonesArray("AliAODJet", 0); - fTCARandomConesOut->SetName(cName.Data()); - AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data()); + if(fJetTypes&kRC){ + if(!AODEvent()->FindListObject(cName.Data())){ + fTCARandomConesOut = new TClonesArray("AliAODJet", 0); + fTCARandomConesOut->SetName(cName.Data()); + AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data()); + } } // create the branch with the random for the random cones on the random event - cName = Form("%sRandomCone_random",fNonStdBranch.Data()); - if(!AODEvent()->FindListObject(cName.Data())){ - fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0); - fTCARandomConesOutRan->SetName(cName.Data()); - AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data()); + if(fJetTypes&kRCRan){ + cName = Form("%sRandomCone_random",fNonStdBranch.Data()); + if(!AODEvent()->FindListObject(cName.Data())){ + fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0); + fTCARandomConesOutRan->SetName(cName.Data()); + AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data()); + } } } @@ -498,35 +655,109 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() nBinEta,binLimitsEta,nBinPt,binLimitsPt); fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}", - nBinEta,binLimitsEta,nBinPt,binLimitsPt); + nBinEta,binLimitsEta,nBinPt,binLimitsPt); fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}", - nBinEta,binLimitsEta,nBinPt,binLimitsPt); + nBinEta,binLimitsEta,nBinPt,binLimitsPt); fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta", - nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); + nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta", nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + + if(fStoreRhoLeadingTrackCorr) { + fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.); + fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.); + fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.); + fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.); + + + fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.); + fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.); + fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.); + fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.); + + fh3CentvsSigmaLeadingTrackPtQ1 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ1","centrality vs background density Q1; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.); + fh3CentvsSigmaLeadingTrackPtQ2 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ2","centrality vs background density Q2; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.); + fh3CentvsSigmaLeadingTrackPtQ3 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ3","centrality vs background density Q3; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.); + fh3CentvsSigmaLeadingTrackPtQ4 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ4","centrality vs background density Q4; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.); + + fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.); + fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.); + fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.); + fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.); + + fh3MultvsSigmaLeadingTrackPtQ1 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.); + fh3MultvsSigmaLeadingTrackPtQ2 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.); + fh3MultvsSigmaLeadingTrackPtQ3 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.); + fh3MultvsSigmaLeadingTrackPtQ4 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.); + + + fh3CentvsDeltaRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + fh3CentvsDeltaRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + fh3CentvsDeltaRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + fh3CentvsDeltaRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + + fHistList->Add(fh3CentvsRhoLeadingTrackPt); + fHistList->Add(fh3CentvsSigmaLeadingTrackPt); + fHistList->Add(fh3MultvsRhoLeadingTrackPt); + fHistList->Add(fh3MultvsSigmaLeadingTrackPt); + + fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1); + fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2); + fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3); + fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4); + + fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1); + fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2); + fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3); + fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4); + + fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1); + fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2); + fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3); + fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4); + + fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1); + fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2); + fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3); + fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4); + + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4); + } + + //Detector level effects histos + fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt); + + fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt); + fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt); + + fHistList->Add(fh2PtGenPtSmeared); + fHistList->Add(fp1Efficiency); + fHistList->Add(fp1PtResolution); if(fNRandomCones>0&&fUseBackgroundCalc){ for(int i = 0;i<3;i++){ @@ -574,12 +805,12 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() fHistList->Add(fh1BiARandomConesRan[i]); } } - for(int i = 0;i < kMaxCent;i++){ - fHistList->Add(fh2JetsLeadingPhiPtC[i]); - fHistList->Add(fh2JetsLeadingPhiPtWC[i]); - fHistList->Add(fh2TracksLeadingJetPhiPtC[i]); - fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]); - } + for(int i = 0;i < kMaxCent;i++){ + fHistList->Add(fh2JetsLeadingPhiPtC[i]); + fHistList->Add(fh2JetsLeadingPhiPtWC[i]); + fHistList->Add(fh2TracksLeadingJetPhiPtC[i]); + fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]); + } fHistList->Add(fh2NRecJetsPt); fHistList->Add(fh2NRecTracksPt); @@ -623,7 +854,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() TH1::AddDirectory(oldStatus); } -void AliAnalysisTaskJetCluster::Init() +void AliAnalysisTaskJetCluster::LocalInit() { // // Initialization @@ -631,19 +862,16 @@ void AliAnalysisTaskJetCluster::Init() if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n"); -} + if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB(); + if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB(); -void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) -{ - if(fUseGlobalSelection){ - // no selection by the service task, we continue - if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__); - PostData(1, fHistList); - return; - } + FitMomentumResolution(); +} +void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) +{ // handle and reset the output jet branch @@ -669,7 +897,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput); return; } - // fethc the header + // fetch the header } else{ // assume that the AOD is in the general output... @@ -682,7 +910,17 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fESD = dynamic_cast (InputEvent()); } } - + + //Check if information is provided detector level effects + if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE; + if( fEfficiencyFixed < 1. ) { + if (!fUseDiceEfficiency) + fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user + } + else { + if(!fhEffH1 || !fhEffH2 || !fhEffH3 ) fUseDiceEfficiency = kFALSE; + } + Bool_t selectEvent = false; Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB; @@ -704,24 +942,49 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh1ZPhySel->Fill(zVtx); } - - if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){ + if(fEventSelection){ + if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){ Float_t yvtx = vtxAOD->GetY(); Float_t xvtx = vtxAOD->GetX(); Float_t r2 = yvtx*yvtx+xvtx*xvtx; - if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on + if(TMath::Abs(zVtx)0){ - if(centfCentCutUp){ - selectEvent = false; } + if(fCentCutUp>0){ + if(centfCentCutUp){ + selectEvent = false; + } + } + }else{ + selectEvent = true; } + } + + Bool_t T0 = false; + Bool_t V0 = false; + const AliAODVZERO *vzero = fAOD->GetVZEROData(); + if(vzero){ + if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){ + V0 = true; + } + } + + const AliAODTZERO *tzero = fAOD->GetTZEROData(); + if(tzero){ + if(TMath::Abs(tzero->GetT0VertexRaw())<100){ + T0 = true; + } } + + if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0; + else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0; + else if(fRequireVZEROAC)selectEvent = selectEvent&&V0; + + if(!selectEvent){ PostData(1, fHistList); return; @@ -744,6 +1007,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec); Float_t nCh = recParticles.GetEntries(); + Float_t nGen=genParticles.GetEntries(); fh1Nch->Fill(nCh); if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries()); nT = GetListOfTracks(&genParticles,fTrackTypeGen); @@ -759,11 +1023,130 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) AliAODJet vTmpRan(1,0,0,1); for(int i = 0; i < recParticles.GetEntries(); i++){ AliVParticle *vp = (AliVParticle*)recParticles.At(i); + // Carefull energy is not well determined in real data, should not matter for p_T scheme? // we take total momentum here - fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P()); - jInp.set_user_index(i); - inputParticlesRec.push_back(jInp); + + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) { + //Add particles to fastjet in case we are not running toy model + fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P()); + jInp.set_user_index(i); + inputParticlesRec.push_back(jInp); + } + else if(fUseDiceEfficiency) { + + // Dice to decide if particle is kept or not - toy model for efficiency + // + Double_t sumEff = 0.; + Double_t pT = 0.; + Double_t eff[3] = {0.}; + Int_t cat[3] = {0}; + + Double_t rnd = fRandom->Uniform(1.); + if( fEfficiencyFixed<1. ) { + sumEff = fEfficiencyFixed; + if (fUseDiceEfficiency == 2) { + sumEff = (nCh - fEfficiencyFixed*nGen) / nCh; // rescale eff; fEfficiencyFixed is wrt to nGen, but dicing is fraction of nCh + } + } else { + + pT = vp->Pt(); + Double_t pTtmp = pT; + if(pT>10.) pTtmp = 10.; + Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp)); + Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp)); + Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp)); + + //Sort efficiencies from large to small + if(eff1>eff2 && eff1>eff3) { + eff[0] = eff1; + cat[0] = 1; + if(eff2>eff3) { + eff[1] = eff2; + eff[2] = eff3; + cat[1] = 2; + cat[2] = 3; + } else { + eff[1] = eff3; + eff[2] = eff2; + cat[1] = 3; + cat[2] = 2; + } + } + else if(eff2>eff1 && eff2>eff3) { + eff[0] = eff2; + cat[0] = 2; + if(eff1>eff3) { + eff[1] = eff1; + eff[2] = eff3; + cat[1] = 1; + cat[2] = 3; + } else { + eff[1] = eff3; + eff[2] = eff1; + cat[1] = 3; + cat[2] = 1; + } + } + else if(eff3>eff1 && eff3>eff2) { + eff[0] = eff3; + cat[0] = 3; + if(eff1>eff2) { + eff[1] = eff1; + eff[2] = eff2; + cat[1] = 1; + cat[2] = 2; + } else { + eff[1] = eff2; + eff[2] = eff1; + cat[1] = 2; + cat[2] = 1; + } + } + + sumEff = eff[0]+eff[1]+eff[2]; + } + fp1Efficiency->Fill(vp->Pt(),sumEff); + if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue; + + if(fUseTrPtResolutionSmearing) { + //Smear momentum of generated particle + Double_t smear = 1.; + //Select hybrid track category + if(rnd<=eff[2]) + smear = GetMomentumSmearing(cat[2],pT); + else if(rnd<=(eff[2]+eff[1])) + smear = GetMomentumSmearing(cat[1],pT); + else + smear = GetMomentumSmearing(cat[0],pT); + + fp1PtResolution->Fill(vp->Pt(),smear); + + Double_t sigma = vp->Pt()*smear; + Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma); + + Double_t phi = vp->Phi(); + Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta()))); + Double_t pX = pTrec * TMath::Cos(phi); + Double_t pY = pTrec * TMath::Sin(phi); + Double_t pZ = pTrec/TMath::Tan(theta); + Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ); + + fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec); + + fastjet::PseudoJet jInp(pX,pY,pZ,p); + jInp.set_user_index(i); + inputParticlesRec.push_back(jInp); + + } + else { + fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P()); + jInp.set_user_index(i); + inputParticlesRec.push_back(jInp); + + } + + } // the randomized input changes eta and phi, but keeps the p_T if(i>=fNSkipLeadingRan){// eventually skip the leading particles @@ -788,9 +1171,11 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(fTCAJetsOut){ if(i == 0){ fRef->Delete(); // make sure to delete before placement new... - new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) { + new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ... + } } - fRef->Add(vp); + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ... } }// recparticles @@ -805,13 +1190,15 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // now create the object that holds info about ghosts - if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){ + /* + if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){ // reduce CPU time... fGhostArea = 0.5; fActiveAreaRepeats = 0; - } - - fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea); + } + */ + + fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea); fastjet::AreaType areaType = fastjet::active_area; fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec); fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy); @@ -832,7 +1219,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh1NJetsRec->Fill(sortedJets.size()); - // loop over all jets an fill information, first on is the leading jet + // loop over all jets an fill information, first on is the leading jet Int_t nRecOver = inclusiveJets.size(); Int_t nRec = inclusiveJets.size(); @@ -864,8 +1251,8 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Float_t pTback = 0; if(externalBackground){ // carefull has to be filled in a task before - // todo, ReArrange to the botom - pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged(); + // todo, ReArrange to the botom + pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged(); } pt = leadingJet.Pt() - pTback; // correlation of leading jet with tracks @@ -911,7 +1298,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Float_t tmpPtBack = 0; if(externalBackground){ // carefull has to be filled in a task before - // todo, ReArrange to the botom + // todo, ReArrange to the botom tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged(); } tmpPt = tmpPt - tmpPtBack; @@ -927,200 +1314,329 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh2NConstPt->Fill(tmpPt,constituents.size()); // loop over constiutents and fill spectrum + AliVParticle *partLead = 0; + Float_t ptLead = -1; + for(unsigned int ic = 0; ic < constituents.size();ic++){ AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); + if(!part) continue; fh1PtJetConstRec->Fill(part->Pt()); if(aodOutJet){ - aodOutJet->AddTrack(fRef->At(constituents[ic].user_index())); + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index())); + if(part->Pt()>fMaxTrackPtInJet){ + aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered); + } + } + if(part->Pt()>ptLead){ + ptLead = part->Pt(); + partLead = part; } if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt()); } - - // correlation - Float_t tmpPhi = tmpRec.Phi(); - Float_t tmpEta = tmpRec.Eta(); - if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; - if(j==0){ - fh1PtJetsLeadingRecIn->Fill(tmpPt); - fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta); - fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt); - fh1NConstLeadingRec->Fill(constituents.size()); - fh2NConstLeadingPt->Fill(tmpPt,constituents.size()); - continue; - } - fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta); - fh2JetEtaPt->Fill(tmpEta,tmpPt); - Float_t dPhi = phi - tmpPhi; - if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); - if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); - Float_t dEta = eta - tmpRec.Eta(); - fh2JetsLeadingPhiEta->Fill(dPhi,dEta); - fh2JetsLeadingPhiPt->Fill(dPhi,pt); - if(cenClass>=0){ - fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt); - fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt); - } - fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt); + + AliAODTrack *aodT = 0; + if(partLead){ + if(aodOutJet){ + //set pT of leading constituent of jet + aodOutJet->SetPtLeading(partLead->Pt()); + aodT = dynamic_cast(partLead); + if(aodT){ + if(aodT->TestFilterBit(fFilterMaskBestPt)){ + aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest); + } + } + } + } + + // correlation + Float_t tmpPhi = tmpRec.Phi(); + Float_t tmpEta = tmpRec.Eta(); + if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; + if(j==0){ + fh1PtJetsLeadingRecIn->Fill(tmpPt); + fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta); + fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt); + fh1NConstLeadingRec->Fill(constituents.size()); + fh2NConstLeadingPt->Fill(tmpPt,constituents.size()); + continue; + } + fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta); + fh2JetEtaPt->Fill(tmpEta,tmpPt); + Float_t dPhi = phi - tmpPhi; + if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); + if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); + Float_t dEta = eta - tmpRec.Eta(); + fh2JetsLeadingPhiEta->Fill(dPhi,dEta); + fh2JetsLeadingPhiPt->Fill(dPhi,pt); + if(cenClass>=0){ + fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt); + fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt); + } + fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt); }// loop over reconstructed jets - delete recIter; - - - - // Add the random cones... - if(fNRandomCones>0&&fTCARandomConesOut){ - // create a random jet within the acceptance - Double_t etaMax = fTrackEtaWindow - fRparam; - Int_t nCone = 0; - Int_t nConeRan = 0; - Double_t pTC = 1; // small number - for(int ir = 0;ir < fNRandomCones;ir++){ - Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax - Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi - // massless jet - Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); - Double_t pZC = pTC/TMath::Tan(thetaC); - Double_t pXC = pTC * TMath::Cos(phiC); - Double_t pYC = pTC * TMath::Sin(phiC); - Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); - AliAODJet tmpRecC (pXC,pYC,pZC, pC); - bool skip = false; - for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets - AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E()); - if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){ - skip = true; - break; - } - } - // test for overlap with previous cones to avoid double counting - for(int iic = 0;iicAt(iic); - if(iicone){ - if(iicone->DeltaR(&tmpRecC)<2.*fRparam){ - skip = true; - break; - } - } - } - if(skip)continue; - tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below - if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC); - if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC); - }// loop over random cones creation + delete recIter; + + + + // Add the random cones... + if(fNRandomCones>0&&fTCARandomConesOut){ + // create a random jet within the acceptance + Double_t etaMax = fTrackEtaWindow - fRparam; + Int_t nCone = 0; + Int_t nConeRan = 0; + Double_t pTC = 1; // small number + for(int ir = 0;ir < fNRandomCones;ir++){ + Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax + Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi + // massless jet + Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); + Double_t pZC = pTC/TMath::Tan(thetaC); + Double_t pXC = pTC * TMath::Cos(phiC); + Double_t pYC = pTC * TMath::Sin(phiC); + Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); + AliAODJet tmpRecC (pXC,pYC,pZC, pC); + bool skip = false; + for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets + AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E()); + if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){ + skip = true; + break; + } + } + // test for overlap with previous cones to avoid double counting + for(int iic = 0;iicAt(iic); + if(iicone){ + if(iicone->DeltaR(&tmpRecC)<2.*fRparam){ + skip = true; + break; + } + } + } + if(skip)continue; + tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below + tmpRecC.SetPtLeading(-1.); + if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC); + if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC); + }// loop over random cones creation - // loop over the reconstructed particles and add up the pT in the random cones - // maybe better to loop over randomized particles not in the real jets... - // but this by definition brings dow average energy in the whole event - AliAODJet vTmpRanR(1,0,0,1); - for(int i = 0; i < recParticles.GetEntries(); i++){ - AliVParticle *vp = (AliVParticle*)recParticles.At(i); - if(fTCARandomConesOut){ - for(int ir = 0;ir < fNRandomCones;ir++){ - AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir); - if(jC&&jC->DeltaR(vp)SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0); - } - } - }// add up energy in cone - - // the randomized input changes eta and phi, but keeps the p_T - if(i>=fNSkipLeadingRan){// eventually skip the leading particles - Double_t pTR = vp->Pt(); - Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow; - Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm(); + // loop over the reconstructed particles and add up the pT in the random cones + // maybe better to loop over randomized particles not in the real jets... + // but this by definition brings dow average energy in the whole event + AliAODJet vTmpRanR(1,0,0,1); + for(int i = 0; i < recParticles.GetEntries(); i++){ + AliVParticle *vp = (AliVParticle*)recParticles.At(i); + if(fTCARandomConesOut){ + for(int ir = 0;ir < fNRandomCones;ir++){ + AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir); + if(jC&&jC->DeltaR(vp)Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered); + jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0); + if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt()); + } + } + }// add up energy in cone + + // the randomized input changes eta and phi, but keeps the p_T + if(i>=fNSkipLeadingRan){// eventually skip the leading particles + Double_t pTR = vp->Pt(); + Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow; + Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm(); - Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR)); - Double_t pZR = pTR/TMath::Tan(thetaR); + Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR)); + Double_t pZR = pTR/TMath::Tan(thetaR); - Double_t pXR = pTR * TMath::Cos(phiR); - Double_t pYR = pTR * TMath::Sin(phiR); - Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR); - vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR); - if(fTCARandomConesOutRan){ - for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ - AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir); - if(jC&&jC->DeltaR(&vTmpRanR)SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0); - } - } - } - } - }// loop over recparticles + Double_t pXR = pTR * TMath::Cos(phiR); + Double_t pYR = pTR * TMath::Sin(phiR); + Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR); + vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR); + if(fTCARandomConesOutRan){ + for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ + AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir); + if(jC&&jC->DeltaR(&vTmpRanR)fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered); + jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0); + if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt()); + } + } + } + } + }// loop over recparticles - Float_t jetArea = fRparam*fRparam*TMath::Pi(); - if(fTCARandomConesOut){ - for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){ - // rescale the momntum vectors for the random cones + Float_t jetArea = fRparam*fRparam*TMath::Pi(); + if(fTCARandomConesOut){ + for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){ + // rescale the momentum vectors for the random cones - AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir); - if(rC){ - Double_t etaC = rC->Eta(); - Double_t phiC = rC->Phi(); - // massless jet, unit vector - pTC = rC->ChargedBgEnergy(); - if(pTC<=0)pTC = 0.001; // for almost empty events - Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); - Double_t pZC = pTC/TMath::Tan(thetaC); - Double_t pXC = pTC * TMath::Cos(phiC); - Double_t pYC = pTC * TMath::Sin(phiC); - Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); - rC->SetPxPyPzE(pXC,pYC,pZC, pC); - rC->SetBgEnergy(0,0); - rC->SetEffArea(jetArea,0); - } - } - } - if(fTCARandomConesOutRan){ - for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ - AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir); - // same wit random - if(rC){ - Double_t etaC = rC->Eta(); - Double_t phiC = rC->Phi(); - // massless jet, unit vector - pTC = rC->ChargedBgEnergy(); - if(pTC<=0)pTC = 0.001;// for almost empty events - Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); - Double_t pZC = pTC/TMath::Tan(thetaC); - Double_t pXC = pTC * TMath::Cos(phiC); - Double_t pYC = pTC * TMath::Sin(phiC); - Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); - rC->SetPxPyPzE(pXC,pYC,pZC, pC); - rC->SetBgEnergy(0,0); - rC->SetEffArea(jetArea,0); - } - } - } - }// if(fNRandomCones + AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir); + if(rC){ + Double_t etaC = rC->Eta(); + Double_t phiC = rC->Phi(); + // massless jet, unit vector + pTC = rC->ChargedBgEnergy(); + if(pTC<=0)pTC = 0.001; // for almost empty events + Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); + Double_t pZC = pTC/TMath::Tan(thetaC); + Double_t pXC = pTC * TMath::Cos(phiC); + Double_t pYC = pTC * TMath::Sin(phiC); + Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); + rC->SetPxPyPzE(pXC,pYC,pZC, pC); + rC->SetBgEnergy(0,0); + rC->SetEffArea(jetArea,0); + } + } + } + if(fTCARandomConesOutRan){ + for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){ + AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir); + // same wit random + if(rC){ + Double_t etaC = rC->Eta(); + Double_t phiC = rC->Phi(); + // massless jet, unit vector + pTC = rC->ChargedBgEnergy(); + if(pTC<=0)pTC = 0.001;// for almost empty events + Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC)); + Double_t pZC = pTC/TMath::Tan(thetaC); + Double_t pXC = pTC * TMath::Cos(phiC); + Double_t pYC = pTC * TMath::Sin(phiC); + Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC); + rC->SetPxPyPzE(pXC,pYC,pZC, pC); + rC->SetBgEnergy(0,0); + rC->SetEffArea(jetArea,0); + } + } + } + }// if(fNRandomCones - //background estimates:all bckg jets(0) & wo the 2 hardest(1) - + //background estimates:all bckg jets(0) & wo the 2 hardest(1) + if(fAODJetBackgroundOut){ + vector jets2=sortedJets; + if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); + Double_t bkg1=0; + Double_t sigma1=0.; + Double_t meanarea1=0.; + Double_t bkg2=0; + Double_t sigma2=0.; + Double_t meanarea2=0.; + Double_t bkg4=0; + Double_t sigma4=0.; + Double_t meanarea4=0.; - if(fAODJetBackgroundOut){ - vector jets2=sortedJets; - if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); - Double_t bkg1=0; - Double_t sigma1=0.; - Double_t meanarea1=0.; - Double_t bkg2=0; - Double_t sigma2=0.; - Double_t meanarea2=0.; + clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true); + fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1); - clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true); - fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1); + // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone)); + // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone)); + clustSeq.get_median_rho_and_sigma(sortedJets, range, true, bkg4, sigma4, meanarea4, true); + fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4); - // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone)); - // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone)); + ////////////////////// - clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true); - fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2); - // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone)); - // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone)); + clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true); + fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2); + // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone)); + // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone)); + + } + } + + if(fStoreRhoLeadingTrackCorr) { + vector jets3=sortedJets; + if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2); + + Double_t bkg2=0; + Double_t sigma2=0.; + Double_t meanarea2=0.; + + clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true); + fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2); + + //Get leading track in event + AliVParticle *leading = (AliVParticle*)recParticles.At(0); + + fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt()); + fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt()); + fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt()); + fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt()); + + + //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event + //exclude 2 leading jets from event + //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track) + //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track) + //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track) + //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track) + + //Assuming R=0.2 for background clusters + + Double_t bkg2Q[4] = {0.}; + Double_t sigma2Q[4] = {0.}; + Double_t meanarea2Q[4] = {0.}; + + //Get phi of leading track in event + Float_t phiLeadingTrack = leading->Phi(); + Float_t phiStep = (TMath::Pi()/2./2. - 0.2); + + //Quadrant1 - near side + phiMin = phiLeadingTrack - phiStep; + phiMax = phiLeadingTrack + phiStep; + fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax); + clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true); + + //Quadrant2 - orthogonal + Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.; + if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi(); + phiMin = phiQ2 - phiStep; + phiMax = phiQ2 + phiStep; + fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax); + clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true); + + //Quadrant3 - away side + Double_t phiQ3 = phiLeadingTrack + TMath::Pi(); + if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi(); + phiMin = phiQ3 - phiStep; + phiMax = phiQ3 + phiStep; + fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax); + clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true); + + //Quadrant4 - othogonal + Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi(); + if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi(); + phiMin = phiQ4 - phiStep; + phiMax = phiQ4 + phiStep; + fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax); + clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true); + + fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt()); + fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt()); + fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt()); + fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt()); + + fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt()); + fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt()); + fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt()); + fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt()); + + fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt()); + fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt()); + fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt()); + fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt()); + + fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt()); + fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt()); + fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt()); + fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt()); + + fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt()); + fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt()); + fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt()); + fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt()); - } } @@ -1132,79 +1648,79 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // do the same for tracks and jets if(nTrackOver>0){ - TIterator *recIter = recParticles.MakeIterator(); - AliVParticle *tmpRec = (AliVParticle*)(recIter->Next()); - Float_t pt = tmpRec->Pt(); - - // Printf("Leading track p_t %3.3E",pt); - for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){ - Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i); - while(ptNext()); - if(tmpRec){ - pt = tmpRec->Pt(); - } - } - if(nTrackOver<=0)break; - fh2NRecTracksPt->Fill(ptCut,nTrackOver); - } + TIterator *recIter = recParticles.MakeIterator(); + AliVParticle *tmpRec = (AliVParticle*)(recIter->Next()); + Float_t pt = tmpRec->Pt(); + + // Printf("Leading track p_t %3.3E",pt); + for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){ + Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i); + while(ptNext()); + if(tmpRec){ + pt = tmpRec->Pt(); + } + } + if(nTrackOver<=0)break; + fh2NRecTracksPt->Fill(ptCut,nTrackOver); + } - recIter->Reset(); - AliVParticle *leading = (AliVParticle*)recParticles.At(0); - Float_t phi = leading->Phi(); - if(phi<0)phi+=TMath::Pi()*2.; - Float_t eta = leading->Eta(); - pt = leading->Pt(); - while((tmpRec = (AliVParticle*)(recIter->Next()))){ - Float_t tmpPt = tmpRec->Pt(); - Float_t tmpEta = tmpRec->Eta(); - fh1PtTracksRecIn->Fill(tmpPt); - fh2TrackEtaPt->Fill(tmpEta,tmpPt); - if(tmpRec==leading){ - fh1PtTracksLeadingRecIn->Fill(tmpPt); - fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt); - continue; - } + recIter->Reset(); + AliVParticle *leading = (AliVParticle*)recParticles.At(0); + Float_t phi = leading->Phi(); + if(phi<0)phi+=TMath::Pi()*2.; + Float_t eta = leading->Eta(); + pt = leading->Pt(); + while((tmpRec = (AliVParticle*)(recIter->Next()))){ + Float_t tmpPt = tmpRec->Pt(); + Float_t tmpEta = tmpRec->Eta(); + fh1PtTracksRecIn->Fill(tmpPt); + fh2TrackEtaPt->Fill(tmpEta,tmpPt); + if(tmpRec==leading){ + fh1PtTracksLeadingRecIn->Fill(tmpPt); + fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt); + continue; + } // correlation - Float_t tmpPhi = tmpRec->Phi(); + Float_t tmpPhi = tmpRec->Phi(); - if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; - Float_t dPhi = phi - tmpPhi; - if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); - if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); - Float_t dEta = eta - tmpRec->Eta(); - fh2TracksLeadingPhiEta->Fill(dPhi,dEta); - fh2TracksLeadingPhiPt->Fill(dPhi,pt); - fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt); - } - delete recIter; - } + if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; + Float_t dPhi = phi - tmpPhi; + if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); + if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); + Float_t dEta = eta - tmpRec->Eta(); + fh2TracksLeadingPhiEta->Fill(dPhi,dEta); + fh2TracksLeadingPhiPt->Fill(dPhi,pt); + fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt); + } + delete recIter; + } - // find the random jets + // find the random jets - fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef); + fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef); - // fill the jet information from random track - const vector &inclusiveJetsRan = clustSeqRan.inclusive_jets(); - const vector &sortedJetsRan = sorted_by_pt(inclusiveJetsRan); + // fill the jet information from random track + const vector &inclusiveJetsRan = clustSeqRan.inclusive_jets(); + const vector &sortedJetsRan = sorted_by_pt(inclusiveJetsRan); fh1NJetsRecRan->Fill(sortedJetsRan.size()); - // loop over all jets an fill information, first on is the leading jet + // loop over all jets an fill information, first on is the leading jet - Int_t nRecOverRan = inclusiveJetsRan.size(); - Int_t nRecRan = inclusiveJetsRan.size(); + Int_t nRecOverRan = inclusiveJetsRan.size(); + Int_t nRecRan = inclusiveJetsRan.size(); - if(inclusiveJetsRan.size()>0){ - AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E()); - Float_t pt = leadingJet.Pt(); + if(inclusiveJetsRan.size()>0){ + AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E()); + Float_t pt = leadingJet.Pt(); - Int_t iCount = 0; - TLorentzVector vecarearanb; + Int_t iCount = 0; + TLorentzVector vecarearanb; - for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){ - Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i); + for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){ + Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i); while(ptFill(nCh,tmpPt,constituents.size()); - if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){ - aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec); - Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); - aodOutJetRan->GetRefTracks()->Clear(); - aodOutJetRan->SetEffArea(arearan,0); - fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]); - vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e()); - aodOutJetRan->SetVectorAreaCharged(&vecarearanb); + if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){ + aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec); + Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); + aodOutJetRan->GetRefTracks()->Clear(); + aodOutJetRan->SetEffArea(arearan,0); + fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]); + vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e()); + aodOutJetRan->SetVectorAreaCharged(&vecarearanb); - } + } // correlation Float_t tmpPhi = tmpRec.Phi(); @@ -1273,79 +1789,90 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(fAODJetBackgroundOut){ - Double_t bkg3=0.; - Double_t sigma3=0.; - Double_t meanarea3=0.; - clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true); - fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3); - // float areaRandomCone = rRandomCone2 *TMath::Pi(); - /* - fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone)); - fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone)); - */ + Double_t bkg3=0.; + Double_t sigma3=0.; + Double_t meanarea3=0.; + clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true); + fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3); + // float areaRandomCone = rRandomCone2 *TMath::Pi(); + /* + fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone)); + fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone)); + */ } - } + } - // do the event selection if activated - if(fJetTriggerPtCut>0){ - bool select = false; - Float_t minPt = fJetTriggerPtCut; - /* - // hard coded for now ... - // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3 - if(cent<10)minPt = 50; - else if(cent<30)minPt = 42; - else if(cent<50)minPt = 28; - else if(cent<80)minPt = 18; - */ - float rho = 0; - if(externalBackground)rho = externalBackground->GetBackground(2); - if(fTCAJetsOut){ - for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){ - AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i); - Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged(); - if(ptSub>=minPt){ - select = true; - break; - } - } - } + // do the event selection if activated + if(fJetTriggerPtCut>0){ + bool select = false; + Float_t minPt = fJetTriggerPtCut; + /* + // hard coded for now ... + // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3 + if(cent<10)minPt = 50; + else if(cent<30)minPt = 42; + else if(cent<50)minPt = 28; + else if(cent<80)minPt = 18; + */ + float rho = 0; + if(externalBackground)rho = externalBackground->GetBackground(2); + if(fTCAJetsOut){ + for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){ + AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i); + Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged(); + if(ptSub>=minPt){ + select = true; + break; + } + } + } - if(select){ - static AliAODHandler *aodH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); - fh1CentralitySelect->Fill(cent); - fh1ZSelect->Fill(zVtx); - aodH->SetFillAOD(kTRUE); - } - } - if (fDebug > 2){ - if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast()); - if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast()); - if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast()); - if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast()); - } - PostData(1, fHistList); + if(select){ + static AliAODHandler *aodH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); + fh1CentralitySelect->Fill(cent); + fh1ZSelect->Fill(zVtx); + aodH->SetFillAOD(kTRUE); + } + } + if (fDebug > 2){ + if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast()); + if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast()); + if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast()); + if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast()); + } + PostData(1, fHistList); } void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/) { -// Terminate analysis -// - if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n"); + // + // Terminate analysis + // + if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n"); + + if(fMomResH1Fit) delete fMomResH1Fit; + if(fMomResH2Fit) delete fMomResH2Fit; + if(fMomResH3Fit) delete fMomResH3Fit; + } Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ - if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type); + // + // get list of tracks/particles for different types + // + + if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type); Int_t iCount = 0; - if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){ - if(type!=kTrackAODextraonly) { + if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly || type==kTrackAODMCextra || type==kTrackAODMCextraonly){ + + if(type!=kTrackAODextraonly && type!=kTrackAODMCextraonly) { AliAODEvent *aod = 0; if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); else aod = AODEvent(); @@ -1353,12 +1880,36 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__); return iCount; } + for(int it = 0;it < aod->GetNumberOfTracks();++it){ AliAODTrack *tr = aod->GetTrack(it); - if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))){ + Bool_t bGood = false; + if(fFilterType == 0)bGood = true; + else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal(); + else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal(); + if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){ if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap()); continue; } + + // heavy flavor jets + if(fFilterMask==528 && fUseHFcuts){ + Double_t ntpcClus = tr->GetTPCNcls(); + Double_t trPt=tr->Pt(); + TFormula NTPCClsCut("f1NClustersTPCLinearPtDep","70.+30./20.*x"); + + if (trPt <= 20. && (ntpcClus < NTPCClsCut.Eval(trPt))) continue; + else if (trPt > 20. && ntpcClus < 100) continue; + + if (AvoidDoubleCountingHF(aod,tr)) continue; + } + // + + if(fRequireITSRefit){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;} + if (fApplySharedClusterCut) { + Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls()); + if (frac > 0.4) continue; + } if(TMath::Abs(tr->Eta())>fTrackEtaWindow){ if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks()); continue; @@ -1388,11 +1939,72 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ AliAODTrack *trackAOD = dynamic_cast (track); if(!trackAOD)continue; + Bool_t bGood = false; + if(fFilterType == 0)bGood = true; + else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal(); + else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal(); + if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue; + if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;} + if (fApplySharedClusterCut) { + Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls()); + if (frac > 0.4) continue; + } + + + if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue; if(trackAOD->Pt()Pt()); list->Add(trackAOD); iCount++; } } + + if(type==kTrackAODMCextra || type==kTrackAODMCextraonly) { //embed generator level particles + AliAODEvent *aod = 0; + if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + if(!aod){ + return iCount; + } + TClonesArray *aodExtraTracks = dynamic_cast(aod->FindListObject("aodExtraMCparticles")); + if(!aodExtraTracks)return iCount; + for(int it =0; itGetEntries(); it++) { + AliVParticle *track = dynamic_cast ((*aodExtraTracks)[it]); + AliAODMCParticle *partmc = dynamic_cast ((*aodExtraTracks)[it]); + if (!track) { + if(fDebug>10) printf("track %d does not exist\n",it); + continue; + } + + if(!partmc) continue; + if(!partmc->IsPhysicalPrimary())continue; + + if (track->Pt()10) printf("track %d has too low pt %.2f\n",it,track->Pt()); + continue; + } + + /* + AliAODTrack *trackAOD = dynamic_cast((*aodExtraTracks)[it]);//(track); + + if(!trackAOD) { + if(fDebug>10) printf("trackAOD %d does not exist\n",it); + continue; + } + + trackAOD->SetFlags(AliESDtrack::kEmbedded); + trackAOD->SetFilterMap(fFilterMask); + */ + if(fDebug>10) printf("pt extra track %.2f \n", track->Pt()); + + if(TMath::Abs(track->Eta())>fTrackEtaWindow) continue; + if(track->Pt()Add(track); + + iCount++; + } + } + } else if (type == kTrackKineAll||type == kTrackKineCharged){ AliMCEvent* mcEvent = MCEvent(); @@ -1443,21 +2055,248 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ } } }// AODMCparticle + else if (type == kTrackAODMCHF){ + + AliAODEvent *aod = 0; + if(fUseAODMCInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + if(!aod)return iCount; + TClonesArray *tca = dynamic_cast(aod->FindListObject(AliAODMCParticle::StdBranchName())); + if(!tca)return iCount; + for(int it = 0;it < tca->GetEntriesFast();++it){ + AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it)); + if(!part) continue; + int partpdg= part->PdgCode(); + if(!part->IsPhysicalPrimary() && !IsBMeson(partpdg) && !IsDMeson(partpdg) )continue; + + if (IsBMeson(partpdg) || IsDMeson(partpdg)) { + iCount+= AddDaughters( list , part,tca); + } + else { + + if(part->Pt()Eta())>fTrackEtaWindow) continue; + if(part->Charge()==0) continue; + + if((part->Pt()>=fTrackPtCut) && (TMath::Abs(part->Eta())<=fTrackEtaWindow) && (part->Charge()!=0))list->Add(part); + iCount++; + } + } + } + list->Sort(); return iCount; } +Int_t AliAnalysisTaskJetCluster::AddDaughters(TList * list, AliAODMCParticle *part, TClonesArray * tca){ + Int_t count=0; + Int_t nDaugthers = part->GetNDaughters(); + for(Int_t i=0;i< nDaugthers;++i){ + AliAODMCParticle *partdaughter = (AliAODMCParticle*)(tca->At(i)); + if(!partdaughter) continue; + if(partdaughter->Pt()Eta())>fTrackEtaWindow)continue; + if(partdaughter->Charge()==0)continue; + + if(!IsDMeson(partdaughter->PdgCode()) && !IsBMeson(partdaughter->PdgCode()) ){ + list->Add(partdaughter); + count++; + } + else count+=AddDaughters(list,part,tca); + } +return count; +} + + +Bool_t AliAnalysisTaskJetCluster::AvoidDoubleCountingHF(AliAODEvent *aod, AliAODTrack *tr1){ + + if(!(tr1->TestFilterBit(BIT(9)))) return kFALSE; + + Int_t idtr1 = tr1->GetID(); + + for(int jt = 0;jt < aod->GetNumberOfTracks();++jt){ + + const AliAODTrack *tr2 = aod->GetTrack(jt); + Int_t idtr2 = tr2->GetID(); + + if (!(tr2->TestFilterBit(BIT(4)))) continue; + if (idtr1==(idtr2+1)*-1.) return kTRUE; + + } + return kFALSE; +} +void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() { + + if (!gGrid) { + AliInfo("Trying to connect to AliEn ..."); + TGrid::Connect("alien://"); + } + + TFile *f = TFile::Open(fPathTrPtResolution.Data()); + + if(!f)return; + + TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt"); + TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS"); + TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD"); + + SetSmearResolution(kTRUE); + SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD); + + +} + +void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() { + + if (!gGrid) { + AliInfo("Trying to connect to AliEn ..."); + TGrid::Connect("alien://"); + } + + TFile *f = TFile::Open(fPathTrEfficiency.Data()); + if(!f)return; + + TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt"); + TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS"); + TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD"); + + SetDiceEfficiency(kTRUE); + + if(fChangeEfficiencyFraction>0.) { + + TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin"); + + for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) { + Double_t content = hEffPosGlobSt->GetBinContent(bin); + hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction); + } + + SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD); + + } + else + SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD); + +} + +void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) { + + // + // set mom res profiles + // + + if(fMomResH1) delete fMomResH1; + if(fMomResH2) delete fMomResH2; + if(fMomResH3) delete fMomResH3; + + fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1"); + fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2"); + fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3"); + +} + +void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) { + // + // set tracking efficiency histos + // + + fhEffH1 = (TH1*)h1->Clone("fhEffH1"); + fhEffH2 = (TH1*)h2->Clone("fhEffH2"); + fhEffH3 = (TH1*)h3->Clone("fhEffH3"); +} + +Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) { + + // + // Get smearing on generated momentum + // + + //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt); + + TProfile *fMomRes = 0x0; + if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes"); + if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes"); + if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes"); + + if(!fMomRes) { + return 0.; + } + + + Double_t smear = 0.; + + if(pt>20.) { + if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt); + if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt); + if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt); + } + else { + + Int_t bin = fMomRes->FindBin(pt); + + smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin)); + + } + + if(fMomRes) delete fMomRes; + + return smear; +} + +void AliAnalysisTaskJetCluster::FitMomentumResolution() { + // + // Fit linear function on momentum resolution at high pT + // + + if(!fMomResH1Fit && fMomResH1) { + fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.); + fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.); + fMomResH1Fit ->SetRange(5.,100.); + } + + if(!fMomResH2Fit && fMomResH2) { + fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.); + fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.); + fMomResH2Fit ->SetRange(5.,100.); + } + + if(!fMomResH3Fit && fMomResH3) { + fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.); + fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.); + fMomResH3Fit ->SetRange(5.,100.); + } + +} + /* -Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector &inputParticles){ + Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector &inputParticles){ for(int i = 0; i < particles.GetEntries(); i++){ - AliVParticle *vp = (AliVParticle*)particles.At(i); - // Carefull energy is not well determined in real data, should not matter for p_T scheme? - fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E()); - jInp.set_user_index(i); - inputParticles.push_back(jInp); + AliVParticle *vp = (AliVParticle*)particles.At(i); + // Carefull energy is not well determined in real data, should not matter for p_T scheme? + fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E()); + jInp.set_user_index(i); + inputParticles.push_back(jInp); } return 0; -} + } */ + + +bool AliAnalysisTaskJetCluster::IsBMeson(int pc){ + int bPdG[] = {511,521,10511,10521,513,523,10513,10523,20513,20523,20513,20523,515,525,531, + 10531,533,10533,20533,535,541,10541,543,10543,20543,545,551,10551,100551, + 110551,200551,210551,553,10553,20553,30553,100553,110553,120553,130553,200553,210553,220553, + 300553,9000533,9010553,555,10555,20555,100555,110555,120555,200555,557,100557}; + for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true; +return false; +} +bool AliAnalysisTaskJetCluster::IsDMeson(int pc){ + int bPdG[] = {411,421,10411,10421,413,423,10413,10423,20431,20423,415, + 425,431,10431,433,10433,20433,435,441,10441,100441,443,10443,20443, + 100443,30443,9000443,9010443,9020443,445,100445}; + for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true; +return false; +} +