From fa7c34ba27bd192f9d929589e9da8580bbb04fd3 Mon Sep 17 00:00:00 2001 From: kleinb Date: Mon, 18 Jun 2012 11:41:44 +0000 Subject: [PATCH] Added data driven estimate for leading particle efficiency + coverity fixes (M. Verweij) --- JETAN/AliAnalysisTaskJetCluster.cxx | 91 +++++++++++++++++---- JETAN/AliAnalysisTaskJetCluster.h | 20 +++-- JETAN/DEV/AliAnalysisTaskJetCluster.cxx | 91 +++++++++++++++++---- JETAN/DEV/AliAnalysisTaskJetCluster.h | 20 +++-- PWGJE/AliAnaChargedJetResponseMaker.cxx | 102 ++++++++++++++++++++---- PWGJE/AliAnaChargedJetResponseMaker.h | 74 ++++++++--------- PWGJE/macros/AddTaskJetCluster.C | 13 ++- 7 files changed, 311 insertions(+), 100 deletions(-) diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx index b088c9cf5a8..3e1bf13c762 100644 --- a/JETAN/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/AliAnalysisTaskJetCluster.cxx @@ -136,8 +136,13 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): fhEffH1(0x0), fhEffH2(0x0), fhEffH3(0x0), - fUseTrMomentumSmearing(kFALSE), + fUseTrPtResolutionSmearing(kFALSE), fUseDiceEfficiency(kFALSE), + fUseTrPtResolutionFromOADB(kFALSE), + fUseTrEfficiencyFromOADB(kFALSE), + fPathTrPtResolution(""), + fPathTrEfficiency(""), + fChangeEfficiencyFraction(0.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -268,8 +273,13 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fhEffH1(0x0), fhEffH2(0x0), fhEffH3(0x0), - fUseTrMomentumSmearing(kFALSE), + fUseTrPtResolutionSmearing(kFALSE), fUseDiceEfficiency(kFALSE), + fUseTrPtResolutionFromOADB(kFALSE), + fUseTrEfficiencyFromOADB(kFALSE), + fPathTrPtResolution(""), + fPathTrEfficiency(""), + fChangeEfficiencyFraction(0.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -404,8 +414,8 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() if(fJetTypes&kJetRan){ fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); - if(fUseDiceEfficiency || fUseTrMomentumSmearing) - fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%d",fNonStdBranch.Data(),"random",fUseTrMomentumSmearing,fUseDiceEfficiency)); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) + fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data()); } @@ -413,16 +423,16 @@ 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 || fUseTrMomentumSmearing) - fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency)); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) + 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 || fUseTrMomentumSmearing) - cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) + cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone); if(fNRandomCones>0){ if(fJetTypes&kRC){ @@ -453,8 +463,6 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() } } - // FitMomentumResolution(); - if(!fHistList)fHistList = new TList(); fHistList->SetOwner(); @@ -705,7 +713,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() TH1::AddDirectory(oldStatus); } -void AliAnalysisTaskJetCluster::Init() +void AliAnalysisTaskJetCluster::LocalInit() { // // Initialization @@ -713,6 +721,10 @@ void AliAnalysisTaskJetCluster::Init() if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n"); + if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB(); + if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB(); + + FitMomentumResolution(); } @@ -759,7 +771,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) } //Check if information is provided detector level effects - if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE; + if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE; if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE; Bool_t selectEvent = false; @@ -846,7 +858,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // Carefull energy is not well determined in real data, should not matter for p_T scheme? // we take total momentum here - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) { + 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); @@ -916,7 +928,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fp1Efficiency->Fill(vp->Pt(),sumEff); if(rnd>sumEff) continue; - if(fUseTrMomentumSmearing) { + if(fUseTrPtResolutionSmearing) { //Smear momentum of generated particle Double_t smear = 1.; //Select hybrid track category @@ -978,11 +990,11 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(fTCAJetsOut){ if(i == 0){ fRef->Delete(); // make sure to delete before placement new... - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) { + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) { new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ... } } - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ... + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ... } }// recparticles @@ -1129,7 +1141,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(!part) continue; fh1PtJetConstRec->Fill(part->Pt()); if(aodOutJet){ - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) 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); } @@ -1690,6 +1702,51 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ return iCount; } +void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() { + + TFile *f = new TFile(fPathTrPtResolution.Data()); + + TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt"); + TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS"); + TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD"); + + SetSmearResolution(kTRUE); + SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD); + + + if(f) delete f; + +} + +void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() { + + TFile *f = new TFile(fPathTrEfficiency.Data()); + + 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); + + if(f) delete f; + +} + void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) { // diff --git a/JETAN/AliAnalysisTaskJetCluster.h b/JETAN/AliAnalysisTaskJetCluster.h index e1f93dae23b..f055d3c2ef0 100644 --- a/JETAN/AliAnalysisTaskJetCluster.h +++ b/JETAN/AliAnalysisTaskJetCluster.h @@ -44,8 +44,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE virtual ~AliAnalysisTaskJetCluster(); // Implementation of interface methods virtual void UserCreateOutputObjects(); - virtual void Init(); - virtual void LocalInit() { Init(); } + virtual void LocalInit(); virtual void UserExec(Option_t *option); virtual void Terminate(Option_t *option); virtual Bool_t Notify(); @@ -83,7 +82,12 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE virtual void SetBackgroundCalc(Bool_t b){fUseBackgroundCalc = b;} //Setters for detector level effects - virtual void SetSmearResolution(Bool_t b){fUseTrMomentumSmearing = b;} + virtual void SetUseTrResolutionFromOADB(Bool_t b=kTRUE, TString path="$ALICE_ROOT/OADB/PWGJE/Resolution/PtResol_LHCh_Cent0-10_v1.root") {fUseTrPtResolutionFromOADB = b; fPathTrPtResolution=path;} + virtual void SetUseTrEfficiencyFromOADB(Bool_t b=kTRUE, TString path="$ALICE_ROOT/OADB/PWGJE/Efficiency/Efficiency_LHC11a2aj_Cent0_v1.root") {fUseTrEfficiencyFromOADB = b; fPathTrEfficiency=path;} + virtual void LoadTrEfficiencyRootFileFromOADB(); + virtual void LoadTrPtResolutionRootFileFromOADB(); + virtual void SetChangeEfficiencyFraction(Double_t p) {fChangeEfficiencyFraction = p;} + virtual void SetSmearResolution(Bool_t b){fUseTrPtResolutionSmearing = b;} virtual void SetDiceEfficiency(Bool_t b){fUseDiceEfficiency = b;} virtual void SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3); virtual void SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3); @@ -180,8 +184,14 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE TH1 *fhEffH1; // Efficiency for Spectra Hybrid Category 1 TH1 *fhEffH2; // Efficiency for Spectra Hybrid Category 2 TH1 *fhEffH3; // Efficiency for Spectra Hybrid Category 3 - Bool_t fUseTrMomentumSmearing; // Apply momentum smearing on track level - Bool_t fUseDiceEfficiency; // Apply efficiency on track level by dicing + Bool_t fUseTrPtResolutionSmearing; // Apply momentum smearing on track level + Bool_t fUseDiceEfficiency; // Apply efficiency on track level by dicing + Bool_t fUseTrPtResolutionFromOADB; // Load track pt resolution root file from OADB path + Bool_t fUseTrEfficiencyFromOADB; // Load tracking efficiency root file from OADB path + TString fPathTrPtResolution; // OADB path to root file + TString fPathTrEfficiency; // OADB path to root file + Double_t fChangeEfficiencyFraction; //change efficiency by fraction + // Fast jet Double_t fRparam; // fastjet distance parameter diff --git a/JETAN/DEV/AliAnalysisTaskJetCluster.cxx b/JETAN/DEV/AliAnalysisTaskJetCluster.cxx index b088c9cf5a8..3e1bf13c762 100644 --- a/JETAN/DEV/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/DEV/AliAnalysisTaskJetCluster.cxx @@ -136,8 +136,13 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): fhEffH1(0x0), fhEffH2(0x0), fhEffH3(0x0), - fUseTrMomentumSmearing(kFALSE), + fUseTrPtResolutionSmearing(kFALSE), fUseDiceEfficiency(kFALSE), + fUseTrPtResolutionFromOADB(kFALSE), + fUseTrEfficiencyFromOADB(kFALSE), + fPathTrPtResolution(""), + fPathTrEfficiency(""), + fChangeEfficiencyFraction(0.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -268,8 +273,13 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fhEffH1(0x0), fhEffH2(0x0), fhEffH3(0x0), - fUseTrMomentumSmearing(kFALSE), + fUseTrPtResolutionSmearing(kFALSE), fUseDiceEfficiency(kFALSE), + fUseTrPtResolutionFromOADB(kFALSE), + fUseTrEfficiencyFromOADB(kFALSE), + fPathTrPtResolution(""), + fPathTrEfficiency(""), + fChangeEfficiencyFraction(0.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -404,8 +414,8 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() if(fJetTypes&kJetRan){ fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); - if(fUseDiceEfficiency || fUseTrMomentumSmearing) - fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%d",fNonStdBranch.Data(),"random",fUseTrMomentumSmearing,fUseDiceEfficiency)); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) + fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.))); AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data()); } @@ -413,16 +423,16 @@ 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 || fUseTrMomentumSmearing) - fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency)); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) + 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 || fUseTrMomentumSmearing) - cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) + cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone); if(fNRandomCones>0){ if(fJetTypes&kRC){ @@ -453,8 +463,6 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() } } - // FitMomentumResolution(); - if(!fHistList)fHistList = new TList(); fHistList->SetOwner(); @@ -705,7 +713,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() TH1::AddDirectory(oldStatus); } -void AliAnalysisTaskJetCluster::Init() +void AliAnalysisTaskJetCluster::LocalInit() { // // Initialization @@ -713,6 +721,10 @@ void AliAnalysisTaskJetCluster::Init() if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n"); + if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB(); + if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB(); + + FitMomentumResolution(); } @@ -759,7 +771,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) } //Check if information is provided detector level effects - if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE; + if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE; if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE; Bool_t selectEvent = false; @@ -846,7 +858,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // Carefull energy is not well determined in real data, should not matter for p_T scheme? // we take total momentum here - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) { + 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); @@ -916,7 +928,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fp1Efficiency->Fill(vp->Pt(),sumEff); if(rnd>sumEff) continue; - if(fUseTrMomentumSmearing) { + if(fUseTrPtResolutionSmearing) { //Smear momentum of generated particle Double_t smear = 1.; //Select hybrid track category @@ -978,11 +990,11 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(fTCAJetsOut){ if(i == 0){ fRef->Delete(); // make sure to delete before placement new... - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) { + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) { new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ... } } - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ... + if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ... } }// recparticles @@ -1129,7 +1141,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) if(!part) continue; fh1PtJetConstRec->Fill(part->Pt()); if(aodOutJet){ - if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) 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); } @@ -1690,6 +1702,51 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ return iCount; } +void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() { + + TFile *f = new TFile(fPathTrPtResolution.Data()); + + TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt"); + TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS"); + TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD"); + + SetSmearResolution(kTRUE); + SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD); + + + if(f) delete f; + +} + +void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() { + + TFile *f = new TFile(fPathTrEfficiency.Data()); + + 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); + + if(f) delete f; + +} + void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) { // diff --git a/JETAN/DEV/AliAnalysisTaskJetCluster.h b/JETAN/DEV/AliAnalysisTaskJetCluster.h index e1f93dae23b..f055d3c2ef0 100644 --- a/JETAN/DEV/AliAnalysisTaskJetCluster.h +++ b/JETAN/DEV/AliAnalysisTaskJetCluster.h @@ -44,8 +44,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE virtual ~AliAnalysisTaskJetCluster(); // Implementation of interface methods virtual void UserCreateOutputObjects(); - virtual void Init(); - virtual void LocalInit() { Init(); } + virtual void LocalInit(); virtual void UserExec(Option_t *option); virtual void Terminate(Option_t *option); virtual Bool_t Notify(); @@ -83,7 +82,12 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE virtual void SetBackgroundCalc(Bool_t b){fUseBackgroundCalc = b;} //Setters for detector level effects - virtual void SetSmearResolution(Bool_t b){fUseTrMomentumSmearing = b;} + virtual void SetUseTrResolutionFromOADB(Bool_t b=kTRUE, TString path="$ALICE_ROOT/OADB/PWGJE/Resolution/PtResol_LHCh_Cent0-10_v1.root") {fUseTrPtResolutionFromOADB = b; fPathTrPtResolution=path;} + virtual void SetUseTrEfficiencyFromOADB(Bool_t b=kTRUE, TString path="$ALICE_ROOT/OADB/PWGJE/Efficiency/Efficiency_LHC11a2aj_Cent0_v1.root") {fUseTrEfficiencyFromOADB = b; fPathTrEfficiency=path;} + virtual void LoadTrEfficiencyRootFileFromOADB(); + virtual void LoadTrPtResolutionRootFileFromOADB(); + virtual void SetChangeEfficiencyFraction(Double_t p) {fChangeEfficiencyFraction = p;} + virtual void SetSmearResolution(Bool_t b){fUseTrPtResolutionSmearing = b;} virtual void SetDiceEfficiency(Bool_t b){fUseDiceEfficiency = b;} virtual void SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3); virtual void SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3); @@ -180,8 +184,14 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE TH1 *fhEffH1; // Efficiency for Spectra Hybrid Category 1 TH1 *fhEffH2; // Efficiency for Spectra Hybrid Category 2 TH1 *fhEffH3; // Efficiency for Spectra Hybrid Category 3 - Bool_t fUseTrMomentumSmearing; // Apply momentum smearing on track level - Bool_t fUseDiceEfficiency; // Apply efficiency on track level by dicing + Bool_t fUseTrPtResolutionSmearing; // Apply momentum smearing on track level + Bool_t fUseDiceEfficiency; // Apply efficiency on track level by dicing + Bool_t fUseTrPtResolutionFromOADB; // Load track pt resolution root file from OADB path + Bool_t fUseTrEfficiencyFromOADB; // Load tracking efficiency root file from OADB path + TString fPathTrPtResolution; // OADB path to root file + TString fPathTrEfficiency; // OADB path to root file + Double_t fChangeEfficiencyFraction; //change efficiency by fraction + // Fast jet Double_t fRparam; // fastjet distance parameter diff --git a/PWGJE/AliAnaChargedJetResponseMaker.cxx b/PWGJE/AliAnaChargedJetResponseMaker.cxx index 535fa32f10d..80e6bb82272 100644 --- a/PWGJE/AliAnaChargedJetResponseMaker.cxx +++ b/PWGJE/AliAnaChargedJetResponseMaker.cxx @@ -45,9 +45,75 @@ AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(): fbCalcErrors(kFALSE) {;} -AliAnaChargedJetResponseMaker::~AliAnaChargedJetResponseMaker() { - //destructor +//-------------------------------------------------------------------------------------------------------------------------------------------------- +AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj): + fDebug(obj.fDebug), + fResolutionType(obj.fResolutionType), + fDeltaPt(obj.fDeltaPt), + fhDeltaPt(obj.fhDeltaPt), + fDimensions(obj.fDimensions), + fDimRec(obj.fDimRec), + fDimGen(obj.fDimGen), + fPtMin(obj.fPtMin), + fPtMax(obj.fPtMax), + fNbins(obj.fNbins), + fBinArrayPtRec(obj.fBinArrayPtRec), + fPtMeasured(obj.fPtMeasured), + fEffFlat(obj.fEffFlat), + fEfficiency(obj.fEfficiency), + fEfficiencyFine(obj.fEfficiencyFine), + fResponseMatrix(obj.fResponseMatrix), + fResponseMatrixFine(obj.fResponseMatrixFine), + fPtMinUnfolded(obj.fPtMinUnfolded), + fPtMaxUnfolded(obj.fPtMaxUnfolded), + fPtMaxUnfoldedHigh(obj.fPtMaxUnfoldedHigh), + fBinWidthFactorUnfolded(obj.fBinWidthFactorUnfolded), + fSkipBinsUnfolded(obj.fSkipBinsUnfolded), + fExtraBinsUnfolded(obj.fExtraBinsUnfolded), + fbVariableBinning(obj.fbVariableBinning), + fPtMaxUnfVarBinning(obj.fPtMaxUnfVarBinning), + f1MergeFunction(obj.f1MergeFunction), + fFineFrac(obj.fFineFrac), + fbCalcErrors(obj.fbCalcErrors) +{;} + +//-------------------------------------------------------------------------------------------------------------------------------------------------- +AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other) +{ +// Assignment + if(&other == this) return *this; + AliAnaChargedJetResponseMaker::operator=(other); + fDebug = other.fDebug; + fResolutionType = other.fResolutionType; + fDeltaPt = other.fDeltaPt; + fhDeltaPt = other.fhDeltaPt; + fDimensions = other.fDimensions; + fDimRec = other.fDimRec; + fDimGen = other.fDimGen; + fPtMin = other.fPtMin; + fPtMax = other.fPtMax; + fNbins = other.fNbins; + fBinArrayPtRec = other.fBinArrayPtRec; + fPtMeasured = other.fPtMeasured; + fEffFlat = other.fEffFlat; + fEfficiency = other.fEfficiency; + fEfficiencyFine = other.fEfficiencyFine; + fResponseMatrix = other.fResponseMatrix; + fResponseMatrixFine = other.fResponseMatrixFine; + fPtMinUnfolded = other.fPtMinUnfolded; + fPtMaxUnfolded = other.fPtMaxUnfolded; + fPtMaxUnfoldedHigh = other.fPtMaxUnfoldedHigh; + fBinWidthFactorUnfolded = other.fBinWidthFactorUnfolded; + fSkipBinsUnfolded = other.fSkipBinsUnfolded; + fExtraBinsUnfolded = other.fExtraBinsUnfolded; + fbVariableBinning = other.fbVariableBinning; + fPtMaxUnfVarBinning = other.fPtMaxUnfVarBinning; + f1MergeFunction = other.f1MergeFunction; + fFineFrac = other.fFineFrac; + fbCalcErrors = other.fbCalcErrors; + + return *this; } //-------------------------------------------------------------------------------------------------------------------------------------------------- @@ -206,6 +272,7 @@ void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() { Int_t nbins[fDimensions*2]; nbins[fDimRec] = fNbins; + nbins[fDimGen] = fNbins; double binWidthMeas = (double)((fPtMax-fPtMin)/fNbins); double binWidthUnf = (double)fBinWidthFactorUnfolded*binWidthMeas; @@ -263,13 +330,15 @@ void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() { binArrayPt0[nbins[fDimRec]]= xmax[fDimRec]; //Define bin limits generated/unfolded axis - for(int i=0; iGetAxis(fDimRec)->GetNbins()*fFineFrac; // dimension 0 is reconstructed value - else nbinsFine[i] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; // dimension 1 is generated value - - xminFine[i] = TMath::Min(fPtMin,0.); - xmaxFine[i] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.; - - pTarrayFine[i] = 0.; - } + nbinsFine[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; + nbinsFine[fDimGen] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; + xminFine[fDimRec] = TMath::Min(fPtMin,0.); + xminFine[fDimGen] = TMath::Min(fPtMin,0.); + xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.; + xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.; + pTarrayFine[fDimRec] = 0.; + pTarrayFine[fDimGen] = 0.; Double_t binWidth[2]; binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1); @@ -443,6 +510,8 @@ void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() { xminFine[fDimRec] = recAxis->GetXmin(); xmaxFine[fDimGen] = genAxis->GetXmax(); xmaxFine[fDimRec] = recAxis->GetXmax(); + pTarrayFine[fDimGen] = 0.; + pTarrayFine[fDimRec] = 0.; double sumyield = 0.; double sumyielderror2 = 0.; @@ -623,7 +692,6 @@ TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *h Double_t sumYield = 0.; const Int_t nbinsRec = hRec->GetNbinsX(); Double_t sumError2[nbinsRec+1]; - for(int i=0; i<=nbinsRec; i++) sumError2[i] = 0.; Double_t eff = 0.; for(int igen=1; igen<=hGen->GetNbinsX(); igen++) { diff --git a/PWGJE/AliAnaChargedJetResponseMaker.h b/PWGJE/AliAnaChargedJetResponseMaker.h index 79f661d06a7..520e6fca5a2 100644 --- a/PWGJE/AliAnaChargedJetResponseMaker.h +++ b/PWGJE/AliAnaChargedJetResponseMaker.h @@ -18,7 +18,9 @@ class TGraphErrors; class AliAnaChargedJetResponseMaker { public: AliAnaChargedJetResponseMaker(); - ~AliAnaChargedJetResponseMaker(); + AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj); // copy constructor + AliAnaChargedJetResponseMaker& operator=(const AliAnaChargedJetResponseMaker& other); // assignment + virtual ~AliAnaChargedJetResponseMaker() {;} // kParam = use parametrization of response // kResiduals = use response as measured, w/o statistical error propagation @@ -26,59 +28,59 @@ class AliAnaChargedJetResponseMaker { enum ResolutionType {kParam,kResiduals,kResidualsErr}; //Setters - void SetDebugMode(Bool_t b) {fDebug=b;} + virtual void SetDebugMode(Bool_t b) {fDebug=b;} - void SetResolutionType(ResolutionType r) {fResolutionType=r;} + virtual void SetResolutionType(ResolutionType r) {fResolutionType=r;} - void SetDeltaPtJetsFunc(TF1 *f1) {fDeltaPt=f1;} - void SetDeltaPtJetsHist(TH1D *h1) {fhDeltaPt=h1;} - void SetNDimensions(Int_t dim) {fDimensions = dim;} - void SetMeasuredSpectrum(TH1D *hPtMeasured); + virtual void SetDeltaPtJetsFunc(TF1 *f1) {fDeltaPt=f1;} + virtual void SetDeltaPtJetsHist(TH1D *h1) {fhDeltaPt=h1;} + virtual void SetNDimensions(Int_t dim) {fDimensions = dim;} + virtual void SetMeasuredSpectrum(TH1D *hPtMeasured); - void SetFlatEfficiency(Double_t eff); - void SetEfficiency(TGraphErrors *grEff); + virtual void SetFlatEfficiency(Double_t eff); + virtual void SetEfficiency(TGraphErrors *grEff); - void SetPtMinUnfolded(Double_t ptmin) {fPtMinUnfolded = ptmin;} - void SetPtMaxUnfolded(Double_t ptmax) {fPtMaxUnfolded = ptmax;} - void SetPtMaxUnfoldedHigh(Double_t ptmaxh) {fPtMaxUnfoldedHigh = ptmaxh;} - void SetBinWidthFactorUnfolded(Int_t fac) {fBinWidthFactorUnfolded = fac;} - void SetSkipBinsUnfolded(Int_t skip) {fSkipBinsUnfolded=skip;} - void SetExtraBinsUnfolded(Int_t extra) {fExtraBinsUnfolded=extra;} - void SetVariableBinning(Bool_t b, double ptmax) { + virtual void SetPtMinUnfolded(Double_t ptmin) {fPtMinUnfolded = ptmin;} + virtual void SetPtMaxUnfolded(Double_t ptmax) {fPtMaxUnfolded = ptmax;} + virtual void SetPtMaxUnfoldedHigh(Double_t ptmaxh) {fPtMaxUnfoldedHigh = ptmaxh;} + virtual void SetBinWidthFactorUnfolded(Int_t fac) {fBinWidthFactorUnfolded = fac;} + virtual void SetSkipBinsUnfolded(Int_t skip) {fSkipBinsUnfolded=skip;} + virtual void SetExtraBinsUnfolded(Int_t extra) {fExtraBinsUnfolded=extra;} + virtual void SetVariableBinning(Bool_t b, double ptmax) { fbVariableBinning=b; fPtMaxUnfVarBinning=ptmax; } - void SetCalcErrors(Bool_t b) {fbCalcErrors=b;} + virtual void SetCalcErrors(Bool_t b) {fbCalcErrors=b;} //Setters for merging fine to normal response matrix - void SetFineFrac(Int_t i = 10) {fFineFrac = i;} - void SetRMMergeWeightFunction(TF1 *f1) {f1MergeFunction = f1;} + virtual void SetFineFrac(Int_t i = 10) {fFineFrac = i;} + virtual void SetRMMergeWeightFunction(TF1 *f1) {f1MergeFunction = f1;} //Getters - TF1 *GetDeltaPtJetsFunc() {return fDeltaPt;} - TH1D *GetDeltaPtJetsHist() {return fhDeltaPt;} - THnSparse *GetMeasuredSpectrum() {return fPtMeasured;} - THnSparse *GetEfficiency() {return fEfficiency;} - THnSparse *GetEfficiencyFine() {return fEfficiencyFine;} - THnSparse *GetResponseMatrix() {return fResponseMatrix;} - THnSparse *GetResponseMatrixFine() {return fResponseMatrixFine;} + virtual TF1 *GetDeltaPtJetsFunc() {return fDeltaPt;} + virtual TH1D *GetDeltaPtJetsHist() {return fhDeltaPt;} + virtual THnSparse *GetMeasuredSpectrum() {return fPtMeasured;} + virtual THnSparse *GetEfficiency() {return fEfficiency;} + virtual THnSparse *GetEfficiencyFine() {return fEfficiencyFine;} + virtual THnSparse *GetResponseMatrix() {return fResponseMatrix;} + virtual THnSparse *GetResponseMatrixFine() {return fResponseMatrixFine;} //Utility functions - Double_t InterpolateFast(TGraph *gr, Double_t x); - Double_t InterpolateFast(TH1 *h, Double_t x); + virtual Double_t InterpolateFast(TGraph *gr, Double_t x); + virtual Double_t InterpolateFast(TH1 *h, Double_t x); - TH1D *MultiplyResponseGenerated(TH1 *hGen=0, TH2 *hResponse=0,TH1 *hEfficiency=0,Bool_t bDrawSlices=kFALSE); - TH1D *MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency); + virtual TH1D *MultiplyResponseGenerated(TH1 *hGen=0, TH2 *hResponse=0,TH1 *hEfficiency=0,Bool_t bDrawSlices=kFALSE); + virtual TH1D *MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency); - void MakeResponseMatrixJetsFineMerged(Int_t skipBins =0, Int_t binWidthFactor = 2, Int_t extraBins = 0, Bool_t bVariableBinning = kFALSE, Double_t ptmin = 0.); + virtual void MakeResponseMatrixJetsFineMerged(Int_t skipBins =0, Int_t binWidthFactor = 2, Int_t extraBins = 0, Bool_t bVariableBinning = kFALSE, Double_t ptmin = 0.); - void InitializeResponseMatrix(); - void InitializeResponseMatrixFine(); + virtual void InitializeResponseMatrix(); + virtual void InitializeResponseMatrixFine(); - void InitializeEfficiency(); - void InitializeEfficiencyFine(); + virtual void InitializeEfficiency(); + virtual void InitializeEfficiencyFine(); - void FillResponseMatrixFineAndMerge(); + virtual void FillResponseMatrixFineAndMerge(); protected: Bool_t fDebug; diff --git a/PWGJE/macros/AddTaskJetCluster.C b/PWGJE/macros/AddTaskJetCluster.C index 984ff61159e..c95170d18e1 100644 --- a/PWGJE/macros/AddTaskJetCluster.C +++ b/PWGJE/macros/AddTaskJetCluster.C @@ -1,4 +1,4 @@ -AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec = "AOD",char* bGen = "",UInt_t filterMask = 16, UInt_t iPhysicsSelectionFlag = AliVEvent::kAny,Char_t *jf = "KT", Float_t radius = 0.4,Int_t nSkip = 0,Int_t kWriteAOD = kFALSE,char* deltaFile = "",Float_t ptTrackCut = 0.15, Float_t etaTrackWindow = 0.9,Float_t vertexWindow = 10.,Int_t nSkipCone = 2,Int_t dice=0,Int_t smear=0); +AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec = "AOD",char* bGen = "",UInt_t filterMask = 16, UInt_t iPhysicsSelectionFlag = AliVEvent::kAny,Char_t *jf = "KT", Float_t radius = 0.4,Int_t nSkip = 0,Int_t kWriteAOD = kFALSE,char* deltaFile = "",Float_t ptTrackCut = 0.15, Float_t etaTrackWindow = 0.9,Float_t vertexWindow = 10.,Int_t nSkipCone = 2,Int_t dice=0,Int_t smear=0,Bool_t useOADB=kFALSE,Double_t changeEfficiencyFraction=0.); Int_t kBackgroundModeCl = 0; Float_t kPtTrackCutCl = 0.15; @@ -30,7 +30,7 @@ AliAnalysisTaskJetCluster *AddTaskJetClusterDelta(UInt_t filterMask = 16,Bool_t } -AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filterMask,UInt_t iPhysicsSelectionFlag,Char_t *jf,Float_t radius,Int_t nSkip,Int_t kWriteAOD,char *deltaFile,Float_t ptTrackCut,Float_t etaTrackWindow,Float_t vertexWindow,Int_t nSkipCone,Int_t dice,Int_t smear) +AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filterMask,UInt_t iPhysicsSelectionFlag,Char_t *jf,Float_t radius,Int_t nSkip,Int_t kWriteAOD,char *deltaFile,Float_t ptTrackCut,Float_t etaTrackWindow,Float_t vertexWindow,Int_t nSkipCone,Int_t dice,Int_t smear,Bool_t useOADB,Double_t changeEfficiencyFraction) { // Creates a jet fider task, configures it and adds it to the analysis manager. kPtTrackCutCl = ptTrackCut; @@ -74,7 +74,8 @@ AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filte cAdd += Form("_Cut%05d",(int)(1000.*kPtTrackCutCl)); cAdd += Form("_Skip%02d",nSkip); if(dice>0 || smear>0) - cAdd += Form("Detector%d%d",dice,smear); + cAdd += Form("Detector%d%dFr%d",dice,smear,(int)(changeEfficiencyFraction*100.)); + Printf("%s %E %d %d",cAdd.Data(),kPtTrackCutCl,dice,smear); AliAnalysisTaskJetCluster* clus = new AliAnalysisTaskJetCluster(Form("JetCluster%s_%s%s",bRec,jf,cAdd.Data())); @@ -155,6 +156,12 @@ AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filte if(iPhysicsSelectionFlag)clus->SelectCollisionCandidates(iPhysicsSelectionFlag); + if(useOADB) { + clus->SetUseTrResolutionFromOADB(); + clus->SetUseTrEfficiencyFromOADB(); + clus->SetChangeEfficiencyFraction(changeEfficiencyFraction); + } + mgr->AddTask(clus); // Create ONLY the output containers for the data produced by the task. -- 2.39.3