X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliAnalysisTaskJetCluster.cxx;h=270b6140abd53242c739aaca74e0cefc4ef904a4;hb=4d019bab422f7f1f1fd4c6a1c472ab41a787f9a0;hp=43b75ddbfb954953e76525ad94ef6ad3edd68d7f;hpb=7bedcd31e780ab24608475655d42474c90352d9e;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx index 43b75ddbfb9..270b6140abd 100644 --- a/JETAN/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/AliAnalysisTaskJetCluster.cxx @@ -37,6 +37,7 @@ #include #include #include "TDatabasePDG.h" +#include #include "AliAnalysisTaskJetCluster.h" #include "AliAnalysisManager.h" @@ -103,6 +104,8 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): fUseAODMCInput(kFALSE), fUseBackgroundCalc(kFALSE), fEventSelection(kFALSE), + fRequireVZEROAC(kFALSE), + fRequireTZEROvtx(kFALSE), fFilterMask(0), fFilterMaskBestPt(0), fFilterType(0), @@ -115,6 +118,7 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): fAvgTrials(1), fExternalWeight(1), fTrackEtaWindow(0.9), + fRequireITSRefit(0), fRecEtaWindow(0.5), fTrackPtCut(0.), fJetOutputMinPt(0.150), @@ -139,12 +143,13 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): fhEffH3(0x0), fUseTrPtResolutionSmearing(kFALSE), fUseDiceEfficiency(kFALSE), - fDiceEfficiencyMinPt(0.), + fDiceEfficiencyMinPt(-1.), fUseTrPtResolutionFromOADB(kFALSE), fUseTrEfficiencyFromOADB(kFALSE), fPathTrPtResolution(""), fPathTrEfficiency(""), fChangeEfficiencyFraction(0.), + fEfficiencyFixed(1.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -213,10 +218,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): fh2PtNchNRan(0x0), fh2TracksLeadingJetPhiPtRan(0x0), fh2TracksLeadingJetPhiPtWRan(0x0), - fh2CentvsRho(0x0), - fh2CentvsSigma(0x0), - fh2MultvsRho(0x0), - fh2MultvsSigma(0x0), + fh3CentvsRhoLeadingTrackPt(0x0), + fh3CentvsSigmaLeadingTrackPt(0x0), + fh3MultvsRhoLeadingTrackPt(0x0), + fh3MultvsSigmaLeadingTrackPt(0x0), fh3CentvsRhoLeadingTrackPtQ1(0x0), fh3CentvsRhoLeadingTrackPtQ2(0x0), fh3CentvsRhoLeadingTrackPtQ3(0x0), @@ -233,6 +238,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): fh3MultvsSigmaLeadingTrackPtQ2(0x0), fh3MultvsSigmaLeadingTrackPtQ3(0x0), fh3MultvsSigmaLeadingTrackPtQ4(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0), fh2PtGenPtSmeared(0x0), fp1Efficiency(0x0), fp1PtResolution(0x0), @@ -262,7 +271,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fUseAODTrackInput(kFALSE), fUseAODMCInput(kFALSE), fUseBackgroundCalc(kFALSE), - fEventSelection(kFALSE), fFilterMask(0), + fEventSelection(kFALSE), + fRequireVZEROAC(kFALSE), + fRequireTZEROvtx(kFALSE), + fFilterMask(0), fFilterMaskBestPt(0), fFilterType(0), fJetTypes(kJet), @@ -274,6 +286,7 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fAvgTrials(1), fExternalWeight(1), fTrackEtaWindow(0.9), + fRequireITSRefit(0), fRecEtaWindow(0.5), fTrackPtCut(0.), fJetOutputMinPt(0.150), @@ -298,12 +311,13 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fhEffH3(0x0), fUseTrPtResolutionSmearing(kFALSE), fUseDiceEfficiency(kFALSE), - fDiceEfficiencyMinPt(0.), + fDiceEfficiencyMinPt(-1.), fUseTrPtResolutionFromOADB(kFALSE), fUseTrEfficiencyFromOADB(kFALSE), fPathTrPtResolution(""), fPathTrEfficiency(""), fChangeEfficiencyFraction(0.), + fEfficiencyFixed(1.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -372,10 +386,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fh2PtNchNRan(0x0), fh2TracksLeadingJetPhiPtRan(0x0), fh2TracksLeadingJetPhiPtWRan(0x0), - fh2CentvsRho(0x0), - fh2CentvsSigma(0x0), - fh2MultvsRho(0x0), - fh2MultvsSigma(0x0), + fh3CentvsRhoLeadingTrackPt(0x0), + fh3CentvsSigmaLeadingTrackPt(0x0), + fh3MultvsRhoLeadingTrackPt(0x0), + fh3MultvsSigmaLeadingTrackPt(0x0), fh3CentvsRhoLeadingTrackPtQ1(0x0), fh3CentvsRhoLeadingTrackPtQ2(0x0), fh3CentvsRhoLeadingTrackPtQ3(0x0), @@ -392,6 +406,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fh3MultvsSigmaLeadingTrackPtQ2(0x0), fh3MultvsSigmaLeadingTrackPtQ3(0x0), fh3MultvsSigmaLeadingTrackPtQ4(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0), fh2PtGenPtSmeared(0x0), fp1Efficiency(0x0), fp1PtResolution(0x0), @@ -411,7 +429,7 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fh2TracksLeadingJetPhiPtC[i] = 0; fh2TracksLeadingJetPhiPtWC[i] = 0; } - DefineOutput(1, TList::Class()); + DefineOutput(1, TList::Class()); } @@ -458,8 +476,12 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() if(fJetTypes&kJetRan){ fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); - if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) - fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); + 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()); } @@ -467,17 +489,23 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() 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) - fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); - + 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) - cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),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(fJetTypes&kRC){ if(!AODEvent()->FindListObject(cName.Data())){ @@ -623,40 +651,40 @@ 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) { - fh2CentvsRho = new TH2F("fh2CentvsRho","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.); - fh2CentvsSigma = new TH2F("fh2CentvsSigma","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.); - fh2MultvsRho = new TH2F("fh2MultvsRho","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.); - fh2MultvsSigma = new TH2F("fh2MultvsSigma","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.); + 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.); @@ -679,10 +707,16 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() 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.); - fHistList->Add(fh2CentvsRho); - fHistList->Add(fh2CentvsSigma); - fHistList->Add(fh2MultvsRho); - fHistList->Add(fh2MultvsSigma); + + 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); @@ -704,6 +738,11 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3); fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4); + } //Detector level effects histos @@ -762,12 +801,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); @@ -870,8 +909,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) //Check if information is provided detector level effects if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE; - if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = 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; @@ -913,7 +958,28 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) 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); @@ -937,6 +1003,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); @@ -966,63 +1033,75 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // Dice to decide if particle is kept or not - toy model for efficiency // - Double_t rnd = fRandom->Uniform(1.); - Double_t pT = vp->Pt(); + Double_t sumEff = 0.; + Double_t pT = 0.; Double_t eff[3] = {0.}; - 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)); Int_t cat[3] = {0}; - //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; + + 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(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; + 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]; } - - Double_t sumEff = eff[0]+eff[1]+eff[2]; fp1Efficiency->Fill(vp->Pt(),sumEff); if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue; @@ -1108,11 +1187,11 @@ 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); @@ -1136,7 +1215,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(); @@ -1169,7 +1248,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(externalBackground){ // carefull has to be filled in a task before // todo, ReArrange to the botom - pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged(); + pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged(); } pt = leadingJet.Pt() - pTback; // correlation of leading jet with tracks @@ -1215,7 +1294,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; @@ -1265,274 +1344,295 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) } } - // 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); + // 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 - tmpRecC.SetPtLeading(-1.); - 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)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(); + // 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)fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered); - jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0); - if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt()); - } - } - } - } - }// 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 momentum 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)); - - if(fStoreRhoLeadingTrackCorr) { - fh2CentvsRho->Fill(cent,bkg2); - fh2CentvsSigma->Fill(cent,sigma2); - fh2MultvsRho->Fill(nCh,bkg2); - fh2MultvsSigma->Fill(nCh,sigma2); + ////////////////////// + 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)); + + } + } - //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 - AliVParticle *leading = (AliVParticle*)recParticles.At(0); - 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(jets2, 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(jets2, 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(jets2, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true); - - //Quadrant4 - othogonal - Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi(); - if(phiQ3>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(jets2, 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()); - - } + 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()); - } } @@ -1544,79 +1644,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(); @@ -1685,62 +1785,62 @@ 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*/) @@ -1748,11 +1848,11 @@ void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/) // // Terminate analysis // - if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n"); + if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n"); - if(fMomResH1Fit) delete fMomResH1Fit; - if(fMomResH2Fit) delete fMomResH2Fit; - if(fMomResH3Fit) delete fMomResH3Fit; + if(fMomResH1Fit) delete fMomResH1Fit; + if(fMomResH2Fit) delete fMomResH2Fit; + if(fMomResH3Fit) delete fMomResH3Fit; } @@ -1763,11 +1863,12 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ // get list of tracks/particles for different types // - if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type); + 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(); @@ -1786,6 +1887,7 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap()); continue; } + if(fRequireITSRefit){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)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; @@ -1820,12 +1922,61 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal(); else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal(); if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue; - if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue; + if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)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(); @@ -1876,12 +2027,67 @@ 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; +} + + + void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() { + if (!gGrid) { + AliInfo("Trying to connect to AliEn ..."); + TGrid::Connect("alien://"); + } + TFile *f = TFile::Open(fPathTrPtResolution.Data()); if(!f)return; @@ -1898,6 +2104,11 @@ void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() { void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() { + if (!gGrid) { + AliInfo("Trying to connect to AliEn ..."); + TGrid::Connect("alien://"); + } + TFile *f = TFile::Open(fPathTrEfficiency.Data()); if(!f)return; @@ -2014,16 +2225,34 @@ void AliAnalysisTaskJetCluster::FitMomentumResolution() { } /* -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; +} +