From: loizides Date: Mon, 24 Jun 2013 10:03:40 +0000 (+0000) Subject: from Marta: X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=56f370b088530d9ef2aacb0323bec2bd1ee47343;p=u%2Fmrichter%2FAliRoot.git from Marta: - delta-pT task picks up automatically the number of centrality bins specified + histos for partial exclusion in pA analysis - AliAnalysisTaskEmcalJet+__AliAnalysisTaskEmcalJetDev+__AliJetContainer: adjust jet acceptance based on run number in case of 2012 and 2013 data (new SM not active/masked) --- diff --git a/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx b/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx index 105115e40ed..bb712e67b7b 100644 --- a/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx +++ b/PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx @@ -41,20 +41,68 @@ AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() : fRandCaloClusters(0), fEmbeddedClusterNIds(0), fEmbeddedTrackNIds(0), - fHistRCPhiEta(0), + fHistRCPhiEta(0), + fHistRCPt(0), + fHistRCPtExLJ(0), + fHistRCPtExPartialLJ(0), + fHistRCPtRand(0), + fHistRhoVSRCPt(0), + fHistDeltaPtRC(0), + fHistDeltaPtRCExLJ(0), + fHistDeltaPtRCExPartialLJ(0), + fHistDeltaPtRCRand(0), + fHistEmbNotFoundPt(0), + fHistEmbNotFoundPhiEta(0), + fHistEmbRejectedJetsPhiEta(0), + fHistEmbRejectedJetsPtArea(0), + fHistEmbJetsPtArea(0), + fHistEmbJetsCorrPtArea(0), + fHistEmbPartPtvsJetPt(0), + fHistEmbPartPtvsJetCorrPt(0), + fHistJetPtvsJetCorrPt(0), + fHistDistLeadPart2JetAxis(0), + fHistEmbBkgArea(0), + fHistRhoVSEmbBkg(0), + fHistDeltaPtEmbArea(0), fHistRCPtExLJVSDPhiLJ(0), + fHistRCPtExPartialLJVSDPhiLJ(0), fHistEmbJetsPhiEta(0), fHistLeadPartPhiEta(0) { // Default constructor. - for (Int_t i = 0; i < 4; i++) { + fHistRCPt = new TH1*[fNcentBins]; + fHistRCPtExLJ = new TH1*[fNcentBins]; + fHistRCPtExPartialLJ = new TH1*[fNcentBins]; + fHistRCPtRand = new TH1*[fNcentBins]; + fHistRhoVSRCPt = new TH2*[fNcentBins]; + fHistDeltaPtRC = new TH1*[fNcentBins]; + fHistDeltaPtRCExLJ = new TH1*[fNcentBins]; + fHistDeltaPtRCExPartialLJ = new TH1*[fNcentBins]; + fHistDeltaPtRCRand = new TH1*[fNcentBins]; + fHistEmbNotFoundPt = new TH1*[fNcentBins]; + fHistEmbNotFoundPhiEta = new TH2*[fNcentBins]; + fHistEmbRejectedJetsPhiEta = new TH2*[fNcentBins]; + fHistEmbRejectedJetsPtArea = new TH1*[fNcentBins]; + fHistEmbJetsPtArea = new TH3*[fNcentBins]; + fHistEmbJetsCorrPtArea = new TH3*[fNcentBins]; + fHistEmbPartPtvsJetPt = new TH2*[fNcentBins]; + fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins]; + fHistJetPtvsJetCorrPt = new TH2*[fNcentBins]; + fHistDistLeadPart2JetAxis = new TH1*[fNcentBins]; + fHistEmbBkgArea = new TH2*[fNcentBins]; + fHistRhoVSEmbBkg = new TH2*[fNcentBins]; + fHistDeltaPtEmbArea = new TH2*[fNcentBins]; + + for (Int_t i = 0; i < fNcentBins; i++) { fHistRCPt[i] = 0; fHistRCPtExLJ[i] = 0; - fHistRCPtRand[i] = 0; + fHistRCPtExPartialLJ[i] = 0; + fHistRCPtRand[i] = 0; fHistRhoVSRCPt[i] = 0; fHistDeltaPtRC[i] = 0; fHistDeltaPtRCExLJ[i] = 0; + fHistDeltaPtRCExPartialLJ[i] = 0; fHistDeltaPtRCRand[i] = 0; fHistEmbRejectedJetsPhiEta[i] = 0; fHistEmbRejectedJetsPtArea[i] = 0; @@ -95,20 +143,68 @@ AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) : fRandCaloClusters(0), fEmbeddedClusterNIds(0), fEmbeddedTrackNIds(0), - fHistRCPhiEta(0), + fHistRCPhiEta(0), + fHistRCPt(0), + fHistRCPtExLJ(0), + fHistRCPtExPartialLJ(0), + fHistRCPtRand(0), + fHistRhoVSRCPt(0), + fHistDeltaPtRC(0), + fHistDeltaPtRCExLJ(0), + fHistDeltaPtRCExPartialLJ(0), + fHistDeltaPtRCRand(0), + fHistEmbNotFoundPt(0), + fHistEmbNotFoundPhiEta(0), + fHistEmbRejectedJetsPhiEta(0), + fHistEmbRejectedJetsPtArea(0), + fHistEmbJetsPtArea(0), + fHistEmbJetsCorrPtArea(0), + fHistEmbPartPtvsJetPt(0), + fHistEmbPartPtvsJetCorrPt(0), + fHistJetPtvsJetCorrPt(0), + fHistDistLeadPart2JetAxis(0), + fHistEmbBkgArea(0), + fHistRhoVSEmbBkg(0), + fHistDeltaPtEmbArea(0), fHistRCPtExLJVSDPhiLJ(0), + fHistRCPtExPartialLJVSDPhiLJ(0), fHistEmbJetsPhiEta(0), fHistLeadPartPhiEta(0) { // Standard constructor. - for (Int_t i = 0; i < 4; i++) { + fHistRCPt = new TH1*[fNcentBins]; + fHistRCPtExLJ = new TH1*[fNcentBins]; + fHistRCPtExPartialLJ = new TH1*[fNcentBins]; + fHistRCPtRand = new TH1*[fNcentBins]; + fHistRhoVSRCPt = new TH2*[fNcentBins]; + fHistDeltaPtRC = new TH1*[fNcentBins]; + fHistDeltaPtRCExLJ = new TH1*[fNcentBins]; + fHistDeltaPtRCExPartialLJ = new TH1*[fNcentBins]; + fHistDeltaPtRCRand = new TH1*[fNcentBins]; + fHistEmbRejectedJetsPhiEta = new TH2*[fNcentBins]; + fHistEmbNotFoundPt = new TH1*[fNcentBins]; + fHistEmbNotFoundPhiEta = new TH2*[fNcentBins]; + fHistEmbRejectedJetsPtArea = new TH1*[fNcentBins]; + fHistEmbJetsPtArea = new TH3*[fNcentBins]; + fHistEmbJetsCorrPtArea = new TH3*[fNcentBins]; + fHistEmbPartPtvsJetPt = new TH2*[fNcentBins]; + fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins]; + fHistJetPtvsJetCorrPt = new TH2*[fNcentBins]; + fHistDistLeadPart2JetAxis = new TH1*[fNcentBins]; + fHistEmbBkgArea = new TH2*[fNcentBins]; + fHistRhoVSEmbBkg = new TH2*[fNcentBins]; + fHistDeltaPtEmbArea = new TH2*[fNcentBins]; + + for (Int_t i = 0; i < fNcentBins; i++) { fHistRCPt[i] = 0; fHistRCPtExLJ[i] = 0; - fHistRCPtRand[i] = 0; + fHistRCPtExPartialLJ[i] = 0; + fHistRCPtRand[i] = 0; fHistRhoVSRCPt[i] = 0; fHistDeltaPtRC[i] = 0; fHistDeltaPtRCExLJ[i] = 0; + fHistDeltaPtRCExPartialLJ[i] = 0; fHistDeltaPtRCRand[i] = 0; fHistEmbRejectedJetsPhiEta[i] = 0; fHistEmbRejectedJetsPtArea[i] = 0; @@ -155,6 +251,11 @@ void AliAnalysisTaskDeltaPt::UserCreateOutputObjects() fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi"); fOutput->Add(fHistRCPtExLJVSDPhiLJ); + + fHistRCPtExPartialLJVSDPhiLJ = new TH2F("fHistRCPtExPartialLJVSDPhiLJ","fHistRCPtExPartialLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8); + fHistRCPtExPartialLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); + fHistRCPtExPartialLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi"); + fOutput->Add(fHistRCPtExPartialLJVSDPhiLJ); } } @@ -216,6 +317,20 @@ void AliAnalysisTaskDeltaPt::UserCreateOutputObjects() fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})"); fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtRCExLJ[i]); + + histname = "fHistRCPtExPartialLJ_"; + histname += i; + fHistRCPtExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2); + fHistRCPtExPartialLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})"); + fHistRCPtExPartialLJ[i]->GetYaxis()->SetTitle("counts"); + fOutput->Add(fHistRCPtExPartialLJ[i]); + + histname = "fHistDeltaPtRCExPartialLJ_"; + histname += i; + fHistDeltaPtRCExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtRCExPartialLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})"); + fHistDeltaPtRCExPartialLJ[i]->GetYaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtRCExPartialLJ[i]); } } @@ -397,6 +512,27 @@ Bool_t AliAnalysisTaskDeltaPt::FillHistograms() fHistRCPtExLJ[fCentBin]->Fill(RCpt); fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal); } + + //partial exclusion + if(fBeamType == kpA) { + + RCpt = 0; + RCeta = 0; + RCphi = 0; + GetRandomCone(RCpt, RCeta, RCphi, jet,0,0,kTRUE); + + if (RCpt > 0) { + if (jet) { + Float_t dphi = RCphi - jet->Phi(); + if (dphi > 4.8) dphi -= TMath::Pi() * 2; + if (dphi < -1.6) dphi += TMath::Pi() * 2; + fHistRCPtExPartialLJVSDPhiLJ->Fill(RCpt, dphi); + } + fHistRCPtExPartialLJ[fCentBin]->Fill(RCpt); + fHistDeltaPtRCExPartialLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal); + } + } + } } } @@ -638,7 +774,8 @@ AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Int_t i) //________________________________________________________________________ void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi, - AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const + AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters, + Bool_t bPartialExclusion) const { // Get rigid cone. @@ -673,12 +810,29 @@ void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &p Float_t dLJ = 0; Int_t repeats = 0; + Bool_t reject = kTRUE; do { eta = gRandom->Rndm() * (maxEta - minEta) + minEta; phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi; dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi)); + + if(bPartialExclusion) { + reject = kFALSE; + + TRandom3 rnd; + rnd.SetSeed(0); + + Double_t ncoll = GetNColl(); + + Double_t prob = 0.; + if(ncoll>0) + prob = 1./ncoll; + + if(rnd.Rndm()<=prob) reject = kTRUE; //reject cone + } + repeats++; - } while (dLJ < fMinRC2LJ && repeats < 999); + } while (dLJ < fMinRC2LJ && repeats < 999 && reject); if (repeats == 999) { AliWarning(Form("%s: Could not get random cone!", GetName())); @@ -825,6 +979,36 @@ void AliAnalysisTaskDeltaPt::ExecOnce() } } +//________________________________________________________________________ +Double_t AliAnalysisTaskDeltaPt::GetNColl() const { + // Get NColl - returns value of corresponding bin + // only works for pA + // values taken from V0A slicing https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Tables_with_centrality_bins_from + + if(fBeamType == kpA) { + + const Int_t nNCollBins = 7; + + Double_t centMin[nNCollBins] = {0.,5.,10.,20.,40.,60.,80.}; + Double_t centMax[nNCollBins] = {5.,10.,20.,40.,60.,80.,100.}; + + Double_t nColl[nNCollBins] = {14.7,13.,11.7,9.38,6.49,3.96,1.52}; + + for(Int_t i = 0; i=centMin[i] && fCentGetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius); - SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius); + + Int_t runno = InputEvent()->GetRunNumber(); + if(runno>=177295 && runno<=197470) //small SM masked in 2012 and 2013 + SetJetPhiLimits(1.4 + fJetRadius, TMath::Pi() - fJetRadius); + else + SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius); + } if (!fRhoName.IsNull() && !fRho) { diff --git a/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.cxx b/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.cxx index 5f54aa50f85..75417027095 100644 --- a/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.cxx +++ b/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.cxx @@ -2,7 +2,7 @@ // // Emcal jet analysis base task. // -// Author: S.Aiola +// Author: S.Aiola, M. Verweij #include "AliAnalysisTaskEmcalJetDev.h" @@ -153,14 +153,14 @@ void AliAnalysisTaskEmcalJetDev::ExecOnce() for(Int_t i =0; i(fJetCollArray.At(i)); - cont->SetJetArray(InputEvent()); + cont->SetRunNumber(InputEvent()->GetRunNumber()); cont->SetEMCALGeometry(); + cont->SetJetArray(InputEvent()); cont->LoadRho(InputEvent()); } //Get Jets, cuts and rho for first jet container AliJetContainer *cont = GetJetContainer(0); - // fJetsName = cont->GetArrayName(); if (fAnaType == kTPC) { cont->SetJetAcceptanceType(AliJetContainer::kTPC); cont->SetJetEtaPhiTPC(); @@ -320,14 +320,18 @@ Bool_t AliAnalysisTaskEmcalJetDev::RetrieveEventObjects() } //________________________________________________________________________ -void AliAnalysisTaskEmcalJetDev::AddJetContainer(const char *n, TString defaultCutType) { +void AliAnalysisTaskEmcalJetDev::AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius) { // Add particle container // will be called in AddTask macro + TString tmp = TString(n); + if(tmp.IsNull()) return; + AliJetContainer *cont = 0x0; cont = new AliJetContainer(); cont->SetArrayName(n); + cont->SetJetRadius(jetRadius); if(!defaultCutType.IsNull()) { if(defaultCutType.EqualTo("TPC")) diff --git a/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.h b/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.h index 17252d84446..26b61a5e8d6 100644 --- a/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.h +++ b/PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetDev.h @@ -39,7 +39,7 @@ class AliAnalysisTaskEmcalJetDev : public AliAnalysisTaskEmcalDev { void SetJetsName(const char *n) { fJetsName = n; AddJetContainer(n); } virtual void SetRhoName(const char *n, Int_t c = 0); - void AddJetContainer(const char *n, TString defaultCutType = ""); + void AddJetContainer(const char *n, TString defaultCutType = "", Float_t jetRadius = 0.4); void RemoveJetContainer(Int_t i) { fJetCollArray.RemoveAt(i);} AliJetContainer *GetJetContainer(Int_t i) const; diff --git a/PWGJE/EMCALJetTasks/AliJetContainer.cxx b/PWGJE/EMCALJetTasks/AliJetContainer.cxx index 34e563869cf..c1ed8d4617a 100644 --- a/PWGJE/EMCALJetTasks/AliJetContainer.cxx +++ b/PWGJE/EMCALJetTasks/AliJetContainer.cxx @@ -42,7 +42,8 @@ AliJetContainer::AliJetContainer(): fNLeadingJets(1), fJetBitMap(0), fRho(0), - fGeom(0) + fGeom(0), + fRunNumber(0) { // Default constructor. @@ -69,7 +70,8 @@ AliJetContainer::AliJetContainer(const char *name): fNLeadingJets(1), fJetBitMap(0), fRho(0), - fGeom(0) + fGeom(0), + fRunNumber(0) { // Standard constructor. @@ -86,15 +88,19 @@ void AliJetContainer::SetJetArray(AliVEvent *event) { // Set jet array - // SetArray(event, fClArrayName.GetName(), "AliEmcalJet"); SetArray(event, "AliEmcalJet"); - if(fJetAcceptanceType==kTPC) + if(fJetAcceptanceType==kTPC) { + AliDebug(2,Form("%s: set TPC acceptance cuts",GetName())); SetJetEtaPhiTPC(); - else if(fJetAcceptanceType==kEMCAL) + } + else if(fJetAcceptanceType==kEMCAL) { + AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName())); SetJetEtaPhiEMCAL(); + } } + //________________________________________________________________________ void AliJetContainer::SetEMCALGeometry() { fGeom = AliEMCALGeometry::GetInstance(); @@ -222,10 +228,15 @@ void AliJetContainer::SetJetEtaPhiEMCAL() if(!fGeom) SetEMCALGeometry(); if(fGeom) { SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius); - SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius); + + if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013 + SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius); + else + SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius); + } else { - AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings"); + AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!"); SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius); SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius); } @@ -235,9 +246,9 @@ void AliJetContainer::SetJetEtaPhiEMCAL() //________________________________________________________________________ void AliJetContainer::SetJetEtaPhiTPC() { - //Set default cuts for full jets + //Set default cuts for charged jets - SetJetEtaLimits(-0.5, 0.5); + SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius); SetJetPhiLimits(-10, 10); } diff --git a/PWGJE/EMCALJetTasks/AliJetContainer.h b/PWGJE/EMCALJetTasks/AliJetContainer.h index 201dc7822f4..f0232f54595 100644 --- a/PWGJE/EMCALJetTasks/AliJetContainer.h +++ b/PWGJE/EMCALJetTasks/AliJetContainer.h @@ -52,6 +52,7 @@ class AliJetContainer : public AliEmcalContainer { void SetJetEtaPhiEMCAL() ; void SetJetEtaPhiTPC() ; + void SetRunNumber(Int_t r) { fRunNumber = r; } void SetJetEtaLimits(Float_t min, Float_t max) { fJetMinEta = min, fJetMaxEta = max ; } void SetJetPhiLimits(Float_t min, Float_t max) { fJetMinPhi = min, fJetMaxPhi = max ; } @@ -81,6 +82,7 @@ class AliJetContainer : public AliEmcalContainer { Double_t GetRhoVal() const {return fRho->GetVal();} const TString& GetRhoName() const {return fRhoName;} Double_t GetJetPtCorr(Int_t i) const; + Float_t GetJetRadius() const {return fJetRadius;} protected: JetAcceptanceType fJetAcceptanceType; // acceptance type @@ -104,6 +106,7 @@ class AliJetContainer : public AliEmcalContainer { AliRhoParameter *fRho; //!event rho for these jets AliEMCALGeometry *fGeom; //!emcal geometry + Int_t fRunNumber; // run number private: AliJetContainer(const AliJetContainer& obj); // copy constructor