From a8ef19993122928e20e8a567f6f679557d7c111d Mon Sep 17 00:00:00 2001 From: rbailhac Date: Tue, 6 Mar 2012 10:05:06 +0000 Subject: [PATCH] Add flow tasks to the compilation, update others --- PWGHF/CMakelibPWGHFhfe.pkg | 4 +- PWGHF/PWGHFhfeLinkDef.h | 3 + PWGHF/hfe/AliAnalysisTaskHFE.cxx | 232 +++-- PWGHF/hfe/AliAnalysisTaskHFE.h | 4 +- PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx | 291 ++++++- PWGHF/hfe/AliAnalysisTaskHFEFlow.h | 83 +- PWGHF/hfe/AliHFEVZEROEventPlane.cxx | 407 +++++++++ PWGHF/hfe/AliHFEVZEROEventPlane.h | 59 ++ PWGHF/hfe/AliHFEcollection.cxx | 32 + PWGHF/hfe/AliHFEcollection.h | 2 + PWGHF/hfe/AliHFEdebugTreeTask.cxx | 30 +- PWGHF/hfe/AliHFEdebugTreeTask.h | 2 + PWGHF/hfe/AliHFEemcalPIDqa.cxx | 8 +- PWGHF/hfe/AliHFEextraCuts.h | 4 +- PWGHF/hfe/AliHFEmcQA.cxx | 539 +++++++++--- PWGHF/hfe/AliHFEmcQA.h | 26 +- PWGHF/hfe/AliHFEpidEMCAL.cxx | 29 +- PWGHF/hfe/AliHFEpidQA.cxx | 2 +- PWGHF/hfe/AliHFEpidQA.h | 2 +- PWGHF/hfe/AliHFEsignalCuts.cxx | 18 + PWGHF/hfe/AliHFEspectrum.cxx | 1171 +++++++++++++++++++++----- PWGHF/hfe/AliHFEspectrum.h | 43 +- PWGHF/hfe/AliHFEtools.cxx | 161 ++++ PWGHF/hfe/AliHFEtools.h | 8 + PWGHF/hfe/AliHFEtrdPIDqa.cxx | 38 +- PWGHF/hfe/AliHFEtrdPIDqa.h | 11 +- 26 files changed, 2644 insertions(+), 565 deletions(-) create mode 100644 PWGHF/hfe/AliHFEVZEROEventPlane.cxx create mode 100644 PWGHF/hfe/AliHFEVZEROEventPlane.h diff --git a/PWGHF/CMakelibPWGHFhfe.pkg b/PWGHF/CMakelibPWGHFhfe.pkg index a8c1964286e..c6f8284c650 100644 --- a/PWGHF/CMakelibPWGHFhfe.pkg +++ b/PWGHF/CMakelibPWGHFhfe.pkg @@ -79,10 +79,12 @@ set (SRCS hfe/AliHFEpidObject.cxx hfe/AliAnalysisTaskElecHadronCorrel.cxx hfe/AliHFEdebugTreeTask.cxx + hfe/AliHFEVZEROEventPlane.cxx + hfe/AliAnalysisTaskHFEFlow.cxx ) string (REPLACE ".cxx" ".h" HDRS "${SRCS}") set ( DHDR PWGHFhfeLinkDef.h) -set ( EINCLUDE PWGHF/hfe TPC CORRFW STEER/STEERBase) +set ( EINCLUDE PWGHF/hfe TPC CORRFW STEER/STEERBase PWG/FLOW/Base PWG/FLOW/Tasks) diff --git a/PWGHF/PWGHFhfeLinkDef.h b/PWGHF/PWGHFhfeLinkDef.h index 5cc2f59f42d..7bdd0f03f8a 100644 --- a/PWGHF/PWGHFhfeLinkDef.h +++ b/PWGHF/PWGHFhfeLinkDef.h @@ -70,4 +70,7 @@ #pragma link C++ class AliHFEdebugTreeTask+; +#pragma link C++ class AliHFEVZEROEventPlane+; +#pragma link C++ class AliAnalysisTaskHFEFlow+; + #endif diff --git a/PWGHF/hfe/AliAnalysisTaskHFE.cxx b/PWGHF/hfe/AliAnalysisTaskHFE.cxx index c5386ab00cb..c9e6eb0b8fd 100644 --- a/PWGHF/hfe/AliAnalysisTaskHFE.cxx +++ b/PWGHF/hfe/AliAnalysisTaskHFE.cxx @@ -13,7 +13,6 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ // -// // The analysis task: // Filling an AliCFContainer with the quantities pt, eta and phi // for tracks which survivied the particle cuts (MC resp. ESD tracks) @@ -432,6 +431,9 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){ fHistMCQA->SetOwner(); if(IsPbPb()) fMCQA->SetPbPb(); if(fisppMultiBin) fMCQA->SetPPMultiBin(); + if(TestBit(kTreeStream)){ + fMCQA->EnableDebugStreamer(); + } fMCQA->CreatDefaultHistograms(fHistMCQA); fMCQA->SetBackgroundWeightFactor(fElecBackgroundFactor[0][0][0],fBinLimit); fQA->Add(fHistMCQA); @@ -711,19 +713,9 @@ void AliAnalysisTaskHFE::ProcessMC(){ fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty); fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm); fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty); - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 0); // no accept cut - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 0); // no accept cut - if (TMath::Abs(mcpart->Eta()) < 0.9) { - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9 - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9 - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9 - } - if (TMath::Abs(AliHFEtools::GetRapidity(mcpart)) < 0.5) { - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5 - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5 - fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5 - } + fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG); // no accept cut + fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG); // no accept cut + fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG); // no accept cut } //fMCQA->EndOfEventAna(AliHFEmcQA::kCharm); //fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty); @@ -830,6 +822,7 @@ void AliAnalysisTaskHFE::ProcessESD(){ memset(container, 0, sizeof(Double_t) * 10); // container for the output THnSparse Double_t dataE[6]; // [pT, eta, Phi, type, 'C' or 'B'] + Double_t dataDca[7]; // [source, pT, dca, dcaSig, centrality] Int_t nElectronCandidates = 0; AliESDtrack *track = NULL, *htrack = NULL; AliMCParticle *mctrack = NULL; @@ -938,10 +931,22 @@ void AliAnalysisTaskHFE::ProcessESD(){ if(HasMCData()){ FillProductionVertex(track); - if(fMCQA){ + if(fMCQA && signal){ fMCQA->SetCentrality(fCentralityF); if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){ Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.}; + Double_t hfeimpactRtmp=0., hfeimpactnsigmaRtmp=0.; + fExtraCuts->GetHFEImpactParameters(track, hfeimpactRtmp, hfeimpactnsigmaRtmp); + UChar_t itsPixel = track->GetITSClusterMap(); + Double_t ilyrhit=0, ilyrstat=0; + for(Int_t ilyr=0; ilyr<6; ilyr++){ + if(TESTBIT(itsPixel, ilyr)) ilyrhit += TMath::Power(2,ilyr); + if(fExtraCuts->CheckITSstatus(fExtraCuts->GetITSstatus(track,ilyr))) ilyrstat += TMath::Power(2,ilyr); + } + fMCQA->SetITSInfo(ilyrhit,ilyrstat); + fMCQA->SetHFEImpactParameters(hfeimpactRtmp, hfeimpactnsigmaRtmp); + fMCQA->SetTrkKine(track->Pt(),track->Eta(), track->Phi()); + fMCQA->SetContainerStep(3); for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE if(!fisNonHFEsystematics)break; @@ -958,8 +963,12 @@ void AliAnalysisTaskHFE::ProcessESD(){ if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue; if(elecSource == iSource){ for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ - if(weightElecBgV0[iLevel]>0){ fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);} - else if(weightElecBgV0[iLevel]<0){ fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);} + if(weightElecBgV0[iLevel]>0){ + fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]); + } + else if(weightElecBgV0[iLevel]<0){ + fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]); + } } break; } @@ -968,15 +977,22 @@ void AliAnalysisTaskHFE::ProcessESD(){ } } //else{ - if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]); - else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]); + if(weightElecBgV0[0]>0) { + fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]); + fVarManager->FillContainer(fContainer, "conversionElecs", 4, kTRUE, weightElecBgV0[0]); + } + else if(weightElecBgV0[0]<0) { + fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]); + fVarManager->FillContainer(fContainer, "mesonElecs", 4, kTRUE, -1*weightElecBgV0[0]); + } //} } } } if(TMath::Abs(track->Eta()) < 0.5){ - fQACollection->Fill("TPCdEdxBeforePID", track->P(), track->GetTPCsignal()); + if(track->GetInnerParam()) + fQACollection->Fill("TPCdEdxBeforePID", track->GetInnerParam()->P(), track->GetTPCsignal()); fQACollection->Fill("TPCnSigmaBeforePID", track->P(), fInputHandler->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron)); } @@ -1002,11 +1018,11 @@ void AliAnalysisTaskHFE::ProcessESD(){ if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls; Double_t itsChi2[7] = {track->Pt(),track->Eta(), track->Phi(), - fCentralityF,track->GetTPCsignalN(), sharebit, itschi2percluster}; + fCentralityF,track->GetTPCsignalN(), sharebit,itschi2percluster}; fQACollection->Fill("fChi2perITScluster", itsChi2); } else{ - + Double_t itschi2percluster = 0.0; Double_t itsnbcls = static_cast(track->GetNcls(0)); if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls; @@ -1024,9 +1040,9 @@ void AliAnalysisTaskHFE::ProcessESD(){ Int_t glabel=TMath::Abs(mctrack->GetMother()); if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==321) - fQACollection->Fill("Ke3Kecorr",mctrackmother->Pt(),mctrack->Pt()); + fQACollection->Fill("Ke3Kecorr",mctrack->Pt(),mctrackmother->Pt()); else if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==130) - fQACollection->Fill("Ke3K0Lecorr",mctrackmother->Pt(),mctrack->Pt()); + fQACollection->Fill("Ke3K0Lecorr",mctrack->Pt(),mctrackmother->Pt()); } } } @@ -1107,42 +1123,41 @@ void AliAnalysisTaskHFE::ProcessESD(){ } // end of electron background analysis - + Int_t sourceDca =-1; if (GetPlugin(kDEstep)) { Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.,}; Int_t elecSource = 0; // minjung for IP QA(temporary ~ 2weeks) Double_t hfeimpactR=0., hfeimpactnsigmaR=0.; fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR); - fQACollection->Fill("dataDca",track->Pt(),hfeimpactR); - fQACollection->Fill("dataDcaSig",track->Pt(),hfeimpactnsigmaR); - fQACollection->Fill("dataDcaSigDca",hfeimpactR,hfeimpactnsigmaR); - if(HasMCData()){ + sourceDca=0; + if(HasMCData()) + { // minjung for IP QA(temporary ~ 2weeks) - if(fSignalCuts->IsCharmElectron(track)){ - fQACollection->Fill("charmDca",track->Pt(),hfeimpactR); - fQACollection->Fill("charmDcaSig",track->Pt(),hfeimpactnsigmaR); + if(fSignalCuts->IsCharmElectron(track)){ + sourceDca=1; } - else if(fSignalCuts->IsBeautyElectron(track)){ - fQACollection->Fill("beautyDca",track->Pt(),hfeimpactR); - fQACollection->Fill("beautyDcaSig",track->Pt(),hfeimpactnsigmaR); + else if(fSignalCuts->IsBeautyElectron(track)){ + sourceDca=2; } - else if(fSignalCuts->IsGammaElectron(track)){ - fQACollection->Fill("conversionDca",track->Pt(),hfeimpactR); - fQACollection->Fill("conversionDcaSig",track->Pt(),hfeimpactnsigmaR); - fQACollection->Fill("conversionDcaSigDca",hfeimpactR,hfeimpactnsigmaR); + else if(fSignalCuts->IsGammaElectron(track)){ + sourceDca=3; } - else if(fSignalCuts->IsNonHFElectron(track)){ - fQACollection->Fill("nonhfeDca",track->Pt(),hfeimpactR); - fQACollection->Fill("nonhfeDcaSig",track->Pt(),hfeimpactnsigmaR); + else if(fSignalCuts->IsNonHFElectron(track)){ + sourceDca=4; } - - if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){ + else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){ + sourceDca=5; fQACollection->Fill("hadronsBeforeIPcut",track->Pt()); fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt()); } - if(fMCQA) { + else { + sourceDca=6; + } + + if(fMCQA && signal) { + fMCQA->SetContainerStep(0); for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE if(!fisNonHFEsystematics)break; @@ -1168,19 +1183,45 @@ void AliAnalysisTaskHFE::ProcessESD(){ } } //else{ - if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]); - else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]); + if(weightElecBgV0[0]>0) { + fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]); + fVarManager->FillContainer(fContainer, "conversionElecs", 5, kTRUE, weightElecBgV0[0]); + } + else if(weightElecBgV0[0]<0) { + fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]); + fVarManager->FillContainer(fContainer, "mesonElecs", 5, kTRUE, -1*weightElecBgV0[0]); + } //} if(bTagged){ // bg estimation for the secondary vertex tagged signals if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 2, kFALSE, weightElecBgV0[0]); else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 2, kFALSE, -1*weightElecBgV0[0]); } } - } + } // end of MC + + dataDca[0]=sourceDca; + dataDca[1]=track->Pt(); + dataDca[2]=hfeimpactR; + dataDca[3]=hfeimpactnsigmaR; + dataDca[4]=fCentralityF; + dataDca[5] = 49; + Double_t xr[3]={49,49,49}; + if(HasMCData()) { + mctrack->XvYvZv(xr); + dataDca[5] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); + } + dataDca[6] = v0pid; + + // printf("Entries dca: [%.3f|%.3f|%.3f|%f|%f]\n", dataDca[0],dataDca[1],dataDca[2],dataDca[3],dataDca[4]); + if (!HasMCData()) fQACollection->Fill("Dca", dataDca); + else if(signal) fQACollection->Fill("Dca", dataDca); + + // Fill Containers for impact parameter analysis if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue; if(HasMCData()){ - if(fMCQA) { + if(fMCQA && signal) { + fMCQA->SetContainerStep(1); for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE if(!fisNonHFEsystematics)break; @@ -1205,8 +1246,14 @@ void AliAnalysisTaskHFE::ProcessESD(){ } } // else{ - if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]); - else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]); + if(weightElecBgV0[0]>0) { + fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]); + fVarManager->FillContainer(fContainer, "conversionElecs", 6, kTRUE, weightElecBgV0[0]); + } + else if(weightElecBgV0[0]<0) { + fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]); + fVarManager->FillContainer(fContainer, "mesonElecs", 6, kTRUE, -1*weightElecBgV0[0]); + } //} } } @@ -1517,8 +1564,8 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){ fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1); if(HasMCData()){ - fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",4); - fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",4); + fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",7); + fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",7); fContainer->Sumw2("conversionElecs"); fContainer->Sumw2("mesonElecs"); @@ -1588,36 +1635,55 @@ void AliAnalysisTaskHFE::InitContaminationQA(){ // // Add QA for Impact Parameter cut // - const Double_t kPtbound[2] = {0.1, 20.}; - Int_t iBin[1]; - iBin[0] = 44; // bins in pt - fQACollection->CreateTH1F("hadronsBeforeIPcut", "Hadrons before IP cut", iBin[0], kPtbound[0], kPtbound[1], 1); - fQACollection->CreateTH1F("hadronsAfterIPcut", "Hadrons after IP cut", iBin[0], kPtbound[0], kPtbound[1], 1); - fQACollection->CreateTH1F("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", iBin[0], kPtbound[0], kPtbound[1], 1); - fQACollection->CreateTH1F("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1); - - fQACollection->CreateTH2F("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1); - fQACollection->CreateTH2F("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ",20,0.,20.,iBin[0],kPtbound[0],kPtbound[1], 1); - fQACollection->CreateTH1F("Kptspectra", "Charged Kaons: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1); - fQACollection->CreateTH1F("K0Lptspectra", "K0L: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1); - - const Double_t kDCAbound[2] = {-5., 5.}; - const Double_t kDCAsigbound[2] = {-50., 50.}; - - fQACollection->CreateTH2F("dataDcaSig", "data dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0); - fQACollection->CreateTH2F("charmDcaSig", "charm dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0); - fQACollection->CreateTH2F("beautyDcaSig", "beauty dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0); - fQACollection->CreateTH2F("conversionDcaSig", "conversion dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0); - fQACollection->CreateTH2F("nonhfeDcaSig", "nonhfe dca significance: dca sig ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAsigbound[0], kDCAsigbound[1], 0); - - fQACollection->CreateTH2F("dataDca", "data dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0); - fQACollection->CreateTH2F("charmDca", "charm dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0); - fQACollection->CreateTH2F("beautyDca", "beauty dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0); - fQACollection->CreateTH2F("conversionDca", "conversion dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0); - fQACollection->CreateTH2F("nonhfeDca", "nonhfe dca : dca ",iBin[0], kPtbound[0], kPtbound[1], 2000,kDCAbound[0], kDCAbound[1], 0); - - fQACollection->CreateTH2F("dataDcaSigDca", "data dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0); - fQACollection->CreateTH2F("conversionDcaSigDca", "conversion dca significance and dca correlation; dca; dca sig ", 2000, kDCAbound[0], kDCAbound[1], 2000, kDCAsigbound[0], kDCAsigbound[1], 0); + + TObjArray *array = fVarManager->GetVariables(); + Int_t nvars = array->GetEntriesFast(); + for(Int_t v = 0; v < nvars; v++) { + AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v); + if(!variable) continue; + TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName()); + if(!name.CompareTo("pt")) { + const Int_t nBinPt = variable->GetNumberOfBins(); + const Double_t *kPtRange = variable->GetBinning(); + + fQACollection->CreateTH1Farray("hadronsBeforeIPcut", "Hadrons before IP cut", nBinPt, kPtRange); + fQACollection->CreateTH1Farray("hadronsAfterIPcut", "Hadrons after IP cut", nBinPt, kPtRange); + fQACollection->CreateTH1Farray("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", nBinPt, kPtRange); + fQACollection->CreateTH1Farray("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", nBinPt, kPtRange); + + fQACollection->CreateTH2Farray("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.); + fQACollection->CreateTH2Farray("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.); + fQACollection->CreateTH1Farray("Kptspectra", "Charged Kaons: MC p_{t} ", nBinPt, kPtRange); + fQACollection->CreateTH1Farray("K0Lptspectra", "K0L: MC p_{t} ", nBinPt, kPtRange); + + const Double_t kDCAbound[2] = {-5., 5.}; + const Double_t kDCAsigbound[2] = {-50., 50.}; + + const Int_t nDimDca=7; + const Int_t nBinDca[nDimDca] = { 8, nBinPt, 2000, 2000, 12, 500, 6}; + Double_t minimaDca[nDimDca] = { -1., 0., kDCAbound[0], kDCAsigbound[0], -1., 0, -1}; + Double_t maximaDca[nDimDca] = { 7., 20., kDCAbound[1], kDCAsigbound[1], 11., 50, 5}; + + Double_t *sourceBins = AliHFEtools::MakeLinearBinning(nBinDca[0], minimaDca[0], maximaDca[0]); + Double_t *dcaBins = AliHFEtools::MakeLinearBinning(nBinDca[2], minimaDca[2], maximaDca[2]); + Double_t *dcaSigBins = AliHFEtools::MakeLinearBinning(nBinDca[3], minimaDca[3], maximaDca[3]); + Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBinDca[4], minimaDca[4], maximaDca[4]); + Double_t *eProdRBins = AliHFEtools::MakeLinearBinning(nBinDca[5], minimaDca[5], maximaDca[5]); + Double_t *v0PIDBins = AliHFEtools::MakeLinearBinning(nBinDca[6], minimaDca[6], maximaDca[6]); + + fQACollection->CreateTHnSparseNoLimits("Dca", "Dca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; dcasig; centrality bin; eProdR; v0pid", nDimDca, nBinDca); + ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(0, sourceBins); + ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(1, kPtRange); + ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(2, dcaBins); + ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(3, dcaSigBins); + ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(4, centralityBins); + ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(5, eProdRBins); + ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(6, v0PIDBins); + + break; + } + } + } //____________________________________________________________ diff --git a/PWGHF/hfe/AliAnalysisTaskHFE.h b/PWGHF/hfe/AliAnalysisTaskHFE.h index 9021a302a36..e10ee959028 100644 --- a/PWGHF/hfe/AliAnalysisTaskHFE.h +++ b/PWGHF/hfe/AliAnalysisTaskHFE.h @@ -133,13 +133,15 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{ Bool_t ReadCentrality(); void RejectionPileUpVertexRangeEventCut(); void SelectSpecialTrigger(const Char_t *trgclust, Int_t runMin = 0, Int_t runMax = 999999999); + void SetDebugStreaming() {SetBit(kTreeStream);}; private: enum{ kHasMCdata = BIT(19), kAODanalysis = BIT(20), kBeamType = BIT(21), - kBackgroundInitialized = BIT(22) + kBackgroundInitialized = BIT(22), + kTreeStream = BIT(23) }; Bool_t FillProductionVertex(const AliVParticle * const track) const; diff --git a/PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx b/PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx index 437424dbc1c..e11df38329f 100644 --- a/PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx +++ b/PWGHF/hfe/AliAnalysisTaskHFEFlow.cxx @@ -46,13 +46,11 @@ #include "AliFlowVector.h" #include "AliFlowCommonConstants.h" - #include "AliHFEcuts.h" #include "AliHFEpid.h" #include "AliHFEpidQAmanager.h" #include "AliHFEtools.h" - - +#include "AliHFEVZEROEventPlane.h" #include "AliCentrality.h" #include "AliEventplane.h" @@ -62,7 +60,7 @@ //____________________________________________________________________ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() : AliAnalysisTaskSE(), - fListHist(), + fListHist(0x0), fVZEROEventPlane(kFALSE), fVZEROEventPlaneA(kFALSE), fVZEROEventPlaneC(kFALSE), @@ -83,6 +81,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() : fPrecisionPhi(0.001), fUseMCReactionPlane(kFALSE), fMCPID(kFALSE), + fNoPID(kFALSE), fDebugLevel(0), fcutsRP(0), fcutsPOI(0), @@ -90,6 +89,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() : fPID(0), fPIDqa(0), fflowEvent(0x0), + fHFEVZEROEventPlane(0x0), fHistEV(0), fEventPlane(0x0), fEventPlaneaftersubtraction(0x0), @@ -100,10 +100,12 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() : fSin2phiep(0x0), fSin2phiephiep(0x0), fCosResabc(0x0), + fSinResabc(0x0), fProfileCosResab(0x0), fProfileCosResac(0x0), fProfileCosResbc(0x0), fCosRes(0x0), + fSinRes(0x0), fProfileCosRes(0x0), fDeltaPhiMaps(0x0), fCosPhiMaps(0x0), @@ -119,7 +121,7 @@ AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow() : //______________________________________________________________________________ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) : AliAnalysisTaskSE(name), - fListHist(), + fListHist(0x0), fVZEROEventPlane(kFALSE), fVZEROEventPlaneA(kFALSE), fVZEROEventPlaneC(kFALSE), @@ -140,6 +142,7 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) : fPrecisionPhi(0.001), fUseMCReactionPlane(kFALSE), fMCPID(kFALSE), + fNoPID(kFALSE), fDebugLevel(0), fcutsRP(0), fcutsPOI(0), @@ -147,6 +150,7 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) : fPID(0), fPIDqa(0), fflowEvent(0x0), + fHFEVZEROEventPlane(0x0), fHistEV(0), fEventPlane(0x0), fEventPlaneaftersubtraction(0x0), @@ -157,10 +161,12 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) : fSin2phiep(0x0), fSin2phiephiep(0x0), fCosResabc(0x0), + fSinResabc(0x0), fProfileCosResab(0x0), fProfileCosResac(0x0), fProfileCosResbc(0x0), fCosRes(0x0), + fSinRes(0x0), fProfileCosRes(0x0), fDeltaPhiMaps(0x0), fCosPhiMaps(0x0), @@ -188,6 +194,147 @@ AliAnalysisTaskHFEFlow:: AliAnalysisTaskHFEFlow(const char *name) : DefineOutput(bincless+2,AliFlowEventSimple::Class()); } +} +//____________________________________________________________ +AliAnalysisTaskHFEFlow::AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref): + AliAnalysisTaskSE(ref), + fListHist(0x0), + fVZEROEventPlane(ref.fVZEROEventPlane), + fVZEROEventPlaneA(ref.fVZEROEventPlaneA), + fVZEROEventPlaneC(ref.fVZEROEventPlaneC), + fSubEtaGapTPC(ref.fSubEtaGapTPC), + fEtaGap(ref.fEtaGap), + fNbBinsCentralityQCumulant(ref.fNbBinsCentralityQCumulant), + fNbBinsPtQCumulant(ref.fNbBinsPtQCumulant), + fMinPtQCumulant(ref.fMinPtQCumulant), + fMaxPtQCumulant(ref.fMaxPtQCumulant), + fAfterBurnerOn(ref.fAfterBurnerOn), + fNonFlowNumberOfTrackClones(ref.fNonFlowNumberOfTrackClones), + fV1(ref.fV1), + fV2(ref.fV2), + fV3(ref.fV3), + fV4(ref.fV4), + fV5(ref.fV5), + fMaxNumberOfIterations(ref.fMaxNumberOfIterations), + fPrecisionPhi(ref.fPrecisionPhi), + fUseMCReactionPlane(ref.fUseMCReactionPlane), + fMCPID(ref.fMCPID), + fNoPID(ref.fNoPID), + fDebugLevel(ref.fDebugLevel), + fcutsRP(0), + fcutsPOI(0), + fHFECuts(0), + fPID(0), + fPIDqa(0), + fflowEvent(0x0), + fHFEVZEROEventPlane(0x0), + fHistEV(0), + fEventPlane(0x0), + fEventPlaneaftersubtraction(0x0), + fCosSin2phiep(0x0), + fCos2phie(0x0), + fSin2phie(0x0), + fCos2phiep(0x0), + fSin2phiep(0x0), + fSin2phiephiep(0x0), + fCosResabc(0x0), + fSinResabc(0x0), + fProfileCosResab(0x0), + fProfileCosResac(0x0), + fProfileCosResbc(0x0), + fCosRes(0x0), + fSinRes(0x0), + fProfileCosRes(0x0), + fDeltaPhiMaps(0x0), + fCosPhiMaps(0x0), + fProfileCosPhiMaps(0x0) +{ + // + // Copy Constructor + // + ref.Copy(*this); +} + +//____________________________________________________________ +AliAnalysisTaskHFEFlow &AliAnalysisTaskHFEFlow::operator=(const AliAnalysisTaskHFEFlow &ref){ + // + // Assignment operator + // + if(this == &ref) + ref.Copy(*this); + return *this; +} + +//____________________________________________________________ +void AliAnalysisTaskHFEFlow::Copy(TObject &o) const { + // + // Copy into object o + // + AliAnalysisTaskHFEFlow &target = dynamic_cast(o); + target.fVZEROEventPlane = fVZEROEventPlane; + target.fVZEROEventPlaneA = fVZEROEventPlaneA; + target.fVZEROEventPlaneC = fVZEROEventPlaneC; + target.fSubEtaGapTPC = fSubEtaGapTPC; + target.fEtaGap = fEtaGap; + target.fNbBinsCentralityQCumulant = fNbBinsCentralityQCumulant; + target.fNbBinsPtQCumulant = fNbBinsPtQCumulant; + target.fMinPtQCumulant = fMinPtQCumulant; + target.fMaxPtQCumulant = fMaxPtQCumulant; + target.fAfterBurnerOn = fAfterBurnerOn; + target.fNonFlowNumberOfTrackClones = fNonFlowNumberOfTrackClones; + target.fV1 = fV1; + target.fV2 = fV2; + target.fV3 = fV3; + target.fV4 = fV4; + target.fV5 = fV5; + target.fMaxNumberOfIterations = fMaxNumberOfIterations; + target.fPrecisionPhi = fPrecisionPhi; + target.fUseMCReactionPlane = fUseMCReactionPlane; + target.fMCPID = fMCPID; + target.fNoPID = fNoPID; + target.fDebugLevel = fDebugLevel; + target.fcutsRP = fcutsRP; + target.fcutsPOI = fcutsPOI; + target.fHFECuts = fHFECuts; + target.fPID = fPID; + target.fPIDqa = fPIDqa; + target.fHFEVZEROEventPlane = fHFEVZEROEventPlane; + +} +//____________________________________________________________ +AliAnalysisTaskHFEFlow::~AliAnalysisTaskHFEFlow(){ + // + // Destructor + // + if(fListHist) delete fListHist; + if(fcutsRP) delete fcutsRP; + if(fcutsPOI) delete fcutsPOI; + if(fHFECuts) delete fHFECuts; + if(fPID) delete fPID; + if(fPIDqa) delete fPIDqa; + if(fflowEvent) delete fflowEvent; + if(fHFEVZEROEventPlane) delete fHFEVZEROEventPlane; + if(fHistEV) delete fHistEV; + if(fEventPlane) delete fEventPlane; + if(fEventPlaneaftersubtraction) delete fEventPlaneaftersubtraction; + if(fCosSin2phiep) delete fCosSin2phiep; + if(fCos2phie) delete fCos2phie; + if(fSin2phie) delete fSin2phie; + if(fCos2phiep) delete fCos2phiep; + if(fSin2phiep) delete fSin2phiep; + if(fSin2phiephiep) delete fSin2phiephiep; + if(fCosResabc) delete fCosResabc; + if(fSinResabc) delete fSinResabc; + if(fProfileCosResab) delete fProfileCosResab; + if(fProfileCosResac) delete fProfileCosResac; + if(fProfileCosResbc) delete fProfileCosResbc; + if(fCosRes) delete fCosRes; + if(fSinRes) delete fSinRes; + if(fProfileCosRes) delete fProfileCosRes; + if(fDeltaPhiMaps) delete fDeltaPhiMaps; + if(fCosPhiMaps) delete fCosPhiMaps; + if(fProfileCosPhiMaps) delete fProfileCosPhiMaps; + } //________________________________________________________________________ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects() @@ -391,6 +538,15 @@ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects() fCosResabc->SetBinEdges(3,binLimCMore); fCosResabc->Sumw2(); + const Int_t nDimfbiss=4; + Int_t nBinfbiss[nDimfbiss] = {nBinsCos,nBinsCos,nBinsCos,nBinsC}; + fSinResabc = new THnSparseF("SinRes_abc","SinRes_abc",nDimfbiss,nBinfbiss); + fSinResabc->SetBinEdges(0,binLimCos); + fSinResabc->SetBinEdges(1,binLimCos); + fSinResabc->SetBinEdges(2,binLimCos); + fSinResabc->SetBinEdges(3,binLimC); + fSinResabc->Sumw2(); + // Profile cosres centrality with 3 subevents fProfileCosResab = new TProfile("ProfileCosRes_a_b","ProfileCosRes_a_b",nBinsCMore,binLimCMore); fProfileCosResab->Sumw2(); @@ -407,6 +563,13 @@ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects() fCosRes->SetBinEdges(1,binLimCMore); fCosRes->Sumw2(); + const Int_t nDimff=2; + Int_t nBinff[nDimff] = {nBinsCos, nBinsC}; + fSinRes = new THnSparseF("SinRes","SinRes",nDimff,nBinff); + fSinRes->SetBinEdges(0,binLimCos); + fSinRes->SetBinEdges(1,binLimC); + fSinRes->Sumw2(); + // Profile cosres centrality fProfileCosRes = new TProfile("ProfileCosRes","ProfileCosRes",nBinsCMore,binLimCMore); fProfileCosRes->Sumw2(); @@ -455,9 +618,13 @@ void AliAnalysisTaskHFEFlow::UserCreateOutputObjects() fListHist->Add(fSin2phiephiep); fListHist->Add(fCosRes); fListHist->Add(fCosResabc); + fListHist->Add(fSinRes); + fListHist->Add(fSinResabc); fListHist->Add(fDeltaPhiMaps); fListHist->Add(fCosPhiMaps); fListHist->Add(fProfileCosPhiMaps); + + if(fHFEVZEROEventPlane && (fDebugLevel > 2)) fListHist->Add(fHFEVZEROEventPlane->GetOutputList()); PostData(1, fListHist); @@ -487,7 +654,9 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) Double_t valuensparseabis[5]; Double_t valuensparsee[4]; Double_t valuensparsef[2]; + Double_t valuensparsefsin[2]; Double_t valuensparsefbis[4]; + Double_t valuensparsefbissin[4]; Double_t valuensparseg[3]; Double_t valuensparseh[3]; Double_t valuensparsehprofile[3]; @@ -542,14 +711,7 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) if((95.0< cntr) && (cntr<100.0)) binctMore = 19.5; binctLess = cntr; - - //for(Int_t bincless = 0; bincless < fNbBinsCentralityQCumulant; bincless++) { - //if((fBinCentralityLess[bincless]< cntr) && (cntr < fBinCentralityLess[bincless+1])) { - // binctLess = bincless+0.5; - //printf("fBinCentralityLess[bincless] %f and binctLess %f\n",fBinCentralityLess[bincless],binctLess); - //} - //} - + if(binct > 11.0) return; @@ -558,7 +720,9 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) valuensparseabis[1] = binct; valuensparsee[1] = binct; valuensparsef[1] = binctMore; + valuensparsefsin[1] = binct; valuensparsefbis[3] = binctMore; + valuensparsefbissin[3] = binct; valuensparseg[1] = binct; valuensparseh[1] = binct; valuensparsehprofile[1] = binctLess; @@ -615,14 +779,41 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) TVector2 *qsub2a = 0x0; // V0 - - eventPlaneV0 = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0", esd,2)); - if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi(); - eventPlaneV0A = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0A", esd,2)); - if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi(); - eventPlaneV0C = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0C", esd,2)); - if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi(); + + if(fHFEVZEROEventPlane){ + + fHFEVZEROEventPlane->ProcessEvent(esd); + + if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0A()+100) < 0.0000001) eventPlaneV0A = -100.0; + else { + eventPlaneV0A = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0A()); + if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi(); + } + + if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0C()+100) < 0.0000001) eventPlaneV0C = -100.0; + else { + eventPlaneV0C = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0C()); + if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi(); + } + + if(TMath::Abs(fHFEVZEROEventPlane->GetEventPlaneV0()+100) < 0.0000001) eventPlaneV0 = -100.0; + else { + eventPlaneV0 = TVector2::Phi_0_2pi(fHFEVZEROEventPlane->GetEventPlaneV0()); + if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi(); + } + + } + else { + + eventPlaneV0 = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0", esd,2)); + if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi(); + eventPlaneV0A = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0A", esd,2)); + if(eventPlaneV0A > TMath::Pi()) eventPlaneV0A = eventPlaneV0A - TMath::Pi(); + eventPlaneV0C = TVector2::Phi_0_2pi(esdEPa->GetEventplane("V0C", esd,2)); + if(eventPlaneV0C > TMath::Pi()) eventPlaneV0C = eventPlaneV0C - TMath::Pi(); + } + // TPC standardQ = esdEPa->GetQVector(); @@ -649,21 +840,32 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) Float_t eventPlanesub1a = -100.0; Float_t eventPlanesub2a = -100.0; Double_t diffsub1sub2a = -100.0; + Double_t diffsub1sub2asin = -100.0; Double_t diffsubasubb = -100.0; Double_t diffsubasubc = -100.0; Double_t diffsubbsubc = -100.0; + Double_t diffsubasubbsin = -100.0; + Double_t diffsubasubcsin = -100.0; + Double_t diffsubbsubcsin = -100.0; //if(fVZEROEventPlane) { diffsubasubb = TMath::Cos(2.*(eventPlaneV0A - eventPlaneV0C)); diffsubasubc = TMath::Cos(2.*(eventPlaneV0A - eventPlaneTPC)); diffsubbsubc = TMath::Cos(2.*(eventPlaneV0C - eventPlaneTPC)); + + diffsubasubbsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneV0C)); + diffsubasubcsin = TMath::Sin(2.*(eventPlaneV0A - eventPlaneTPC)); + diffsubbsubcsin = TMath::Sin(2.*(eventPlaneV0C - eventPlaneTPC)); //} //else { qsub1a = esdEPa->GetQsub1(); qsub2a = esdEPa->GetQsub2(); if(qsub1a) eventPlanesub1a = TVector2::Phi_0_2pi(qsub1a->Phi())/2.; if(qsub2a) eventPlanesub2a = TVector2::Phi_0_2pi(qsub2a->Phi())/2.; - if(qsub1a && qsub2a) diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.)); + if(qsub1a && qsub2a) { + diffsub1sub2a = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.)); + diffsub1sub2asin = TMath::Sin(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.)); + } //} @@ -736,7 +938,9 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) if(!fVZEROEventPlane) { valuensparsef[0] = diffsub1sub2a; fCosRes->Fill(&valuensparsef[0]); - if(fDebugLevel > 0) { + valuensparsefsin[0] = diffsub1sub2asin; + fSinRes->Fill(&valuensparsefsin[0]); + if(fDebugLevel > 1) { fProfileCosRes->Fill(valuensparsef[1],valuensparsef[0]); } } @@ -745,7 +949,11 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) valuensparsefbis[1] = diffsubbsubc; valuensparsefbis[2] = diffsubasubc; fCosResabc->Fill(&valuensparsefbis[0]); - if(fDebugLevel > 0) { + valuensparsefbissin[0] = diffsubasubbsin; + valuensparsefbissin[1] = diffsubbsubcsin; + valuensparsefbissin[2] = diffsubasubcsin; + fSinResabc->Fill(&valuensparsefbissin[0]); + if(fDebugLevel > 1) { fProfileCosResab->Fill(valuensparsefbis[3],valuensparsefbis[0]); fProfileCosResac->Fill(valuensparsefbis[3],valuensparsefbis[1]); fProfileCosResbc->Fill(valuensparsefbis[3],valuensparsefbis[2]); @@ -775,23 +983,26 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) if(!survived) continue; - // Apply PID for Data - if(!fMCPID) { - AliHFEpidObject hfetrack; - hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis); - hfetrack.SetRecTrack(track); - hfetrack.SetCentrality((Int_t)binct); - //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality()); - hfetrack.SetPbPb(); - if(!fPID->IsSelected(&hfetrack,0x0,"",fPIDqa)) { - continue; + // Apply PID + if(!fNoPID) { + // Apply PID for Data + if(!fMCPID) { + AliHFEpidObject hfetrack; + hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis); + hfetrack.SetRecTrack(track); + hfetrack.SetCentrality((Int_t)binct); + //printf("centrality %f and %d\n",binct,hfetrack.GetCentrality()); + hfetrack.SetPbPb(); + if(!fPID->IsSelected(&hfetrack,0x0,"",fPIDqa)) { + continue; + } + } + else { + if(!mcEvent) continue; + if(!(mctrack = dynamic_cast(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue; + //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode())); + if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue; } - } - else { - if(!mcEvent) continue; - if(!(mctrack = dynamic_cast(mcEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue; - //printf("PdgCode %d\n",TMath::Abs(mctrack->Particle()->GetPdgCode())); - if(TMath::Abs(mctrack->Particle()->GetPdgCode())!=11) continue; } @@ -894,7 +1105,7 @@ void AliAnalysisTaskHFEFlow::UserExec(Option_t */*option*/) if(fillEventPlane) fEventPlaneaftersubtraction->Fill(&valuensparseabis[0]); - if(fDebugLevel > 1) + if(fDebugLevel > 3) { // valuensparsee[0] = TMath::Cos(2*phitrack); diff --git a/PWGHF/hfe/AliAnalysisTaskHFEFlow.h b/PWGHF/hfe/AliAnalysisTaskHFEFlow.h index ebf989506d3..7850d00ce07 100644 --- a/PWGHF/hfe/AliAnalysisTaskHFEFlow.h +++ b/PWGHF/hfe/AliAnalysisTaskHFEFlow.h @@ -36,25 +36,31 @@ class TProfile2D; class THnSparse; class AliHFEpidQAmanager; class AliFlowEvent; +class AliHFEVZEROEventPlane; class AliAnalysisTaskHFEFlow: public AliAnalysisTaskSE { public: AliAnalysisTaskHFEFlow(); AliAnalysisTaskHFEFlow(const char *name); - //TODO: if we get the ESDpid object from somewhere else, we should not delete it anymore!!! - virtual ~AliAnalysisTaskHFEFlow(){;} + AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow &ref); + AliAnalysisTaskHFEFlow& operator=(const AliAnalysisTaskHFEFlow &ref); + virtual void Copy(TObject &o) const; + virtual ~AliAnalysisTaskHFEFlow(); virtual void UserExec(Option_t */*option*/); virtual void UserCreateOutputObjects(); AliHFEpid *GetPID() const { return fPID; } + AliHFEpidQAmanager *GetPIDQAManager() const { return fPIDqa; } + void SetHFECuts(AliHFEcuts * const cuts) { fHFECuts = cuts; }; void SetSubEtaGapTPC(Bool_t subEtaGapTPC) { fSubEtaGapTPC = subEtaGapTPC; }; void SetEtaGap(Double_t etaGap) { fEtaGap = etaGap; }; void SetVZEROEventPlane(Bool_t vzeroEventPlane) { fVZEROEventPlane = vzeroEventPlane; }; void SetVZEROEventPlaneA(Bool_t vzeroEventPlaneA) { fVZEROEventPlaneA = vzeroEventPlaneA; }; void SetVZEROEventPlaneC(Bool_t vzeroEventPlaneC) { fVZEROEventPlaneC = vzeroEventPlaneC; }; + void SetHFEVZEROEventPlane(AliHFEVZEROEventPlane *hfeVZEROEventPlane) { fHFEVZEROEventPlane = hfeVZEROEventPlane; }; void SetNbBinsCentralityQCumulant(Int_t nbBinsCentralityQCumulant) { fNbBinsCentralityQCumulant = nbBinsCentralityQCumulant; }; void SetBinCentralityLess(Int_t k, Float_t value) { fBinCentralityLess[k] = value; }; @@ -69,6 +75,7 @@ public: void SetPrecisionPhi(Double_t precisionPhi) { fPrecisionPhi = precisionPhi;}; void SetUseMCReactionPlane(Bool_t useMCReactionPlane) { fUseMCReactionPlane = useMCReactionPlane;}; void SetMCPID(Bool_t mcPID) { fMCPID = mcPID;}; + void SetNoPID(Bool_t noPID) { fNoPID = noPID;}; void SetDebugLevel(Int_t debugLevel) { fDebugLevel = debugLevel;}; @@ -79,21 +86,21 @@ public: Double_t GetPhiAfterAddV2(Double_t phi,Double_t reactionPlaneAngle) const; private: - TList *fListHist; //TH list + TList *fListHist; //! TH list - Bool_t fVZEROEventPlane; // Use Event Planes from VZERO + Bool_t fVZEROEventPlane; // Use Event Planes from VZERO Bool_t fVZEROEventPlaneA; // Use Event Planes from VZERO A Bool_t fVZEROEventPlaneC; // Use Event Planes from VZERO C - Bool_t fSubEtaGapTPC; // bool to fill with eta gap - Double_t fEtaGap; // Value of the eta gap + Bool_t fSubEtaGapTPC; // bool to fill with eta gap + Double_t fEtaGap; // Value of the eta gap Int_t fNbBinsCentralityQCumulant; // Number of Bins Q Cumulant - Double_t fBinCentralityLess[10]; // Centrality Bin lower value - Int_t fNbBinsPtQCumulant; // Nbbinspt QCumulant method - Double_t fMinPtQCumulant; // Min pt QCumulant method - Double_t fMaxPtQCumulant; // Max pt QCumulant method - Bool_t fAfterBurnerOn; // Add flow to all tracks + Double_t fBinCentralityLess[10]; // Centrality Bin lower value + Int_t fNbBinsPtQCumulant; // Nbbinspt QCumulant method + Double_t fMinPtQCumulant; // Min pt QCumulant method + Double_t fMaxPtQCumulant; // Max pt QCumulant method + Bool_t fAfterBurnerOn; // Add flow to all tracks Int_t fNonFlowNumberOfTrackClones; // number of times to clone the particles (nonflow) Double_t fV1; // Add Flow. Must be in range [0,0.5]. Double_t fV2; // Add Flow. Must be in range [0,0.5]. @@ -105,6 +112,7 @@ private: Bool_t fUseMCReactionPlane; // use MC reaction plane Bool_t fMCPID; // MC PID for electrons + Bool_t fNoPID; // No PID for checks Int_t fDebugLevel; // Debug Level @@ -113,53 +121,54 @@ private: AliFlowTrackCuts* fcutsPOI; // Particle Of Interest cut // Cuts for HFE - AliHFEcuts *fHFECuts; // HFE cuts - AliHFEpid *fPID; // PID cuts - AliHFEpidQAmanager *fPIDqa; // QA Manager - AliFlowEvent *fflowEvent; // Flow event - + AliHFEcuts *fHFECuts; // HFE cuts + AliHFEpid *fPID; // PID cuts + AliHFEpidQAmanager *fPIDqa; // QA Manager + AliFlowEvent *fflowEvent; //! Flow event + // VZERO Event plane after calibration 2010 + AliHFEVZEROEventPlane *fHFEVZEROEventPlane; // VZERO event plane calibrated + // Histos - TH2D *fHistEV; // Number of events + TH2D *fHistEV; //! Number of events // A Event plane as function of phiepa, phiepb, phiepc, phiepd centrality // a V0A, b V0C, c TPC, d V0 - THnSparseF *fEventPlane; // Event plane + THnSparseF *fEventPlane; //! Event plane // B Event Plane after subtraction as function of phiep, centrality - THnSparseF *fEventPlaneaftersubtraction; // Event plane + THnSparseF *fEventPlaneaftersubtraction; //! Event plane // Monitoring Event plane: cos2phi, sin2phi, centrality - THnSparseF *fCosSin2phiep; // Cos(2phi), Sin(2phi) + THnSparseF *fCosSin2phiep; //! Cos(2phi), Sin(2phi) // E Monitoring Event plane after subtraction of the track: cos, centrality, pt, eta - THnSparseF *fCos2phie; // Monitoring - THnSparseF *fSin2phie; // Monitoring - THnSparseF *fCos2phiep; // Monitoring - THnSparseF *fSin2phiep; // Monitoring - THnSparseF *fSin2phiephiep; // Monitoring + THnSparseF *fCos2phie; //! Monitoring + THnSparseF *fSin2phie; //! Monitoring + THnSparseF *fCos2phiep; //! Monitoring + THnSparseF *fSin2phiep; //! Monitoring + THnSparseF *fSin2phiephiep; //! Monitoring // Fbis Resolution as function of cosres, cosres, cosres, centrality for three subevents (V0) // a V0A, b V0C, c TPC - THnSparseF *fCosResabc; // Res - TProfile *fProfileCosResab; // Profile Res_a_b - TProfile *fProfileCosResac; // Profile Res_a_c - TProfile *fProfileCosResbc; // Profile Res_b_c + THnSparseF *fCosResabc; //! Res + THnSparseF *fSinResabc; //! Res + TProfile *fProfileCosResab; //! Profile Res_a_b + TProfile *fProfileCosResac; //! Profile Res_a_c + TProfile *fProfileCosResbc; //! Profile Res_b_c // F Resolution as function of cosres, centrality for two subevents (TPC) - THnSparseF *fCosRes; // Res - TProfile *fProfileCosRes; // Profile Res + THnSparseF *fCosRes; //! Res + THnSparseF *fSinRes; //! Res + TProfile *fProfileCosRes; //! Profile Res // G Maps delta phi as function of deltaphi, centrality, pt - THnSparseF *fDeltaPhiMaps; // Delta phi + THnSparseF *fDeltaPhiMaps; //! Delta phi // H Maps cos phi : cos, centrality, pt - THnSparseF *fCosPhiMaps; //Cos - TProfile2D *fProfileCosPhiMaps; // Profile Cos - + THnSparseF *fCosPhiMaps; //! Cos + TProfile2D *fProfileCosPhiMaps; //! Profile Cos - AliAnalysisTaskHFEFlow(const AliAnalysisTaskHFEFlow&); // not implemented - AliAnalysisTaskHFEFlow& operator=(const AliAnalysisTaskHFEFlow&); // not implemented ClassDef(AliAnalysisTaskHFEFlow, 1); // analysisclass }; diff --git a/PWGHF/hfe/AliHFEVZEROEventPlane.cxx b/PWGHF/hfe/AliHFEVZEROEventPlane.cxx new file mode 100644 index 00000000000..3580da9438a --- /dev/null +++ b/PWGHF/hfe/AliHFEVZEROEventPlane.cxx @@ -0,0 +1,407 @@ +#include "AliHFEVZEROEventPlane.h" + +// AliRoot includes +#include "AliAnalysisManager.h" +#include "AliInputEventHandler.h" +#include "AliVEvent.h" +#include "AliESDEvent.h" +#include "AliESDtrack.h" +#include "AliCentrality.h" +#include "AliESDVZERO.h" +#include "TFile.h" +#include "AliOADBContainer.h" +#include "TH2F.h" +#include "TF1.h" +#include "TList.h" +#include "TChain.h" +#include "THnSparse.h" +#include "TString.h" +#include "TVector2.h" +#include "AliEventplane.h" + + +// STL includes +#include +using namespace std; + + +//_____________________________________________________________________________ +AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(): + TNamed(), + fEventPlaneV0A(-100.0), + fEventPlaneV0C(-100.0), + fEventPlaneV0(-100.0), + fESD(0), + fRun(-1), + fMultV0(0x0), + fV0Cpol(100), + fV0Apol(100), + fnamefile(TString("")), + fOutputList(0x0), + fMultV0Before(0x0), + fMultV0After(0x0) +{ + // + // Default constructor (should not be used) + // + + for(Int_t k = 0; k < fgknCentrBin; k++) { + for(Int_t iside = 0; iside < 2; iside++) { + for(Int_t icoord = 0; icoord < 2; icoord++) { + fQBefore[k][iside][icoord] = 0x0; + fQAfter[k][iside][icoord] = 0x0; + } + } + } + +} + +//______________________________________________________________________________ +AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const char *name, const Char_t *title): + TNamed(name,title), + fEventPlaneV0A(-100.0), + fEventPlaneV0C(-100.0), + fEventPlaneV0(-100.0), + fESD(0), + fRun(-1), + fMultV0(0x0), + fV0Cpol(100), + fV0Apol(100), + fnamefile(TString("")), + fOutputList(0x0), + fMultV0Before(0x0), + fMultV0After(0x0) +{ + // + // Constructor + // + + char namelist[100]; + sprintf(namelist,"QA_hist_%s",name); + fOutputList = new TList(); + fOutputList->SetName((const char*)namelist); + fOutputList->SetOwner(kTRUE); + + // Multiplicity + fMultV0Before = new TProfile("MultV0Before","",64,0,64); + fMultV0Before->Sumw2(); + fMultV0After = new TProfile("MultV0After","",64,0,64); + fMultV0After->Sumw2(); + fOutputList->Add(fMultV0Before); + fOutputList->Add(fMultV0After); + + // Recentering + char namecontbefore[100]; + char namecontafter[100]; + for(Int_t k = 0; k < fgknCentrBin; k++) { + for(Int_t iside = 0; iside < 2; iside++) { + for(Int_t icoord = 0; icoord < 2; icoord++) { + + if(iside==0 && icoord==0) + sprintf(namecontbefore,"hQxc2_%i",k); + else if(iside==1 && icoord==0) + sprintf(namecontbefore,"hQxa2_%i",k); + else if(iside==0 && icoord==1) + sprintf(namecontbefore,"hQyc2_%i",k); + else if(iside==1 && icoord==1) + sprintf(namecontbefore,"hQya2_%i",k); + // + if(iside==0 && icoord==0) + sprintf(namecontafter,"hQxc2_%i_after",k); + else if(iside==1 && icoord==0) + sprintf(namecontafter,"hQxa2_%i_after",k); + else if(iside==0 && icoord==1) + sprintf(namecontafter,"hQyc2_%i_after",k); + else if(iside==1 && icoord==1) + sprintf(namecontafter,"hQya2_%i_after",k); + + // + fQBefore[k][iside][icoord] = new TH1F(((const char*)namecontbefore),"",800,-400.0,400.0); + fQBefore[k][iside][icoord]->Sumw2(); + fOutputList->Add(fQBefore[k][iside][icoord]); + fQAfter[k][iside][icoord] = new TH1F(((const char*)namecontafter),"",800,-400.0,400.0); + fQAfter[k][iside][icoord]->Sumw2(); + fOutputList->Add(fQAfter[k][iside][icoord]); + // + } + } + } + +} +//____________________________________________________________ +AliHFEVZEROEventPlane::AliHFEVZEROEventPlane(const AliHFEVZEROEventPlane &ref): + TNamed(ref), + fEventPlaneV0A(ref.fEventPlaneV0A), + fEventPlaneV0C(ref.fEventPlaneV0C), + fEventPlaneV0(ref.fEventPlaneV0), + fESD(0), + fRun(-1), + fMultV0(0x0), + fV0Cpol(100), + fV0Apol(100), + fnamefile(ref.fnamefile), + fOutputList(0x0), + fMultV0Before(0x0), + fMultV0After(0x0) +{ + // + // Copy Constructor + // + ref.Copy(*this); +} + +//____________________________________________________________ +AliHFEVZEROEventPlane &AliHFEVZEROEventPlane::operator=(const AliHFEVZEROEventPlane &ref){ + // + // Assignment operator + // + if(this == &ref) + ref.Copy(*this); + return *this; +} + +//____________________________________________________________ +void AliHFEVZEROEventPlane::Copy(TObject &o) const { + // + // Copy into object o + // + AliHFEVZEROEventPlane &target = dynamic_cast(o); + + target.fEventPlaneV0A = fEventPlaneV0A; + target.fEventPlaneV0C = fEventPlaneV0C; + target.fEventPlaneV0 = fEventPlaneV0; + target.fnamefile = fnamefile; + +} + +//__________________________________________________________________ +AliHFEVZEROEventPlane::~AliHFEVZEROEventPlane(){ + // + // Destruktor + // + if(fOutputList){ + fOutputList->Delete(); + delete fOutputList; + } + fOutputList = 0x0; + if(fMultV0Before) delete fMultV0Before; + if(fMultV0After) delete fMultV0After; + + for(Int_t k = 0; k < fgknCentrBin; k++) { + for(Int_t iside = 0; iside < 2; iside++) { + for(Int_t icoord = 0; icoord < 2; icoord++) { + if(fQBefore[k][iside][icoord]) delete fQBefore[k][iside][icoord]; + if(fQAfter[k][iside][icoord]) delete fQAfter[k][iside][icoord]; + } + } + } + +} +//______________________________________________________________________________ +void AliHFEVZEROEventPlane::ProcessEvent(AliESDEvent *event) +{ + // + // Process the event + // + + // Reset + fEventPlaneV0A = -100.0; + fEventPlaneV0C = -100.0; + fEventPlaneV0 = -100.0; + + // + Int_t run = event->GetRunNumber(); + Bool_t ok = kTRUE; + if(run != fRun){ + if(fMultV0) delete fMultV0; + fMultV0 = 0x0; + // Load the calibrations run dependent + if(!OpenInfoCalbration(run)) ok = kFALSE; + fRun=run; + } + + if(ok) Analyze(event); + +} + +//________________________________________________________________________ +void AliHFEVZEROEventPlane::Analyze(AliESDEvent* esdEvent) +{ + // + // Do VZERO calibration + centering + // + + if(!fMultV0) { + //printf("Did not find calibration VZERO\n"); + return; + } + //printf("HERE!!!\n"); + + //Centrality + Float_t v0Centr = -10.; + AliCentrality* centrality = esdEvent->GetCentrality(); + if (centrality){ + v0Centr = centrality->GetCentralityPercentile("V0M"); + } + + // Analyse only for 0-80% PbPb collisions + Int_t iC = -1; + if (v0Centr >0 && v0Centr < 80){ + if(v0Centr < 5) iC = 0; + else if(v0Centr < 10) iC = 1; + else if(v0Centr < 20) iC = 2; + else if(v0Centr < 30) iC = 3; + else if(v0Centr < 40) iC = 4; + else if(v0Centr < 50) iC = 5; + else if(v0Centr < 60) iC = 6; + else if(v0Centr < 70) iC = 7; + else iC = 8; + + //general + Double_t qxa2 = 0, qya2 = 0; + Double_t qxc2 = 0, qyc2 = 0; + + //V0 info + AliESDVZERO* esdV0 = esdEvent->GetVZEROData(); + + for (Int_t iv0 = 0; iv0 < 64; iv0++) { + Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8); + Float_t multv0 = esdV0->GetMultiplicity(iv0); + if (iv0 < 32){ + //printf("Value %f\n",fMultV0->GetBinContent(iv0+1)); + if(fMultV0->GetBinContent(iv0+1)>0.0) { + qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1); + qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1); + fMultV0After->Fill(iv0,multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1)); + } + } else { + if(fMultV0->GetBinContent(iv0+1)>0.0) { + qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1); + qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1); + fMultV0After->Fill(iv0,multv0*fV0Apol/fMultV0->GetBinContent(iv0+1)); + } + } + fMultV0Before->Fill(iv0,multv0); + } + + // Fill histos + fQBefore[iC][1][0] -> Fill(qxa2); + fQBefore[iC][1][1] -> Fill(qya2); + fQBefore[iC][0][0] -> Fill(qxc2); + fQBefore[iC][0][1] -> Fill(qyc2); + + //grab for each centrality the proper histo with the Qx and Qy to do the recentering + Double_t qxamean2 = fMeanQ[iC][1][0]; + Double_t qxarms2 = fWidthQ[iC][1][0]; + Double_t qyamean2 = fMeanQ[iC][1][1]; + Double_t qyarms2 = fWidthQ[iC][1][1]; + + Double_t qxcmean2 = fMeanQ[iC][0][0]; + Double_t qxcrms2 = fWidthQ[iC][0][0]; + Double_t qycmean2 = fMeanQ[iC][0][1]; + Double_t qycrms2 = fWidthQ[iC][0][1]; + + if((TMath::Abs(qxarms2) < 0.00000001) || (TMath::Abs(qyarms2) < 0.00000001) || (TMath::Abs(qxcrms2) < 0.00000001) || (TMath::Abs(qycrms2) < 0.00000001)) return; + + Double_t qxaCor2 = (qxa2 - qxamean2)/qxarms2; + Double_t qyaCor2 = (qya2 - qyamean2)/qyarms2; + Double_t qxcCor2 = (qxc2 - qxcmean2)/qxcrms2; + Double_t qycCor2 = (qyc2 - qycmean2)/qycrms2; + + Double_t qxCor2 = qxaCor2 + qxcCor2; + Double_t qyCor2 = qyaCor2 + qycCor2; + + // Fill histos + fQAfter[iC][1][0] -> Fill(qxaCor2); + fQAfter[iC][1][1] -> Fill(qyaCor2); + fQAfter[iC][0][0] -> Fill(qxcCor2); + fQAfter[iC][0][1] -> Fill(qycCor2); + + Double_t evPlAngV0ACor2 = TVector2::Phi_0_2pi(TMath::ATan2(qyaCor2, qxaCor2))/2.; + Double_t evPlAngV0CCor2 = TVector2::Phi_0_2pi(TMath::ATan2(qycCor2, qxcCor2))/2.; + Double_t evPlAngV0Cor2 = TVector2::Phi_0_2pi(TMath::ATan2(qyCor2, qxCor2))/2.; + + fEventPlaneV0A = evPlAngV0ACor2; + fEventPlaneV0C = evPlAngV0CCor2; + fEventPlaneV0 = evPlAngV0Cor2; + + //printf("Eventplane V0A %f, V0C %f\n",fEventPlaneV0A,fEventPlaneV0C); + + } +} +//_____________________________________________________________________________ +Bool_t AliHFEVZEROEventPlane::OpenInfoCalbration(Int_t run) +{ + // + // Take the calibration coefficients + // + + //printf("Name of the file %s\n",(const char*)fnamefile); + + TFile *foadb = TFile::Open(fnamefile.Data()); + if(!foadb){ + printf("OADB file %s cannot be opened\n",fnamefile.Data()); + return kFALSE; + } + + //printf("test\n"); + + AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr"); + if(!cont){ + printf("OADB object hMultV0BefCorr is not available in the file\n"); + return kFALSE; + } + + if(!(cont->GetObject(run))){ + printf("OADB object hMultV0BefCorr is not available for run %i\n",run); + return kFALSE; + } + //TProfile *multV0 = ((TH2F *) cont->GetObject(run))->ProfileX(); + TProfile *multV0 = ((TProfile *) cont->GetObject(run)); + if(!multV0) return kFALSE; + fMultV0 = (TProfile *) multV0->Clone(); + fMultV0->SetDirectory(0); + + TF1 *fpol0 = new TF1("fpol0","pol0"); + fMultV0->Fit(fpol0,"Q0","",0,31); + fV0Cpol = fpol0->GetParameter(0); + fMultV0->Fit(fpol0,"Q0","",32,64); + fV0Apol = fpol0->GetParameter(0); + + for(Int_t iside=0;iside<2;iside++){ + for(Int_t icoord=0;icoord<2;icoord++){ + for(Int_t i=0;i < fgknCentrBin;i++){ + char namecont[100]; + if(iside==0 && icoord==0) + sprintf(namecont,"hQxc2_%i",i); + else if(iside==1 && icoord==0) + sprintf(namecont,"hQxa2_%i",i); + else if(iside==0 && icoord==1) + sprintf(namecont,"hQyc2_%i",i); + else if(iside==1 && icoord==1) + sprintf(namecont,"hQya2_%i",i); + + cont = (AliOADBContainer*) foadb->Get(namecont); + if(!cont){ + printf("OADB object %s is not available in the file\n",namecont); + delete fpol0; + return kFALSE; + } + + if(!(cont->GetObject(run))){ + printf("OADB object %s is not available for run %i\n",namecont,run); + delete fpol0; + return kFALSE; + } + fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean(); + fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS(); + } + } + } + + delete fpol0; + foadb->Close(); + + return kTRUE; + +} diff --git a/PWGHF/hfe/AliHFEVZEROEventPlane.h b/PWGHF/hfe/AliHFEVZEROEventPlane.h new file mode 100644 index 00000000000..e0f1ac3e17c --- /dev/null +++ b/PWGHF/hfe/AliHFEVZEROEventPlane.h @@ -0,0 +1,59 @@ +// +// VZERO event plane task for 2010 +// +#ifndef AliHFEVZEROEVENTPLANE_H +#define AliHFEVZEROEVENTPLANE_H + +#include "TProfile.h" +#include +#include "TString.h" +#include "TH1F.h" +#include "TList.h" + + +class AliHFEVZEROEventPlane : public TNamed { + public: + AliHFEVZEROEventPlane(); + AliHFEVZEROEventPlane(const char *name, const Char_t *title); + AliHFEVZEROEventPlane(const AliHFEVZEROEventPlane &ref); + AliHFEVZEROEventPlane& operator=(const AliHFEVZEROEventPlane &ref); + virtual void Copy(TObject &o) const; + ~AliHFEVZEROEventPlane(); + + void ProcessEvent(AliESDEvent *event); + + void SetNameFile(TString namefile) {fnamefile = namefile;}; + Bool_t OpenInfoCalbration(Int_t run); + + Double_t GetEventPlaneV0A() const {return fEventPlaneV0A;}; + Double_t GetEventPlaneV0C() const {return fEventPlaneV0C;}; + Double_t GetEventPlaneV0() const {return fEventPlaneV0;}; + + TList *GetOutputList() const {return fOutputList;}; + + private: + virtual void Analyze(AliESDEvent* esdEvent); + + Double_t fEventPlaneV0A; // Corrected event plane V0A + Double_t fEventPlaneV0C; // Corrected event plane V0C + Double_t fEventPlaneV0; // Corrected event plane V0 + + AliESDEvent* fESD; //! ESD object + Int_t fRun; // Run number + TProfile *fMultV0; //! fMultiplicityV0 + Float_t fV0Cpol,fV0Apol; // fV0Cpol, fV0Apol + static const Int_t fgknCentrBin = 9; // Centrality bins + Float_t fMeanQ[fgknCentrBin][2][2]; // mean for centering + Float_t fWidthQ[fgknCentrBin][2][2]; // rms for centering + TString fnamefile; // name of the file with the coefficient + TList *fOutputList; //! Output list + TProfile *fMultV0Before; //! fMultiplicityV0 Before + TProfile *fMultV0After; //! fMultiplicityV0 After + TH1F *fQBefore[fgknCentrBin][2][2]; //! Q centering before + TH1F *fQAfter[fgknCentrBin][2][2]; //! Q centering after + + + ClassDef(AliHFEVZEROEventPlane, 1); //Analysis task for high pt analysis +}; + +#endif diff --git a/PWGHF/hfe/AliHFEcollection.cxx b/PWGHF/hfe/AliHFEcollection.cxx index 1f2a4814de7..75f93b72295 100644 --- a/PWGHF/hfe/AliHFEcollection.cxx +++ b/PWGHF/hfe/AliHFEcollection.cxx @@ -149,6 +149,23 @@ Bool_t AliHFEcollection::CreateTH1Farray(const char* name, const char* title, In } } +//___________________________________________________________________ +Bool_t AliHFEcollection::CreateTH2Farray(const char* name, const char* title, Int_t nBin, const Double_t* xbins, Int_t nBinY, Float_t nMinY, Float_t nMaxY){ + + // + // Creates a TH1F histogram for the collection 2nd version + // + + if(!fList){ + AliError("No TList pointer ! "); + return kFALSE; + } + else{ + fList->Add(new TH2F(name, title, nBin, xbins, nBinY, nMinY, nMaxY)); + return CheckObject(name); + } +} + //___________________________________________________________________ Bool_t AliHFEcollection::CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t logAxis){ @@ -301,6 +318,21 @@ Bool_t AliHFEcollection::CreateTHnSparse(const char* name, const char* title, In fList->Add(new THnSparseF(name, title, dim, nbins, xmin, xmax)); return CheckObject(name); +} +//___________________________________________________________________ +Bool_t AliHFEcollection::CreateTHnSparseNoLimits(const char* name, const char* title, Int_t dim, const Int_t* nbins){ + + // + // create 'dim' dimensional THnSparse without limits + // + + if(!fList){ + AliError("No TList pointer ! "); + return kFALSE; + } + fList->Add(new THnSparseF(name, title, dim, nbins)); + return CheckObject(name); + } //___________________________________________________________________ TObject* AliHFEcollection::Get(const char* name){ diff --git a/PWGHF/hfe/AliHFEcollection.h b/PWGHF/hfe/AliHFEcollection.h index 5de56cef86e..1464f0f93d5 100644 --- a/PWGHF/hfe/AliHFEcollection.h +++ b/PWGHF/hfe/AliHFEcollection.h @@ -55,6 +55,7 @@ class AliHFEcollection : public TNamed{ Bool_t CreateTH1Farray(const char* name, const char* title, Int_t nBin, const Double_t* xbins); Bool_t CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t logAxis = -1); + Bool_t CreateTH2Farray(const char* name, const char* title, Int_t nBin, const Double_t* xbins, Int_t nBinY, Float_t nMinY, Float_t nMaxY); Bool_t CreateTH3F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t nBinZ, Float_t minZ, Float_t maxZ, Int_t logAxis = -1); Bool_t CreateTH1Fvector1(Int_t X, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax, Int_t logAxis = -1); @@ -62,6 +63,7 @@ class AliHFEcollection : public TNamed{ Bool_t CreateTH2Fvector1(Int_t X, const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY, Int_t logAxis = -1); Bool_t CreateProfile(const char* name, const char* title, Int_t nbins, Double_t xmin, Double_t xmax); Bool_t CreateTHnSparse(const char* name, const char* title, Int_t dim, const Int_t* nbins, const Double_t* xmin, const Double_t* xmax); + Bool_t CreateTHnSparseNoLimits(const char* name, const char* title, Int_t dim, const Int_t* nbins); Bool_t BinLogAxis(const char* name, Int_t dim); Bool_t Sumw2(const char*name); diff --git a/PWGHF/hfe/AliHFEdebugTreeTask.cxx b/PWGHF/hfe/AliHFEdebugTreeTask.cxx index 9edd39d1925..234cf902afc 100644 --- a/PWGHF/hfe/AliHFEdebugTreeTask.cxx +++ b/PWGHF/hfe/AliHFEdebugTreeTask.cxx @@ -18,7 +18,6 @@ // Authors: // Markus Fasel // -// #include #include @@ -37,6 +36,7 @@ #include "AliMCParticle.h" #include "AliPIDResponse.h" #include "AliVEvent.h" +#include "AliHFEpidTRD.h" #include "TTreeStream.h" #include "AliHFEdebugTreeTask.h" @@ -47,6 +47,7 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask(): AliAnalysisTaskSE(), fTrackCuts(NULL), fSignalCuts(NULL), + fTRDpid(NULL), fNclustersTPC(70), fNclustersTPCPID(0), fNclustersITS(2), @@ -60,6 +61,7 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name): AliAnalysisTaskSE(name), fTrackCuts(NULL), fSignalCuts(NULL), + fTRDpid(NULL), fNclustersTPC(70), fNclustersTPCPID(0), fNclustersITS(2), @@ -70,7 +72,8 @@ AliHFEdebugTreeTask::AliHFEdebugTreeTask(const char *name): } AliHFEdebugTreeTask::~AliHFEdebugTreeTask(){ - if(fDebugTree) delete fDebugTree; + if(fDebugTree) delete fDebugTree; + if(fTRDpid) delete fTRDpid; } void AliHFEdebugTreeTask::UserCreateOutputObjects(){ @@ -93,6 +96,9 @@ void AliHFEdebugTreeTask::UserCreateOutputObjects(){ fTrackCuts->SetUseMixedVertex(kTRUE); fTrackCuts->SetVertexRange(10.); fTrackCuts->Initialize(); + + fTRDpid = new AliHFEpidTRD("QAtrdPID"); + } void AliHFEdebugTreeTask::UserExec(Option_t *){ @@ -281,13 +287,16 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){ for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++; // TRD clusters and tracklets UChar_t nclustersTRD = track->GetTRDncls(); - UChar_t ntracklets = track->GetTRDntrackletsPID(); + UChar_t ntrackletsTRDPID = track->GetTRDntrackletsPID(); // ITS and TRD acceptance maps UChar_t hasClusterITS[6], hasTrackletTRD[6]; UChar_t itsPixel = track->GetITSClusterMap(); for(Int_t icl = 0; icl < 6; icl++) hasClusterITS[icl] = TESTBIT(itsPixel, icl) ? 1 : 0; + Double_t trddEdxSum[6]; + for(Int_t a=0;a<6;a++) { trddEdxSum[a]= 0.;} for(Int_t itl = 0; itl < 6; itl++){ - Int_t nSliceNonZero = 0; + Int_t nSliceNonZero = 0; + trddEdxSum[itl] = track->GetTRDslice(itl, 0); // in new reconstruction slice 0 contains the total charge for(Int_t islice = 0; islice < 8; islice++){ if(track->GetTRDslice(itl, islice) > 0.001) nSliceNonZero++; } @@ -298,6 +307,9 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){ track->GetTRDpid(pidprobs); Double_t likeEleTRD = pidprobs[0]; Double_t likeEleTRDn = likeEleTRD/(likeEleTRD + pidprobs[2]); + Double_t trdtruncmean1 = fTRDpid->GetTRDSignalV1(track, 0.6); + Double_t trdtruncmean2 = fTRDpid->GetTRDSignalV2(track, 0.6); + // DCA Float_t b[2] = {0.,0.}; Float_t bCov[3] = {0.,0.,0.}; @@ -326,7 +338,7 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){ << "pt=" << transversemomentum << "eta=" << eta << "phi=" << phi - << "ntracklets=" << ntracklets + << "ntracklets=" << ntrackletsTRDPID << "nclustersTPC=" << nclustersTPC << "nclustersTPCPID=" << nclustersTPCPID << "nclustersTPCshared=" << nclustersTPCshared @@ -348,11 +360,19 @@ void AliHFEdebugTreeTask::UserExec(Option_t *){ << "trd3=" << hasTrackletTRD[3] << "trd4=" << hasTrackletTRD[4] << "trd5=" << hasTrackletTRD[5] + << "TRDdEdxl0=" << trddEdxSum[0] + << "TRDdEdxl1=" << trddEdxSum[1] + << "TRDdEdxl2=" << trddEdxSum[2] + << "TRDdEdxl3=" << trddEdxSum[3] + << "TRDdEdxl4=" << trddEdxSum[4] + << "TRDdEdxl5=" << trddEdxSum[5] << "TOFsigmaEl=" << nSigmaTOF << "TPCsigmaEl=" << nSigmaTPC << "TPCdEdx=" << tPCdEdx << "TRDlikeEl=" << likeEleTRD << "TRDlikeEln=" << likeEleTRDn + << "trdtruncmean1=" << trdtruncmean1 + << "trdtruncmean2=" << trdtruncmean2 << "dcaR=" << b[0] << "dcaZ=" << b[1] << "dca=" << dca diff --git a/PWGHF/hfe/AliHFEdebugTreeTask.h b/PWGHF/hfe/AliHFEdebugTreeTask.h index 07f578c95d9..47a790797dd 100644 --- a/PWGHF/hfe/AliHFEdebugTreeTask.h +++ b/PWGHF/hfe/AliHFEdebugTreeTask.h @@ -25,6 +25,7 @@ class AliHFEcuts; class AliHFEsignalCuts; class TString; class TTreeSRedirector; +class AliHFEpidTRD; class AliHFEdebugTreeTask : public AliAnalysisTaskSE{ public: @@ -48,6 +49,7 @@ class AliHFEdebugTreeTask : public AliAnalysisTaskSE{ AliHFEcuts *fTrackCuts; // Track AliHFEsignalCuts *fSignalCuts; // Signal Cuts + AliHFEpidTRD *fTRDpid; // TRD PID Int_t fNclustersTPC; // Min Number of clusters in TPC Int_t fNclustersTPCPID; // Min Number of clusters for TPC PID Int_t fNclustersITS; // Min Number of clusters in ITS diff --git a/PWGHF/hfe/AliHFEemcalPIDqa.cxx b/PWGHF/hfe/AliHFEemcalPIDqa.cxx index 46ea73a14f2..bb34feec7c7 100644 --- a/PWGHF/hfe/AliHFEemcalPIDqa.cxx +++ b/PWGHF/hfe/AliHFEemcalPIDqa.cxx @@ -158,15 +158,15 @@ void AliHFEemcalPIDqa::Initialize(){ const Double_t kTPCSigMax = 140.; // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, Centrality, select) - Int_t nBins[12] = {AliPID::kSPECIES + 1, 500, 500, 400, 630, 200, 600, 300, 125, kCentralityBins, 6, 2}; + Int_t nBins[12] = {AliPID::kSPECIES + 1, 500, 500, 400, 630, 200, 600, 300, 100, kCentralityBins, 6, 2}; Double_t min[12] = {-1, kMinP, kMinP, kTPCSigMim, 0., -1.0, -8.0, 0, 0, 0, -3.0, 0.}; - Double_t max[12] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax, 6.3, 1.0, 4.0, 3.0, 0.5, 11., 3.0, 2.}; + Double_t max[12] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax, 6.3, 1.0, 4.0, 3.0, 0.1, 11., 3.0, 2.}; fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; phi ; eta ; nSig ; E/p ; Rmatch ;Centrality; PID Step; ", 12, nBins, min, max); //2nd histogram: EMCAL signal - E/p - Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 125, 600, 2}; + Int_t nBins2[7] = {AliPID::kSPECIES + 1, 500, 500, 500, 100, 600, 2}; Double_t min2[7] = {-1, kMinP, kMinP, 0, 0, -8., 0.}; - Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5, 0.5, 4., 2.}; + Double_t max2[7] = {AliPID::kSPECIES, kMaxP, kMaxP, 5, 0.1, 4., 2.}; fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; TPCnsigma; PID Step", 7, nBins2, min2, max2); } diff --git a/PWGHF/hfe/AliHFEextraCuts.h b/PWGHF/hfe/AliHFEextraCuts.h index 5ee7b15f181..00c14772383 100644 --- a/PWGHF/hfe/AliHFEextraCuts.h +++ b/PWGHF/hfe/AliHFEextraCuts.h @@ -82,19 +82,19 @@ class AliHFEextraCuts: public AliCFCutBase{ Bool_t GetCheckITSstatus() const { return fCheck; }; Int_t GetDebugLevel() const { return fDebugLevel; }; void GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy); // temporary moved from protected to publich for IP QA + Int_t GetITSstatus(AliVTrack *track, Int_t layer); + Bool_t CheckITSstatus(Int_t itsStatus) const; protected: virtual void AddQAHistograms(TList *qaList); Bool_t CheckRecCuts(AliVTrack *track); Bool_t CheckMCCuts(AliVParticle * /*track*/) const; - Bool_t CheckITSstatus(Int_t itsStatus) const; void FillQAhistosRec(AliVTrack *track, UInt_t when); // void FillQAhistosMC(AliMCParticle *track, UInt_t when); void FillCutCorrelation(ULong64_t survivedCut); void PrintBitMap(Int_t bitmap); // Getter Functions for ESD/AOD compatible mode - Int_t GetITSstatus(AliVTrack *track, Int_t layer); UInt_t GetTPCncls(AliVTrack *track); UInt_t GetTPCnclusdEdx(AliVTrack *track); Bool_t GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track); diff --git a/PWGHF/hfe/AliHFEmcQA.cxx b/PWGHF/hfe/AliHFEmcQA.cxx index 418eccc121e..6ddce290bf0 100644 --- a/PWGHF/hfe/AliHFEmcQA.cxx +++ b/PWGHF/hfe/AliHFEmcQA.cxx @@ -32,6 +32,7 @@ #include #include #include +#include "TTreeStream.h" #include #include @@ -56,6 +57,16 @@ ClassImp(AliHFEmcQA) ,fCentrality(0) ,fIsPbPb(kFALSE) ,fIsppMultiBin(kFALSE) + ,fContainerStep(0) + ,fIsDebugStreamerON(kFALSE) + ,fRecPt(-999) + ,fRecEta(-999) + ,fRecPhi(-999) + ,fLyrhit(0) + ,fLyrstat(0) + ,fHfeImpactR(-999) + ,fHfeImpactnsigmaR(-999) + ,fTreeStream(NULL) { // Default constructor for(Int_t mom = 0; mom < 9; mom++){ @@ -83,6 +94,16 @@ TObject(p) ,fCentrality(0) ,fIsPbPb(kFALSE) ,fIsppMultiBin(kFALSE) + ,fContainerStep(0) + ,fIsDebugStreamerON(kFALSE) + ,fRecPt(-999) + ,fRecEta(-999) + ,fRecPhi(-999) + ,fLyrhit(0) + ,fLyrstat(0) + ,fHfeImpactR(0) + ,fHfeImpactnsigmaR(0) + ,fTreeStream(NULL) { // Copy constructor for(Int_t mom = 0; mom < 9; mom++){ @@ -113,6 +134,7 @@ AliHFEmcQA::~AliHFEmcQA() { // Destructor + if(fTreeStream && fIsDebugStreamerON) delete fTreeStream; AliInfo("Analysis Done."); } @@ -147,15 +169,9 @@ void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList) fQAhistos = qaList; fQAhistos->SetName("MCqa"); - CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm - CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty - CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_"); // create histograms for beauty - CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm - CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty - CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_"); // create histograms for beauty - CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm - CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty - CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_"); // create histograms for beauty + CreateHistograms(AliHFEmcQA::kCharm); // create histograms for charm + CreateHistograms(AliHFEmcQA::kBeauty); // create histograms for beauty + CreateHistograms(AliHFEmcQA::kOthers); // create histograms for beauty // prepare 2D(pt vs Y) histogram for D spectra, we consider following 9 particles const Int_t nbspecies = 9; @@ -229,10 +245,14 @@ void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList) fQAhistos->Add(fMCQACollection->GetList()); + if(!fTreeStream && fIsDebugStreamerON){ + fTreeStream = new TTreeSRedirector(Form("HFEmcqadebugTree%s.root", GetName())); + } + } //__________________________________________ -void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) +void AliHFEmcQA::CreateHistograms(const Int_t kquark) { // create histograms @@ -266,6 +286,11 @@ void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) kqTypeLabel[kMisID-4]="miside"; } + TString kqEtaRangeLabel[fgkEtaRanges]; + kqEtaRangeLabel[0] = "mcqa_"; + kqEtaRangeLabel[1] = "mcqa_barrel_"; + kqEtaRangeLabel[2] = "mcqa_unitY_"; + const Double_t kPtbound[2] = {0.1, 20.}; //bin taken for considering inclusive e analysis binning const Int_t nptbinning1 = 35; Int_t iBin[2]; @@ -285,64 +310,62 @@ void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) TString hname; if(kquark == kOthers){ + for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){ for (Int_t iqType = 0; iqType < 4; iqType++ ){ - hname = hnopt+"Pt_"+kqTypeLabel[iqType]; - //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL + hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType]; fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); - hname = hnopt+"Y_"+kqTypeLabel[iqType]; + hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType]; fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5); - hname = hnopt+"Eta_"+kqTypeLabel[iqType]; + hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType]; fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5); // Fill List if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos); } - return; + } + return; } - for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){ + for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){ + for (Int_t iqType = 0; iqType < fgkqType; iqType++ ){ if (iqType < keHadron && icut > 0) continue; // don't duplicate histogram for quark and hadron - hname = hnopt+"PdgCode_"+kqTypeLabel[iqType]; + hname = kqEtaRangeLabel[icut]+"PdgCode_"+kqTypeLabel[iqType]; fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5); - hname = hnopt+"Pt_"+kqTypeLabel[iqType]; + hname = kqEtaRangeLabel[icut]+"Pt_"+kqTypeLabel[iqType]; fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[1],kPtbinning1); // new binning - //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); // old log binning - //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.,30.); // mj to compare with FONLL - //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL - hname = hnopt+"Y_"+kqTypeLabel[iqType]; + hname = kqEtaRangeLabel[icut]+"Y_"+kqTypeLabel[iqType]; fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5); - hname = hnopt+"Eta_"+kqTypeLabel[iqType]; + hname = kqEtaRangeLabel[icut]+"Eta_"+kqTypeLabel[iqType]; fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5); // Fill List if(fQAhistos) fHist[iq][iqType][icut].FillList(fQAhistos); + } } - if (icut == 0){ - hname = hnopt+"Nq_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5); - hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5); + for (Int_t icut = 0; icut < fgkEtaRanges; icut++ ){ + hname = kqEtaRangeLabel[icut]+"PtCorr_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); + hname = kqEtaRangeLabel[icut]+"PtCorrDp_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); + hname = kqEtaRangeLabel[icut]+"PtCorrD0_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); + hname = kqEtaRangeLabel[icut]+"PtCorrDrest_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); + + hname = kqEtaRangeLabel[icut]+"ePtRatio_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1); + hname = kqEtaRangeLabel[icut]+"DePtRatio_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1); + hname = kqEtaRangeLabel[icut]+"eDistance_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2); + hname = kqEtaRangeLabel[icut]+"DeDistance_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2); + if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos); } - hname = hnopt+"PtCorr_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning - //fHistComm[iq][icut].fPtCorr = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]); - hname = hnopt+"PtCorrDp_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning - //fHistComm[iq][icut].fPtCorrDp= new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]); - hname = hnopt+"PtCorrD0_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning - //fHistComm[iq][icut].fPtCorrD0 = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]); - hname = hnopt+"PtCorrDrest_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[1],kPtbinning1); // new binning - //fHistComm[iq][icut].fPtCorrDrest = new TH2F(hname,hname+";D p_{T} (GeV/c);e p_{T} (GeV/c)",ndptbins,xcorrbin,iBin[0],binEdges[0]); + + hname = kqEtaRangeLabel[0]+"Nq_"+kqTypeLabel[kQuark]; + fHistComm[iq][0].fNq = new TH1F(hname,hname,50,-0.5,49.5); + hname = kqEtaRangeLabel[0]+"ProcessID_"+kqTypeLabel[kQuark]; + fHistComm[iq][0].fProcessID = new TH1F(hname,hname,21,-10.5,10.5); - hname = hnopt+"ePtRatio_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",200,0,20,100,0,1); - hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,20,100,0,1); - hname = hnopt+"eDistance_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2); - hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark]; - fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,20,200,0,2); - if(fQAhistos) fHistComm[iq][icut].FillList(fQAhistos); } //__________________________________________ @@ -791,7 +814,7 @@ void AliHFEmcQA::GetHadronKine(TParticle* mcpart, const Int_t kquark) } //__________________________________________ -void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut) +void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kdecayed) { // decay electron kinematics @@ -807,14 +830,27 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde return; } + Double_t eabsEta = TMath::Abs(mcpart->Eta()); + Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart)); + if(kquark==kOthers){ Int_t esource = -1; if ( abs(mcpart->GetPdgCode()) != kdecayed ) esource = kMisID-4; else esource =GetSource(mcpart)-4; // return for the cases kGamma=4, kPi0=5, kElse=6 if(esource==0|| esource==1 || esource==2 || esource==3){ - fHist[iq][esource][icut].fPt->Fill(mcpart->Pt()); - fHist[iq][esource][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart)); - fHist[iq][esource][icut].fEta->Fill(mcpart->Eta()); + fHist[iq][esource][0].fPt->Fill(mcpart->Pt()); + fHist[iq][esource][0].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][esource][0].fEta->Fill(mcpart->Eta()); + if(eabsEta<0.9){ + fHist[iq][esource][1].fPt->Fill(mcpart->Pt()); + fHist[iq][esource][1].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][esource][1].fEta->Fill(mcpart->Eta()); + } + if(eabsY<0.5){ + fHist[iq][esource][2].fPt->Fill(mcpart->Pt()); + fHist[iq][esource][2].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][esource][2].fEta->Fill(mcpart->Eta()); + } return; } else { @@ -886,22 +922,62 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde if (kquark == kCharm) return; // fill electron kinematics - fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode()); - fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt()); - fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart)); - fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta()); + fHist[iq][kElectron2nd][0].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][0].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][0].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron2nd][0].fEta->Fill(mcpart->Eta()); // fill mother hadron kinematics - fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); - fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt()); - fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa)); - fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta()); - + fHist[iq][kDeHadron][0].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][0].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][0].fY->Fill(AliHFEtools::GetRapidity(grandMa)); + fHist[iq][kDeHadron][0].fEta->Fill(grandMa->Eta()); + + if(eabsEta<0.9){ + fHist[iq][kElectron2nd][1].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][1].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][1].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron2nd][1].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][kDeHadron][1].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][1].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][1].fY->Fill(AliHFEtools::GetRapidity(grandMa)); + fHist[iq][kDeHadron][1].fEta->Fill(grandMa->Eta()); + } + + if(eabsY<0.5){ + fHist[iq][kElectron2nd][2].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][2].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][2].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron2nd][2].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][kDeHadron][2].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][2].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][2].fY->Fill(AliHFEtools::GetRapidity(grandMa)); + fHist[iq][kDeHadron][2].fEta->Fill(grandMa->Eta()); + } + // ratio between pT of electron and pT of mother B hadron - if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt()); + if(grandMa->Pt()) { + fHistComm[iq][0].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt()); + if(eabsEta<0.9){ + fHistComm[iq][1].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt()); + } + if(eabsY<0.5){ + fHistComm[iq][2].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt()); + } + } // distance between electron production point and primary vertex - fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy); + fHistComm[iq][0].fDeDistance->Fill(grandMa->Pt(),decayLxy); + if(eabsEta<0.9){ + fHistComm[iq][1].fDeDistance->Fill(grandMa->Pt(),decayLxy); + } + if(eabsY<0.5){ + fHistComm[iq][2].fDeDistance->Fill(grandMa->Pt(),decayLxy); + } return; } } @@ -913,51 +989,103 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde for (Int_t i=0; iPt(); - Double_t wfactor=2*(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225)); // 2* considering in pythia I used particle + antiparticle differently from the measurement - //Double_t wfactor=(703.681*mpt/TMath::Power((1+TMath::Power(mpt/1.73926,2)),2.34821))/(368.608*mpt/TMath::Power((1+TMath::Power(mpt/2.74868,2)),2.34225)); - // fill electron kinematics - if(iq==0){ - fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode(),wfactor); - fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt(),wfactor); - fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart),wfactor); - fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta(),wfactor); - - fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode()); - fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt()); - fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart)); - fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta()); - } - else{ - fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode()); - fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt()); - fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart)); - fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta()); - } -//-------------- + fHist[iq][kElectron][0].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron][0].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][0].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron][0].fEta->Fill(mcpart->Eta()); + // fill mother hadron kinematics - fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); - fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt()); - fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); - fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta()); + fHist[iq][keHadron][0].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][0].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][0].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); + fHist[iq][keHadron][0].fEta->Fill(partMotherCopy->Eta()); + + if(eabsEta<0.9){ + fHist[iq][kElectron][1].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron][1].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][1].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron][1].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][keHadron][1].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][1].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][1].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); + fHist[iq][keHadron][1].fEta->Fill(partMotherCopy->Eta()); + } + + if(eabsY<0.5){ + fHist[iq][kElectron][2].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron][2].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][2].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron][2].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][keHadron][2].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][2].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][2].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); + fHist[iq][keHadron][2].fEta->Fill(partMotherCopy->Eta()); + } // ratio between pT of electron and pT of mother B or direct D hadron - if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt()); - fHistComm[iq][icut].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt()); - if(TMath::Abs(partMotherCopy->GetPdgCode())==411) fHistComm[iq][icut].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt()); - else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) fHistComm[iq][icut].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt()); - else fHistComm[iq][icut].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt()); + if(partMotherCopy->Pt()) { + fHistComm[iq][0].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt()); + if(eabsEta<0.9){ + fHistComm[iq][1].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt()); + } + if(eabsY<0.5){ + fHistComm[iq][2].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt()); + } + } + fHistComm[iq][0].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt()); + if(eabsEta<0.9){ + fHistComm[iq][1].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + if(eabsY<0.5){ + fHistComm[iq][2].fPtCorr->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + if(TMath::Abs(partMotherCopy->GetPdgCode())==411) { + fHistComm[iq][0].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt()); + if(eabsEta<0.9){ + fHistComm[iq][1].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + if(eabsY<0.5){ + fHistComm[iq][2].fPtCorrDp->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + } + else if(TMath::Abs(partMotherCopy->GetPdgCode())==421) { + fHistComm[iq][0].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt()); + if(eabsEta<0.9){ + fHistComm[iq][1].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + if(eabsY<0.5){ + fHistComm[iq][2].fPtCorrD0->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + } + else { + fHistComm[iq][0].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt()); + if(eabsEta<0.9){ + fHistComm[iq][1].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + if(eabsY<0.5){ + fHistComm[iq][2].fPtCorrDrest->Fill(partMotherCopy->Pt(),mcpart->Pt()); + } + } // distance between electron production point and primary vertex - fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy); + fHistComm[iq][0].feDistance->Fill(partMotherCopy->Pt(),decayLxy); + if(eabsEta<0.9){ + fHistComm[iq][1].feDistance->Fill(partMotherCopy->Pt(),decayLxy); + } + if(eabsY<0.5){ + fHistComm[iq][2].feDistance->Fill(partMotherCopy->Pt(),decayLxy); + } } } } // end of if } //____________________________________________________________________ -void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut) +void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed) { // decay electron kinematics @@ -976,6 +1104,9 @@ void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I if ( abs(mcpart->GetPdgCode()) != kdecayed ) return; + Double_t eabsEta = TMath::Abs(mcpart->Eta()); + Double_t eabsY = TMath::Abs(AliHFEtools::GetRapidity(mcpart)); + // mother Int_t iLabel = mcpart->GetMother(); if (iLabel<0){ @@ -1019,16 +1150,43 @@ void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I if (kquark == kCharm) return; // fill electron kinematics - fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode()); - fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt()); - fHist[iq][kElectron2nd][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart)); - fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta()); + fHist[iq][kElectron2nd][0].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][0].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][0].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron2nd][0].fEta->Fill(mcpart->Eta()); // fill mother hadron kinematics - fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); - fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt()); - fHist[iq][kDeHadron][icut].fY->Fill(AliHFEtools::GetRapidity(grandMa)); - fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta()); + fHist[iq][kDeHadron][0].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][0].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][0].fY->Fill(AliHFEtools::GetRapidity(grandMa)); + fHist[iq][kDeHadron][0].fEta->Fill(grandMa->Eta()); + + if(eabsEta<0.9){ + // fill electron kinematics + fHist[iq][kElectron2nd][1].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][1].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][1].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron2nd][1].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][kDeHadron][1].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][1].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][1].fY->Fill(AliHFEtools::GetRapidity(grandMa)); + fHist[iq][kDeHadron][1].fEta->Fill(grandMa->Eta()); + } + if(eabsY<0.5){ + // fill electron kinematics + fHist[iq][kElectron2nd][2].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][2].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][2].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron2nd][2].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][kDeHadron][2].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][2].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][2].fY->Fill(AliHFEtools::GetRapidity(grandMa)); + fHist[iq][kDeHadron][2].fEta->Fill(grandMa->Eta()); + } return; } @@ -1042,16 +1200,43 @@ void AliHFEmcQA::GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, I if (abs(maPdgcodeCopy)==fParentSelect[iq][i]){ // fill electron kinematics - fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode()); - fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt()); - fHist[iq][kElectron][icut].fY->Fill(AliHFEtools::GetRapidity(mcpart)); - fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta()); + fHist[iq][kElectron][0].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron][0].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][0].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron][0].fEta->Fill(mcpart->Eta()); // fill mother hadron kinematics - fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); - fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt()); - fHist[iq][keHadron][icut].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); - fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta()); + fHist[iq][keHadron][0].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][0].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][0].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); + fHist[iq][keHadron][0].fEta->Fill(partMotherCopy->Eta()); + + if(eabsEta<0.9){ + // fill electron kinematics + fHist[iq][kElectron][1].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron][1].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][1].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron][1].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][keHadron][1].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][1].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][1].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); + fHist[iq][keHadron][1].fEta->Fill(partMotherCopy->Eta()); + } + if(eabsY<0.5){ + // fill electron kinematics + fHist[iq][kElectron][2].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron][2].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][2].fY->Fill(AliHFEtools::GetRapidity(mcpart)); + fHist[iq][kElectron][2].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][keHadron][2].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][2].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][2].fY->Fill(AliHFEtools::GetRapidity(partMotherCopy)); + fHist[iq][keHadron][2].fEta->Fill(partMotherCopy->Eta()); + } } } @@ -1486,6 +1671,21 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve else if(mesonID==kGammaEtaPrime || mesonID==kEtaPrime) mArr=4; //etaprime else if(mesonID==kGammaRho0 || mesonID==kRho0) mArr=5; //rho + Double_t datamc[24]={-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999,-999, -999, -999, -999, -999, -999, -999, -999, -999}; + Double_t xr[3]={-999,-999,-999}; + datamc[0] = mesonID; + datamc[17] = mctrack->Pt(); //electron pt + datamc[18] = mctrack->Eta(); //electron eta + + mctrack->XvYvZv(xr); + datamc[9] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); + datamc[10] = xr[2]; + + TParticle *mcpart = mctrack->Particle(); + if(mcpart){ + datamc[14] = mcpart->GetUniqueID(); + } + if(!(mArr<0)){ if(mesonID>=kGammaPi0) { // conversion electron, be careful with the enum odering Int_t glabel=TMath::Abs(mctrack->GetMother()); // gamma label @@ -1494,6 +1694,37 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ mesonPt = mctrackmother->Pt(); //meson pt bgcategory = 1.; + datamc[1] = bgcategory; + datamc[2] = mesonPt; + mctrackmother->XvYvZv(xr); + datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); + datamc[12] = xr[2]; + + mcpart = mctrackmother->Particle(); + if(mcpart){ + datamc[15] = mcpart->GetUniqueID(); + } + if(glabel>fMCEvent->GetNumberOfPrimaries()) { + bgcategory = 2.; + datamc[1] = bgcategory; + //printf("I should be gamma meson = %d mesonlabel= %d NumberOfPrimaries= %d \n",mctrackmother->PdgCode(),glabel,fMCEvent->GetNumberOfPrimaries()); + glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother + if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ + datamc[3]=mctrackmother->PdgCode(); + datamc[4]=mctrackmother->Pt(); + if(TMath::Abs(mctrackmother->PdgCode())==310){ + bgcategory = 3.; + glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother + if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ + datamc[5]=mctrackmother->PdgCode(); + glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother + if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ + datamc[6]=mctrackmother->PdgCode(); + } + } + } + } + } } } } @@ -1502,9 +1733,40 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ mesonPt=mctrackmother->Pt(); //meson pt bgcategory = -1.; + datamc[1] = bgcategory; + datamc[2] = mesonPt; + mctrackmother->XvYvZv(xr); + datamc[11] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); + datamc[12] = xr[2]; + + mcpart = mctrackmother->Particle(); + if(mcpart){ + datamc[15] = mcpart->GetUniqueID(); + } + if(glabel>fMCEvent->GetNumberOfPrimaries()) { + bgcategory = -2.; + datamc[1] = bgcategory; + glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother + if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ + datamc[3]=mctrackmother->PdgCode(); + datamc[4]=mctrackmother->Pt(); + if(TMath::Abs(mctrackmother->PdgCode())==310){ + bgcategory = -3.; + glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother's mother + if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ + datamc[5]=mctrackmother->PdgCode(); + glabel=TMath::Abs(mctrackmother->GetMother()); // gamma's mother's mother + if((mctrackmother = dynamic_cast(fMCEvent->GetTrack(glabel)))){ + datamc[6]=mctrackmother->PdgCode(); + } + } + } + } + } } } + if(fIsPbPb){ if(fCentrality < 0)return 0.; weightElecBg=fElecBackgroundFactor[iBgLevel][fCentrality][mArr][kBgPtBins-1]; @@ -1525,7 +1787,54 @@ Double_t AliHFEmcQA::GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLeve } } } - return bgcategory*weightElecBg; + + datamc[13] = weightElecBg; + datamc[16] = Double_t(fContainerStep); + + datamc[7] = fHfeImpactR; + datamc[8] = fHfeImpactnsigmaR; + + datamc[19] = fRecPt; + datamc[20] = fRecEta; + datamc[21] = fRecPhi; + datamc[22] = fLyrhit; + datamc[23] = fLyrstat; + + if(fIsDebugStreamerON && fTreeStream){ + if(!iBgLevel){ + (*fTreeStream)<<"nonhfeQA"<< + "mesonID="<1) returnval = bgcategory/2.; + + return returnval; } //__________________________________________ diff --git a/PWGHF/hfe/AliHFEmcQA.h b/PWGHF/hfe/AliHFEmcQA.h index b8cf76d2997..55d2c7aceb0 100644 --- a/PWGHF/hfe/AliHFEmcQA.h +++ b/PWGHF/hfe/AliHFEmcQA.h @@ -40,6 +40,7 @@ class AliGenEventHeader; class AliMCParticle; class AliAODMCParticle; class AliHFEcollection; +class TTreeSRedirector; //________________________________________________________________ class AliHFEmcQA: public TObject { @@ -68,7 +69,7 @@ class AliHFEmcQA: public TObject { TList *GetList() const { return fQAhistos; }; void PostAnalyze() const; void CreatDefaultHistograms(TList * const qaList); // create default histograms - void CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt=""); // create histograms for mc qa analysis + void CreateHistograms(const Int_t kquark); // create histograms for mc qa analysis void SetMCEvent(AliMCEvent* const mcEvent){fMCEvent = mcEvent;} void SetGenEventHeader(AliGenEventHeader* const mcHeader){fMCHeader=mcHeader;} // set stack pointer void SetMCArray(TClonesArray* const mcarry){fMCArray=mcarry;} // set mcarray pointer @@ -76,16 +77,21 @@ class AliHFEmcQA: public TObject { void GetQuarkKine(TParticle *part, Int_t iTrack, const Int_t kquark); // get heavy quark kinematics distribution void GetHadronKine(TParticle *part, const Int_t kquark); // get heavy hadron kinematics distribution - void GetDecayedKine(TParticle *part, const Int_t kquark, const Int_t kdecayed, Int_t icut); // get decay electron kinematics distribution - void GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut); // get decay electron kinematics for AOD + void GetDecayedKine(TParticle *part, const Int_t kquark, const Int_t kdecayed); // get decay electron kinematics distribution + void GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed); // get decay electron kinematics for AOD void GetMesonKine(); // get meson and its decay electron pt spectra void EndOfEventAna(const Int_t kquark); // run analysis which should be done at the end of the event loop Int_t GetSource(const TParticle * const mcpart); // return source id Int_t GetElecSource(TParticle * const mcpart); // return electron source id Int_t GetSource(const AliAODMCParticle * const mcpart); // return electron source id for AOD Double_t GetWeightFactor(AliMCParticle *mctrack, const Int_t iBgLevel); // return best/lower/upper weighting factor for electron's mother meson + void EnableDebugStreamer() { fIsDebugStreamerON = kTRUE;}; void SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit); + void SetContainerStep(Int_t containerStep) { fContainerStep = containerStep;}; + void SetHFEImpactParameters(Double_t hfeimpactR, Double_t hfeimpactnsigmaR) {fHfeImpactR = hfeimpactR; fHfeImpactnsigmaR = hfeimpactnsigmaR; }; + void SetTrkKine(Double_t pt, Double_t eta, Double_t phi) {fRecPt = pt; fRecEta = eta; fRecPhi = phi;}; + void SetITSInfo(Double_t ilyrhit, Double_t ilyrstat) { fLyrhit = ilyrhit; fLyrstat = ilyrstat;}; void SetCentrality(Int_t centrality) { fCentrality = centrality; }; void SetPbPb() { fIsPbPb = kTRUE; }; @@ -110,6 +116,7 @@ class AliHFEmcQA: public TObject { static const Int_t fgkMaxGener=10; // ancester level wanted to be checked static const Int_t fgkMaxIter=100; // number of iteration to find out matching particle static const Int_t fgkqType=7; // number of particle type to be checked + static const Int_t fgkEtaRanges=3; // cuts for different eta ranges struct AliHists{ TH1F *fPdgCode; // histogram to store particle pdg code @@ -206,6 +213,19 @@ private: Int_t fCentrality; // Centrality Bool_t fIsPbPb; // Analysis Type: pp or PbPb Bool_t fIsppMultiBin; // pp multiplicity bin analysis + Int_t fContainerStep; // step the weighting factor called + Bool_t fIsDebugStreamerON; // check if the debugstreamer is on + + Double_t fRecPt; //reconstructed pt + Double_t fRecEta; //reconstructed eta + Double_t fRecPhi; //reconstructed phi + Double_t fLyrhit; //its layer hit + Double_t fLyrstat; // its layer status + + Double_t fHfeImpactR; // absolute impact parameter R + Double_t fHfeImpactnsigmaR; // absolute impact parameter sigma R + + TTreeSRedirector *fTreeStream; //! TreeStream ClassDef(AliHFEmcQA,1); }; diff --git a/PWGHF/hfe/AliHFEpidEMCAL.cxx b/PWGHF/hfe/AliHFEpidEMCAL.cxx index 8453be93312..027c2bd1f5c 100644 --- a/PWGHF/hfe/AliHFEpidEMCAL.cxx +++ b/PWGHF/hfe/AliHFEpidEMCAL.cxx @@ -159,8 +159,8 @@ Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanage feopMaxCut = CalEopCutMax(track->GetRecTrack(),1); } - //printf("eop cuts org; %g; %g \n",feopMim,feopMax); - //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut); + //printf("eop cuts org; %g; %g \n",feopMim,feopMax); + //printf("eop cuts ; %g; %g \n",feopMimCut,feopMaxCut); Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?) AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop)); @@ -198,20 +198,23 @@ Double_t AliHFEpidEMCAL::CalEopCutMax(const AliVParticle *const track, Int_t fla if(flageop<0.5) { //<--- new para for updated non-liniarity - double meanP[3] = {1.02902,1.22293,1.67287e-01}; - double sigP[3] = {5.36355e-02,-3.25999e-02,4.42442e-01}; + double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02}; + double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04}; double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); maxCut = mean+3.0*sig; } else { - double meanP[3] = {0.99,1.299,0.35}; - double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; + double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01}; + double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01}; double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); maxCut = mean+3.0*sig; } + + //printf("eop cuts org; %g; %g ; %g\n",feopMim,feopMax,pt); + return maxCut; } @@ -226,22 +229,20 @@ Double_t AliHFEpidEMCAL::CalEopCutMim(const AliVParticle *const track, Int_t fla if(flageop<0.5) // real { //printf("real"); - //double meanP[3] = {1.056,1.169,0.1339}; - //double sigP[3] = {3.37e-05,1.298e-04,1.403e-04}; - double meanP[3] = {1.02902,1.22293,1.67287e-01}; - double sigP[3] = {5.36355e-02,-3.25999e-02,4.42442e-01}; + double meanP[3] = {1.08833e+00,1.17837e+00,7.47923e-02}; + double sigP[3] = {9.93208e-05,7.13074e-04,2.45145e-04}; double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); - mimCut = mean+0.0*sig; + mimCut = mean-3.0*sig; } else // MC { //printf("MC"); - double meanP[3] = {0.99,1.299,0.35}; - double sigP[3] = {4.11e-02,1.588e-01,2.664e-01}; + double meanP[3] = {9.94183e-01,1.23952e+00,3.55571e-01}; + double sigP[3] = {4.61207e-02,-3.39978e-02,4.26198e-01}; double mean = meanP[0]*tanh(meanP[1]+meanP[2]*pt); double sig = sigP[0]/tanh(sigP[1]+sigP[2]*pt); - mimCut = mean+0.0*sig; + mimCut = mean-3.0*sig; } return mimCut; } diff --git a/PWGHF/hfe/AliHFEpidQA.cxx b/PWGHF/hfe/AliHFEpidQA.cxx index 9a96f698d2d..f1f823440be 100644 --- a/PWGHF/hfe/AliHFEpidQA.cxx +++ b/PWGHF/hfe/AliHFEpidQA.cxx @@ -1309,7 +1309,7 @@ Int_t AliHFEpidQA::GetCentrality(AliVEvent* const fInputEvent){ } //___________________________________________________ -Int_t AliHFEpidQA::GetMultiplicityITS(AliVEvent* const fInputEvent) +Int_t AliHFEpidQA::GetMultiplicityITS(AliVEvent* const fInputEvent) const { // // Definition of the Multiplicity according to the JPSI group (F. Kramer) diff --git a/PWGHF/hfe/AliHFEpidQA.h b/PWGHF/hfe/AliHFEpidQA.h index 721b5d9951e..fcd242fe2dd 100644 --- a/PWGHF/hfe/AliHFEpidQA.h +++ b/PWGHF/hfe/AliHFEpidQA.h @@ -80,7 +80,7 @@ class AliHFEpidQA : public TObject{ AliHFEtrdPIDqa *GetTRDQA() const { return fTRDpidQA; } Int_t GetCentrality(AliVEvent* const fInputEvent); - Int_t GetMultiplicityITS(AliVEvent* const fInputEvent); + Int_t GetMultiplicityITS(AliVEvent* const fInputEvent) const; void SetPbPb() { fIsPbPb = kTRUE; }; void SetPP() { fIsPbPb = kFALSE; }; void SetPPMultiBin() { fIsppMultiBin = kFALSE; }; diff --git a/PWGHF/hfe/AliHFEsignalCuts.cxx b/PWGHF/hfe/AliHFEsignalCuts.cxx index bddaf85037e..02990c71a1d 100644 --- a/PWGHF/hfe/AliHFEsignalCuts.cxx +++ b/PWGHF/hfe/AliHFEsignalCuts.cxx @@ -320,6 +320,7 @@ Int_t AliHFEsignalCuts::GetElecSource(const AliVParticle * const track) const { TClass *tracktype; const AliVParticle *mctrack = NULL; TParticle *mcpart = NULL; + //AliMCParticle *esdmcmother = NULL; if((tracktype = track->IsA()) == AliESDtrack::Class() || tracktype == AliAODTrack::Class()){ mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())); } else mctrack = track; @@ -331,6 +332,23 @@ Int_t AliHFEsignalCuts::GetElecSource(const AliVParticle * const track) const { if(esdmc){ mcpart = esdmc->Particle(); eSource=fMCQA->GetElecSource(mcpart); +/* // considering secondary pions + if(eSource>=AliHFEmcQA::kGammaPi0) { // conversion electron, be careful with the enum odering + Int_t glabel=TMath::Abs(esdmc->GetMother()); // gamma label + if((esdmcmother= dynamic_cast(fMC->GetTrack(glabel)))){ + glabel=TMath::Abs(esdmcmother->GetMother()); // gamma's mother's label + if((esdmcmother= dynamic_cast(fMC->GetTrack(glabel)))){ + if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kElse; + } + } + } + else if(eSource==AliHFEmcQA::kPi0 || (eSource>=AliHFEmcQA::kEta && eSource<=AliHFEmcQA::kRho0) ){ // nonHFE except for the conversion electron + Int_t glabel=TMath::Abs(esdmc->GetMother()); + if((esdmcmother= dynamic_cast(fMC->GetTrack(glabel)))){ + if(glabel>fMC->GetNumberOfPrimaries()) eSource=AliHFEmcQA::kElse; + } + } +*/ } } else { return -1; diff --git a/PWGHF/hfe/AliHFEspectrum.cxx b/PWGHF/hfe/AliHFEspectrum.cxx index 2c59c9d446a..61af42b2a48 100644 --- a/PWGHF/hfe/AliHFEspectrum.cxx +++ b/PWGHF/hfe/AliHFEspectrum.cxx @@ -1,3 +1,4 @@ + /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -76,11 +77,11 @@ AliHFEspectrum::AliHFEspectrum(const char *name): fIPanaCharmBgSubtract(kFALSE), fIPanaConversionBgSubtract(kFALSE), fIPanaNonHFEBgSubtract(kFALSE), - fNonHFEbgMethod2(kFALSE), + fIPParameterizedEff(kFALSE), fNonHFEsyst(0), + fBeauty2ndMethod(kFALSE), fNbDimensions(1), fNMCEvents(0), - fNMCbgEvents(0), fStepMC(-1), fStepTrue(-1), fStepData(-1), @@ -91,8 +92,8 @@ AliHFEspectrum::AliHFEspectrum(const char *name): fChargeChoosen(kAllCharge), fNCentralityBinAtTheEnd(0), fHadronEffbyIPcut(NULL), - fConversionEff(NULL), - fNonHFEEff(NULL), + fBSpectrum2ndMethod(NULL), + fkBeauty2ndMethodfilename(""), fBeamType(0), fDebugLevel(0), fWriteToFile(kFALSE) @@ -102,9 +103,25 @@ AliHFEspectrum::AliHFEspectrum(const char *name): // for(Int_t k = 0; k < 20; k++){ - fNEvents[k] = 0; - fLowBoundaryCentralityBinAtTheEnd[k] = 0; - fHighBoundaryCentralityBinAtTheEnd[k] = 0; + fNEvents[k] = 0; + fNMCbgEvents[k] = 0; + fLowBoundaryCentralityBinAtTheEnd[k] = 0; + fHighBoundaryCentralityBinAtTheEnd[k] = 0; + if(kGetCFContainer("recTrackContReco"); } else{ - datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco"); + + datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco"); } AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground"); if((!datacontainer) || (!contaminationcontainer)) return kFALSE; @@ -219,10 +248,21 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){ for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel])); - fNonHFESourceContainer[iSource][iLevel] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen); - convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel])); - fConvSourceContainer[iSource][iLevel] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen); - if((!fConvSourceContainer[iSource][iLevel])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE; + convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel])); + for(Int_t icentr=0;icentrCalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step # - fConversionEff = (TH1D*)efficiencyConv->Project(0); - - source = 3; //non HFE except for the conversion - mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen); - AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD); - efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. be careful with step # - fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0); + SetParameterizedEff(mccontainermc, mccontainermcbg, dims); if(!fNonHFEsyst){ AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen); @@ -346,6 +376,17 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF } +//____________________________________________________________ +void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){ + // + // get spectrum for beauty 2nd method + // + // + TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename); + fBSpectrum2ndMethod = new TH1D(*static_cast(inputfile2ndmethod->Get("BSpectrum"))); +} + + //____________________________________________________________ Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){ // @@ -708,12 +749,15 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ AliCFDataGrid *unnormalizedRawSpectrum = 0x0; TGraphErrors *gNormalizedRawSpectrum = 0x0; if(subtractcontamination) { - dataspectrumaftersubstraction = SubtractBackground(kTRUE); - unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone(); - dataGridAfterFirstSteps = dataspectrumaftersubstraction; - gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum); + if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE); + else dataspectrumaftersubstraction = GetRawBspectra2ndMethod(); + unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone(); + dataGridAfterFirstSteps = dataspectrumaftersubstraction; + gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum); } + printf("after normalize getting IP \n"); + ///////////////////////////////////////////////////////////////////////////////////////// // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds ///////////////////////////////////////////////////////////////////////////////////////// @@ -724,7 +768,7 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ if(fEfficiencyFunction){ dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps); - dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection; + dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection; } else if(dataContainerV0){ dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction); @@ -732,6 +776,7 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ } + /////////////// // Unfold ////////////// @@ -754,13 +799,19 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ ///////////////////// // Simply correct //////////////////// + AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps); ////////// // Plot ////////// + if(fDebugLevel > 0.0) { + + Int_t ptpr = 0; + if(fBeamType==0) ptpr=0; + if(fBeamType==1) ptpr=1; TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700); ccorrected->Divide(2,1); @@ -787,8 +838,8 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p"); legcorrected->Draw("same"); ccorrected->cd(2); - TH1D *correctedTH1D = correctedspectrum->Projection(0); - TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0); + TH1D *correctedTH1D = correctedspectrum->Projection(ptpr); + TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr); TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone(); ratiocorrected->SetName("ratiocorrected"); ratiocorrected->SetTitle(""); @@ -799,15 +850,17 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ ratiocorrected->Draw(); if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps"); - if(fNonHFEsyst){ - CalculateNonHFEsyst(); + if(fBeamType == 0){ + if(fNonHFEsyst){ + CalculateNonHFEsyst(0); + } } - // Dump to file if needed if(fDumpToFile) { - + // to do centrality dependent + TFile *out; out = new TFile("finalSpectrum.root","recreate"); out->cd(); @@ -824,16 +877,73 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){ alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum"); alltogetherCorrection->Write(); // - if(unnormalizedRawSpectrum) { - unnormalizedRawSpectrum->SetName("beautyAfterIP"); - unnormalizedRawSpectrum->Write(); - } + unnormalizedRawSpectrum->SetName("beautyAfterIP"); + unnormalizedRawSpectrum->Write(); if(gNormalizedRawSpectrum){ gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP"); gNormalizedRawSpectrum->Write(); } + if(fBeamType==1) { + + TGraphErrors* correctedspectrumDc[kCentrality]; + TGraphErrors* alltogetherspectrumDc[kCentrality]; + for(Int_t i=0;iGetAxis(0)->SetRange(i+1,i+1); + correctedspectrumDc[i] = Normalize(correctedspectrum,i); + correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i)); + correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i)); + alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1); + alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i); + alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i)); + alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i)); + correctedspectrumDc[i]->Write(); + alltogetherspectrumDc[i]->Write(); + + + TH1D *centrcrosscheck = correctedspectrum->Projection(0); + centrcrosscheck->SetTitle(Form("centrality_%i",i)); + centrcrosscheck->SetName(Form("centrality_%i",i)); + centrcrosscheck->Write(); + + TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr); + TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr); + + TH1D *centrcrosscheck2 = (TH1D *) alltogetherCorrection->Project(0); + centrcrosscheck2->SetTitle(Form("centrality2_%i",i)); + centrcrosscheck2->SetName(Form("centrality2_%i",i)); + centrcrosscheck2->Write(); + + TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone(); + ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1); + ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i)); + ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i)); + ratiocorrectedc->Write(); + + fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i)); + fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i)); + fEfficiencyCharmSigD[i]->Write(); + fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i)); + fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i)); + fEfficiencyBeautySigD[i]->Write(); + fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i)); + fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i)); + fCharmEff[i]->Write(); + fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i)); + fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i)); + fBeautyEff[i]->Write(); + fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i)); + fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i)); + fConversionEff[i]->Write(); + fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i)); + fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i)); + fNonHFEEff[i]->Write(); + } + + } + out->Close(); delete out; } @@ -847,7 +957,20 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ // // Apply background subtraction // - + + Int_t ptpr = 0; + Int_t nbins=1; + if(fBeamType==0) + { + ptpr=0; + nbins=1; + } + if(fBeamType==1) + { + ptpr=1; + nbins=2; + } + // Raw spectrum AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ @@ -872,27 +995,39 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ if(!fInclusiveSpectrum){ //Background subtraction for IP analysis - TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(0); + TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr); CorrectFromTheWidth(measuredTH1Draw); TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500); rawspectra->cd(); - rawspectra->SetLogx(); rawspectra->SetLogy(); TLegend *lRaw = new TLegend(0.65,0.65,0.95,0.95); measuredTH1Draw->SetMarkerStyle(20); measuredTH1Draw->Draw(); lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum"); + TH1D* htemp; + Int_t* bins=new Int_t[1]; if(fIPanaHadronBgSubtract){ // Hadron background - printf("Hadron background for IP analysis subtracted!\n"); - Int_t* bins=new Int_t[1]; - TH1D* htemp = (TH1D *) fHadronEffbyIPcut->Projection(0); - bins[0]=htemp->GetNbinsX(); - AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",1,bins); + printf("Hadron background for IP analysis subtracted!\n"); + if(fBeamType==0) + { + + htemp = (TH1D *) fHadronEffbyIPcut->Projection(0); + bins[0]=htemp->GetNbinsX(); + } + if(fBeamType==1) + { + bins=new Int_t[2]; + htemp = (TH1D *) fHadronEffbyIPcut->Projection(0); + bins[0]=htemp->GetNbinsX(); + htemp = (TH1D *) fHadronEffbyIPcut->Projection(1); + bins[1]=htemp->GetNbinsX(); + } + AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins); hbgContainer->SetGrid(fHadronEffbyIPcut); backgroundGrid->Multiply(hbgContainer,1); // draw raw hadron bg spectra - TH1D *hadronbg= (TH1D *) backgroundGrid->Project(0); + TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr); CorrectFromTheWidth(hadronbg); hadronbg->SetMarkerColor(7); hadronbg->SetMarkerStyle(20); @@ -907,7 +1042,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ printf("Charm background for IP analysis subtracted!\n"); AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground(); // draw charm bg spectra - TH1D *charmbg= (TH1D *) charmbgContainer->Project(0); + TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr); CorrectFromTheWidth(charmbg); charmbg->SetMarkerColor(3); charmbg->SetMarkerStyle(20); @@ -919,9 +1054,9 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ } if(fIPanaConversionBgSubtract){ // Conversion background - AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground(); + AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground(); // draw conversion bg spectra - TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(0); + TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr); CorrectFromTheWidth(conversionbg); conversionbg->SetMarkerColor(4); conversionbg->SetMarkerStyle(20); @@ -936,7 +1071,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ // NonHFE background AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground(); // draw Dalitz/dielectron bg spectra - TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(0); + TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr); CorrectFromTheWidth(nonhfebg); nonhfebg->SetMarkerColor(6); nonhfebg->SetMarkerStyle(20); @@ -947,7 +1082,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ spectrumSubtracted->Add(nonHFEbgContainer,-1.0); printf("Non HFE background subtraction is preliminary!\n"); } - TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(0); + TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr); CorrectFromTheWidth(rawbgsubtracted); rawbgsubtracted->SetMarkerStyle(24); rawspectra->cd(); @@ -968,16 +1103,16 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ if(fDebugLevel > 0) { - Int_t ptpr; - if(fBeamType==0) ptpr=0; - if(fBeamType==1) ptpr=1; + Int_t ptprd; + if(fBeamType==0) ptprd=0; + if(fBeamType==1) ptprd=1; TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700); cbackgroundsubtraction->Divide(3,1); cbackgroundsubtraction->cd(1); gPad->SetLogy(); - TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptpr); - TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr); + TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd); + TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd); CorrectFromTheWidth(measuredTH1Daftersubstraction); CorrectFromTheWidth(measuredTH1Dbeforesubstraction); measuredTH1Daftersubstraction->SetStats(0); @@ -1019,7 +1154,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){ } ratiomeasuredcontamination->Draw("P"); cbackgroundsubtraction->cd(3); - TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptpr); + TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd); CorrectFromTheWidth(measuredTH1background); measuredTH1background->SetStats(0); measuredTH1background->SetTitle(""); @@ -1133,9 +1268,20 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){ // // calculate charm background // + Int_t ptpr = 0; + Int_t nDim = 1; + if(fBeamType==0) + { + ptpr=0; + } + if(fBeamType==1) + { + ptpr=1; + nDim=2; + } Double_t evtnorm=0; - if(fNMCbgEvents) evtnorm= double(fNEvents[0])/double(fNMCbgEvents); + if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]); AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC); if(!mcContainer){ @@ -1150,29 +1296,58 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){ AliCFDataGrid *charmBackgroundGrid= 0x0; charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID - TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0); + TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0); Int_t* bins=new Int_t[1]; bins[0]=charmbgaftertofpid->GetNbinsX(); - AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",1,bins); - ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency - TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0); - charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm); - TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0); + if(fBeamType==1) + { + charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0); + bins=new Int_t[2]; + bins[0]=charmbgaftertofpid->GetNbinsX(); + charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1); + bins[1]=charmbgaftertofpid->GetNbinsX(); - AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",1,bins); + } + + AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins); + ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency + TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr); + + if(fBeamType==0)charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm); + Double_t contents[2]; + AliCFDataGrid *eventTemp = new AliCFDataGrid("eventTemp","eventTemp",nDim,bins); + if(fBeamType==1) + { + for(Int_t kCentr=0;kCentrFill(contents,evtnormPbPb); + } + } + charmBackgroundGrid->Multiply(eventTemp,1); + } + TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr); + + AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins); weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors - TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0); + TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr); charmBackgroundGrid->Multiply(weightedCharmContainer,1.); - TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0); + TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr); // Efficiency (set efficiency to 1 for only folding) AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer); efficiencyD->CalculateEfficiency(0,0); - // Folding - Int_t nDim = 1; + // Folding + if(fBeamType==0)nDim = 1; + if(fBeamType==1)nDim = 2; AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid()); folding.SetMaxNumberOfIterations(1); folding.Unfold(); @@ -1181,7 +1356,7 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){ THnSparse* result1= folding.GetEstMeasured(); // folded spectra THnSparse* result=(THnSparse*)result1->Clone(); charmBackgroundGrid->SetGrid(result); - TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0); + TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr); //Charm background evaluation plots @@ -1189,7 +1364,18 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){ cCharmBgEval->Divide(3,1); cCharmBgEval->cd(1); - charmbgaftertofpid->Scale(evtnorm); + + if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm); + if(fBeamType==1) + { + Double_t evtnormPbPb=0; + for(Int_t kCentr=0;kCentrScale(evtnormPbPb); + } + CorrectFromTheWidth(charmbgaftertofpid); charmbgaftertofpid->SetMarkerStyle(25); charmbgaftertofpid->Draw("p"); @@ -1248,15 +1434,15 @@ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){ // Double_t evtnorm[1] = {0.0}; - if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents); + if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]); printf("check event!!! %lf \n",evtnorm[0]); AliCFContainer *backgroundContainer = 0x0; if(fNonHFEsyst){ - backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0]->Clone(); + backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone(); for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){ - backgroundContainer->Add(fConvSourceContainer[iSource][0]); + backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent } } else{ @@ -1269,22 +1455,46 @@ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){ } Int_t stepbackground = 3; + if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4; + AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground); - //backgroundGrid->Scale(evtnorm); - //workaround for statistical errors: "Scale" method does not work for our errors, since Sumw2 is not defined for AliCFContainer - Int_t nBin[1]; - - Int_t bins[1]; - bins[0]=fConversionEff->GetNbinsX(); - - for(Long_t iBin=1; iBin<= bins[0];iBin++){ - nBin[0]=iBin; - backgroundGrid->SetElementError(nBin, backgroundGrid->GetElementError(nBin)*evtnorm[0]); - backgroundGrid->SetElement(nBin,backgroundGrid->GetElement(nBin)*evtnorm[0]); + Int_t *nBinpp=new Int_t[1]; + Int_t* binspp=new Int_t[1]; + binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins + + Int_t *nBinPbPb=new Int_t[2]; + Int_t* binsPbPb=new Int_t[2]; + binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins + binsPbPb[0]=12; + + Int_t looppt=binspp[0]; + if(fBeamType==1) looppt=binsPbPb[1]; + + for(Long_t iBin=1; iBin<= looppt;iBin++){ + if(fBeamType==0) + { + nBinpp[0]=iBin; + backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]); + backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]); + } + if(fBeamType==1) + { + // loop over centrality + for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){ + nBinPbPb[0]=iiBin; + nBinPbPb[1]=iBin; + Double_t evtnormPbPb=0; + if(fNMCbgEvents[iiBin]) evtnormPbPb= double(fNEvents[iiBin])/double(fNMCbgEvents[iiBin]); + backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb); + backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb); + } + } } //end of workaround for statistical errors - - AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,&bins[0]); + + AliCFDataGrid *weightedConversionContainer; + if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp); + else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb); weightedConversionContainer->SetGrid(GetPIDxIPEff(2)); backgroundGrid->Multiply(weightedConversionContainer,1.0); @@ -1299,14 +1509,14 @@ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){ // Double_t evtnorm[1] = {0.0}; - if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents); + if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]); printf("check event!!! %lf \n",evtnorm[0]); AliCFContainer *backgroundContainer = 0x0; if(fNonHFEsyst){ - backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0]->Clone(); + backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone(); for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){ - backgroundContainer->Add(fNonHFESourceContainer[iSource][0]); + backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]); } } else{ @@ -1320,22 +1530,45 @@ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){ Int_t stepbackground = 3; + if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4; + AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground); - //backgroundGrid->Scale(evtnorm); - //workaround for statistical errors: "Scale" method does not work for our errors, since Sumw2 is not defined for AliCFContainer - Int_t nBin[1]; - - Int_t bins[1]; - bins[0]=fConversionEff->GetNbinsX(); - - for(Long_t iBin=1; iBin<= bins[0];iBin++){ - nBin[0]=iBin; - backgroundGrid->SetElementError(nBin, backgroundGrid->GetElementError(nBin)*evtnorm[0]); - backgroundGrid->SetElement(nBin,backgroundGrid->GetElement(nBin)*evtnorm[0]); + Int_t *nBinpp=new Int_t[1]; + Int_t* binspp=new Int_t[1]; + binspp[0]=fConversionEff[0]->GetNbinsX(); // number of pt bins + + Int_t *nBinPbPb=new Int_t[2]; + Int_t* binsPbPb=new Int_t[2]; + binsPbPb[1]=fConversionEff[0]->GetNbinsX(); // number of pt bins + binsPbPb[0]=12; + + Int_t looppt=binspp[0]; + if(fBeamType==1) looppt=binsPbPb[1]; + + + for(Long_t iBin=1; iBin<= looppt;iBin++){ + if(fBeamType==0) + { + nBinpp[0]=iBin; + backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]); + backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]); + } + if(fBeamType==1) + { + for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){ + nBinPbPb[0]=iiBin; + nBinPbPb[1]=iBin; + Double_t evtnormPbPb=0; + if(fNMCbgEvents[iiBin]) evtnormPbPb= double(fNEvents[iiBin])/double(fNMCbgEvents[iiBin]); + backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb); + backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb); + } + } } //end of workaround for statistical errors - - AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,&bins[0]); + AliCFDataGrid *weightedNonHFEContainer; + if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp); + else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb); weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3)); backgroundGrid->Multiply(weightedNonHFEContainer,1.0); @@ -1361,16 +1594,13 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons AliError("Data Container not available"); return NULL; } - dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData); } - AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone(); result->SetName("ParametrizedEfficiencyBefore"); THnSparse *h = result->GetGrid(); Int_t nbdimensions = h->GetNdimensions(); //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions); - AliCFContainer *dataContainer = GetContainer(kDataContainer); if(!dataContainer){ AliError("Data Container not available"); @@ -1384,7 +1614,6 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons memset(coord, 0, sizeof(Int_t) * nbdimensions); Double_t* points = new Double_t[nbdimensions]; - ULong64_t nEntries = h->GetNbins(); for (ULong64_t i = 0; i < nEntries; ++i) { @@ -1421,7 +1650,7 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData); if(fDebugLevel > 0) { - + TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700); cEfficiencyParametrized->Divide(2,1); cEfficiencyParametrized->cd(1); @@ -1458,9 +1687,9 @@ AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* cons //ratioefficiency->Draw(); if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps"); + } - return resultt; } @@ -1673,15 +1902,18 @@ TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){ // Efficiency AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer); efficiencyD->CalculateEfficiency(fStepMC,fStepTrue); - - // Consider parameterized IP cut efficiency - if(!fInclusiveSpectrum){ - Int_t* bins=new Int_t[1]; - bins[0]=35; - AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins); - beffContainer->SetGrid(GetBeautyIPEff()); - efficiencyD->Multiply(beffContainer,1); + if(!fBeauty2ndMethod) + { + // Consider parameterized IP cut efficiency + if(!fInclusiveSpectrum){ + Int_t* bins=new Int_t[1]; + bins[0]=35; + + AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins); + beffContainer->SetGrid(GetBeautyIPEff()); + efficiencyD->Multiply(beffContainer,1); + } } @@ -1774,14 +2006,18 @@ AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpe AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer); efficiencyD->CalculateEfficiency(fStepMC,fStepTrue); - // Consider parameterized IP cut efficiency - if(!fInclusiveSpectrum){ - Int_t* bins=new Int_t[1]; - bins[0]=35; - - AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins); - beffContainer->SetGrid(GetBeautyIPEff()); - efficiencyD->Multiply(beffContainer,1); + + if(!fBeauty2ndMethod) + { + // Consider parameterized IP cut efficiency + if(!fInclusiveSpectrum){ + Int_t* bins=new Int_t[1]; + bins[0]=35; + + AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins); + beffContainer->SetGrid(GetBeautyIPEff()); + efficiencyD->Multiply(beffContainer,1); + } } // Data in the right format @@ -1958,6 +2194,7 @@ TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) TH1D* projection = (TH1D *) spectrum->Project(ptpr); CorrectFromTheWidth(projection); TGraphErrors *graphError = NormalizeTH1(projection,i); + return graphError; } @@ -1990,7 +2227,6 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const { dp = input->GetXaxis()->GetBinWidth(ibin)/2.; n = input->GetBinContent(ibin); dN = input->GetBinError(ibin); - // New point nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n; errdN = 1./(2. * TMath::Pi() * p); @@ -2005,7 +2241,7 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const { spectrumNormalized->SetMarkerStyle(22); spectrumNormalized->SetMarkerColor(kBlue); spectrumNormalized->SetLineColor(kBlue); - + return spectrumNormalized; } @@ -2038,7 +2274,6 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) cons dp = input->GetXaxis()->GetBinWidth(ibin)/2.; n = input->GetBinContent(ibin); dN = input->GetBinError(ibin); - // New point nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n; errdN = 1./(2. * TMath::Pi() * p); @@ -2080,13 +2315,14 @@ AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type) return dynamic_cast(fCFContainers->At(type)); } //____________________________________________________________ -AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge) { +AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centrality) { // // Slice bin for a given source of electron // nDim is the number of dimension the corrections are done // dimensions are the definition of the dimensions // source is if we want to keep only one MC source (-1 means we don't cut on the MC source) // positivenegative if we want to keep positive (1) or negative (0) or both (-1) + // centrality (-1 means we do not cut on centrality) // Double_t *varMin = new Double_t[container->GetNVar()], @@ -2124,7 +2360,12 @@ AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, In fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1])); AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0])); } - + if(ivar == 5){ + if((centrality>= 0) && (centralityGetNBins(ivar))) { + varMin[ivar] = binLimits[centrality]; + varMax[ivar] = binLimits[centrality]; + } + } delete[] binLimits; } @@ -2255,86 +2496,353 @@ THnSparse* AliHFEspectrum::GetCharmWeights(){ // Measured D->e based weighting factors // - const Int_t nDim=1; - Int_t nBin[nDim] = {35}; - + const Int_t nDimpp=1; + Int_t nBinpp[nDimpp] = {35}; Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.}; + const Int_t nDimPbPb=2; + Int_t nBinPbPb[nDimPbPb] = {11,35}; + Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.}; + Int_t loopcentr=1; + Int_t looppt=nBinpp[0]; + if(fBeamType==0) + { + fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp); + fWeightCharm->SetBinEdges(0, ptbinning1); + } + if(fBeamType==1) + { + fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb); + fWeightCharm->SetBinEdges(1, ptbinning1); + fWeightCharm->SetBinEdges(0, kCentralityRange); + loopcentr=nBinPbPb[0]; + } + - fWeightCharm = new THnSparseF("weightHisto", "weiting factor; pt[GeV/c]", nDim, nBin); - for(Int_t idim = 0; idim < nDim; idim++) - fWeightCharm->SetBinEdges(idim, ptbinning1); Double_t weight[35]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385}; //{1.11865,1.11444,1.13395,1.15751,1.13412,1.14446,1.15372,1.09458,1.13446,1.11707,1.1352,1.10979,1.08887,1.11164,1.10772,1.10727,1.07027,1.07016,1.04826,1.03248,1.0093,0.973534,0.932574,0.869838,0.799316,0.734015,0.643001,0.584319,0.527147,0.495661,0.470002,0.441129,0.431156,0.417242,0.39246,0.379611,0.390845,0.368857,0.34959,0.387186,0.376364,0.384437,0.349585,0.383225}; //points Double_t pt[1]; - for(int i=0; iFill(pt,weight[i]); + Double_t contents[2]; + + for(int icentr=0; icentrFill(contents,weight[i]); + } + } + + Int_t nDimSparse = fWeightCharm->GetNdimensions(); + Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable + Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse + + for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); + nBins *= binsvar[iVar]; + } + + Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) + // loop that sets 0 error in each bin + for (Long_t iBin=0; iBinSetBinError(binfill,0.); // put 0 everywhere } - Int_t* ibins[nDim]; - for(Int_t ivar = 0; ivar < nDim; ivar++) - ibins[ivar] = new Int_t[nBin[ivar] + 1]; - fWeightCharm->SetBinError(ibins[0],0); return fWeightCharm; } //____________________________________________________________________________ -THnSparse* AliHFEspectrum::GetBeautyIPEff(){ - // - // Return beauty electron IP cut efficiency - // +void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, Int_t *dimensions){ + + // TOF PID efficiencies + Int_t ptpr; + if(fBeamType==0) ptpr=0; + if(fBeamType==1) ptpr=1; + + Int_t loopcentr=1; + const Int_t nCentralitybinning=11; //number of centrality bins + if(fBeamType==1) + { + loopcentr=nCentralitybinning; + } - const Int_t nDim=1; - Int_t nBin[nDim] = {35}; + TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,20.); + TF1 *fipfit2 = new TF1("fipfit2","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,10.0); + TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",1.0,7.); + + TCanvas * cefficiencyParam = new TCanvas("efficiencyParam","efficiencyParam",1200,600); + cefficiencyParam->Divide(2,1); + cefficiencyParam->cd(1); + + AliCFContainer *mccontainermcD = 0x0; + TH1D* efficiencysigTOFPIDD[nCentralitybinning]; + TH1D* efficiencyTOFPIDD[nCentralitybinning]; + Int_t source = -1; //get parameterized TOF PID efficiencies + + for(int icentr=0; icentrCalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies + + // mb sample for double check + if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); + if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1); + AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD); + efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies + + efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr); + efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr); + efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr)); + efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr)); + + fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); + fittofpid->SetLineColor(3); + efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R"); + efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency"); + fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid"); + } - Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.}; + // draw (for PbPb, only 1st bin) + efficiencysigTOFPIDD[0]->SetTitle(""); + efficiencysigTOFPIDD[0]->SetStats(0); + efficiencysigTOFPIDD[0]->SetMarkerStyle(25); + efficiencysigTOFPIDD[0]->SetMarkerColor(2); + efficiencysigTOFPIDD[0]->SetLineColor(2); + efficiencysigTOFPIDD[0]->Draw(); - THnSparseF *ipcut = new THnSparseF("beff", "b eff; pt[GeV/c]", nDim, nBin); - for(Int_t idim = 0; idim < nDim; idim++) - ipcut->SetBinEdges(idim, ptbinning1); - Double_t pt[1]; - Double_t weight; - for(int i=0; iFill(pt,weight); - } - Int_t* ibins[nDim]; - for(Int_t ivar = 0; ivar < nDim; ivar++) - ibins[ivar] = new Int_t[nBin[ivar] + 1]; - ipcut->SetBinError(ibins[0],0); + fEfficiencyTOFPIDD[0]->Draw("same"); - return ipcut; + efficiencyTOFPIDD[0]->SetTitle(""); + efficiencyTOFPIDD[0]->SetStats(0); + efficiencyTOFPIDD[0]->SetMarkerStyle(24); + efficiencyTOFPIDD[0]->Draw("same"); + + + cefficiencyParam->cd(2); + + // IP cut efficiencies + + for(int icentr=0; icentrCalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. + + source = 1; //beauty + if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen); + if(fBeamType==1) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); + AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD); + efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. + + source = 0; //charm + if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); + if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); + AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD); + efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. + + source = 1; //beauty + if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); + if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); + AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD); + efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. + + source = 2; //conversion + if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); + if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); + AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD); + efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. + + source = 3; //non HFE except for the conversion + if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen); + if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1); + AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD); + efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. + + fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); + fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); + fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); + fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); + fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); + fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); + } + + for(int icentr=0; icentrSetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); + fipfit->SetLineColor(2); + fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R"); + fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency"); + if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); + else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit"); + + if(fIPParameterizedEff){ + fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); + fipfit->SetLineColor(1); + fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R"); + fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency"); + fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit"); + + fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); + fipfit2->SetLineColor(3); + fConversionEff[icentr]->Fit(fipfit2,"R"); + fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency"); + fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfit2"); + + fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08); + fipfit2->SetLineColor(4); + fNonHFEEff[icentr]->Fit(fipfit2,"R"); + fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency"); + fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfit2"); + } + } + + // draw (for PbPb, only 1st bin) + + fEfficiencyCharmSigD[0]->SetMarkerStyle(21); + fEfficiencyCharmSigD[0]->SetMarkerColor(1); + fEfficiencyCharmSigD[0]->SetLineColor(1); + fEfficiencyBeautySigD[0]->SetMarkerStyle(21); + fEfficiencyBeautySigD[0]->SetMarkerColor(2); + fEfficiencyBeautySigD[0]->SetLineColor(2); + fCharmEff[0]->SetMarkerStyle(24); + fCharmEff[0]->SetMarkerColor(1); + fCharmEff[0]->SetLineColor(1); + fBeautyEff[0]->SetMarkerStyle(24); + fBeautyEff[0]->SetMarkerColor(2); + fBeautyEff[0]->SetLineColor(2); + fConversionEff[0]->SetMarkerStyle(24); + fConversionEff[0]->SetMarkerColor(3); + fConversionEff[0]->SetLineColor(3); + fNonHFEEff[0]->SetMarkerStyle(24); + fNonHFEEff[0]->SetMarkerColor(4); + fNonHFEEff[0]->SetLineColor(4); + + fEfficiencyCharmSigD[0]->Draw(); + fEfficiencyBeautySigD[0]->Draw("same"); + fCharmEff[0]->Draw("same"); + fBeautyEff[0]->Draw("same"); + fConversionEff[0]->Draw("same"); + fNonHFEEff[0]->Draw("same"); + + fEfficiencyIPBeautyD[0]->Draw("same"); + if(fIPParameterizedEff){ + fEfficiencyIPCharmD[0]->Draw("same"); + fEfficiencyIPConversionD[0]->Draw("same"); + fEfficiencyIPNonhfeD[0]->Draw("same"); + } } //____________________________________________________________________________ -THnSparse* AliHFEspectrum::GetCharmEff(){ +THnSparse* AliHFEspectrum::GetBeautyIPEff(){ // - // Return charm electron IP cut efficiency + // Return beauty electron IP cut efficiency // - const Int_t nDim=1; - Int_t nBin[nDim] = {35}; + const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning + const Int_t nCentralitybinning=11;//number of centrality bins + Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits + Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.}; + Int_t ptpr = 0; + Int_t nDim=1; //dimensions of the efficiency weighting grid + if(fBeamType==1) + { + ptpr=1; + nDim=2; //dimensions of the efficiency weighting grid + } + Int_t nBin[1] = {nPtbinning1}; + Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1}; - Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.}; - THnSparseF *ipcut = new THnSparseF("ceff", "c eff; pt[GeV/c]", nDim, nBin); + THnSparseF *ipcut; + if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin); + else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb); + for(Int_t idim = 0; idim < nDim; idim++) - ipcut->SetBinEdges(idim, ptbinning1); + { + if(nDim==1) ipcut->SetBinEdges(idim, kPtRange); + if(nDim==2) + { + ipcut->SetBinEdges(0, kCentralityRange); + ipcut->SetBinEdges(1, kPtRange); + } + } Double_t pt[1]; - Double_t weight; - for(int i=0; iFill(pt,weight); - } - Int_t* ibins[nDim]; - for(Int_t ivar = 0; ivar < nDim; ivar++) - ibins[ivar] = new Int_t[nBin[ivar] + 1]; - ipcut->SetBinError(ibins[0],0); + Double_t weight[nCentralitybinning]; + Double_t contents[2]; + + for(Int_t a=0;a<11;a++) + { + weight[a] = 1.0; + } + + + Int_t looppt=nBin[0]; + Int_t loopcentr=1; + if(fBeamType==1) + { + loopcentr=nBinPbPb[0]; + } + + + for(int icentr=0; icentrEval(pt[0]); + else{ + printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr); + weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1); + } + + if(fBeamType==1){ + contents[0]=icentr; + contents[1]=pt[0]; + } + if(fBeamType==0){ + contents[0]=pt[0]; + } + ipcut->Fill(contents,weight[icentr]); + } + } + + + Int_t nDimSparse = ipcut->GetNdimensions(); + Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable + Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse + + for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); + nBins *= binsvar[iVar]; + } + + Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) + // loop that sets 0 error in each bin + for (Long_t iBin=0; iBinSetBinError(binfill,0.); // put 0 everywhere + } return ipcut; } @@ -2344,54 +2852,255 @@ THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){ // // Return PID x IP cut efficiency // + const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning + const Int_t nCentralitybinning=11;//number of centrality bins + Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits + Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.}; + Int_t ptpr = 0; + Int_t nDim=1; //dimensions of the efficiency weighting grid + if(fBeamType==1) + { + ptpr=1; + nDim=2; //dimensions of the efficiency weighting grid + } + Int_t nBin[1] = {nPtbinning1}; + Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1}; + + THnSparseF *pideff; + if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin); + else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb); + for(Int_t idim = 0; idim < nDim; idim++) + { + + if(nDim==1) pideff->SetBinEdges(idim, kPtRange); + if(nDim==2) + { + pideff->SetBinEdges(0, kCentralityRange); + pideff->SetBinEdges(1, kPtRange); + } + } - const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning - - const Int_t nDim=1;//dimensions of the efficiency weighting grid - Int_t nBin[nDim] = {nPtbinning1}; + Double_t pt[1]; + Double_t weight[nCentralitybinning]; + Double_t contents[2]; + + for(Int_t a=0;a<11;a++) + { + weight[a] = 1.0; + } + + Int_t looppt=nBin[0]; + Int_t loopcentr=1; + if(fBeamType==1) + { + loopcentr=nBinPbPb[0]; + } + + for(int icentr=0; icentrEval(0); // assume we have constant TRD+TPC PID efficiency + for(int i=0; iEval(pt[0]); + else{ + printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr); + } + + //tof pid eff x tpc pid eff x ip cut eff + if(fIPParameterizedEff){ + if(source==0) { + if(fEfficiencyIPCharmD[icentr]) + weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]); + else{ + printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr); + weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1); + } + } + else if(source==2) { + if(fEfficiencyIPConversionD[icentr]) + weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]); + else{ + printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr); + weight[icentr] = tofpideff*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); + } + } + else if(source==3) { + if(fEfficiencyIPNonhfeD[icentr]) + weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]); + else{ + printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr); + weight[icentr] = tofpideff*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); + } + } + } + else{ + if(source==0) weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1); //charm + else if(source==2) weight[icentr] = tofpideff*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion + else if(source==3) weight[icentr] = tofpideff*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz + } + + if(fBeamType==1){ + contents[0]=icentr; + contents[1]=pt[0]; + } + if(fBeamType==0){ + contents[0]=pt[0]; + } + + pideff->Fill(contents,weight[icentr]); + } + } + + Int_t nDimSparse = pideff->GetNdimensions(); + Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable + Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse + + for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); + nBins *= binsvar[iVar]; + } + + Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) + // loop that sets 0 error in each bin + for (Long_t iBin=0; iBinSetBinError(binfill,0.); // put 0 everywhere + } - Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits - THnSparseF *pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin); - for(Int_t idim = 0; idim < nDim; idim++) - pideff->SetBinEdges(idim, kPtRange); - Double_t pt[1]; - Double_t weight = 1.0; - Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency - for(int i=0; iGetBinContent(i+1); // multiply ip cut eff for conversion electrons - if(source==3) weight = tofeff*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // multiply ip cut eff for Dalitz/dielectrons - //printf("tofeff= %lf trdtpcPidEfficiency= %lf conversion= %lf nonhfe= %lf\n",tofeff,trdtpcPidEfficiency,fConversionEff->GetBinContent(i+1),fNonHFEEff->GetBinContent(i+1)); - - pideff->Fill(pt,weight); - } - Int_t* ibins[nDim]; - for(Int_t ivar = 0; ivar < nDim; ivar++) - ibins[ivar] = new Int_t[nBin[ivar] + 1]; - pideff->SetBinError(ibins[0],0); return pideff; } + +//__________________________________________________________________________ +AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){ + // + // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method + // + Int_t ptpr = 0; + Int_t nDim = 1; + if(fBeamType==0) + { + ptpr=0; + } + if(fBeamType==1) + { + ptpr=1; + nDim=2; + } + + const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning + const Int_t nCentralitybinning=11;//number of centrality bins + Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits + Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.}; + Int_t nBin[1] = {nPtbinning1}; + Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1}; + //fBSpectrum2ndMethod->GetNbinsX() + + AliCFDataGrid *rawBeautyContainer; + if(fBeamType==0) rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin); + else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb); + // printf("number of bins= %d\n",bins[0]); + + + + + THnSparseF *brawspectra; + if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin); + else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb); + if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange); + if(fBeamType==1) + { + // brawspectra->SetBinEdges(0, centralityBins); + brawspectra->SetBinEdges(0, kCentralityRange); + brawspectra->SetBinEdges(1, kPtRange); + } + + Double_t pt[1]; + Double_t yields= 0.; + Double_t valuesb[2]; + + //Int_t looppt=nBin[0]; + Int_t loopcentr=1; + if(fBeamType==1) + { + loopcentr=nBinPbPb[0]; + } + + for(int icentr=0; icentrGetNbinsX(); i++){ + pt[0]=(kPtRange[i]+kPtRange[i+1])/2.; + + yields = fBSpectrum2ndMethod->GetBinContent(i+1); + + if(fBeamType==1) + { + valuesb[0]=icentr; + valuesb[1]=pt[0]; + } + if(fBeamType==0) valuesb[0]=pt[0]; + brawspectra->Fill(valuesb,yields); + } + } + + + + Int_t nDimSparse = brawspectra->GetNdimensions(); + Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable + Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse + + for (Int_t iVar=0; iVarGetAxis(iVar)->GetNbins(); + nBins *= binsvar[iVar]; + } + + Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates) + // loop that sets 0 error in each bin + for (Long_t iBin=0; iBinSetBinError(binfill,0.); // put 0 everywhere + } + + + rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency + TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr); + + new TCanvas; + fBSpectrum2ndMethod->SetMarkerStyle(24); + fBSpectrum2ndMethod->Draw("p"); + hRawBeautySpectra->SetMarkerStyle(25); + hRawBeautySpectra->Draw("samep"); + + return rawBeautyContainer; +} + //__________________________________________________________________________ -void AliHFEspectrum::CalculateNonHFEsyst(){ +void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){ + // + // Calculate non HFE sys // - //Calculate systematic errors for non-HFE and conversion background, - //based on correlated errors from pi0 input and uncorrelated errors - //from m_t scaling // + if(!fNonHFEsyst) return; Double_t evtnorm[1] = {0.0}; - if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents); + if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]); AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels]; AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels]; @@ -2402,22 +3111,23 @@ void AliHFEspectrum::CalculateNonHFEsyst(){ Int_t stepbackground = 3; Int_t* bins=new Int_t[1]; - bins[0]=fConversionEff->GetNbinsX(); + + bins[0]=fConversionEff[centrality]->GetNbinsX(); + AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins); AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins); - for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ - + for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){ for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){ - convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground); + convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),*fConvSourceContainer[iSource][iLevel][centrality],stepbackground); weightedConversionContainer->SetGrid(GetPIDxIPEff(2)); convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0); - - nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground); + + nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),*fNonHFESourceContainer[iSource][iLevel][centrality],stepbackground); weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3)); nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0); } - + bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone(); for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){ bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]); @@ -2427,11 +3137,12 @@ void AliHFEspectrum::CalculateNonHFEsyst(){ for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){//add other sources to get overall background from all meson decays bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]); } - + bgLevelGrid[iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone(); bgLevelGrid[iLevel]->Add(bgNonHFEGrid[iLevel]); } + //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated) AliCFDataGrid *bgErrorGrid[2]; bgErrorGrid[0] = (AliCFDataGrid*)bgLevelGrid[1]->Clone(); @@ -2588,5 +3299,7 @@ void AliHFEspectrum::CalculateNonHFEsyst(){ gOverallSystErrLow->Write(); output->Close(); - delete output; + delete output; + } + diff --git a/PWGHF/hfe/AliHFEspectrum.h b/PWGHF/hfe/AliHFEspectrum.h index 8a1b70c2f87..15f2c96a563 100644 --- a/PWGHF/hfe/AliHFEspectrum.h +++ b/PWGHF/hfe/AliHFEspectrum.h @@ -52,7 +52,8 @@ class AliHFEspectrum : public TNamed{ enum{ kElecBgSources = 6, kBgLevels = 3, - kBgPtBins = 44 + kBgPtBins = 44, + kCentrality = 12 }; enum Chargetype_t{ @@ -87,10 +88,12 @@ class AliHFEspectrum : public TNamed{ void SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type); void SetEfficiencyFunction(TF1 *efficiencyFunction) { fEfficiencyFunction = efficiencyFunction; }; void SetPbPbAnalysis(Bool_t isPbPb = kFALSE) { fBeamType=(Char_t) isPbPb; }; + + void SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, Int_t *dimensions); void SetNumberOfEvents(Int_t nEvents,Int_t i = 0) { fNEvents[i] = nEvents; }; void SetNumberOfMCEvents(Int_t nEvents) { fNMCEvents = nEvents; }; - void SetNumberOfMC2Events(Int_t nEvents) { fNMCbgEvents = nEvents; }; + void SetNumberOfMC2Events(Int_t nEvents,Int_t i = 0) { fNMCbgEvents[i] = nEvents; }; void SetMCEffStep(Int_t step) { fStepMC = step; }; void SetMCTruthStep(Int_t step) { fStepTrue = step; }; void SetStepToCorrect(Int_t step) { fStepData = step; }; @@ -106,8 +109,10 @@ class AliHFEspectrum : public TNamed{ void SetLowHighBoundaryCentralityBinAtTheEnd(Int_t low, Int_t high, Int_t i) { fLowBoundaryCentralityBinAtTheEnd[i] = low; fHighBoundaryCentralityBinAtTheEnd[i] = high;}; void SetBeautyAnalysis() { fInclusiveSpectrum = kFALSE; }; + void CallInputFileForBeauty2ndMethod(); + void SetInputFileForBeauty2ndMethod(const char *filenameb = "BSpectrum2ndmethod.root"){fkBeauty2ndMethodfilename = filenameb; }; + void SetBeautyAnalysis2ndMethod(Bool_t beauty2ndmethod) { fBeauty2ndMethod = beauty2ndmethod; } void SetHadronEffbyIPcut(THnSparseF* hsHadronEffbyIPcut) { fHadronEffbyIPcut = hsHadronEffbyIPcut;}; - void SetNonHFEBackground2ndMethod() { fNonHFEbgMethod2 = kTRUE; }; void SetNonHFEsyst(Bool_t syst){ fNonHFEsyst = syst; }; void SetStepGuessedUnfolding(Int_t stepGuessedUnfolding) { fStepGuessedUnfolding = stepGuessedUnfolding; }; @@ -117,24 +122,26 @@ class AliHFEspectrum : public TNamed{ void SetDebugLevel(Int_t debugLevel, Bool_t writeToFile = kFALSE) { fDebugLevel = debugLevel; fWriteToFile = writeToFile; }; + + AliCFDataGrid* GetRawBspectra2ndMethod(); AliCFDataGrid* GetCharmBackground(); AliCFDataGrid* GetConversionBackground(); AliCFDataGrid* GetNonHFEBackground(); THnSparse* GetCharmWeights(); THnSparse* GetBeautyIPEff(); - THnSparse* GetCharmEff(); THnSparse* GetPIDxIPEff(Int_t source); - void CalculateNonHFEsyst(); + void CalculateNonHFEsyst(Int_t centrality = 0); void EnableIPanaHadronBgSubtract() { fIPanaHadronBgSubtract = kTRUE; }; void EnableIPanaCharmBgSubtract() { fIPanaCharmBgSubtract = kTRUE; }; void EnableIPanaConversionBgSubtract() { fIPanaConversionBgSubtract = kTRUE; }; void EnableIPanaNonHFEBgSubtract() { fIPanaNonHFEBgSubtract = kTRUE; }; + void EnableIPParameterizedEff() { fIPParameterizedEff = kTRUE; }; protected: AliCFContainer *GetContainer(AliHFEspectrum::CFContainer_t contt); - AliCFContainer *GetSlicedContainer(AliCFContainer *cont, Int_t ndim, Int_t *dimensions,Int_t source=-1,Chargetype_t charge=kAllCharge); + AliCFContainer *GetSlicedContainer(AliCFContainer *cont, Int_t ndim, Int_t *dimensions,Int_t source=-1,Chargetype_t charge=kAllCharge,Int_t centrality=-1); THnSparseF *GetSlicedCorrelation(THnSparseF *correlationmatrix,Int_t nDim, Int_t *dimensions) const; TObject* GetSpectrum(const AliCFContainer * const c, Int_t step); TObject* GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0); @@ -154,11 +161,16 @@ class AliHFEspectrum : public TNamed{ THnSparseF *fCorrelation; // Correlation Matrices AliCFDataGrid *fBackground; // Background Grid TF1 *fEfficiencyFunction; // Efficiency Function + TF1 *fEfficiencyTOFPIDD[kCentrality]; // TOF PID efficiency parameterized + TF1 *fEfficiencyIPCharmD[kCentrality]; // IP efficiency parameterized for charm + TF1 *fEfficiencyIPBeautyD[kCentrality]; // IP efficiency parameterized for beauty + TF1 *fEfficiencyIPConversionD[kCentrality]; // IP efficiency parameterized for conversion + TF1 *fEfficiencyIPNonhfeD[kCentrality]; // IP efficiency parameterized for nonhfe THnSparseF *fWeightCharm; // Weight for charm bg - AliCFContainer *fConvSourceContainer[kElecBgSources][kBgLevels]; //container for conversion electrons, divided into different photon sources - AliCFContainer *fNonHFESourceContainer[kElecBgSources][kBgLevels]; //container for non-HF electrons, divided into different sources + AliCFContainer *fConvSourceContainer[kElecBgSources][kBgLevels][kCentrality]; //container for conversion electrons, divided into different photon sources + AliCFContainer *fNonHFESourceContainer[kElecBgSources][kBgLevels][kCentrality]; //container for non-HF electrons, divided into different sources Bool_t fInclusiveSpectrum; // Inclusive Spectrum Bool_t fDumpToFile; // Write Result in a file @@ -171,13 +183,14 @@ class AliHFEspectrum : public TNamed{ Bool_t fIPanaCharmBgSubtract; // Charm background subtraction Bool_t fIPanaConversionBgSubtract; // Conversion background subtraction Bool_t fIPanaNonHFEBgSubtract; // nonHFE except for conversion background subtraction - Bool_t fNonHFEbgMethod2; // switch for 2nd method to subtract non HFE background + Bool_t fIPParameterizedEff; // switch to use parameterized efficiency for ip analysis Bool_t fNonHFEsyst; // choose NonHFE background level (upper, lower, central) + Bool_t fBeauty2ndMethod; // 2nd method to get beauty spectrum Int_t fNbDimensions; // Number of dimensions for the correction Int_t fNEvents[20]; // Number of Events Int_t fNMCEvents; // Number of MC Events - Int_t fNMCbgEvents; // Number of BG MC Events + Int_t fNMCbgEvents[20]; // Number of BG MC Events Int_t fStepMC; // MC step (for unfolding) Int_t fStepTrue; // MC step of the final spectrum Int_t fStepData; // Data Step (various applications) @@ -195,8 +208,14 @@ class AliHFEspectrum : public TNamed{ Int_t fHighBoundaryCentralityBinAtTheEnd[20]; // Boundary of the bins THnSparseF *fHadronEffbyIPcut;// container for hadron efficiency by IP cut - TH1D *fConversionEff; // conversion IP cut eff - TH1D *fNonHFEEff; // nonhfe IP cut eff + TH1D *fEfficiencyCharmSigD[kCentrality]; // charm IP cut eff from signal enhanced MC + TH1D *fEfficiencyBeautySigD[kCentrality]; // beauty IP cut eff from signal enhanced MC + TH1D *fConversionEff[kCentrality]; // conversion IP cut eff + TH1D *fNonHFEEff[kCentrality]; // nonhfe IP cut eff + TH1D *fCharmEff[kCentrality]; // charm IP cut eff + TH1D *fBeautyEff[kCentrality]; // beauty IP cut eff + TH1D *fBSpectrum2ndMethod; // beauty spectrum for 2nd method + const char *fkBeauty2ndMethodfilename; // name of file, which contains beauty spectrum for 2ndmethod Char_t fBeamType; // beamtype; default -1; pp =0; PbPb=1 diff --git a/PWGHF/hfe/AliHFEtools.cxx b/PWGHF/hfe/AliHFEtools.cxx index c2a3d090312..ec20ab39a7c 100644 --- a/PWGHF/hfe/AliHFEtools.cxx +++ b/PWGHF/hfe/AliHFEtools.cxx @@ -23,9 +23,14 @@ #include #include #include "TH1.h" +#include "TH1D.h" #include "TH2.h" +#include "TGraph.h" +#include "TGraphErrors.h" #include "THnSparse.h" #include "TAxis.h" +#include "TMath.h" +#include "TString.h" #include "AliAODMCParticle.h" #include "AliAODpidUtil.h" @@ -242,3 +247,159 @@ Int_t AliHFEtools::PDG2AliPID(Int_t pdg){ }; return pid; } + +//__________________________________________ +TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize) +{ + // Creates a TH1D from TGraph g. The binwidth of the first bin has to + // specified. The others bins are calculated automatically. Supports also Graphs + // with non constant x steps. The axis of the Graph can be exchanged if + // exchange=kTRUE (modifies the Graph). + + TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize); + if( result == 0) return result; + + //------------------------------------------ + // setup the final hist + Int_t nBinX = g->GetN(); + Double_t* err = g->GetEY(); // error y is still ok even if exchanged + for(Int_t i = 0; i < nBinX; i ++){ + result->SetBinError(i + 1, err[i]); + } + if(exchange){ + AliHFEtools::ExchangeXYGraph(g); // undo what has been done in GraphToHist + AliHFEtools::ExchangeXYGraphErrors(g); // now exchange including errors + } + + return result; +} + +//__________________________________________ +Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g) +{ + // exchanges x-values and y-values. + if(g==0) return kFALSE; + Int_t nbin=g->GetN(); + Double_t x,y; + for(Int_t i = 0; i < nbin; i ++) + { + g->GetPoint(i,x,y); + g->SetPoint(i,y,x); + } + + return kTRUE; +} + +//__________________________________________ +Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g) +{ + // exchanges x-values and y-values and + // corresponding errors. + if(g==0) return kFALSE; + Int_t nbin=g->GetN(); + Double_t x,y; + Double_t ex,ey; + for(Int_t i = 0; i < nbin; i ++) + { + g->GetPoint(i,x,y); + ex = g->GetErrorX(i); + ey = g->GetErrorY(i); + g->SetPoint(i,y,x); + g->SetPointError(i,ey,ex); + } + + return kTRUE; + +} + +//__________________________________________ +TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize) +{ + // Creates a TH1D from TGraph g. The binwidth of the first bin has to + // specified. The others bins are calculated automatically. Supports also Graphs + // with non constant x steps. The axis of the Graph can be exchanged if + // exchange=kTRUE (modifies the Graph). + + + TH1D* result = 0; + if(g == 0) return result; + if(firstBinWidth == -1) return result; + TString myname=""; + myname = g->GetName(); + myname += "_graph"; + if(exchange) AliHFEtools::ExchangeXYGraph(g); + + Int_t nBinX = g->GetN(); + Double_t* x = g->GetX(); + Double_t* y = g->GetY(); + + if(nBinX < 1) return result; + + //------------------------------------------ + // create the Matrix for the equation system + // and init the values + + Int_t nDim = nBinX - 1; + TMatrixD a(nDim,nDim); + TMatrixD b(nDim,1); + + Double_t* aA = a.GetMatrixArray(); + Double_t* aB = b.GetMatrixArray(); + memset(aA,0,nDim * nDim * sizeof(Double_t)); + memset(aB,0,nDim * sizeof(Double_t)); + //------------------------------------------ + + //------------------------------------------ + // setup equation system + // width for 1st bin is given therefore + // we shift bin parameter (column) by one to the left + // to reduce the matrix size + + Double_t* xAxis = new Double_t [nBinX + 1]; + Double_t* binW = new Double_t [nBinX ]; + binW[0] = firstBinWidth; + + aB[0] = x[1] - x[0] - 0.5 * binW[0]; + aA[0] = 0.5; + + for(Int_t col = 1; col < nDim ; col ++) + { + Int_t row = col; + aB[col] = x[col + 1] - x[ col ]; + aA[row * nDim + col - 1 ] = 0.5; + aA[row * nDim + col ] = 0.5; + } + //------------------------------------------ + + //------------------------------------------ + // solve the equations + a.Invert(); + TMatrixD c = a * b; + //------------------------------------------ + + //------------------------------------------ + // calculate the bin boundaries + xAxis[0] = x[0] - 0.5 * binW[0]; + memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t)); + for(Int_t col = 0; col < nBinX ; col ++) { + xAxis[col + 1] = x[col] + 0.5 * binW[col]; + } + //------------------------------------------ + + //------------------------------------------ + // setup the final hist + result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis); + for(Int_t i = 0; i < nBinX; i ++){ + result->SetBinContent(i + 1, y[i]); + } + result->SetMarkerColor(markercolor); + result->SetMarkerStyle(markerstyle); + result->SetMarkerSize(markersize); + //------------------------------------------ + + delete [] xAxis; + delete [] binW; + + + return result; +} diff --git a/PWGHF/hfe/AliHFEtools.h b/PWGHF/hfe/AliHFEtools.h index b88469f76e5..f548bcb3109 100644 --- a/PWGHF/hfe/AliHFEtools.h +++ b/PWGHF/hfe/AliHFEtools.h @@ -26,6 +26,10 @@ class TParticle; class AliAODMCParticle; class AliPIDResponse; class AliVParticle; +class TGraph; +class TGraphErrors; +class TH1D; +class TString; class AliHFEtools : public TObject{ public: @@ -42,6 +46,10 @@ class AliHFEtools : public TObject{ static AliPIDResponse *GetDefaultPID(Bool_t isMC = kTRUE, Bool_t isESD = kTRUE); static void DestroyDefaultPID(); static void SetLogLevel(Int_t loglevel) { fgLogLevel = loglevel ;} + static TH1D* GraphErrorsToHist(TGraphErrors* g = 0, Double_t firstBinWidth = -1, Bool_t exchange=kFALSE, Int_t markerstyle=8, Int_t markercolor=2, Float_t markersize=0.7); + static TH1D* GraphToHist(TGraph* g = 0, Double_t firstBinWidth = -1, Bool_t exchange=kFALSE, Int_t markerstyle=8, Int_t markercolor=2, Float_t markersize=0.7); + static Bool_t ExchangeXYGraph(TGraph* g = 0); + static Bool_t ExchangeXYGraphErrors(TGraphErrors* g = 0); private: AliHFEtools(const AliHFEtools &); diff --git a/PWGHF/hfe/AliHFEtrdPIDqa.cxx b/PWGHF/hfe/AliHFEtrdPIDqa.cxx index cbfda4151f1..e402dfafb76 100644 --- a/PWGHF/hfe/AliHFEtrdPIDqa.cxx +++ b/PWGHF/hfe/AliHFEtrdPIDqa.cxx @@ -224,15 +224,15 @@ void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){ Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon); nbins[kElectronLike] = 100; nbins[kNClustersLike] = 200; - nbins[kCentralityBin] = 12; + nbins[kCentralityBinLike] = 12; Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon); Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon); binMin[kElectronLike] = 0.; binMin[kNClustersLike] = 0.; - binMin[kCentralityBin] = -1.; + binMin[kCentralityBinLike] = -1.; binMax[kElectronLike] = 1.; binMax[kNClustersLike] = 200.; - binMax[kCentralityBin] = 11.; + binMax[kCentralityBinLike] = 11.; fHistos->CreateTHnSparse("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax); fHistos->BinLogAxis("fLikeTRD", kP); @@ -247,13 +247,15 @@ void AliHFEtrdPIDqa::CreateQAHistogram(){ Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon); nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1; nbins[kNClusters] = 200; + nbins[kCentralityBinQA] = 12; Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon); binMin[kNonZeroTrackletCharge] = 0.; - binMin[kNClusters] = 0.; + binMin[kNClusters] = 0.; + binMin[kCentralityBinQA] = -1.; Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon); binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.; binMax[kNClusters] = 200.; - + binMax[kCentralityBinQA] = 11.; fHistos->CreateTHnSparse("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax); fHistos->BinLogAxis("fQAtrack", kP); } @@ -267,15 +269,18 @@ void AliHFEtrdPIDqa::CreatedEdxHistogram(){ Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon); nbins[kdEdx] = 100; nbins[kNclusters] = 261; - nbins[kNonZeroSlices] = 9; + nbins[kNonZeroSlices] = 9; + nbins[kCentralityBindEdx] = 12; Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon); binMin[kdEdx] = 0.; binMin[kNclusters] = 0; binMin[kNonZeroSlices] = 0.; + binMin[kCentralityBindEdx] = -1.; Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon); binMax[kdEdx] = 10000.; binMax[kNclusters] = 260.; binMax[kNonZeroSlices] = 8.; + binMax[kCentralityBindEdx] = 11.; fHistos->CreateTHnSparse("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax); fHistos->BinLogAxis("fQAdEdx", kP); @@ -292,14 +297,17 @@ void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){ nbins[kTPCdEdx] = 600; nbins[kTRDdEdxMethod1] = 1000; nbins[kTRDdEdxMethod2] = 1000; + nbins[kCentralityBinTruncMean] = 12; Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon); binMin[kTPCdEdx] = 0.; binMin[kTRDdEdxMethod1] = 0.; - binMin[kTRDdEdxMethod2] = 0.; + binMin[kTRDdEdxMethod2] = 0.; + binMin[kCentralityBinTruncMean] = -1.; Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon); binMax[kTPCdEdx] = 600; binMax[kTRDdEdxMethod1] = 20000.; binMax[kTRDdEdxMethod2] = 20000.; + binMax[kCentralityBinTruncMean] = 11.; fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax); fHistos->BinLogAxis("fTRDtruncMean", kP); @@ -375,12 +383,13 @@ void AliHFEtrdPIDqa::FillTRDLikelihoods(const AliESDtrack * const track, Int_t s // p // ntracklets // Electron Likelihood + // Centrality quantities[kSpecies] = species; quantities[kP] = outerPars ? outerPars->P() : track->P(); quantities[kNTracklets] = track->GetTRDntrackletsPID(); quantities[kElectronLike] = likeEle; quantities[kNClustersLike] = track->GetTRDncls(); - quantities[kCentralityBin] = fCentralityBin; + quantities[kCentralityBinLike] = fCentralityBin; fHistos->Fill("fLikeTRD", quantities); } @@ -455,14 +464,17 @@ void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t speci quantitiesdEdx[kdEdx] = fTotalChargeInSlice0 ? track->GetTRDslice(iplane, 0) : dEdxSum; // hack by mfasel: In the new reconstruction, the total charge is stored in the first slice, in the old reconstruction it has to be calculated from the slice charges. if(dEdxSum) ntrackletsNonZero++; // Fill dEdx histogram + quantitiesdEdx[kCentralityBindEdx] = fCentralityBin; if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries } quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero; + quantitiesQA[kCentralityBinQA] = fCentralityBin; fHistos->Fill("fQAtrack", quantitiesQA); quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal(); quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, 0.6); quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, 0.6); + quantitiesTruncMean[kCentralityBinTruncMean] = fCentralityBin; fHistos->Fill("fTRDtruncMean", quantitiesTruncMean); } @@ -609,8 +621,8 @@ void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets, Int_t nCentrality, Bool hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, isGreaterEqual ? 7 : binTracklets); if(isPbPb){ - Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBin)->FindBin(nCentrality); - hLikeTRD->GetAxis(kCentralityBin)->SetRange(binCentrality, isGreaterEqual ? 11 : binCentrality); + Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(nCentrality); + hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, isGreaterEqual ? 11 : binCentrality); /* new TCanvas; TH2 *test = hLikeTRD->Projection(kCentralityBin, kP); @@ -637,7 +649,7 @@ void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets, Int_t nCentrality, Bool // Undo ranges hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins()); hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins()); - hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kCentralityBin)->GetNbins()); + hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kCentralityBinLike)->GetNbins()); // Prepare List for output Char_t *listname=Form("%dTracklets", nTracklets); @@ -842,8 +854,8 @@ Double_t AliHFEtrdPIDqa::CalculateIntegratedPionEfficiency(UInt_t nTracklets, Do hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets); if(icentrality!=1){ - Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBin)->FindBin(icentrality); - hLikeTRD->GetAxis(kCentralityBin)->SetRange(binCentrality, binCentrality); + Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(icentrality); + hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, binCentrality); } Int_t pbinMin = hLikeTRD->GetAxis(kP)->FindBin(pmax), diff --git a/PWGHF/hfe/AliHFEtrdPIDqa.h b/PWGHF/hfe/AliHFEtrdPIDqa.h index 656deb22463..7c5ef340a42 100644 --- a/PWGHF/hfe/AliHFEtrdPIDqa.h +++ b/PWGHF/hfe/AliHFEtrdPIDqa.h @@ -119,25 +119,28 @@ class AliHFEtrdPIDqa : public TNamed{ enum QuantitiesLike_t{ kElectronLike = 3, kNClustersLike = 4, - kCentralityBin = 5, + kCentralityBinLike = 5, kQuantitiesLike = 6 }; enum QuantitiesQAtrack_t{ kNonZeroTrackletCharge = 3, kNClusters = 4, - kQuantitiesQA = 5 + kCentralityBinQA = 5, + kQuantitiesQA = 6 }; enum QuantitiesdEdx_t{ kNclusters = 3, kNonZeroSlices = 4, kdEdx = 5, - kQuantitiesdEdx = 6 + kCentralityBindEdx = 6, + kQuantitiesdEdx = 7 }; enum QuantitiesTruncMean_t{ kTPCdEdx = 3, kTRDdEdxMethod1 = 4, kTRDdEdxMethod2 = 5, - kQuantitiesTruncMean = 6 + kCentralityBinTruncMean = 6, + kQuantitiesTruncMean = 7 }; void ProcessTrackESD(const AliESDtrack * const track, Int_t species); -- 2.43.0