From dbe3abbee80a0f97526838620af3c6720fe7dffd Mon Sep 17 00:00:00 2001 From: dainese Date: Tue, 23 Jun 2009 12:40:35 +0000 Subject: [PATCH] Fix coding rules violations --- PWG3/hfe/AliAnalysisTaskHFE.cxx | 589 +++++++++++++++++++++++--------- PWG3/hfe/AliAnalysisTaskHFE.h | 56 +-- PWG3/hfe/AliHFEcuts.cxx | 111 ++++-- PWG3/hfe/AliHFEcuts.h | 29 +- PWG3/hfe/AliHFEextraCuts.cxx | 71 +++- PWG3/hfe/AliHFEextraCuts.h | 7 + PWG3/hfe/AliHFEmcQA.cxx | 288 ++++++++++------ PWG3/hfe/AliHFEmcQA.h | 74 ++-- PWG3/hfe/AliHFEpidTRD.cxx | 160 +++------ PWG3/hfe/AliHFEpidTRD.h | 22 +- PWG3/hfe/AliHFEpriVtx.h | 6 +- PWG3/hfe/AliHFEsecVtx.cxx | 454 +++++++++++++++++++----- PWG3/hfe/AliHFEsecVtx.h | 90 +++-- 13 files changed, 1358 insertions(+), 599 deletions(-) diff --git a/PWG3/hfe/AliAnalysisTaskHFE.cxx b/PWG3/hfe/AliAnalysisTaskHFE.cxx index d2d02279fc4..c8cd573d17b 100644 --- a/PWG3/hfe/AliAnalysisTaskHFE.cxx +++ b/PWG3/hfe/AliAnalysisTaskHFE.cxx @@ -30,7 +30,6 @@ #include #include #include -#include #include #include #include @@ -67,19 +66,23 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(): , fMC(0x0) , fCFM(0x0) , fCorrelation(0x0) + , fFakeElectrons(0x0) , fPID(0x0) , fCuts(0x0) , fSecVtx(0x0) - , fAnalysisMCQA(0x0) + , fMCQA(0x0) , fNEvents(0x0) , fQA(0x0) , fOutput(0x0) , fHistMCQA(0x0) , fHistSECVTX(0x0) { - DefineInput(0, TChain::Class()); - DefineOutput(0, TH1I::Class()); - DefineOutput(1, TList::Class()); + // + // Default constructor + // + DefineInput(0, TChain::Class()); + DefineOutput(0, TH1I::Class()); + DefineOutput(1, TList::Class()); DefineOutput(2, TList::Class()); // Initialize cuts @@ -87,12 +90,55 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(): fPID = new AliHFEpid; } +//____________________________________________________________ +AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref): + AliAnalysisTask("PID efficiency Analysis COPY", "") + , fESD(0x0) + , fMC(0x0) + , fCFM(0x0) + , fCorrelation(0x0) + , fFakeElectrons(0x0) + , fPID(0x0) + , fCuts(0x0) + , fSecVtx(0x0) + , fMCQA(0x0) + , fNEvents(0x0) + , fQA(0x0) + , fOutput(0x0) + , fHistMCQA(0x0) + , fHistSECVTX(0x0) +{ + // + // Copy Constructor + // + DefineInput(0, TChain::Class()); + DefineOutput(0, TH1I::Class()); + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); + + ref.Copy(*this); +} + +//____________________________________________________________ +AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){ + // + // Assignment operator + // + if(this != &ref) + ref.Copy(*this); + + return *this; +} + //____________________________________________________________ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){ - if(fESD) delete fESD; - if(fMC) delete fMC; - if(fPID) delete fPID; - if(fQA){ + // + // Destructor + // + if(fESD) delete fESD; + if(fMC) delete fMC; + if(fPID) delete fPID; + if(fQA){ fQA->Clear(); delete fQA; } @@ -110,49 +156,65 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){ } if(fCuts) delete fCuts; if(fSecVtx) delete fSecVtx; - if(fAnalysisMCQA) delete fAnalysisMCQA; - if(fNEvents) delete fNEvents; + if(fMCQA) delete fMCQA; + if(fNEvents) delete fNEvents; if(fCorrelation) delete fCorrelation; if(fFakeElectrons) delete fFakeElectrons; } +//____________________________________________________________ +void AliAnalysisTaskHFE::Copy(TObject &o) const { + // + // Make a copy of the task objecy + // + + AliAnalysisTaskHFE &target = dynamic_cast(o); + + // copy objects + if(fPID) target.fPID = new AliHFEpid(*fPID); + if(fCuts) target.fCuts = new AliHFEcuts(*fCuts); + + if(fESD || fMC) target.ConnectInputData(0x0); + if(fOutput || fQA) target.CreateOutputObjects(); +} + //____________________________________________________________ void AliAnalysisTaskHFE::ConnectInputData(Option_t *){ - TTree *esdchain = dynamic_cast(GetInputData(0)); - if(!esdchain){ - AliError("ESD chain empty"); - return; - } else { - esdchain->SetBranchStatus("Tracks", 1); - } - AliESDInputHandler *esdH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); - if(!esdH){ - AliError("No ESD input handler"); - return; - } else { - fESD = esdH->GetEvent(); - } - AliMCEventHandler *mcH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); - if(!mcH){ - AliError("No MC truth handler"); - return; - } else { - fMC = mcH->MCEvent(); - } + TTree *esdchain = dynamic_cast(GetInputData(0)); + if(!esdchain){ + AliError("ESD chain empty"); + return; + } else { + esdchain->SetBranchStatus("Tracks", 1); + } + AliESDInputHandler *esdH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + if(!esdH){ + AliError("No ESD input handler"); + return; + } else { + fESD = esdH->GetEvent(); + } + AliMCEventHandler *mcH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); + if(!mcH){ + AliError("No MC truth handler"); + return; + } else { + fMC = mcH->MCEvent(); + } } //____________________________________________________________ void AliAnalysisTaskHFE::CreateOutputObjects(){ - fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram - // First Step: TRD alone - if(!fQA) fQA = new TList; - fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0); - fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1); - fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2); - fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3); - fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4); - fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5); - fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6); + fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram + // First Step: TRD alone + if(!fQA) fQA = new TList; + fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0); + fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1); + fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2); + fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3); + fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4); + fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5); + fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6); if(!fOutput) fOutput = new TList; // Initialize correction Framework and Cuts @@ -174,183 +236,393 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){ // Initialize PID if(IsQAOn()){ + AliInfo("QA on for Cuts and PID"); fPID->SetQAOn(); fQA->Add(fPID->GetQAhistograms()); } fPID->SetHasMCData(kTRUE); fPID->InitializePID("TRD"); - // [mj] mcQA---------------------------------- + // mcQA---------------------------------- if (IsMCQAOn()) { + AliInfo("MC QA on"); + if(!fMCQA) fMCQA = new AliHFEmcQA; if(!fHistMCQA) fHistMCQA = new TList(); + fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_"); // create histograms for charm + fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty + fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm + fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty + TIter next_(gDirectory->GetList()); + TObject *obj_; + int counter_ = 0; + TString objname; + while ((obj_ = next_.Next())) { + objname = obj_->GetName(); + TObjArray *toks = objname.Tokenize("_"); + if (toks->GetEntriesFast()){ + TObjString *fpart = (TObjString *)(toks->UncheckedAt(0)); + if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj_, counter_++); + } + } fQA->Add(fHistMCQA); - fAnalysisMCQA->CreateHistograms(AliHFEmcQA::fkCharm,"mcqa_"); // create histograms for charm - fAnalysisMCQA->CreateHistograms(AliHFEmcQA::fkBeauty,"mcqa_"); // create histograms for beauty - //Rossella -// fAnalysisMCQA->CreateHistosHadrons(); // create histograms for hadrons + } - } - // [mj] secvtx---------------------------- + // secvtx---------------------------------- if (IsSecVtxOn()) { - fOutput->Add(fHistSECVTX); + AliInfo("Secondary Vertex Analysis on"); + fSecVtx = new AliHFEsecVtx; + if(!fHistSECVTX) fHistSECVTX = new TList(); fSecVtx->CreateHistograms("secvtx_"); - } - - TIter next_(gDirectory->GetList()); - TObject *obj_; - int counter_ = 0; - int counter__ = 0; - TString objname; - while ((obj_ = next_.Next())) { - objname = obj_->GetName(); - TObjArray *toks = objname.Tokenize("_"); - if (toks->GetEntriesFast()){ - TObjString *fpart = (TObjString *)(toks->UncheckedAt(0)); - if ((fpart->String()).CompareTo("secvtx") == 0){ - fHistSECVTX->AddAt(obj_, counter_++); - } - else if ((fpart->String()).CompareTo("mcqa") == 0){ - fHistMCQA->AddAt(obj_, counter__++); + TIter next__(gDirectory->GetList()); + TObject *obj__; + int counter__ = 0; + TString objname; + while ((obj__ = next__.Next())) { + objname = obj__->GetName(); + TObjArray *toks = objname.Tokenize("_"); + if (toks->GetEntriesFast()){ + TObjString *fpart = (TObjString *)(toks->UncheckedAt(0)); + if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj__, counter__++); } } - } - //--------------------------------------- + fOutput->Add(fHistSECVTX); + } } //____________________________________________________________ void AliAnalysisTaskHFE::Exec(Option_t *){ - // - // Run the analysis - // - if(!fESD){ - AliError("No ESD Event"); - return; - } - if(!fMC){ - AliError("No MC Event"); - return; - } - fCFM->SetEventInfo(fMC); + // + // Run the analysis + // + if(!fESD){ + AliError("No ESD Event"); + return; + } + if(!fMC){ + AliError("No MC Event"); + return; + } + fCFM->SetEventInfo(fMC); fPID->SetMCEvent(fMC); - //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC); + //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC); + Double_t container[6]; - Double_t container[6]; + // Loop over the Monte Carlo tracks to see whether we have overlooked any track + AliMCParticle *mctrack = 0x0; + Int_t nElectrons = 0; - // Loop over the Monte Carlo tracks to see whether we have overlooked any track - AliMCParticle *mctrack = 0x0; - Int_t nElectrons = 0; + if (IsSecVtxOn()) { + fSecVtx->SetEvent(fESD); + fSecVtx->SetStack(fMC->Stack()); + } - // [mj] run MC QA ------------------------------------------------ + // run mc QA ------------------------------------------------ if (IsMCQAOn()) { + AliDebug(2, "Running MC QA"); - fAnalysisMCQA->SetStack(fMC->Stack()); - fAnalysisMCQA->Init(); + fMCQA->SetStack(fMC->Stack()); + fMCQA->Init(); Int_t nPrims = fMC->Stack()->GetNprimary(); Int_t nMCTracks = fMC->Stack()->GetNtrack(); // loop over primary particles for quark and heavy hadrons for (Int_t igen = 0; igen < nPrims; igen++){ - fAnalysisMCQA->GetQuarkKine(igen, AliHFEmcQA::fkCharm); - fAnalysisMCQA->GetQuarkKine(igen, AliHFEmcQA::fkBeauty); + fMCQA->GetQuarkKine(igen, AliHFEmcQA::kCharm); + fMCQA->GetQuarkKine(igen, AliHFEmcQA::kBeauty); + fMCQA->GetHadronKine(igen, AliHFEmcQA::kCharm); + fMCQA->GetHadronKine(igen, AliHFEmcQA::kBeauty); } - fAnalysisMCQA->EndOfEventAna(AliHFEmcQA::fkCharm); - fAnalysisMCQA->EndOfEventAna(AliHFEmcQA::fkBeauty); + fMCQA->EndOfEventAna(AliHFEmcQA::kCharm); + fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty); // loop over all tracks for decayed electrons for (Int_t igen = 0; igen < nMCTracks; igen++){ - fAnalysisMCQA->GetDecayedKine(igen, AliHFEmcQA::fkCharm, AliHFEmcQA::fkElectron); - fAnalysisMCQA->GetDecayedKine(igen, AliHFEmcQA::fkBeauty, AliHFEmcQA::fkElectron); + fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 0); + fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); + fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel + fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel } } // end of MC QA loop // ----------------------------------------------------------------- - for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){ - mctrack = fMC->GetTrack(imc); - container[0] = mctrack->Pt(); - container[1] = mctrack->Eta(); + // + // Loop MC + // + for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){ + mctrack = fMC->GetTrack(imc); + + container[0] = mctrack->Pt(); + container[1] = mctrack->Eta(); container[2] = mctrack->Phi(); - if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue; - fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated); - (dynamic_cast(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi()); - if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue; - // find the label in the vector - fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance); - nElectrons++; - } - (dynamic_cast(fQA->At(3)))->Fill(nElectrons); - - // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD); - AliESDtrack *track = 0x0, *htrack = 0x0; + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue; + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated); + (dynamic_cast(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi()); + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue; + // find the label in the vector + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance); + nElectrons++; + } + (dynamic_cast(fQA->At(3)))->Fill(nElectrons); + + // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD); + + + // + // Loop ESD + // + + Int_t nbtrackrec = fESD->GetNumberOfTracks(); + AliESDtrack *track = 0x0, *htrack = 0x0; Int_t pid = 0; - for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){ - track = fESD->GetTrack(itrack); - container[0] = track->Pt(); - container[1] = track->Eta(); + // For double counted tracks + Int_t *indexRecKineTPC = new Int_t[nbtrackrec]; + Int_t *indexRecKineITS = new Int_t[nbtrackrec]; + Int_t *indexRecPrim = new Int_t[nbtrackrec]; + Int_t *indexHFEcutsITS = new Int_t[nbtrackrec]; + Int_t *indexHFEcutsTPC = new Int_t[nbtrackrec]; + Int_t *indexHFEcutsTRD = new Int_t[nbtrackrec]; + Int_t *indexPid = new Int_t[nbtrackrec]; + + for(Int_t nb = 0; nb < nbtrackrec; nb++){ + indexRecKineTPC[nb] = -1; + indexRecKineITS[nb] = -1; + indexRecPrim[nb] = -1; + indexHFEcutsITS[nb] = -1; + indexHFEcutsTPC[nb] = -1; + indexHFEcutsTRD[nb] = -1; + indexPid[nb] = -1; + } + + Bool_t alreadyseen = kFALSE; + + for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){ + + track = fESD->GetTrack(itrack); + + container[0] = track->Pt(); + container[1] = track->Eta(); container[2] = track->Phi(); - if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKine, track)) continue; - fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKine); - - // Check TRD criterions (outside the correction framework) - if(track->GetTRDncls()){ - (dynamic_cast(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls()); - (dynamic_cast(fQA->At(1)))->Fill(track->GetAlpha()); // Check the acceptance without tight cuts - (dynamic_cast(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality()); - (dynamic_cast(fQA->At(5)))->Fill(container[0], track->GetTRDncls()); - } - if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue; + + Bool_t signal = kTRUE; + + // RecKine: TPC cuts + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineTPC, track)) continue; + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineTPC); + + + // Check if it is signal electrons + if(!(mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())))) continue; + + container[3] = mctrack->Pt(); + container[4] = mctrack->Eta(); + container[5] = mctrack->Phi(); + + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE; + + if(signal) { + alreadyseen = kFALSE; + for(Int_t nb = 0; nb < itrack; nb++){ + if(indexRecKineTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE; + } + indexRecKineTPC[itrack] = TMath::Abs(track->GetLabel()); + + fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineTPC + AliHFEcuts::kNcutESDSteps + 1)); + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 2*(AliHFEcuts::kNcutESDSteps + 1))); + if(alreadyseen) { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 4*(AliHFEcuts::kNcutESDSteps + 1))); + } + else { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 3*(AliHFEcuts::kNcutESDSteps + 1))); + } + } + + + // RecKine: ITS cuts + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITS, track)) continue; + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineITS); + if(signal) { + alreadyseen = kFALSE; + for(Int_t nb = 0; nb < itrack; nb++){ + if(indexRecKineITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE; + } + indexRecKineITS[itrack] = TMath::Abs(track->GetLabel()); + fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineITS + AliHFEcuts::kNcutESDSteps + 1)); + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 2*(AliHFEcuts::kNcutESDSteps + 1))); + if(alreadyseen) { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 4*(AliHFEcuts::kNcutESDSteps + 1))); + } + else { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 3*(AliHFEcuts::kNcutESDSteps + 1))); + } + } + + + // Check TRD criterions (outside the correction framework) + if(track->GetTRDncls()){ + (dynamic_cast(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls()); + (dynamic_cast(fQA->At(1)))->Fill(track->GetAlpha()); // Check the acceptance without tight cuts + (dynamic_cast(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality()); + (dynamic_cast(fQA->At(5)))->Fill(container[0], track->GetTRDncls()); + } + + + // RecPrim + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue; fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim); - if(fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcuts, track)) continue; - fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts); + if(signal) { + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1); + fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1)); + alreadyseen = kFALSE; + for(Int_t nb = 0; nb < itrack; nb++){ + if(indexRecPrim[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE; + } + indexRecPrim[itrack] = TMath::Abs(track->GetLabel()); + if(alreadyseen) { + fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 4*(AliHFEcuts::kNcutESDSteps + 1)); + } + else { + fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 3*(AliHFEcuts::kNcutESDSteps + 1)); + } + } + + + // HFEcuts: ITS layers cuts + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue; + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS); + if(signal) { + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1); + fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1)); + alreadyseen = kFALSE; + for(Int_t nb = 0; nb < itrack; nb++){ + if(indexHFEcutsITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE; + } + indexHFEcutsITS[itrack] = TMath::Abs(track->GetLabel()); + if(alreadyseen) { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 4*(AliHFEcuts::kNcutESDSteps + 1))); + } + else { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 3*(AliHFEcuts::kNcutESDSteps + 1))); + } + } + + // HFEcuts: TPC ratio clusters + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTPC, track)) continue; + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC); + if(signal) { + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC + AliHFEcuts::kNcutESDSteps + 1); + fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTPC + 2*(AliHFEcuts::kNcutESDSteps + 1)); + alreadyseen = kFALSE; + for(Int_t nb = 0; nb < itrack; nb++){ + if(indexHFEcutsTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE; + } + indexHFEcutsTPC[itrack] = TMath::Abs(track->GetLabel()); + if(alreadyseen) { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 4*(AliHFEcuts::kNcutESDSteps + 1))); + } + else { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 3*(AliHFEcuts::kNcutESDSteps + 1))); + } + } + + // HFEcuts: Nb of tracklets TRD + if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue; + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD); + if(signal) { + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 1); + fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1)); + alreadyseen = kFALSE; + for(Int_t nb = 0; nb < itrack; nb++){ + if(indexHFEcutsTRD[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE; + } + indexHFEcutsTRD[itrack] = TMath::Abs(track->GetLabel()); + if(alreadyseen) { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 4*(AliHFEcuts::kNcutESDSteps + 1))); + } + else { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 3*(AliHFEcuts::kNcutESDSteps + 1))); + } + } + + // track accepted, do PID - if(!fPID->IsSelected(track)) continue; - // Track selected: distinguish between true and fake - if(!(mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())))) continue; - if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){ - fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts + 1); + if(!fPID->IsSelected(track)) continue; + + + // Fill Containers + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 1); + if(signal) { + fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 2); + fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1)+1); + alreadyseen = kFALSE; + for(Int_t nb = 0; nb < itrack; nb++){ + if(indexPid[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE; + } + indexPid[itrack] = TMath::Abs(track->GetLabel()); + if(alreadyseen) { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1))); + } + else { + fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 3*(AliHFEcuts::kNcutESDSteps + 1))); + } // dimensions 3&4&5 : pt,eta,phi (MC) - container[3] = mctrack->Pt(); - container[4] = mctrack->Eta(); - container[5] = mctrack->Phi(); fCorrelation->Fill(container); + } + + + // Track selected: distinguish between true and fake + if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){ // pair analysis [mj] if (IsSecVtxOn()) { + AliDebug(2, "Running Secondary Vertex Analysis"); + fSecVtx->InitAnaPair(); for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){ htrack = fESD->GetTrack(jtrack); if ( itrack == jtrack ) continue; //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue; if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue; - fSecVtx->AnaPair(track, htrack, fSecVtx->PairCode(track,htrack), jtrack); + fSecVtx->AnaPair(track, htrack, jtrack); } - } + // based on the partner of e info, you run secandary vertexing function + fSecVtx->RunSECVTX(track); + } // end of pair analysis } else { // Fill THnSparse with the information for Fake Electrons switch(pid){ - case 13: container[3] = AliPID::kMuon; - case 211: container[3] = AliPID::kPion; - case 321: container[3] = AliPID::kKaon; - case 2212: container[3] = AliPID::kProton; + case 13: container[3] = AliPID::kMuon; + case 211: container[3] = AliPID::kPion; + case 321: container[3] = AliPID::kKaon; + case 2212: container[3] = AliPID::kProton; } fFakeElectrons->Fill(container); } - } - fNEvents->Fill(1); - - // Done!!! - PostData(0, fNEvents); - PostData(1, fOutput); - PostData(2, fQA); + } + fNEvents->Fill(1); + + // release memory + delete[] indexRecKineTPC; + delete[] indexRecKineITS; + delete[] indexRecPrim; + delete[] indexHFEcutsITS; + delete[] indexHFEcutsTPC; + delete[] indexHFEcutsTRD; + delete[] indexPid; + + // Done!!! + PostData(0, fNEvents); + PostData(1, fOutput); + PostData(2, fQA); } //____________________________________________________________ void AliAnalysisTaskHFE::Terminate(Option_t *){ - // - // Terminate not implemented at the moment - // + // + // Terminate not implemented at the moment + // } //____________________________________________________________ @@ -360,13 +632,13 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){ // link it // const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi - const Double_t kPtmin = 0., kPtmax = 10.; + const Double_t kPtmin = 0.1, kPtmax = 10.; const Double_t kEtamin = -0.9, kEtamax = 0.9; const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi(); //arrays for the number of bins in each dimension Int_t iBin[kNvar]; - iBin[0] = 20; //bins in pt + iBin[0] = 40; //bins in pt iBin[1] = 8; //bins in eta iBin[2] = 18; // bins in phi @@ -376,14 +648,16 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){ binEdges[ivar] = new Double_t[iBin[ivar] + 1]; //values for bin lower bounds - for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/iBin[0]*(Double_t)i; + for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i); for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i; for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i; //one "container" for MC - AliCFContainer* container = new AliCFContainer("container","container for tracks", AliHFEcuts::kNcutSteps + 1, kNvar, iBin); + AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutSteps + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)), kNvar, iBin); + //setting the bin limits - for(Int_t ivar = 0; ivar < kNvar; ivar++) container -> SetBinLimits(0, binEdges[ivar]); + for(Int_t ivar = 0; ivar < kNvar; ivar++) + container -> SetBinLimits(ivar, binEdges[ivar]); fCFM->SetParticleContainer(container); //create correlation matrix for unfolding @@ -408,3 +682,4 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){ for(Int_t idim = 0; idim < kNvar; idim++) fFakeElectrons->SetBinEdges(idim, binEdges[idim]); } + diff --git a/PWG3/hfe/AliAnalysisTaskHFE.h b/PWG3/hfe/AliAnalysisTaskHFE.h index f158dbb16cd..56af78a002e 100644 --- a/PWG3/hfe/AliAnalysisTaskHFE.h +++ b/PWG3/hfe/AliAnalysisTaskHFE.h @@ -15,7 +15,13 @@ #ifndef ALIANALYSISTASKHFE_H #define ALIANALYSISTASKHFE_H +#ifndef ALIANALYSISTASK_H #include "AliAnalysisTask.h" +#endif + +#ifndef ROOT_THnSparse +#include +#endif class AliHFEpid; class AliHFEcuts; @@ -26,49 +32,57 @@ class AliESDEvent; class AliESDtrackCuts; class AliMCEvent; class TH1I; -class THnSparse; class TList; class AliAnalysisTaskHFE : public AliAnalysisTask{ enum{ - kIsQAOn = BIT(14), - kIsMCQAOn = BIT(15), - kIsSecVtxOn = BIT(16) + kIsQAOn = BIT(18), + kIsMCQAOn = BIT(19), + kIsSecVtxOn = BIT(20), + kIsPriVtxOn = BIT(21) }; - public: - AliAnalysisTaskHFE(); - ~AliAnalysisTaskHFE(); + public: + AliAnalysisTaskHFE(); + AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref); + AliAnalysisTaskHFE& operator=(const AliAnalysisTaskHFE &ref); + ~AliAnalysisTaskHFE(); - virtual void ConnectInputData(Option_t *); - virtual void CreateOutputObjects(); - virtual void Exec(Option_t *); - virtual void Terminate(Option_t *); + virtual void ConnectInputData(Option_t *); + virtual void CreateOutputObjects(); + virtual void Exec(Option_t *); + virtual void Terminate(Option_t *); Bool_t IsQAOn() const { return TestBit(kIsQAOn); }; - Bool_t IsMCQAOn() const { return TestBit(kIsMCQAOn); } + Bool_t IsMCQAOn() const { return TestBit(kIsMCQAOn); }; Bool_t IsSecVtxOn() const { return TestBit(kIsSecVtxOn); }; + Bool_t IsPriVtxOn() const { return TestBit(kIsPriVtxOn); }; void SetQAOn() { SetBit(kIsQAOn, kTRUE); }; void SetMCQAOn() { SetBit(kIsMCQAOn, kTRUE); }; + void SetPriVtxOn() { SetBit(kIsPriVtxOn, kTRUE); }; void SetSecVtxOn() { SetBit(kIsSecVtxOn, kTRUE); }; + + protected: + void Copy(TObject &o) const; - private: + private: void MakeParticleContainer(); - AliESDEvent *fESD; //! The ESD Event - AliMCEvent *fMC; //! The MC Event - AliCFManager *fCFM; //! Correction Framework Manager + AliESDEvent *fESD; //! The ESD Event + AliMCEvent *fMC; //! The MC Event + AliCFManager *fCFM; //! Correction Framework Manager THnSparseF *fCorrelation; //! response matrix for unfolding THnSparseF *fFakeElectrons; //! Contamination from Fake Electrons - AliHFEpid *fPID; //! PID + AliHFEpid *fPID; //! PID AliHFEcuts *fCuts; //! Cut Collection AliHFEsecVtx *fSecVtx; //! Secondary Vertex Analysis - AliHFEmcQA *fAnalysisMCQA; //! MC QA - TH1I *fNEvents; //! counter for the number of Events - TList *fQA; //! QA histos for the cuts + AliHFEmcQA *fMCQA; //! MC QA + TH1I *fNEvents; //! counter for the number of Events + TList *fQA; //! QA histos for the cuts TList *fOutput; //! Container for Task Output TList *fHistMCQA; //! Output container for MC QA histograms TList *fHistSECVTX; //! Output container for sec. vertexing results - ClassDef(AliAnalysisTaskHFE, 1) // The electron Analysis Task + ClassDef(AliAnalysisTaskHFE, 1) // The electron Analysis Task }; #endif + diff --git a/PWG3/hfe/AliHFEcuts.cxx b/PWG3/hfe/AliHFEcuts.cxx index 5dfa40a66e6..608c05b8e5b 100644 --- a/PWG3/hfe/AliHFEcuts.cxx +++ b/PWG3/hfe/AliHFEcuts.cxx @@ -18,6 +18,7 @@ * ALICE Heavy Flavour Electron Group * * * * Authors: * + * Raphaelle Bailhache * * Markus Fasel * * Markus Heide * * Matus Kalisky * @@ -42,8 +43,8 @@ ClassImp(AliHFEcuts) -const Int_t AliHFEcuts::kNcutSteps = 5; - +const Int_t AliHFEcuts::kNcutSteps = 8; +const Int_t AliHFEcuts::kNcutESDSteps = 6; //__________________________________________________________________ AliHFEcuts::AliHFEcuts(): fRequirements(0), @@ -54,7 +55,8 @@ AliHFEcuts::AliHFEcuts(): fMinClusterRatioTPC(0.), fSigmaToVtx(0.), fHistQA(0x0), - fCutList(0x0) + fCutList(0x0), + fDebugLevel(0) { // // Default Constructor @@ -77,7 +79,8 @@ AliHFEcuts::AliHFEcuts(const AliHFEcuts &c): fMinClusterRatioTPC(c.fMinClusterRatioTPC), fSigmaToVtx(c.fSigmaToVtx), fHistQA(0x0), - fCutList(0x0) + fCutList(0x0), + fDebugLevel(0) { // // Copy Constructor @@ -105,19 +108,28 @@ AliHFEcuts::~AliHFEcuts(){ //__________________________________________________________________ void AliHFEcuts::Initialize(AliCFManager *cfm){ + // Call all the setters for the cuts SetParticleGenCutList(); SetAcceptanceCutList(); - SetRecKineCutList(); + SetRecKineTPCCutList(); + SetRecKineITSCutList(); SetRecPrimaryCutList(); - SetHFElectronCuts(); + SetHFElectronITSCuts(); + SetHFElectronTPCCuts(); + SetHFElectronTRDCuts(); + // Connect the cuts cfm->SetParticleCutsList(kStepMCGenerated, dynamic_cast(fCutList->FindObject("fPartGenCuts"))); cfm->SetParticleCutsList(kStepMCInAcceptance, dynamic_cast(fCutList->FindObject("fPartAccCuts"))); - cfm->SetParticleCutsList(kStepRecKine, dynamic_cast(fCutList->FindObject("fPartRecCuts"))); + cfm->SetParticleCutsList(kStepRecKineTPC, dynamic_cast(fCutList->FindObject("fPartRecKineTPCCuts"))); + cfm->SetParticleCutsList(kStepRecKineITS, dynamic_cast(fCutList->FindObject("fPartRecKineITSCuts"))); cfm->SetParticleCutsList(kStepRecPrim, dynamic_cast(fCutList->FindObject("fPartPrimCuts"))); - cfm->SetParticleCutsList(kStepHFEcuts, dynamic_cast(fCutList->FindObject("fPartHFECuts"))); + cfm->SetParticleCutsList(kStepHFEcutsITS, dynamic_cast(fCutList->FindObject("fPartHFECutsITS"))); + cfm->SetParticleCutsList(kStepHFEcutsTPC, dynamic_cast(fCutList->FindObject("fPartHFECutsTPC"))); + cfm->SetParticleCutsList(kStepHFEcutsTRD, dynamic_cast(fCutList->FindObject("fPartHFECutsTRD"))); + } //__________________________________________________________________ @@ -125,9 +137,13 @@ void AliHFEcuts::Initialize(){ // Call all the setters for the cuts SetParticleGenCutList(); SetAcceptanceCutList(); - SetRecKineCutList(); + SetRecKineTPCCutList(); + SetRecKineITSCutList(); SetRecPrimaryCutList(); - SetHFElectronCuts(); + SetHFElectronITSCuts(); + SetHFElectronTPCCuts(); + SetHFElectronTRDCuts(); + } //__________________________________________________________________ @@ -186,7 +202,7 @@ void AliHFEcuts::SetAcceptanceCutList(){ } //__________________________________________________________________ -void AliHFEcuts::SetRecKineCutList(){ +void AliHFEcuts::SetRecKineTPCCutList(){ // // Track Kinematics and Quality cuts (Based on the Standard cuts from PWG0) // @@ -199,7 +215,6 @@ void AliHFEcuts::SetRecKineCutList(){ // Min. Number of Clusters: // TPC: 50 // RefitRequired: - // ITS // TPC // Chi2 per TPC cluster: 3.5 // @@ -210,7 +225,7 @@ void AliHFEcuts::SetRecKineCutList(){ AliCFTrackQualityCuts *trackQuality = new AliCFTrackQualityCuts("fCutsQualityRec","REC Track Quality Cuts"); trackQuality->SetMinNClusterTPC(fMinClustersTPC); trackQuality->SetMaxChi2PerClusterTPC(fMaxChi2clusterTPC); - trackQuality->SetStatus(AliESDtrack::kTPCrefit & AliESDtrack::kITSrefit); + trackQuality->SetStatus(AliESDtrack::kTPCrefit); trackQuality->SetMaxCovDiagonalElements(2., 2., 0.5, 0.5, 2); AliCFTrackKineCuts *kineCuts = new AliCFTrackKineCuts("fCutsKineRec", "REC Kine Cuts"); @@ -223,12 +238,33 @@ void AliHFEcuts::SetRecKineCutList(){ } TObjArray *recCuts = new TObjArray; - recCuts->SetName("fPartRecCuts"); + recCuts->SetName("fPartRecKineTPCCuts"); recCuts->AddLast(trackQuality); recCuts->AddLast(kineCuts); fCutList->AddLast(recCuts); } +//__________________________________________________________________ +void AliHFEcuts::SetRecKineITSCutList(){ + // + // Track Kinematics and Quality cuts (Based on the Standard cuts from PWG0) + // + // RefitRequired: + // ITS + // + AliCFTrackQualityCuts *trackQuality = new AliCFTrackQualityCuts("fCutsQualityRec","REC Track Quality Cuts"); + trackQuality->SetStatus(AliESDtrack::kITSrefit); + + if(IsInDebugMode()){ + trackQuality->SetQAOn(fHistQA); + } + + TObjArray *recCuts = new TObjArray; + recCuts->SetName("fPartRecKineITSCuts"); + recCuts->AddLast(trackQuality); + fCutList->AddLast(recCuts); +} + //__________________________________________________________________ void AliHFEcuts::SetRecPrimaryCutList(){ // @@ -239,11 +275,13 @@ void AliHFEcuts::SetRecPrimaryCutList(){ // No Kink daughters // AliCFTrackIsPrimaryCuts *primaryCut = new AliCFTrackIsPrimaryCuts("fCutsPrimaryCuts", "REC Primary Cuts"); +#ifdef V4_17_00 if(IsRequireDCAToVertex()){ primaryCut->SetDCAToVertex2D(kTRUE); primaryCut->SetMaxDCAToVertexXY(fDCAtoVtx[0]); primaryCut->SetMaxDCAToVertexZ(fDCAtoVtx[1]); } +#endif if(IsRequireSigmaToVertex()){ primaryCut->SetRequireSigmaToVertex(kTRUE); primaryCut->SetMaxNSigmaToVertex(fSigmaToVtx); @@ -258,24 +296,55 @@ void AliHFEcuts::SetRecPrimaryCutList(){ } //__________________________________________________________________ -void AliHFEcuts::SetHFElectronCuts(){ +void AliHFEcuts::SetHFElectronITSCuts(){ // - // Special Cuts introduced by the HFElectron Group + // Special Cuts introduced by the HFElectron Group: ITS // - AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroup","Extra cuts from the HFE group"); + AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupPixels","Extra cuts from the HFE group"); if(IsRequireITSpixel()){ hfecuts->SetRequireITSpixel(AliHFEextraCuts::ITSPixel_t(fCutITSPixel)); } -/* if(IsRequireDCAToVertex()){ +#ifndef V4_17_00 + if(IsRequireDCAToVertex()){ hfecuts->SetMaxImpactParamR(fDCAtoVtx[0]); hfecuts->SetMaxImpactParamZ(fDCAtoVtx[1]); - }*/ - if(fMinTrackletsTRD) hfecuts->SetMinTrackletsTRD(fMinTrackletsTRD); + } +#endif + + if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA); + + TObjArray *hfeCuts = new TObjArray; + hfeCuts->SetName("fPartHFECutsITS"); + hfeCuts->AddLast(hfecuts); + fCutList->AddLast(hfeCuts); +} + +//__________________________________________________________________ +void AliHFEcuts::SetHFElectronTPCCuts(){ + // + // Special Cuts introduced by the HFElectron Group: TPC + // + AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupTPC","Extra cuts from the HFE group"); if(fMinClusterRatioTPC > 0.) hfecuts->SetClusterRatioTPC(fMinClusterRatioTPC); if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA); TObjArray *hfeCuts = new TObjArray; - hfeCuts->SetName("fPartHFECuts"); + hfeCuts->SetName("fPartHFECutsTPC"); + hfeCuts->AddLast(hfecuts); + fCutList->AddLast(hfeCuts); +} + +//__________________________________________________________________ +void AliHFEcuts::SetHFElectronTRDCuts(){ + // + // Special Cuts introduced by the HFElectron Group: TRD + // + AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupTRD","Extra cuts from the HFE group"); + if(fMinTrackletsTRD > 0.) hfecuts->SetMinTrackletsTRD(fMinTrackletsTRD); + if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA); + + TObjArray *hfeCuts = new TObjArray; + hfeCuts->SetName("fPartHFECutsTRD"); hfeCuts->AddLast(hfecuts); fCutList->AddLast(hfeCuts); } diff --git a/PWG3/hfe/AliHFEcuts.h b/PWG3/hfe/AliHFEcuts.h index b0bce923a8a..9bf43f5caf9 100644 --- a/PWG3/hfe/AliHFEcuts.h +++ b/PWG3/hfe/AliHFEcuts.h @@ -47,12 +47,16 @@ class AliHFEcuts : public TObject{ typedef enum{ kStepMCGenerated = 0, kStepMCInAcceptance = 1, - kStepRecKine = 2, - kStepRecPrim = 3, - kStepHFEcuts = 4 - } CutStep_t; + kStepRecKineTPC = 2, + kStepRecKineITS = 3, + kStepRecPrim = 4, + kStepHFEcutsITS = 5, + kStepHFEcutsTPC = 6, + kStepHFEcutsTRD = 7 + } CutStep_t; static const Int_t kNcutSteps; + static const Int_t kNcutESDSteps; AliHFEcuts(); AliHFEcuts(const AliHFEcuts &c); @@ -97,14 +101,20 @@ class AliHFEcuts : public TObject{ void SetRequireITSPixel() { SETBIT(fRequirements, kITSPixel); } void SetRequireProdVertex() { SETBIT(fRequirements, kProductionVertex); }; void SetRequireSigmaToVertex() { SETBIT(fRequirements, kSigmaToVertex); }; + + void SetDebugLevel(Int_t level) { fDebugLevel = level; }; + Int_t GetDebugLevel() const { return fDebugLevel; }; private: void SetParticleGenCutList(); void SetAcceptanceCutList(); - void SetRecKineCutList(); + void SetRecKineTPCCutList(); + void SetRecKineITSCutList(); void SetRecPrimaryCutList(); - void SetHFElectronCuts(); - + void SetHFElectronITSCuts(); + void SetHFElectronTPCCuts(); + void SetHFElectronTRDCuts(); + ULong64_t fRequirements; // Bitmap for requirements Double_t fDCAtoVtx[2]; // DCA to Vertex Double_t fProdVtx[4]; // Production Vertex @@ -118,6 +128,8 @@ class AliHFEcuts : public TObject{ TList *fHistQA; //! QA Histograms TObjArray *fCutList; //! List of cut objects(Correction Framework Manager) + + Int_t fDebugLevel; // Debug Level ClassDef(AliHFEcuts, 1) // Container for HFE cuts }; @@ -166,7 +178,8 @@ void AliHFEcuts::CreateStandardCuts(){ fDCAtoVtx[1] = 10.; fMinClustersTPC = 50; fMinTrackletsTRD = 6; - fCutITSPixel = AliHFEextraCuts::kAny; + SetRequireITSPixel(); + fCutITSPixel = AliHFEextraCuts::kBoth; fMaxChi2clusterTPC = 3.5; fMinClusterRatioTPC = 0.6; fPtRange[0] = 0.1; diff --git a/PWG3/hfe/AliHFEextraCuts.cxx b/PWG3/hfe/AliHFEextraCuts.cxx index 00438915f73..dedbd795a1d 100644 --- a/PWG3/hfe/AliHFEextraCuts.cxx +++ b/PWG3/hfe/AliHFEextraCuts.cxx @@ -47,7 +47,8 @@ AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title): fClusterRatioTPC(0.), fMinTrackletsTRD(0), fPixelITS(0), - fQAlist(0x0) + fQAlist(0x0), + fDebugLevel(0) { // // Default Constructor @@ -63,7 +64,8 @@ AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c): fClusterRatioTPC(c.fClusterRatioTPC), fMinTrackletsTRD(c.fMinTrackletsTRD), fPixelITS(c.fPixelITS), - fQAlist(0x0) + fQAlist(0x0), + fDebugLevel(c.fDebugLevel) { // // Copy constructor @@ -89,6 +91,7 @@ AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){ fClusterRatioTPC = c.fClusterRatioTPC; fMinTrackletsTRD = c.fMinTrackletsTRD; fPixelITS = c.fPixelITS; + fDebugLevel = c.fDebugLevel; memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4); if(IsQAOn()){ @@ -144,6 +147,10 @@ Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){ trdTracklets = track->GetTRDpidQuality(); #endif UChar_t itsPixel = track->GetITSClusterMap(); + Int_t det, status1, status2; + Float_t xloc, zloc; + track->GetITSModuleIndexInfo(0, det, status1, xloc, zloc); + track->GetITSModuleIndexInfo(1, det, status2, xloc, zloc); if(TESTBIT(fRequirements, kMinImpactParamR)){ // cut on min. Impact Parameter in Radial direction if(impact_r >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR); @@ -170,18 +177,38 @@ Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){ } if(TESTBIT(fRequirements, kPixelITS)){ // cut on ITS pixel layers - if(fPixelITS == kFirst) + if(fDebugLevel > 0){ + printf("ITS cluster Map: "); + PrintBitMap(itsPixel); + } switch(fPixelITS){ - case kFirst: if(itsPixel & BIT(0)) SETBIT(survivedCut, kPixelITS); - break; - case kSecond: if(itsPixel & BIT(1)) SETBIT(survivedCut, kPixelITS); - break; - case kBoth: if((itsPixel & BIT(0)) && (itsPixel & BIT(1))) SETBIT(survivedCut, kPixelITS); - break; - case kAny: if((itsPixel & BIT(0)) || (itsPixel & BIT(1))) SETBIT(survivedCut, kPixelITS); - break; + case kFirst: + if(itsPixel & BIT(0)) + SETBIT(survivedCut, kPixelITS); + else if(!CheckITSstatus(status1)) + SETBIT(survivedCut, kPixelITS); + break; + case kSecond: + if(itsPixel & BIT(1)) + SETBIT(survivedCut, kPixelITS); + else if(!CheckITSstatus(status2)) + SETBIT(survivedCut, kPixelITS); + break; + case kBoth: + if((itsPixel & BIT(0)) && (itsPixel & BIT(1))) + SETBIT(survivedCut, kPixelITS); + else if(!CheckITSstatus(status1) && !CheckITSstatus(status2)) + SETBIT(survivedCut, kPixelITS); + break; + case kAny: + if((itsPixel & BIT(0)) || (itsPixel & BIT(1))) + SETBIT(survivedCut, kPixelITS); + else if(!CheckITSstatus(status1) || !CheckITSstatus(status2)) + SETBIT(survivedCut, kPixelITS); + break; default: break; } + if(fDebugLevel > 0)printf("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"); } if(fRequirements == survivedCut){ // @@ -321,3 +348,25 @@ void AliHFEextraCuts::AddQAHistograms(TList *qaList){ } qaList->AddLast(fQAlist); } + +//______________________________________________________ +void AliHFEextraCuts::PrintBitMap(Int_t bitmap){ + for(Int_t ibit = 32; ibit--; ) + printf("%d", bitmap & BIT(ibit) ? 1 : 0); + printf("\n"); +} + +//______________________________________________________ +Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus){ + // + // Check whether ITS area is dead + // + Bool_t status; + switch(itsStatus){ + case 2: status = kFALSE; break; + case 3: status = kFALSE; break; + case 7: status = kFALSE; break; + default: status = kTRUE; + } + return status; +} diff --git a/PWG3/hfe/AliHFEextraCuts.h b/PWG3/hfe/AliHFEextraCuts.h index 8a22fff8003..6aadae9cc4b 100644 --- a/PWG3/hfe/AliHFEextraCuts.h +++ b/PWG3/hfe/AliHFEextraCuts.h @@ -48,14 +48,19 @@ class AliHFEextraCuts : public AliCFCutBase{ inline void SetMinImpactParamZ(Double_t impactParam); inline void SetMaxImpactParamZ(Double_t impactParam); inline void SetMinTrackletsTRD(Int_t minTracklets); + + void SetDebugLevel(Int_t level) { fDebugLevel = level; }; + Int_t GetDebugLevel() const { return fDebugLevel; }; protected: virtual void AddQAHistograms(TList *qaList); Bool_t CheckESDCuts(AliESDtrack *track); Bool_t CheckMCCuts(AliMCParticle */*track*/); + Bool_t CheckITSstatus(Int_t itsStatus); void FillQAhistosESD(AliESDtrack *track, UInt_t when); // void FillQAhistosMC(AliMCParticle *track, UInt_t when); void FillCutCorrelation(ULong64_t survivedCut); + void PrintBitMap(Int_t bitmap); private: typedef enum{ @@ -84,6 +89,8 @@ class AliHFEextraCuts : public AliCFCutBase{ TList *fQAlist; //! Directory for QA histograms + Int_t fDebugLevel; // Debug Level + ClassDef(AliHFEextraCuts, 1) // Additional cuts implemented by the ALICE HFE group }; diff --git a/PWG3/hfe/AliHFEmcQA.cxx b/PWG3/hfe/AliHFEmcQA.cxx index bde5e4e4ba3..e359c93a254 100644 --- a/PWG3/hfe/AliHFEmcQA.cxx +++ b/PWG3/hfe/AliHFEmcQA.cxx @@ -15,6 +15,14 @@ /************************************************************************** * * * QA class of Heavy Flavor quark and fragmeted/decayed particles * + * -Check kinematics of Heavy Quarks/hadrons, and decayed leptons * + * pT, rapidity * + * decay lepton kinematics w/wo acceptance * + * heavy hadron decay length, electron pT fraction carried from decay * + * -Check yield of Heavy Quarks/hadrons * + * Number of produced heavy quark * + * Number of produced hadron of given pdg code * + * * * * * Authors: * * MinJung Kweon * @@ -45,26 +53,23 @@ const Int_t AliHFEmcQA::fgkGluon=21; const Int_t AliHFEmcQA::fgkMaxGener=10; -const Int_t AliHFEmcQA::fgkMaxIter=1000; -const Int_t AliHFEmcQA::fgkqType = 6; // number of species waiting for QA done +const Int_t AliHFEmcQA::fgkMaxIter=100; +const Int_t AliHFEmcQA::fgkqType = 7; // number of species waiting for QA done ClassImp(AliHFEmcQA) //_______________________________________________________________________________________________ AliHFEmcQA::AliHFEmcQA() : - fVerbos(kFALSE) - ,fStack(0x0) + fStack(0x0) ,fNparents(0) { // Default constructor - if (fVerbos) AliInfo("***** Warning! fVerbos is set to TRUE! ******"); } //_______________________________________________________________________________________________ AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p): TObject(p) - ,fVerbos(p.fVerbos) ,fStack(0x0) ,fNparents(p.fNparents) { @@ -90,63 +95,68 @@ AliHFEmcQA::~AliHFEmcQA() } //_______________________________________________________________________________________________ -void AliHFEmcQA::PostAnalyze() +const void AliHFEmcQA::PostAnalyze() { } //__________________________________________ -void AliHFEmcQA::CreateHistograms(const Int_t kquark, TString hnopt) +void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) { // create histograms - if (kquark != fkCharm && kquark != fkBeauty) { + if (kquark != kCharm && kquark != kBeauty) { AliDebug(1, "This task is only for heavy quark QA, return\n"); return; } - Int_t iq = kquark - fkCharm; + Int_t iq = kquark - kCharm; TString kqTypeLabel[fgkqType]; - if (kquark == fkCharm){ - kqTypeLabel[fkQuark]="c"; - kqTypeLabel[fkantiQuark]="cbar"; - kqTypeLabel[fkElectron]="ce"; - kqTypeLabel[fkElectron2nd]="nulle"; - kqTypeLabel[fkeHadron]="ceHadron"; - kqTypeLabel[fkDeHadron]="nullHadron"; - } else if (kquark == fkBeauty){ - kqTypeLabel[fkQuark]="b"; - kqTypeLabel[fkantiQuark]="bbar"; - kqTypeLabel[fkElectron]="be"; - kqTypeLabel[fkElectron2nd]="bce"; - kqTypeLabel[fkeHadron]="beHadron"; - kqTypeLabel[fkDeHadron]="bDeHadron"; + if (kquark == kCharm){ + kqTypeLabel[kQuark]="c"; + kqTypeLabel[kantiQuark]="cbar"; + kqTypeLabel[kHadron]="cHadron"; + kqTypeLabel[keHadron]="ceHadron"; + kqTypeLabel[kDeHadron]="nullHadron"; + kqTypeLabel[kElectron]="ce"; + kqTypeLabel[kElectron2nd]="nulle"; + } else if (kquark == kBeauty){ + kqTypeLabel[kQuark]="b"; + kqTypeLabel[kantiQuark]="bbar"; + kqTypeLabel[kHadron]="bHadron"; + kqTypeLabel[keHadron]="beHadron"; + kqTypeLabel[kDeHadron]="bDeHadron"; + kqTypeLabel[kElectron]="be"; + kqTypeLabel[kElectron2nd]="bce"; } TString hname; 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]; - fHist[iq][iqType].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5); + fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5); hname = hnopt+"Pt_"+kqTypeLabel[iqType]; - fHist[iq][iqType].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",150,0,30); + fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",150,0,30); hname = hnopt+"Y_"+kqTypeLabel[iqType]; - fHist[iq][iqType].fY = new TH1F(hname,hname,150,-7.5,7.5); + fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5); hname = hnopt+"Eta_"+kqTypeLabel[iqType]; - fHist[iq][iqType].fEta = new TH1F(hname,hname,150,-7.5,7.5); + fHist[iq][iqType][icut].fEta = new TH1F(hname,hname,150,-7.5,7.5); } - hname = hnopt+"Nq_"+kqTypeLabel[fkQuark]; - fHistComm[iq].fNq = new TH1F(hname,hname,10,-0.5,9.5); - hname = hnopt+"ProcessID_"+kqTypeLabel[fkQuark]; - fHistComm[iq].fProcessID = new TH1F(hname,hname,21,-10.5,10.5); - hname = hnopt+"ePtRatio_"+kqTypeLabel[fkQuark]; - fHistComm[iq].fePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1); - hname = hnopt+"DePtRatio_"+kqTypeLabel[fkQuark]; - fHistComm[iq].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1); - hname = hnopt+"eDistance_"+kqTypeLabel[fkQuark]; - fHistComm[iq].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2); - hname = hnopt+"DeDistance_"+kqTypeLabel[fkQuark]; - fHistComm[iq].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2); + if (icut == 0){ + hname = hnopt+"Nq_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5); + hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].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",100,0,30,100,0,1); + hname = hnopt+"DePtRatio_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fDePtRatio = new TH2F(hname,hname+";p_{T} (GeV/c);momentum fraction",100,0,30,100,0,1); + hname = hnopt+"eDistance_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].feDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2); + hname = hnopt+"DeDistance_"+kqTypeLabel[kQuark]; + fHistComm[iq][icut].fDeDistance= new TH2F(hname,hname+";p_{T} (GeV/c);distance (cm)",100,0,30,200,0,2); } @@ -184,11 +194,11 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark) { // get heavy quark kinematics - if (kquark != fkCharm && kquark != fkBeauty) { + if (kquark != kCharm && kquark != kBeauty) { AliDebug(1, "This task is only for heavy quark QA, return\n"); return; } - Int_t iq = kquark - fkCharm; + Int_t iq = kquark - kCharm; if (iTrack < 0) { @@ -226,7 +236,7 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark) if ( abs(partMother->GetPdgCode()) != kquark ){ // search fragmented partons in the same string - Bool_t IsSameString = kTRUE; + Bool_t isSameString = kTRUE; for (Int_t i=1; i-1) { partMother = fStack->Particle(iLabel); } @@ -235,8 +245,8 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark) return; } if ( abs(partMother->GetPdgCode()) == kquark ) break; - if ( partMother->GetStatusCode() != 12 ) IsSameString = kFALSE; - if (!IsSameString) return; + if ( partMother->GetStatusCode() != 12 ) isSameString = kFALSE; + if (!isSameString) return; } } AliDebug(1, "Can not find heavy parton of this heavy hadron in the string, return\n"); @@ -248,15 +258,15 @@ void AliHFEmcQA::GetQuarkKine(Int_t iTrack, const Int_t kquark) // fill kinematics for heavy parton if (partMother->GetPdgCode() > 0) { // quark - fHist[iq][fkQuark].fPdgCode->Fill(partMother->GetPdgCode()); - fHist[iq][fkQuark].fPt->Fill(partMother->Pt()); - fHist[iq][fkQuark].fY->Fill(GetRapidity(partMother)); - fHist[iq][fkQuark].fEta->Fill(partMother->Eta()); + fHist[iq][kQuark][0].fPdgCode->Fill(partMother->GetPdgCode()); + fHist[iq][kQuark][0].fPt->Fill(partMother->Pt()); + fHist[iq][kQuark][0].fY->Fill(GetRapidity(partMother)); + fHist[iq][kQuark][0].fEta->Fill(partMother->Eta()); } else{ // antiquark - fHist[iq][fkantiQuark].fPdgCode->Fill(partMother->GetPdgCode()); - fHist[iq][fkantiQuark].fPt->Fill(partMother->Pt()); - fHist[iq][fkantiQuark].fY->Fill(GetRapidity(partMother)); - fHist[iq][fkantiQuark].fEta->Fill(partMother->Eta()); + fHist[iq][kantiQuark][0].fPdgCode->Fill(partMother->GetPdgCode()); + fHist[iq][kantiQuark][0].fPt->Fill(partMother->Pt()); + fHist[iq][kantiQuark][0].fY->Fill(GetRapidity(partMother)); + fHist[iq][kantiQuark][0].fEta->Fill(partMother->Eta()); } } // end of heavy parton slection loop @@ -270,16 +280,16 @@ void AliHFEmcQA::EndOfEventAna(const Int_t kquark) { // end of event analysis - if (kquark != fkCharm && kquark != fkBeauty) { + if (kquark != kCharm && kquark != kBeauty) { AliDebug(1, "This task is only for heavy quark QA, return\n"); return; } - Int_t iq = kquark - fkCharm; + Int_t iq = kquark - kCharm; // # of heavy quark per event AliDebug(1,Form("Number of heavy quark in this event = %d \n",fIsHeavy[iq])); - fHistComm[iq].fNq->Fill(fIsHeavy[iq]); + fHistComm[iq][0].fNq->Fill(fIsHeavy[iq]); Int_t motherID[fgkMaxGener]; Int_t motherType[fgkMaxGener]; @@ -318,7 +328,7 @@ void AliHFEmcQA::EndOfEventAna(const Int_t kquark) if (IsFromInitialShower(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break; if (IsFromFinalParton(ancestorLabel[ig],motherID[i],motherType[i],motherLabel[i])) break; // if it is not the above case, something is strange - reportStrangeness(motherID[i],motherType[i],motherLabel[i]); + ReportStrangeness(motherID[i],motherType[i],motherLabel[i]); break; } if (ancestorLabel[ig] == -1){ // from hard scattering @@ -342,47 +352,117 @@ void AliHFEmcQA::EndOfEventAna(const Int_t kquark) Int_t id2 = ipair*2 + 1; if (motherType[id1] == 2 && motherType[id2] == 2){ - if (motherLabel[id1] == motherLabel[id2]) processID = fkGluonSplitting; // gluon spliting + if (motherLabel[id1] == motherLabel[id2]) processID = kGluonSplitting; // gluon spliting else processID = -9; } else if (motherType[id1] == -1 && motherType[id2] == -1) { if (motherLabel[id1] == -1 && motherLabel[id2] == -1) { - if (motherID[id1] == fgkGluon) processID = fkPairCreationFromg; // gluon fusion - else processID = fkPairCreationFromq; // q-qbar pair creation + if (motherID[id1] == fgkGluon) processID = kPairCreationFromg; // gluon fusion + else processID = kPairCreationFromq; // q-qbar pair creation } else processID = -8; } else if (motherType[id1] == -1 || motherType[id2] == -1) { if ((motherLabel[id1] == -1 || motherLabel[id2] == -1) && (motherLabel[id1]*motherLabel[id2] == -2 || motherLabel[id1]*motherLabel[id2] == -3)) { - if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = fkFlavourExitation; // flavour exitation - else processID = fkLightQuarkShower; + if(motherID[id1]*motherID[id2] == kquark*fgkGluon) processID = kFlavourExitation; // flavour exitation + else processID = kLightQuarkShower; } else processID = -7; } else if (motherType[id1] == -2 || motherType[id2] == -2) { - if (motherLabel[id1] == motherLabel[id2]) processID = fkInitialPartonShower; // initial parton shower + if (motherLabel[id1] == motherLabel[id2]) processID = kInitialPartonShower; // initial parton shower else processID = -6; } else processID = -5; if (nheavypair >1) AliDebug(1,Form("Multi pair found : process ID = %d\n",processID)); - else fHistComm[iq].fProcessID->Fill(processID); + else fHistComm[iq][0].fProcessID->Fill(processID); AliDebug(1,Form("Process ID = %d\n",processID)); } // end of # heavy quark pair loop } //__________________________________________ -void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed, Bool_t isbarrel) +void AliHFEmcQA::GetHadronKine(Int_t iTrack, const Int_t kquark) +{ + // decay electron kinematics + + if (kquark != kCharm && kquark != kBeauty) { + AliDebug(1, "This task is only for heavy quark QA, return\n"); + return; + } + Int_t iq = kquark - kCharm; + + if (iTrack < 0) { + AliDebug(1, "Stack label is negative, return\n"); + return; + } + + TParticle* mcpart = fStack->Particle(iTrack); + + Int_t iLabel = mcpart->GetFirstMother(); + if (iLabel<0){ + AliDebug(1, "Stack label is negative, return\n"); + return; + } + + TParticle *partCopy = mcpart; + Int_t pdgcode = mcpart->GetPdgCode(); + Int_t pdgcodeCopy = pdgcode; + + // if the mother is charmed hadron + Bool_t IsDirectCharm = kFALSE; + if ( int(abs(pdgcode)/100.) == kCharm || int(abs(pdgcode)/1000.) == kCharm ) { + + // iterate until you find B hadron as a mother or become top ancester + for (Int_t i=1; iGetFirstMother(); + if (jLabel == -1){ + IsDirectCharm = kTRUE; + break; // if there is no ancester + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return; + } + // if there is an ancester + TParticle* mother = fStack->Particle(jLabel); + Int_t motherPDG = mother->GetPdgCode(); + + for (Int_t j=0; jFill(pdgcodeCopy); + fHist[iq][kHadron][0].fPt->Fill(partCopy->Pt()); + fHist[iq][kHadron][0].fY->Fill(GetRapidity(partCopy)); + fHist[iq][kHadron][0].fEta->Fill(partCopy->Eta()); + + } + } + } // end of if +} + +//__________________________________________ +void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed, Int_t icut, Bool_t isbarrel) { // decay electron kinematics - if (kquark != fkCharm && kquark != fkBeauty) { + if (kquark != kCharm && kquark != kBeauty) { AliDebug(1, "This task is only for heavy quark QA, return\n"); return; } - Int_t iq = kquark - fkCharm; + Int_t iq = kquark - kCharm; if (iTrack < 0) { AliDebug(1, "Stack label is negative, return\n"); @@ -401,25 +481,24 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed } TParticle *partMother = fStack->Particle(iLabel); + TParticle *partMotherCopy = partMother; Int_t maPdgcode = partMother->GetPdgCode(); - TParticle *partMotherCopy = fStack->Particle(iLabel); - Int_t maPdgcodeCopy = partMotherCopy->GetPdgCode(); + Int_t maPdgcodeCopy = maPdgcode; // get electron production vertex TLorentzVector ePoint; mcpart->ProductionVertex(ePoint); - // if the mother is charmed hadron - Bool_t IsMotherDirectCharm = kFALSE; - if ( int(abs(maPdgcode)/100.) == fkCharm || int(abs(maPdgcode)/1000.) == fkCharm ) { + Bool_t isMotherDirectCharm = kFALSE; + if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { // iterate until you find B hadron as a mother or become top ancester for (Int_t i=1; iGetFirstMother(); if (jLabel == -1){ - IsMotherDirectCharm = kTRUE; + isMotherDirectCharm = kTRUE; break; // if there is no ancester } if (jLabel < 0){ // safety protection @@ -434,27 +513,27 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed for (Int_t j=0; jFill(mcpart->GetPdgCode()); - fHist[iq][fkElectron2nd].fPt->Fill(mcpart->Pt()); - fHist[iq][fkElectron2nd].fY->Fill(GetRapidity(mcpart)); - fHist[iq][fkElectron2nd].fEta->Fill(mcpart->Eta()); + fHist[iq][kElectron2nd][icut].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart)); + fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta()); // fill mother hadron kinematics - fHist[iq][fkDeHadron].fPdgCode->Fill(grandMaPDG); - fHist[iq][fkDeHadron].fPt->Fill(grandMa->Pt()); - fHist[iq][fkDeHadron].fY->Fill(GetRapidity(grandMa)); - fHist[iq][fkDeHadron].fEta->Fill(grandMa->Eta()); + fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa)); + fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta()); // ratio between pT of electron and pT of mother B hadron - if(grandMa->Pt()) fHistComm[iq].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt()); + if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt()); // distance between electron production point and mother hadron production point TLorentzVector debPoint; grandMa->ProductionVertex(debPoint); TLorentzVector dedistance = ePoint - debPoint; - fHistComm[iq].fDeDistance->Fill(grandMa->Pt(),dedistance.P()); + fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),dedistance.P()); return; } } @@ -462,30 +541,30 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed partMother = grandMa; } // end of iteration } // end of if - if((IsMotherDirectCharm == kTRUE && kquark == fkCharm) || kquark == fkBeauty) { + if((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) { for (Int_t i=0; iFill(mcpart->GetPdgCode()); - fHist[iq][fkElectron].fPt->Fill(mcpart->Pt()); - fHist[iq][fkElectron].fY->Fill(GetRapidity(mcpart)); - fHist[iq][fkElectron].fEta->Fill(mcpart->Eta()); + fHist[iq][kElectron][icut].fPdgCode->Fill(mcpart->GetPdgCode()); + fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart)); + fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta()); // fill mother hadron kinematics - fHist[iq][fkeHadron].fPdgCode->Fill(maPdgcodeCopy); - fHist[iq][fkeHadron].fPt->Fill(partMotherCopy->Pt()); - fHist[iq][fkeHadron].fY->Fill(GetRapidity(partMotherCopy)); - fHist[iq][fkeHadron].fEta->Fill(partMotherCopy->Eta()); + fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy)); + fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta()); // ratio between pT of electron and pT of mother B hadron - if(partMotherCopy->Pt()) fHistComm[iq].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt()); + if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt()); // distance between electron production point and mother hadron production point TLorentzVector ebPoint; partMotherCopy->ProductionVertex(ebPoint); TLorentzVector edistance = ePoint - ebPoint; - fHistComm[iq].feDistance->Fill(partMotherCopy->Pt(),edistance.P()); + fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),edistance.P()); } } } // end of if @@ -495,6 +574,8 @@ void AliHFEmcQA::GetDecayedKine(Int_t iTrack, const Int_t kquark, Int_t kdecayed //__________________________________________ void AliHFEmcQA::IdentifyMother(Int_t mother_label, Int_t &mother_pdg, Int_t &grandmother_label) { + // find mother pdg code and label + if (mother_label < 0) { AliDebug(1, "Stack label is negative, return\n"); return; @@ -506,9 +587,10 @@ void AliHFEmcQA::IdentifyMother(Int_t mother_label, Int_t &mother_pdg, Int_t &gr } //__________________________________________ -// mothertype -1 means this heavy quark coming from hard vertex void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel) { + // mothertype -1 means this heavy quark coming from hard vertex + TParticle *afterinitialrad1 = fStack->Particle(4); TParticle *afterinitialrad2 = fStack->Particle(5); @@ -538,9 +620,10 @@ void AliHFEmcQA::HardScattering(const Int_t kquark, Int_t &motherID, Int_t &moth } //__________________________________________ -// mothertype -2 means this heavy quark coming from initial state Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel) { + // mothertype -2 means this heavy quark coming from initial state + if (inputmotherlabel==2 || inputmotherlabel==3){ // mother exist before initial state radiation TParticle *heavysMother = fStack->Particle(inputmotherlabel); motherID = heavysMother->GetPdgCode(); @@ -555,9 +638,10 @@ Bool_t AliHFEmcQA::IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, } //__________________________________________ -// mothertype 2 means this heavy quark coming from final state Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel) { + // mothertype 2 means this heavy quark coming from final state + if (inputmotherlabel > 5){ // mother exist after hard scattering TParticle *heavysMother = fStack->Particle(inputmotherlabel); motherID = heavysMother->GetPdgCode(); @@ -571,8 +655,10 @@ Bool_t AliHFEmcQA::IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, In } //__________________________________________ -void AliHFEmcQA::reportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel) +void AliHFEmcQA::ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel) { + // mark strange behavior + mothertype = -888; motherlabel = -888; motherID = -888; @@ -580,8 +666,10 @@ void AliHFEmcQA::reportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &mo } //__________________________________________ -Float_t AliHFEmcQA::GetRapidity(TParticle *part) +const Float_t AliHFEmcQA::GetRapidity(TParticle *part) { + // return rapidity + Float_t rapidity; if(part->Energy() - part->Pz() == 0 || part->Energy() + part->Pz() == 0) rapidity=-999; else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); diff --git a/PWG3/hfe/AliHFEmcQA.h b/PWG3/hfe/AliHFEmcQA.h index df12d8b8d8f..de8a2437911 100644 --- a/PWG3/hfe/AliHFEmcQA.h +++ b/PWG3/hfe/AliHFEmcQA.h @@ -12,15 +12,21 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - /************************************************************************** * * * QA class of Heavy Flavor quark and fragmeted/decayed particles * + * -Check kinematics of Heavy Quarks/hadrons, and decayed leptons * + * pT, rapidity * + * decay lepton kinematics w/wo acceptance * + * heavy hadron decay length, electron pT fraction carried from decay * + * -Check yield of Heavy Quarks/hadrons * + * Number of produced heavy quark * + * Number of produced hadron of given pdg code * * * **************************************************************************/ -#ifndef _ALIHFEMCQA_H_ -#define _ALIHFEMCQA_H_ +#ifndef ALIHFEMCQA_H +#define ALIHFEMCQA_H #ifndef ROOT_TObject #include @@ -36,49 +42,45 @@ class AliStack; class AliHFEmcQA: public TObject { public: - enum heavyType {fkCharm=4, fkBeauty=5}; - enum qType {fkQuark, fkantiQuark, fkElectron, fkElectron2nd, fkeHadron, fkDeHadron}; + enum heavyType {kCharm=4, kBeauty=5, kElectronPDG=11}; + enum qType {kQuark, kantiQuark, kHadron, keHadron, kDeHadron, kElectron, kElectron2nd}; AliHFEmcQA(); AliHFEmcQA(const AliHFEmcQA &p); // copy constructor AliHFEmcQA &operator=(const AliHFEmcQA &); // assignment operator virtual ~AliHFEmcQA(); - void PostAnalyze(); - void SetVerbosity(Bool_t verb=kTRUE) {fVerbos=verb;} - - protected: - void SetStack(AliStack* stack){fStack=stack;} // set stack pointer - void CreateHistograms(const Int_t kquark, TString hnopt=""); // create histograms for mc qa analysis + const void PostAnalyze(); + void CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt=""); // create histograms for mc qa analysis + void SetStack(AliStack* const stack){fStack=stack;} // set stack pointer void Init(); + void GetQuarkKine(Int_t iTrack, const Int_t kquark); // get heavy quark kinematics distribution - void GetDecayedKine(Int_t iTrack, const Int_t kquark, const Int_t kdecayed, Bool_t isbarrel=kFALSE); // get decay electron kinematics distribution + void GetHadronKine(Int_t iTrack, const Int_t kquark); // get heavy hadron kinematics distribution + void GetDecayedKine(Int_t iTrack, const Int_t kquark, const Int_t kdecayed, Int_t icut, Bool_t isbarrel=kFALSE); // get decay electron kinematics distribution void EndOfEventAna(const Int_t kquark); // run analysis which should be done at the end of the event loop + + protected: void IdentifyMother(Int_t mother_label, Int_t &mother_pdg, Int_t &grandmother_label); // void HardScattering(const Int_t kquark, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // check if the quark is produced from hard scattering - void reportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // report if the quark production process is unknown + void ReportStrangeness(Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // report if the quark production process is unknown Bool_t IsFromInitialShower(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // check if the quark is produced from initial parton shower Bool_t IsFromFinalParton(Int_t inputmotherlabel, Int_t &motherID, Int_t &mothertype, Int_t &motherlabel); // check if the quark is produced from final parton shower - Float_t GetRapidity(TParticle *part); // return rapidity + const Float_t GetRapidity(TParticle *part); // return rapidity - Bool_t fVerbos; // + AliStack* fStack; // stack pointer - AliStack* fStack; // - - static const Int_t fgkGluon; // - static const Int_t fgkPDGInterested; // - static const Int_t fgkMaxGener; // - static const Int_t fgkMaxIter; // - static const Int_t fgkqType; // + static const Int_t fgkGluon; // gluon pdg code + static const Int_t fgkMaxGener; // ancester level wanted to be checked + static const Int_t fgkMaxIter; // number of iteration to find out matching particle + static const Int_t fgkqType; // number of particle type to be checked enum ProcessType_t { - fkPairCreationFromq, fkPairCreationFromg, fkFlavourExitation, fkGluonSplitting, fkInitialPartonShower, fkLightQuarkShower + kPairCreationFromq, kPairCreationFromg, kFlavourExitation, kGluonSplitting, kInitialPartonShower, kLightQuarkShower }; - TString fkqTypeLabel[6]; // - struct hists{ TH1F *fPdgCode; // histogram to store particle pdg code TH1F *fPt; // histogram to store pt @@ -88,22 +90,22 @@ class AliHFEmcQA: public TObject { struct histsComm { TH1F *fNq; // histogram to store number of quark TH1F *fProcessID; // histogram to store process id - TH2F *fePtRatio; // - TH2F *fDePtRatio; // - TH2F *feDistance; // - TH2F *fDeDistance; // + TH2F *fePtRatio; // fraction of electron pT from D or B hadron + TH2F *fDePtRatio; // fraction of D electron pT from B hadron + TH2F *feDistance; // distance between electron production point to mother particle + TH2F *fDeDistance; // distance between D electron production point to mother particle }; - hists fHist[2][6]; // - histsComm fHistComm[2]; // + hists fHist[2][7][5]; // struct of histograms to store kinematics of given particles + histsComm fHistComm[2][5]; // struct of additional histograms of given particles - TParticle *fHeavyQuark[50]; // - Int_t fIsHeavy[2]; // - Int_t fNparents; // - Int_t fParentSelect[2][7]; // + TParticle *fHeavyQuark[50]; // store pointer of heavy flavour quark + Int_t fIsHeavy[2]; // count of heavy flavour + Int_t fNparents; // number of heavy hadrons to be considered + Int_t fParentSelect[2][7]; // heavy hadron species - ClassDef(AliHFEmcQA,0); // HFE Monte Carlo quality assurance + ClassDef(AliHFEmcQA,0); }; #endif diff --git a/PWG3/hfe/AliHFEpidTRD.cxx b/PWG3/hfe/AliHFEpidTRD.cxx index 89973d59689..68e9d292795 100644 --- a/PWG3/hfe/AliHFEpidTRD.cxx +++ b/PWG3/hfe/AliHFEpidTRD.cxx @@ -44,28 +44,23 @@ ClassImp(AliHFEpidTRD) //___________________________________________________________________ AliHFEpidTRD::AliHFEpidTRD(const char* name) : AliHFEpidBase(name) - , fThresholdFile("TRD.PIDthresholds.root") , fPIDMethod(kNN) - , fTRDthresholds(0x0) - , fTRDelectronEfficiencies(0x0) { // // default constructor // - fTRDthresholds = new TMap(); + memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams); } //___________________________________________________________________ AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref): AliHFEpidBase("") - , fThresholdFile("") , fPIDMethod(kLQ) - , fTRDthresholds(0x0) - , fTRDelectronEfficiencies(0x0) { // // Copy constructor // + memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams); ref.Copy(*this); } @@ -87,14 +82,8 @@ void AliHFEpidTRD::Copy(TObject &ref) const { // AliHFEpidTRD &target = dynamic_cast(ref); - target.fThresholdFile = fThresholdFile; target.fPIDMethod = fPIDMethod; - if(fTRDthresholds){ - target.fTRDthresholds = dynamic_cast(fTRDthresholds->Clone()); - } - if(fTRDelectronEfficiencies){ - target.fTRDelectronEfficiencies = new TAxis(*fTRDelectronEfficiencies); - } + memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams); AliHFEpidBase::Copy(ref); } @@ -103,11 +92,6 @@ AliHFEpidTRD::~AliHFEpidTRD(){ // // Destructor // - if(fTRDthresholds){ - fTRDthresholds->Delete(); - delete fTRDthresholds; - } - if(fTRDelectronEfficiencies) delete fTRDelectronEfficiencies; } //______________________________________________________ @@ -116,44 +100,7 @@ Bool_t AliHFEpidTRD::InitializePID(){ // InitializePID: Load TRD thresholds and create the electron efficiency axis // to navigate // - LoadTRDthresholds(); - // Fill the electron efficiencies axis for the TRD alone PID - Int_t nEffs = fTRDthresholds->GetEntries() + 1; - TIterator* thres =fTRDthresholds->MakeIterator(); - TObjString *key = 0x0; - Double_t *tmp = new Double_t[nEffs], *titer = tmp; - while((key = dynamic_cast((*thres)()))){ - (*titer++) = static_cast(key->String().Atoi())/100.; - } - delete thres; - *titer = 1.; - // Sort the electron efficiencies and put them into the TAxis for later navigation - Int_t *ind = new Int_t[nEffs], *iiter = ind; - TMath::Sort(nEffs, tmp, ind, kFALSE); - Double_t *eleffs = new Double_t[nEffs], *eiter = eleffs; - while(eiter < &eleffs[nEffs]) *(eiter++) = tmp[*(iiter++)]; - // print the content - Int_t cnt = 0; - if(GetDebugLevel() > 1){ - printf("Printing electron efficiency bins:\n"); - eiter = eleffs; - thres = fTRDthresholds->MakeIterator(); - TObject *object = 0x0; - while(eiter < &eleffs[nEffs - 1]){ - printf("eleffs[%d] = %f", cnt++, *(eiter++)); - key = dynamic_cast(thres->Next()); - object = fTRDthresholds->GetValue(key->String().Data()); - printf(", Content: %p\n", (void *)object); - } - } - delete[] tmp; delete[] ind; - fTRDelectronEfficiencies = new TAxis(nEffs - 1, eleffs); - if(GetDebugLevel() > 1){ - printf("Printing axis content:\n"); - for(Int_t ibin = fTRDelectronEfficiencies->GetFirst(); ibin <= fTRDelectronEfficiencies->GetLast(); ibin++) - printf("%d.) minimum: %f, maximum? %f\n", ibin, fTRDelectronEfficiencies->GetBinLowEdge(ibin), fTRDelectronEfficiencies->GetBinUpEdge(ibin)); - } - delete[] eleffs; + InitParameters(); return kTRUE; } @@ -169,71 +116,60 @@ Int_t AliHFEpidTRD::IsSelected(AliVParticle *track){ Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P(); if(p < 0.6) return 0; - // Get the Histograms - TH1 *threshist = GetTRDthresholds(0.91); - Int_t bin = 0; - if(p > threshist->GetXaxis()->GetXmax()) - bin = threshist->GetXaxis()->GetLast(); - else if(p < threshist->GetXaxis()->GetXmin()) - bin = threshist->GetXaxis()->GetFirst(); - else - bin = threshist->GetXaxis()->FindBin(p); - Double_t pidProbs[AliPID::kSPECIES]; esdTrack->GetTRDpid(pidProbs); - if(pidProbs[AliPID::kElectron] > threshist->GetBinContent(bin)) return 11; + if(pidProbs[AliPID::kElectron] > GetTRDthresholds(0.91, p)) return 11; return 0; } //___________________________________________________________________ -void AliHFEpidTRD::LoadTRDthresholds(){ - // - // Load TRD threshold histograms from File - // - TFile *mythresholds = TFile::Open(fThresholdFile); - TKey *object = 0x0; - TString electron_eff; - TObjArray *histos = 0x0; - Float_t eff; - TH1F *refhist = 0x0; - TIterator *keyIterator = mythresholds->GetListOfKeys()->MakeIterator(); - TString histnames[2] = {"fHistThreshLQ", "fHistThreshNN"}; - gROOT->cd(); - while((object = dynamic_cast((*keyIterator)()))){ - // Get the electron efficiency bin this histogram was taken with - electron_eff = object->GetName(); - electron_eff = electron_eff.Remove(0,3); - eff = static_cast(electron_eff.Atoi())/100.; - - // Get the threshold according to the selected - histos = dynamic_cast(object->ReadObj()); - refhist = dynamic_cast(histos->FindObject(histnames[fPIDMethod].Data())); - SetTRDthresholds(refhist, eff); - histos->Delete(); - delete histos; - } - delete keyIterator; - mythresholds->Close(); - delete mythresholds; +Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){ + // + // Return momentum dependent and electron efficiency dependent TRD thresholds + // + Double_t params[4]; + GetParameters(electronEff, params); + return 1. - params[0] * p - params[1] * TMath::Exp(-params[2] * p) - params[3]; } //___________________________________________________________________ -void AliHFEpidTRD::SetTRDthresholds(TH1F *thresholds, Float_t electronEff){ - // - // Set the threshold histogram for the TRD pid - // - fTRDthresholds->Add(new TObjString(Form("%d", TMath::Nint(electronEff * 100.))), new TH1F(*thresholds)); +void AliHFEpidTRD::InitParameters(){ + // + // Fill the Parameters into an array + // + + // Parameters for 6 Layers + fThreshParams[0] = -0.001839; // 0.7 electron eff + fThreshParams[1] = 0.000276; + fThreshParams[2] = 0.044902; + fThreshParams[3] = 1.726751; + fThreshParams[4] = -0.002405; // 0.75 electron eff + fThreshParams[5] = 0.000372; + fThreshParams[6] = 0.061775; + fThreshParams[7] = 1.739371; + fThreshParams[8] = -0.003178; // 0.8 electron eff + fThreshParams[9] = 0.000521; + fThreshParams[10] = 0.087585; + fThreshParams[11] = 1.749154; + fThreshParams[12] = -0.004058; // 0.85 electron eff + fThreshParams[13] = 0.000748; + fThreshParams[14] = 0.129583; + fThreshParams[15] = 1.782323; + fThreshParams[16] = -0.004967; // 0.9 electron eff + fThreshParams[17] = 0.001216; + fThreshParams[18] = 0.210128; + fThreshParams[19] = 1.807665; + fThreshParams[20] = -0.000996; // 0.95 electron eff + fThreshParams[21] = 0.002627; + fThreshParams[22] = 0.409099; + fThreshParams[23] = 1.787076; } //___________________________________________________________________ -TH1F *AliHFEpidTRD::GetTRDthresholds(Float_t electronEff){ - Int_t bin = 0; - if(electronEff < fTRDelectronEfficiencies->GetXmin()) bin = fTRDelectronEfficiencies->GetFirst(); - else if(electronEff > fTRDelectronEfficiencies->GetXmax()) bin = fTRDelectronEfficiencies->GetLast(); - else bin = fTRDelectronEfficiencies->FindBin(electronEff); - TObjString keyname = Form("%d", TMath::Nint(fTRDelectronEfficiencies->GetBinLowEdge(bin)* 100.)); -/* printf("Key: %s\n", keyname.String().Data());*/ - TH1F *thresholds = dynamic_cast((dynamic_cast(fTRDthresholds->FindObject(&keyname)))->Value()); -/* printf("thresholds: %p\n", thresholds);*/ - return thresholds; +void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){ + // + // return parameter set for the given efficiency bin + // + Int_t effbin = static_cast((electronEff - 0.7)/0.3 * 6.); + memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4); } diff --git a/PWG3/hfe/AliHFEpidTRD.h b/PWG3/hfe/AliHFEpidTRD.h index e9757d46e12..b2e4ee4ae41 100644 --- a/PWG3/hfe/AliHFEpidTRD.h +++ b/PWG3/hfe/AliHFEpidTRD.h @@ -19,10 +19,6 @@ #include "AliHFEpidBase.h" #endif -class TAxis; -class TH1F; -class TMap; -class TString; class AliVParticle; class AliHFEpidTRD : public AliHFEpidBase{ @@ -31,6 +27,9 @@ class AliHFEpidTRD : public AliHFEpidBase{ kLQ = 0, kNN = 1 } PIDMethodTRD_t; + enum{ + kThreshParams = 24 + }; AliHFEpidTRD(const Char_t *name); AliHFEpidTRD(const AliHFEpidTRD &ref); AliHFEpidTRD& operator=(const AliHFEpidTRD &ref); @@ -40,22 +39,15 @@ class AliHFEpidTRD : public AliHFEpidBase{ virtual Int_t IsSelected(AliVParticle *track); virtual Bool_t HasQAhistos() const { return kFALSE; }; - void LoadTRDthresholds(); - void SetPIDMethod(PIDMethodTRD_t method) { fPIDMethod = method; }; - void SetThresholdFile(Char_t *thresholdFile) { fThresholdFile = thresholdFile; }; - protected: void Copy(TObject &ref) const; - TH1F *GetTRDthresholds(Float_t electronEff); - void SetTRDthresholds(TH1F *thresholds, Float_t electronEff); - + Double_t GetTRDthresholds(Double_t electronEff, Double_t p); + void InitParameters(); + void GetParameters(Double_t electronEff, Double_t *parameters); private: - TString fThresholdFile; // Threshold file name PIDMethodTRD_t fPIDMethod; // PID Method: 2D Likelihood or Neural Network - TMap *fTRDthresholds; //! TRD Thresholds - TAxis *fTRDelectronEfficiencies; //! Electron Efficiencies corresponding to reference histos - + Double_t fThreshParams[kThreshParams]; // Threshold parametrisation ClassDef(AliHFEpidTRD, 1) // TRD electron ID class }; diff --git a/PWG3/hfe/AliHFEpriVtx.h b/PWG3/hfe/AliHFEpriVtx.h index c01be6ea613..ef255c4bea2 100644 --- a/PWG3/hfe/AliHFEpriVtx.h +++ b/PWG3/hfe/AliHFEpriVtx.h @@ -43,11 +43,7 @@ class AliHFEpriVtx : public TObject { AliHFEpriVtx &operator=(const AliHFEpriVtx &); // assignment operator virtual ~AliHFEpriVtx(); - - - protected: void CreateHistograms(TString hnopt=""); // create histograms - void Init(); void SetEvent(AliESDEvent* ESD){fESD1=ESD;}; // set ESD pointer void SetStack(AliStack* stack){fStack=stack;} // set stack pointer @@ -87,7 +83,7 @@ class AliHFEpriVtx : public TObject { TH2F *fDiffDCAvsNt; // histogram to fill DCA difference as a function of pT - ClassDef(AliHFEpriVtx,0); // primary vertex QA for HFE + ClassDef(AliHFEpriVtx,0); }; #endif diff --git a/PWG3/hfe/AliHFEsecVtx.cxx b/PWG3/hfe/AliHFEsecVtx.cxx index e06460c60fb..85d33d88be7 100644 --- a/PWG3/hfe/AliHFEsecVtx.cxx +++ b/PWG3/hfe/AliHFEsecVtx.cxx @@ -15,6 +15,8 @@ /************************************************************************** * * * Secondary vertexing construction Class * + * Construct secondary vertex from Beauty hadron with electron and * + * hadrons, then apply selection criteria * * * * Authors: * * MinJung Kweon * @@ -52,6 +54,7 @@ AliHFEsecVtx::AliHFEsecVtx(): ,finvmass(-1) ,finvmassSigma(-1) ,fBTagged(0) + ,fBElec(0) { // Default constructor @@ -73,6 +76,7 @@ AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p): ,finvmass(p.finvmass) ,finvmassSigma(p.finvmassSigma) ,fBTagged(p.fBTagged) + ,fBElec(p.fBElec) { // Copy constructor } @@ -99,6 +103,8 @@ AliHFEsecVtx::~AliHFEsecVtx() void AliHFEsecVtx::Init() { + // set pdg code and index + fNparents = 7; fParentSelect[0][0] = 411; @@ -161,26 +167,27 @@ void AliHFEsecVtx::Init() fia[3][2][0] = 2; fia[3][2][1] = 3; - - //fid[4][3] = {0,1,2, 0,1,3, 0,2,3, 1,2,3}; - //fia[4][3][2] = {0,1, 0,2, 1,2, 0,1, 0,3, 1,3, 0,2, 0,3, 2,3, 1,2, 1,3, 2,3}; - } //__________________________________________ void AliHFEsecVtx::ResetTagVar() { + // reset tag variables + fdistTwoSecVtx = -1; fcosPhi = -1; fsignedLxy = -1; finvmass = -1; finvmassSigma = -1; fBTagged = kFALSE; + fBElec = kFALSE; } //__________________________________________ void AliHFEsecVtx::InitAnaPair() { + // initialize pair tagging variables + fPairTagged = 0; for (Int_t i=0; i<20; i++){ fpairedTrackID[i] = -1; @@ -196,16 +203,16 @@ void AliHFEsecVtx::CreateHistograms(TString hnopt) { // create histograms - fkSourceLabel[fkAll]="all"; - fkSourceLabel[fkDirectCharm]="directCharm"; - fkSourceLabel[fkDirectBeauty]="directBeauty"; - fkSourceLabel[fkBeautyCharm]="beauty2charm"; - fkSourceLabel[fkGamma]="gamma"; - fkSourceLabel[fkPi0]="pi0"; - fkSourceLabel[fkElse]="others"; - fkSourceLabel[fkBeautyGamma]="beauty22gamma"; - fkSourceLabel[fkBeautyPi0]="beauty22pi0"; - fkSourceLabel[fkBeautyElse]="beauty22others"; + fkSourceLabel[kAll]="all"; + fkSourceLabel[kDirectCharm]="directCharm"; + fkSourceLabel[kDirectBeauty]="directBeauty"; + fkSourceLabel[kBeautyCharm]="beauty2charm"; + fkSourceLabel[kGamma]="gamma"; + fkSourceLabel[kPi0]="pi0"; + fkSourceLabel[kElse]="others"; + fkSourceLabel[kBeautyGamma]="beauty22gamma"; + fkSourceLabel[kBeautyPi0]="beauty22pi0"; + fkSourceLabel[kBeautyElse]="beauty22others"; TString hname; @@ -219,12 +226,16 @@ void AliHFEsecVtx::CreateHistograms(TString hnopt) fHistPair[isource].fInvMassCut2 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10); hname=hnopt+"KFChi2_"+fkSourceLabel[isource]; fHistPair[isource].fKFChi2 = new TH1F(hname,hname,200,0,20); + hname=hnopt+"OpenAngle_"+fkSourceLabel[isource]; + fHistPair[isource].fOpenAngle = new TH1F(hname,hname,100,0,3.14); hname=hnopt+"CosOpenAngle_"+fkSourceLabel[isource]; fHistPair[isource].fCosOpenAngle = new TH1F(hname,hname,100,-1.1,1.1); hname=hnopt+"SignedLxy_"+fkSourceLabel[isource]; fHistPair[isource].fSignedLxy = new TH2F(hname,hname,1000,-5,5,120,-2,10); hname=hnopt+"KFIP_"+fkSourceLabel[isource]; fHistPair[isource].fKFIP = new TH1F(hname,hname,1000,-5,5); + hname=hnopt+"IPMax_"+fkSourceLabel[isource]; + fHistPair[isource].fIPMax= new TH1F(hname,hname,500,0,5); } @@ -237,11 +248,53 @@ void AliHFEsecVtx::CreateHistograms(TString hnopt) hname=hnopt+"pt_wrongtaggedelec"; fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,150,0,30); + hname=hnopt+"InvmassBeautyElecSecVtx"; + fHistTagged.fInvmassBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10); + hname=hnopt+"InvmassNotBeautyElecSecVtx"; + fHistTagged.fInvmassNotBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10); + + hname=hnopt+"SignedLxyBeautyElecSecVtx"; + fHistTagged.fSignedLxyBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5); + hname=hnopt+"SignedLxyNotBeautyElecSecVtx"; + fHistTagged.fSignedLxyNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5); + + hname=hnopt+"DistTwoVtxBeautyElecSecVtx"; + fHistTagged.fDistTwoVtxBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5); + hname=hnopt+"DistTwoVtxNotBeautyElecSecVtx"; + fHistTagged.fDistTwoVtxNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5); + + hname=hnopt+"CosPhiBeautyElecSecVtx"; + fHistTagged.fCosPhiBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1); + hname=hnopt+"CosPhiNotBeautyElecSecVtx"; + fHistTagged.fCosPhiNotBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1); + + hname=hnopt+"Chi2BeautyElecSecVtx"; + fHistTagged.fChi2BeautyElecSecVtx= new TH1F(hname,hname,200,0,20); + hname=hnopt+"Chi2NotBeautyElecSecVtx"; + fHistTagged.fChi2NotBeautyElecSecVtx= new TH1F(hname,hname,200,0,20); + + hname=hnopt+"InvmassBeautyElec2trkVtx"; + fHistTagged.fInvmassBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10); + hname=hnopt+"InvmassNotBeautyElec2trkVtx"; + fHistTagged.fInvmassNotBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10); + + hname=hnopt+"SignedLxyBeautyElec2trkVtx"; + fHistTagged.fSignedLxyBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5); + hname=hnopt+"SignedLxyNotBeautyElec2trkVtx"; + fHistTagged.fSignedLxyNotBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5); + + hname=hnopt+"Chi2BeautyElec2trkVtx"; + fHistTagged.fChi2BeautyElec2trkVtx= new TH1F(hname,hname,200,0,20); + hname=hnopt+"Chi2NotBeautyElec2trkVtx"; + fHistTagged.fChi2NotBeautyElec2trkVtx= new TH1F(hname,hname,200,0,20); } //_______________________________________________________________________________________________ -void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourcePart, Int_t index2) +void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t index2) { + // calculate e-h pair characteristics and tag pair + + Int_t sourcePart = PairCode(track1,track2); // get KF particle input pid Int_t pdg1 = GetMCPID(track1); @@ -291,7 +344,8 @@ void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourc Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF())); // opening angle between two particles in XY plane - Double_t cosphi = TMath::Cos(kfTrack1.GetAngleXY(kfTrack2)); + Double_t phi = kfTrack1.GetAngleXY(kfTrack2); + Double_t cosphi = TMath::Cos(phi); // projection of kf vertex vector to the kf momentum direction Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy); @@ -302,32 +356,35 @@ void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourc Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx); + Float_t dcaR=-1; + Float_t dcaR1=-1, dcaR2=-1; + Float_t dcaZ1=-1, dcaZ2=-1; + track1->GetImpactParameters(dcaR1,dcaZ1); + track2->GetImpactParameters(dcaR2,dcaZ2); - - // pair cuts - if( kfchi2 >2. ) return; + if (TMath::Abs(dcaR1) >= TMath::Abs(dcaR2)) dcaR=dcaR1; + else dcaR=dcaR2; // fill histograms fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma); fHistPair[sourcePart].fKFChi2->Fill(kfchi2); + fHistPair[sourcePart].fOpenAngle->Fill(phi); fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi); fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass); fHistPair[sourcePart].fKFIP->Fill(kfip); - - if ( signedLxy > 0.05 && cosphi > 0.5) { - fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma); - } + fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR)); // pair cuts - if (signedLxy < 0.05) return; - if (cosphi < 0.0) return; + if( kfchi2 >2. ) return; + + if ( signedLxy > 0.05 && cosphi > 0.5) fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma); // pair tagging if it passed the above cuts - fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma); + if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma); // pair tagging condition - if ( signedLxy > 0.0 && cosphi > 0) { + if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut //if ( signedLxy > 0.06 && cosphi > 0) { fpairedTrackID[fPairTagged] = index2; fpairedChi2[fPairTagged] = kfchi2; @@ -341,9 +398,20 @@ void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourc //_______________________________________________________________________________________________ void AliHFEsecVtx::RunSECVTX(AliESDtrack *track) { + // run secondary vertexing algorithm and do tagging ResetTagVar(); + Int_t imclabel = TMath::Abs(track->GetLabel()); + if(imclabel<0) return; + TParticle* mcpart = fStack->Particle(imclabel); + Int_t esource = GetElectronSource(imclabel); + if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){ + fHistTagged.fPtBeautyElec->Fill(mcpart->Pt()); + fBElec = kTRUE; + } + + if (fPairTagged >= 4) { FindSECVTXCandid4Tracks(track); } @@ -357,19 +425,10 @@ void AliHFEsecVtx::RunSECVTX(AliESDtrack *track) ApplyPairTagCut(); } - - Int_t imclabel = TMath::Abs(track->GetLabel()); - if(imclabel<0) return; - TParticle* mcpart = fStack->Particle(imclabel); - Int_t esource = GetElectronSource(imclabel); - if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){ - fHistTagged.fPtBeautyElec->Fill(mcpart->Pt()); - } - if (fBTagged) { fHistTagged.fPtTaggedElec->Fill(mcpart->Pt()); - if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){ + if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){ fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt()); } else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt()); @@ -380,22 +439,34 @@ void AliHFEsecVtx::RunSECVTX(AliESDtrack *track) //_______________________________________________________________________________________________ void AliHFEsecVtx::ApplyPairTagCut() { + // apply tagging cut for e-h pair + if (fBElec){ + fHistTagged.fInvmassBeautyElec2trkVtx->Fill(fpairedInvMass[0]); + fHistTagged.fSignedLxyBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]); + fHistTagged.fChi2BeautyElec2trkVtx->Fill(fpairedChi2[0]); + } + else { + fHistTagged.fInvmassNotBeautyElec2trkVtx->Fill(fpairedInvMass[0]); + fHistTagged.fSignedLxyNotBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]); + fHistTagged.fChi2NotBeautyElec2trkVtx->Fill(fpairedChi2[0]); + } + if (fpairedChi2[0] > 2.0) return; if (fpairedInvMass[0] < 1.5) return; - if (fpairedSignedLxy[0] < 0.5) return; + if (fpairedSignedLxy[0] < 0.05) return; fBTagged = kTRUE; - } //_______________________________________________________________________________________________ Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id) { + // apply tagging cut for e-h pair of given indexed hadron if (fpairedChi2[id] > 2.0) return kFALSE; if (fpairedInvMass[id] < 1.5) return kFALSE; - if (fpairedSignedLxy[id] < 0.5) return kFALSE; + if (fpairedSignedLxy[id] < 0.05) return kFALSE; fBTagged = kTRUE; return kTRUE; @@ -405,10 +476,11 @@ Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id) //_______________________________________________________________________________________________ Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2) { + // apply tagging cut for secondary vertex if (chi2 > 2.0) return kFALSE; if (finvmass < 1.5) return kFALSE; - if (fsignedLxy < 0.5) return kFALSE; + if (fsignedLxy < 0.05) return kFALSE; if (fcosPhi < 0.90) return kFALSE; if (fdistTwoSecVtx > 0.1) return kFALSE; @@ -417,12 +489,12 @@ Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2) } //_______________________________________________________________________________________________ -//void AliHFEsecVtx::ApplyTagCut(Double_t chi2, AliESDtrack *track, AliESDtrack *htrack1, AliESDtrack *htrack2) -void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2) +void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2) { + // apply tagging cut for e-h pair of given indexed hadron if(!ApplyTagCut(chi2)){ - if(!ApplyPairTagCut(0)) ApplyPairTagCut(1); + if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2); } } @@ -430,6 +502,7 @@ void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2) //_______________________________________________________________________________________________ void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track) { + // find secondary vertex for >= 4 e-h pairs // sort pair in increasing order (kFALSE-increasing order) Int_t index[20]; @@ -466,14 +539,21 @@ void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track) // calculate secondary vertex quality variables with 1 electron and 2 hadrons CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]); - ApplyVtxTagCut(sevchi2[indexB[0]]); - //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]); + if (fBElec){ + fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]); + } + else { + fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]); + } + + ApplyVtxTagCut(sevchi2[indexB[0]],index[fia[indexA[0]][indexB[0]][0]],index[fia[indexA[0]][indexB[0]][1]]); } //_______________________________________________________________________________________________ void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track) { + // find secondary vertex for 3 e-h pairs // sort pair in increasing order (kFALSE-increasing order) Int_t indexA[1] = { 0 }; @@ -498,13 +578,20 @@ void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track) // calculate secondary vertex quality variables with 1 electron and 2 hadrons CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]); - ApplyVtxTagCut(sevchi2[indexB[0]]); - //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]); + if (fBElec){ + fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]); + } + else { + fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]); + } + + ApplyVtxTagCut(sevchi2[indexB[0]],fia[indexA[0]][indexB[0]][0],fia[indexA[0]][indexB[0]][1]); } //_______________________________________________________________________________________________ void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track) { + // find secondary vertex for 2 e-h pairs Double_t sevchi2[1]; AliESDtrack *htrack[2]; @@ -518,13 +605,20 @@ void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track) // calculate secondary vertex quality variables with 1 electron and 2 hadrons CalcSECVTXProperty(track,htrack[0],htrack[1]); - ApplyVtxTagCut(sevchi2[0]); - //ApplyTagCut(sevchi2[0],track,htrack[0],htrack[1]); + if (fBElec){ + fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[0]); + } + else { + fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[0]); + } + + ApplyVtxTagCut(sevchi2[0],0,1); } //_______________________________________________________________________________________________ void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3) { + // calculate secondary vertex properties Int_t pdg1 = GetMCPID(track1); Int_t pdg2 = GetMCPID(track2); @@ -568,11 +662,26 @@ void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, // calculate invariant mass of the kf particle kfSecondary.GetMass(finvmass,finvmassSigma); + + if (fBElec){ + fHistTagged.fInvmassBeautyElecSecVtx->Fill(finvmass); + fHistTagged.fSignedLxyBeautyElecSecVtx->Fill(fsignedLxy); + fHistTagged.fDistTwoVtxBeautyElecSecVtx->Fill(fdistTwoSecVtx); + fHistTagged.fCosPhiBeautyElecSecVtx->Fill(fcosPhi); + } + else { + fHistTagged.fInvmassNotBeautyElecSecVtx->Fill(finvmass); + fHistTagged.fSignedLxyNotBeautyElecSecVtx->Fill(fsignedLxy); + fHistTagged.fDistTwoVtxNotBeautyElecSecVtx->Fill(fdistTwoSecVtx); + fHistTagged.fCosPhiNotBeautyElecSecVtx->Fill(fcosPhi); + } + } //_______________________________________________________________________________________________ Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track) { + // return mc pid Int_t label = TMath::Abs(track->GetLabel()); TParticle* mcpart = fStack->Particle(label); @@ -585,6 +694,7 @@ Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track) //_______________________________________________________________________________________________ Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4) { + // return 4 track secondary vertex chi2 Int_t pdg1 = GetMCPID(track1); Int_t pdg2 = GetMCPID(track2); @@ -605,6 +715,7 @@ Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, A //_______________________________________________________________________________________________ Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3) { + // return 3 track secondary vertex chi2 Int_t pdg1 = GetMCPID(track1); Int_t pdg2 = GetMCPID(track2); @@ -623,6 +734,7 @@ Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, A //_______________________________________________________________________________________________ Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2) { + // return 2 track secondary vertex chi2 Int_t pdg1 = GetMCPID(track1); Int_t pdg2 = GetMCPID(track2); @@ -723,44 +835,226 @@ Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2) { // - // return pair code which is predefind as: - // fkDirectCharm, fkDirectBeauty, fkBeautyCharm, fkGamma, fkPi0, fkElse, fkBeautyGamma, fkBeautyPi0, fkBeautyElse + // return pair code which is predefinded as: + // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse // Int_t pairOriginsPDG = PairOrigin(track1,track2); - Int_t sourcePart = fkElse; + Int_t sourcePart = kElse; if (pairOriginsPDG < 0) { - sourcePart = fkBeautyElse; + sourcePart = kBeautyElse; } for (Int_t i=0; i0) sourcePart = fkDirectCharm; + if (pairOriginsPDG>0) sourcePart = kDirectCharm; if (pairOriginsPDG<0) { - sourcePart = fkBeautyCharm; + sourcePart = kBeautyCharm; } } if (abs(pairOriginsPDG)==fParentSelect[1][i]) { if (pairOriginsPDG>0) { - sourcePart = fkDirectBeauty; + sourcePart = kDirectBeauty; } - if (pairOriginsPDG<0) return fkElse; + if (pairOriginsPDG<0) return kElse; } } - if (pairOriginsPDG == 22) sourcePart = fkGamma; + if (pairOriginsPDG == 22) sourcePart = kGamma; if (pairOriginsPDG == -22) { - sourcePart = fkBeautyGamma; + sourcePart = kBeautyGamma; } - if (pairOriginsPDG == 111) sourcePart = fkPi0; + if (pairOriginsPDG == 111) sourcePart = kPi0; if (pairOriginsPDG == -111) { - sourcePart = fkBeautyPi0; + sourcePart = kBeautyPi0; } return sourcePart; } +//_______________________________________________________________________________________________ +Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack) +{ + + // return decay electron's origin + + if (iTrack < 0) { + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + TParticle* mcpart = fStack->Particle(iTrack); + + if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !!!!!!!!!!!!!!!!! + + Int_t iLabel = mcpart->GetFirstMother(); + if (iLabel<0){ + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + Int_t origin = -1; + Bool_t isFinalOpenCharm = kFALSE; + + TParticle *partMother = fStack->Particle(iLabel); + Int_t maPdgcode = partMother->GetPdgCode(); + + // if the mother is charmed hadron + if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { + + for (Int_t i=0; iGetFirstMother(); + if (jLabel == -1){ + origin = kDirectCharm; + return origin; + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + // if there is an ancester + TParticle* grandMa = fStack->Particle(jLabel); + Int_t grandMaPDG = grandMa->GetPdgCode(); + + for (Int_t i=0; iGetFirstMother(); + if (jLabel == -1){ + origin = kGamma; + return origin; + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + // if there is an ancester + TParticle* grandMa = fStack->Particle(jLabel); + Int_t grandMaPDG = grandMa->GetPdgCode(); + + for (Int_t i=0; iGetFirstMother(); + if (jLabel == -1){ + origin = kPi0; + return origin; + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + // if there is an ancester + TParticle* grandMa = fStack->Particle(jLabel); + Int_t grandMaPDG = grandMa->GetPdgCode(); + + for (Int_t i=0; iGetFirstMother(); + if (jLabel == -1){ + origin = kElse; + return origin; + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + // if there is an ancester + TParticle* grandMa = fStack->Particle(jLabel); + Int_t grandMaPDG = grandMa->GetPdgCode(); + + for (Int_t i=0; iParticle(iLabel); Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code //beauty -------------------------- - if ( int(abs(maPdgcode)/100.) == fkBeauty || int(abs(maPdgcode)/1000.) == fkBeauty ) { + if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) { for (Int_t i=0; i excited beauties @@ -799,20 +1093,20 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) } // end of if //charm -------------------------- - else if ( int(abs(maPdgcode)/100.) == fkCharm || int(abs(maPdgcode)/1000.) == fkCharm ) { + else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { for (Int_t i=0; i excited charms + if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms // to prevent any possible double counting for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester Int_t jLabel = partMother->GetFirstMother(); if (jLabel == -1){ - origin = fkDirectCharm; + origin = kDirectCharm; return origin; } if (jLabel < 0){ // safety protection even though not really necessary here @@ -826,7 +1120,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) for (Int_t j=0; jGetFirstMother(); if (jLabel == -1){ - origin = fkGamma; + origin = kGamma; return origin; } if (jLabel < 0){ // safety protection @@ -858,7 +1152,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) for (Int_t j=0; jGetFirstMother(); if (jLabel == -1){ - origin = fkPi0; + origin = kPi0; return origin; } if (jLabel < 0){ // safety protection @@ -892,7 +1186,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) for (Int_t j=0; jGetFirstMother(); if (jLabel == -1){ - origin = fkElse; + origin = kElse; return origin; } if (jLabel < 0){ // safety protection @@ -927,7 +1221,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) for (Int_t j=0; jPt() < 1.0) return kFALSE; //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE; diff --git a/PWG3/hfe/AliHFEsecVtx.h b/PWG3/hfe/AliHFEsecVtx.h index 2edeb143ba1..484f3114a9f 100644 --- a/PWG3/hfe/AliHFEsecVtx.h +++ b/PWG3/hfe/AliHFEsecVtx.h @@ -12,16 +12,17 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - /************************************************************************** * * * Secondary vertexing construction Class * + * Construct secondary vertex from Beauty hadron with electron and * + * hadrons, then apply selection criteria * * * **************************************************************************/ -#ifndef _ALIHFESECVTX_H_ -#define _ALIHFESECVTX_H_ +#ifndef ALIHFESECVTX_H +#define ALIHFESECVTX_H #ifndef ROOT_TObject #include @@ -43,33 +44,35 @@ class AliHFEsecVtx : public TObject { AliHFEsecVtx &operator=(const AliHFEsecVtx &); // assignment operator virtual ~AliHFEsecVtx(); - - protected: void CreateHistograms(TString hnopt=""); - void Init(); + void SetEvent(AliESDEvent* const ESD){fESD1=ESD;}; // set ESD pointer + void SetStack(AliStack* const stack){fStack=stack;} // set stack pointer void InitAnaPair(); // initialize default parameters - void SetEvent(AliESDEvent* ESD){fESD1=ESD;}; // set ESD pointer - void SetStack(AliStack* stack){fStack=stack;} // set stack pointer - void AnaPair(AliESDtrack* ESDtrack1, AliESDtrack* ESDtrack2, Int_t sourcePart, Int_t index2); // do e-h analysis + void AnaPair(AliESDtrack* ESDtrack1, AliESDtrack* ESDtrack2, Int_t index2); // do e-h analysis void RunSECVTX(AliESDtrack *track); // run secondary vertexing algorithm + + Int_t GetMCPID(AliESDtrack *track); // return MC pid + Int_t PairOrigin(AliESDtrack* track1, AliESDtrack* track2); // return pair origin as a pdg code + Int_t PairCode(AliESDtrack* track1,AliESDtrack* track2); // return corresponding pair code to pdg code + Int_t GetElectronSource(Int_t mclabel); // return origin of the electron + + Bool_t SingleTrackCut(AliESDtrack* track1); // single track cut + + protected: + + void Init(); void FindSECVTXCandid4Tracks(AliESDtrack* track1); // find secondary vertex for 4 tracks void FindSECVTXCandid3Tracks(AliESDtrack* track1); // find secondary vertex for 3 tracks void FindSECVTXCandid2Tracks(AliESDtrack* track1); // find secondary vertex for 2 tracks void CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // calculated distinctive variables void ApplyPairTagCut(); // apply strict tagging condition for 1 pair secondary vertex - void ApplyVtxTagCut(Double_t chi2); // apply tagging condition for 3 track secondary vertex - //void ApplyTagCut(Double_t chi2, AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // apply tagging condition for 3 track secondary vertex + void ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2); // apply tagging condition for 3 track secondary vertex void ResetTagVar(); // reset tagging variables Bool_t ApplyPairTagCut(Int_t id); // apply strict tagging condition for 1 pair secondary vertex Bool_t ApplyTagCut(Double_t chi2); // apply tagging condition - Bool_t SingleTrackCut(AliESDtrack* track1); // single track cut - Int_t GetMCPID(AliESDtrack *track); // return MC pid - Int_t PairOrigin(AliESDtrack* track1, AliESDtrack* track2); // return pair origin as a pdg code - Int_t PairCode(AliESDtrack* track1,AliESDtrack* track2); // return corresponding pair code to pdg code - Int_t GetElectronSource(Int_t mclabel); // return origin of the electron Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4); // get secondary vertex chi2 for 4 tracks Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // get secondary vertex chi2 for 3 tracks @@ -81,22 +84,24 @@ class AliHFEsecVtx : public TObject { AliESDEvent* fESD1; // ESD pointer AliStack* fStack; // stack pointer - Int_t fParentSelect[2][7]; // - Int_t fNparents; // + Int_t fParentSelect[2][7]; // heavy hadron species + Int_t fNparents; // number of heavy hadrons to be considered - TString fkSourceLabel[10]; // + TString fkSourceLabel[10]; // electron source label - enum {fkAll, fkDirectCharm, fkDirectBeauty, fkBeautyCharm, fkGamma, fkPi0, fkElse, fkBeautyGamma, fkBeautyPi0, fkBeautyElse}; - enum {fkCharm=4, fkBeauty=5}; + enum {kAll, kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse}; + enum {kCharm=4, kBeauty=5}; struct histspair{ TH2F *fInvMass; // histogram to store pair invariant mass TH2F *fInvMassCut1; // histogram to store pair invariant mass after cut1 TH2F *fInvMassCut2; // histogram to store pair invariant mass after cut2 TH1F *fKFChi2; // histogram to store pair vertex chi2 - TH1F *fCosOpenAngle; // histogram to store pair opening angle + TH1F *fOpenAngle; // histogram to store pair opening angle + TH1F *fCosOpenAngle; // histogram to store cosine of the pair opening angle TH2F *fSignedLxy; // histogram to store signed Lxy TH1F *fKFIP; // histogram to store pair impact parameter + TH1F *fIPMax; // histogram to store larger impact parameter among pair tracks }; struct histstag{ @@ -104,29 +109,46 @@ class AliHFEsecVtx : public TObject { TH1F *fPtTaggedElec; // histogram for total b tagged electrons TH1F *fPtRightTaggedElec; // histogram for right b tagged electrons TH1F *fPtWrongTaggedElec; // histogram for wrong b tagged electrons + TH1F *fInvmassBeautyElecSecVtx; // histogram for right-tagged b invariant mass + TH1F *fInvmassNotBeautyElecSecVtx; // histogram for mis-tagged b invariant mass + TH1F *fSignedLxyBeautyElecSecVtx; // histogram for right-tagged b signed Lxy + TH1F *fSignedLxyNotBeautyElecSecVtx; // histogram for mis-tagged b signed Lxy + TH1F *fDistTwoVtxBeautyElecSecVtx; // histogram for right-tagged b two vertex distance + TH1F *fDistTwoVtxNotBeautyElecSecVtx; // histogram for mis-tagged b two vertex distance + TH1F *fCosPhiBeautyElecSecVtx; // histogram for right-tagged b cos of opening angle + TH1F *fCosPhiNotBeautyElecSecVtx; // histogram for mis-tagged b cos of opening angle + TH1F *fChi2BeautyElecSecVtx; // histogram for right-tagged b chi2 of secondary vertex + TH1F *fChi2NotBeautyElecSecVtx; // histogram for mis-tagged b chi2 of secondary vertex + TH1F *fInvmassBeautyElec2trkVtx; // histogram for right-tagged b invariant mass of two track vertex + TH1F *fInvmassNotBeautyElec2trkVtx; // histogram for mis-tagged b invariant mass of two track vertex + TH1F *fSignedLxyBeautyElec2trkVtx; // histogram for right-tagged b signed Lxy of two track vertex + TH1F *fSignedLxyNotBeautyElec2trkVtx; // histogram for mis-tagged b signed Lxy of two track vertex + TH1F *fChi2BeautyElec2trkVtx; // histogram for right-tagged b chi2 of two track vertex + TH1F *fChi2NotBeautyElec2trkVtx; // histogram for mis-tagged b chi2 of two track vertex }; - histspair fHistPair[10]; // - histstag fHistTagged; // + histspair fHistPair[10]; // struct of above given histogram for different electron sources + histstag fHistTagged; // struct of above given histogram - Int_t fPairTagged; // - Int_t fpairedTrackID[20]; // - Double_t fpairedChi2[20]; // - Double_t fpairedInvMass[20]; // - Double_t fpairedSignedLxy[20]; // + Int_t fPairTagged; // number of tagged e-h pairs + Int_t fpairedTrackID[20]; // paird hadron track track + Double_t fpairedChi2[20]; // pair chi2 + Double_t fpairedInvMass[20]; // pair invariant mass + Double_t fpairedSignedLxy[20]; // pair signed Lxy - Int_t fid[4][3]; // - Int_t fia[4][3][2]; // + Int_t fid[4][3]; // index to store sorting result + Int_t fia[4][3][2]; // index to store sorting result - Double_t fdistTwoSecVtx; // - Double_t fcosPhi; // + Double_t fdistTwoSecVtx; // distance between two pair vertex + Double_t fcosPhi; // cos of opening angle of two pair vertex Double_t fsignedLxy; // signed Lxy of secondary vertex Double_t finvmass; // invariant mass of secondary vertex Double_t finvmassSigma; // invariant mass sigma of secondary vertex Bool_t fBTagged; // if b tagged, set true + Bool_t fBElec; // if set true for b electron, set true - ClassDef(AliHFEsecVtx,0); // secondary vertex for electrons with AliKFParticle package + ClassDef(AliHFEsecVtx,0); }; #endif -- 2.43.0