From 601736df18ce7d884b9bbe88dd5c1920ed46482f Mon Sep 17 00:00:00 2001 From: dainese Date: Tue, 22 Mar 2011 00:40:24 +0000 Subject: [PATCH] Analysis code for D0->4prong (Fabio) --- .../AliAnalysisTaskSESelectHF4Prong.cxx | 1091 +++++++++++++++++ .../AliAnalysisTaskSESelectHF4Prong.h | 141 +++ PWG3/vertexingHF/AliRDHFCuts.cxx | 8 +- PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.cxx | 262 +++- PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.h | 21 +- .../macros/AddTaskSelectHF4Prong.C | 63 + .../vertexingHF/macros/MakeCuts4Charm4Prong.C | 369 ++++++ 7 files changed, 1937 insertions(+), 18 deletions(-) create mode 100644 PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.cxx create mode 100644 PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.h create mode 100644 PWG3/vertexingHF/macros/AddTaskSelectHF4Prong.C create mode 100644 PWG3/vertexingHF/macros/MakeCuts4Charm4Prong.C diff --git a/PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.cxx b/PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.cxx new file mode 100644 index 00000000000..12760f4151c --- /dev/null +++ b/PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.cxx @@ -0,0 +1,1091 @@ +/************************************************************************** + * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +///////////////////////////////////////////////////////////// +// +// AliAnalysisTaskSE for the selection of heavy flavor +// decay candidates and creation a stand-alone AOD for +// 4prong D0 decay. +// +// Author: A.Dainese, andrea.dainese@lnl.infn.it +// F.Colamaria, fabio.colamaria@ba.infn.it +///////////////////////////////////////////////////////////// + + +#include "Riostream.h" +#include "TFile.h" +#include "TList.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TClonesArray.h" +#include "TDatabasePDG.h" +#include "TROOT.h" +#include "TCanvas.h" +#include "TBits.h" +#include "TNtuple.h" + +#include "AliAnalysisDataSlot.h" +#include "AliAnalysisDataContainer.h" +#include "AliAnalysisManager.h" +#include "AliAODHandler.h" +#include "AliAODEvent.h" +#include "AliAODVertex.h" +#include "AliAODTrack.h" +#include "AliAODRecoDecayHF4Prong.h" +#include "AliAnalysisVertexingHF.h" +#include "AliAnalysisTaskSE.h" +#include "AliAnalysisTaskSESelectHF4Prong.h" +#include "AliAODPidHF.h" +#include "AliRDHFCuts.h" + +ClassImp(AliAnalysisTaskSESelectHF4Prong) + + +//________________________________________________________________________ +AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong(): +AliAnalysisTaskSE(), +fVerticesHFTClArr(0), +fCharm4ProngTClArr(0), +fOutput(0), +fOutput2(0), +fOutput3(0), +fOutput4(0), +fOutput5(0), +fOutputC(0), +fhInvMassD0Sum_10Mev_Bin1(0), +fhInvMassD0barSum_10Mev_Bin1(0), +fhInvMassSumAll_10Mev_Bin1(0), +fhInvMassD0Sum_5Mev_Bin1(0), +fhInvMassD0barSum_5Mev_Bin1(0), +fhInvMassSumAll_5Mev_Bin1(0), +fhInvMassD0Sum_10Mev_Bin2(0), +fhInvMassD0barSum_10Mev_Bin2(0), +fhInvMassSumAll_10Mev_Bin2(0), +fhInvMassD0Sum_5Mev_Bin2(0), +fhInvMassD0barSum_5Mev_Bin2(0), +fhInvMassSumAll_5Mev_Bin2(0), +fhInvMassD0Sum_10Mev_Bin3(0), +fhInvMassD0barSum_10Mev_Bin3(0), +fhInvMassSumAll_10Mev_Bin3(0), +fhInvMassD0Sum_5Mev_Bin3(0), +fhInvMassD0barSum_5Mev_Bin3(0), +fhInvMassSumAll_5Mev_Bin3(0), +fhInvMassD0Sum_10Mev_Bin4(0), +fhInvMassD0barSum_10Mev_Bin4(0), +fhInvMassSumAll_10Mev_Bin4(0), +fhInvMassD0Sum_5Mev_Bin4(0), +fhInvMassD0barSum_5Mev_Bin4(0), +fhInvMassSumAll_5Mev_Bin4(0), +fhInvMassD0Sum_10Mev_Bin5(0), +fhInvMassD0barSum_10Mev_Bin5(0), +fhInvMassSumAll_10Mev_Bin5(0), +fhInvMassD0Sum_5Mev_Bin5(0), +fhInvMassD0barSum_5Mev_Bin5(0), +fhInvMassSumAll_5Mev_Bin5(0), +fhInvMassMultipleOnly_Bin1(0), +fhInvMassMultipleOnly_Bin2(0), +fhInvMassMultipleOnly_Bin3(0), +fhInvMassMultipleOnly_Bin4(0), +fhInvMassMultipleOnly_Bin5(0), +fScatterP4PID(0), +fPtVsY(0), +fPtVsYAll(0), +fEventCounter(0), +fCutDCA(0), +fCutDCA3(0), +fCutDCA2(0), +fCutDCA5(0), +fCutVertexDist2(0), +fCutVertexDist3(0), +fCutVertexDist4(0), +fCutCosinePoint(0), +fCutPt(0), +fCutY(0), +fPIDSel(0), +fPIDSel_Bin1(0), +fPIDSel_Bin2(0), +fPIDSel_Bin3(0), +fPIDSel_Bin4(0), +fPIDSel_Bin5(0), +fMultipleHyps(0), +fMultipleHypsType(0), +fPtSel(0), +fCuts(0), +fVHF(0) +{ + // Default constructor +} + +//________________________________________________________________________ +AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong(const char *name,AliRDHFCutsD0toKpipipi* cuts): +AliAnalysisTaskSE(name), +fVerticesHFTClArr(0), +fCharm4ProngTClArr(0), +fOutput(0), +fOutput2(0), +fOutput3(0), +fOutput4(0), +fOutput5(0), +fOutputC(0), +fhInvMassD0Sum_10Mev_Bin1(0), +fhInvMassD0barSum_10Mev_Bin1(0), +fhInvMassSumAll_10Mev_Bin1(0), +fhInvMassD0Sum_5Mev_Bin1(0), +fhInvMassD0barSum_5Mev_Bin1(0), +fhInvMassSumAll_5Mev_Bin1(0), +fhInvMassD0Sum_10Mev_Bin2(0), +fhInvMassD0barSum_10Mev_Bin2(0), +fhInvMassSumAll_10Mev_Bin2(0), +fhInvMassD0Sum_5Mev_Bin2(0), +fhInvMassD0barSum_5Mev_Bin2(0), +fhInvMassSumAll_5Mev_Bin2(0), +fhInvMassD0Sum_10Mev_Bin3(0), +fhInvMassD0barSum_10Mev_Bin3(0), +fhInvMassSumAll_10Mev_Bin3(0), +fhInvMassD0Sum_5Mev_Bin3(0), +fhInvMassD0barSum_5Mev_Bin3(0), +fhInvMassSumAll_5Mev_Bin3(0), +fhInvMassD0Sum_10Mev_Bin4(0), +fhInvMassD0barSum_10Mev_Bin4(0), +fhInvMassSumAll_10Mev_Bin4(0), +fhInvMassD0Sum_5Mev_Bin4(0), +fhInvMassD0barSum_5Mev_Bin4(0), +fhInvMassSumAll_5Mev_Bin4(0), +fhInvMassD0Sum_10Mev_Bin5(0), +fhInvMassD0barSum_10Mev_Bin5(0), +fhInvMassSumAll_10Mev_Bin5(0), +fhInvMassD0Sum_5Mev_Bin5(0), +fhInvMassD0barSum_5Mev_Bin5(0), +fhInvMassSumAll_5Mev_Bin5(0), +fhInvMassMultipleOnly_Bin1(0), +fhInvMassMultipleOnly_Bin2(0), +fhInvMassMultipleOnly_Bin3(0), +fhInvMassMultipleOnly_Bin4(0), +fhInvMassMultipleOnly_Bin5(0), +fScatterP4PID(0), +fPtVsY(0), +fPtVsYAll(0), +fEventCounter(0), +fCutDCA(0), +fCutDCA3(0), +fCutDCA2(0), +fCutDCA5(0), +fCutVertexDist2(0), +fCutVertexDist3(0), +fCutVertexDist4(0), +fCutCosinePoint(0), +fCutPt(0), +fCutY(0), +fPIDSel(0), +fPIDSel_Bin1(0), +fPIDSel_Bin2(0), +fPIDSel_Bin3(0), +fPIDSel_Bin4(0), +fPIDSel_Bin5(0), +fMultipleHyps(0), +fMultipleHypsType(0), +fPtSel(0), +fCuts(0), +fVHF(0) +{ + // Standard constructor + + fCuts=cuts; + + // Input slot #0 works with an Ntuple +// DefineInput(0, TTree::Class()); + + // Output slot #0 writes into a TTree container + // Output slots #1-6 writes into a TList container + DefineOutput(0, TTree::Class()); //default + DefineOutput(1, TList::Class()); //histos inv. mass bin1 + DefineOutput(2, TList::Class()); //histos inv. mass bin2 + DefineOutput(3, TList::Class()); //histos inv. mass bin3 + DefineOutput(4, TList::Class()); //histos inv. mass bin4 + DefineOutput(5, TList::Class()); //histos inv. mass bin5 + DefineOutput(6, TList::Class()); //histos of cuts + DefineOutput(7, AliRDHFCutsD0toKpipipi::Class()); //cuts +} + +//________________________________________________________________________ +AliAnalysisTaskSESelectHF4Prong::~AliAnalysisTaskSESelectHF4Prong() +{ + // Destructor + + if (fOutput) { + delete fOutput; + fOutput = 0; + } + if (fOutput2) { + delete fOutput2; + fOutput2 = 0; + } + if (fOutput3) { + delete fOutput3; + fOutput3 = 0; + } + if (fOutput4) { + delete fOutput4; + fOutput4 = 0; + } + if (fOutput5) { + delete fOutput5; + fOutput5 = 0; + } + if (fOutputC) { + delete fOutputC; + fOutputC = 0; + } + if (fCuts) { + delete fCuts; + fCuts = 0; + } + if (fVHF) { + delete fVHF; + fVHF = 0; + } + +} + +//________________________________________________________________________ +void AliAnalysisTaskSESelectHF4Prong::Init() +{ + // Initialization + + if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::Init() \n"); + + return; +} + +//________________________________________________________________________ +void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects() +{ + // Create the output container + // + if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::UserCreateOutputObjects() \n"); + + fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0); + fVerticesHFTClArr->SetName("VerticesHF"); + AddAODBranch("TClonesArray", &fVerticesHFTClArr); + + fCharm4ProngTClArr = new TClonesArray("AliAODRecoDecayHF4Prong", 0); + fCharm4ProngTClArr->SetName("Charm4Prong"); + AddAODBranch("TClonesArray", &fCharm4ProngTClArr); + + fOutput = new TList(); + fOutput->SetOwner(); + + fOutput2 = new TList(); + fOutput2->SetOwner(); + + fOutput3 = new TList(); + fOutput3->SetOwner(); + + fOutput4 = new TList(); + fOutput4->SetOwner(); + + fOutput5 = new TList(); + fOutput5->SetOwner(); + + fOutputC = new TList(); + fOutputC->SetOwner(); + + fhInvMassD0Sum_10Mev_Bin1 = new TH1F("fhInvMassD0Sum_10Mev_Bin1", "D0 invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0Sum_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_10Mev_Bin1->SetMinimum(0); + fOutput->Add(fhInvMassD0Sum_10Mev_Bin1); + + fhInvMassD0barSum_10Mev_Bin1 = new TH1F("fhInvMassD0barSum_10Mev_Bin1", "D0bar invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0barSum_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_10Mev_Bin1->SetMinimum(0); + fOutput->Add(fhInvMassD0barSum_10Mev_Bin1); + + fhInvMassSumAll_10Mev_Bin1 = new TH1F("fhInvMassSumAll_10Mev_Bin1", "D0/D0bar invariant mass Bin1 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassSumAll_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_10Mev_Bin1->SetMinimum(0); + fOutput->Add(fhInvMassSumAll_10Mev_Bin1); + + fhInvMassD0Sum_5Mev_Bin1 = new TH1F("fhInvMassD0Sum_5Mev_Bin1", "D0 invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0Sum_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_5Mev_Bin1->SetMinimum(0); + fOutput->Add(fhInvMassD0Sum_5Mev_Bin1); + + fhInvMassD0barSum_5Mev_Bin1 = new TH1F("fhInvMassD0barSum_5Mev_Bin1", "D0bar invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0barSum_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_5Mev_Bin1->SetMinimum(0); + fOutput->Add(fhInvMassD0barSum_5Mev_Bin1); + + fhInvMassSumAll_5Mev_Bin1 = new TH1F("fhInvMassSumAll_5Mev_Bin1", "D0/D0bar invariant mass Bin1 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassSumAll_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_5Mev_Bin1->SetMinimum(0); + fOutput->Add(fhInvMassSumAll_5Mev_Bin1); + + fhInvMassD0Sum_10Mev_Bin2 = new TH1F("fhInvMassD0Sum_10Mev_Bin2", "D0 invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0Sum_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_10Mev_Bin2->SetMinimum(0); + fOutput2->Add(fhInvMassD0Sum_10Mev_Bin2); + + fhInvMassD0barSum_10Mev_Bin2 = new TH1F("fhInvMassD0barSum_10Mev_Bin2", "D0bar invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0barSum_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_10Mev_Bin2->SetMinimum(0); + fOutput2->Add(fhInvMassD0barSum_10Mev_Bin2); + + fhInvMassSumAll_10Mev_Bin2 = new TH1F("fhInvMassSumAll_10Mev_Bin2", "D0/D0bar invariant mass Bin2 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassSumAll_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_10Mev_Bin2->SetMinimum(0); + fOutput2->Add(fhInvMassSumAll_10Mev_Bin2); + + fhInvMassD0Sum_5Mev_Bin2 = new TH1F("fhInvMassD0Sum_5Mev_Bin2", "D0 invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0Sum_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_5Mev_Bin2->SetMinimum(0); + fOutput2->Add(fhInvMassD0Sum_5Mev_Bin2); + + fhInvMassD0barSum_5Mev_Bin2 = new TH1F("fhInvMassD0barSum_5Mev_Bin2", "D0bar invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0barSum_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_5Mev_Bin2->SetMinimum(0); + fOutput2->Add(fhInvMassD0barSum_5Mev_Bin2); + + fhInvMassSumAll_5Mev_Bin2 = new TH1F("fhInvMassSumAll_5Mev_Bin2", "D0/D0bar invariant mass Bin2 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassSumAll_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_5Mev_Bin2->SetMinimum(0); + fOutput2->Add(fhInvMassSumAll_5Mev_Bin2); + + fhInvMassD0Sum_10Mev_Bin3 = new TH1F("fhInvMassD0Sum_10Mev_Bin3", "D0 invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0Sum_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_10Mev_Bin3->SetMinimum(0); + fOutput3->Add(fhInvMassD0Sum_10Mev_Bin3); + + fhInvMassD0barSum_10Mev_Bin3 = new TH1F("fhInvMassD0barSum_10Mev_Bin3", "D0bar invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0barSum_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_10Mev_Bin3->SetMinimum(0); + fOutput3->Add(fhInvMassD0barSum_10Mev_Bin3); + + fhInvMassSumAll_10Mev_Bin3 = new TH1F("fhInvMassSumAll_10Mev_Bin3", "D0/D0bar invariant mass Bin3 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassSumAll_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_10Mev_Bin3->SetMinimum(0); + fOutput3->Add(fhInvMassSumAll_10Mev_Bin3); + + fhInvMassD0Sum_5Mev_Bin3 = new TH1F("fhInvMassD0Sum_5Mev_Bin3", "D0 invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0Sum_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_5Mev_Bin3->SetMinimum(0); + fOutput3->Add(fhInvMassD0Sum_5Mev_Bin3); + + fhInvMassD0barSum_5Mev_Bin3 = new TH1F("fhInvMassD0barSum_5Mev_Bin3", "D0bar invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0barSum_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_5Mev_Bin3->SetMinimum(0); + fOutput3->Add(fhInvMassD0barSum_5Mev_Bin3); + + fhInvMassSumAll_5Mev_Bin3 = new TH1F("fhInvMassSumAll_5Mev_Bin3", "D0/D0bar invariant mass Bin3 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassSumAll_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_5Mev_Bin3->SetMinimum(0); + fOutput3->Add(fhInvMassSumAll_5Mev_Bin3); + + fhInvMassD0Sum_10Mev_Bin4 = new TH1F("fhInvMassD0Sum_10Mev_Bin4", "D0 invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0Sum_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_10Mev_Bin4->SetMinimum(0); + fOutput4->Add(fhInvMassD0Sum_10Mev_Bin4); + + fhInvMassD0barSum_10Mev_Bin4 = new TH1F("fhInvMassD0barSum_10Mev_Bin4", "D0bar invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0barSum_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_10Mev_Bin4->SetMinimum(0); + fOutput4->Add(fhInvMassD0barSum_10Mev_Bin4); + + fhInvMassSumAll_10Mev_Bin4 = new TH1F("fhInvMassSumAll_10Mev_Bin4", "D0/D0bar invariant mass Bin4 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassSumAll_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_10Mev_Bin4->SetMinimum(0); + fOutput4->Add(fhInvMassSumAll_10Mev_Bin4); + + fhInvMassD0Sum_5Mev_Bin4 = new TH1F("fhInvMassD0Sum_5Mev_Bin4", "D0 invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0Sum_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_5Mev_Bin4->SetMinimum(0); + fOutput4->Add(fhInvMassD0Sum_5Mev_Bin4); + + fhInvMassD0barSum_5Mev_Bin4 = new TH1F("fhInvMassD0barSum_5Mev_Bin4", "D0bar invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0barSum_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_5Mev_Bin4->SetMinimum(0); + fOutput4->Add(fhInvMassD0barSum_5Mev_Bin4); + + fhInvMassSumAll_5Mev_Bin4 = new TH1F("fhInvMassSumAll_5Mev_Bin4", "D0/D0bar invariant mass Bin4 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassSumAll_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_5Mev_Bin4->SetMinimum(0); + fOutput4->Add(fhInvMassSumAll_5Mev_Bin4); + + fhInvMassD0Sum_10Mev_Bin5 = new TH1F("fhInvMassD0Sum_10Mev_Bin5", "D0 invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0Sum_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_10Mev_Bin5->SetMinimum(0); + fOutput5->Add(fhInvMassD0Sum_10Mev_Bin5); + + fhInvMassD0barSum_10Mev_Bin5 = new TH1F("fhInvMassD0barSum_10Mev_Bin5", "D0bar invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassD0barSum_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_10Mev_Bin5->SetMinimum(0); + fOutput5->Add(fhInvMassD0barSum_10Mev_Bin5); + + fhInvMassSumAll_10Mev_Bin5 = new TH1F("fhInvMassSumAll_10Mev_Bin5", "D0/D0bar invariant mass Bin5 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2); + fhInvMassSumAll_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_10Mev_Bin5->SetMinimum(0); + fOutput5->Add(fhInvMassSumAll_10Mev_Bin5); + + fhInvMassD0Sum_5Mev_Bin5 = new TH1F("fhInvMassD0Sum_5Mev_Bin5", "D0 invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0Sum_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0Sum_5Mev_Bin5->SetMinimum(0); + fOutput5->Add(fhInvMassD0Sum_5Mev_Bin5); + + fhInvMassD0barSum_5Mev_Bin5 = new TH1F("fhInvMassD0barSum_5Mev_Bin5", "D0bar invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassD0barSum_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassD0barSum_5Mev_Bin5->SetMinimum(0); + fOutput5->Add(fhInvMassD0barSum_5Mev_Bin5); + + fhInvMassSumAll_5Mev_Bin5 = new TH1F("fhInvMassSumAll_5Mev_Bin5", "D0/D0bar invariant mass Bin5 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2); + fhInvMassSumAll_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassSumAll_5Mev_Bin5->SetMinimum(0); + fOutput5->Add(fhInvMassSumAll_5Mev_Bin5); + + fhInvMassMultipleOnly_Bin1 = new TH1F("fhInvMassMultipleOnly_Bin1", "D0/D0bar invariant mass Bin1 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2); + fhInvMassMultipleOnly_Bin1->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassMultipleOnly_Bin1->SetMinimum(0); + fOutput->Add(fhInvMassMultipleOnly_Bin1); + + fhInvMassMultipleOnly_Bin2 = new TH1F("fhInvMassMultipleOnly_Bin2", "D0/D0bar invariant mass Bin2 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2); + fhInvMassMultipleOnly_Bin2->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassMultipleOnly_Bin2->SetMinimum(0); + fOutput2->Add(fhInvMassMultipleOnly_Bin2); + + fhInvMassMultipleOnly_Bin3 = new TH1F("fhInvMassMultipleOnly_Bin3", "D0/D0bar invariant mass Bin3 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2); + fhInvMassMultipleOnly_Bin3->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassMultipleOnly_Bin3->SetMinimum(0); + fOutput3->Add(fhInvMassMultipleOnly_Bin3); + + fhInvMassMultipleOnly_Bin4 = new TH1F("fhInvMassMultipleOnly_Bin4", "D0/D0bar invariant mass Bin4 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2); + fhInvMassMultipleOnly_Bin4->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassMultipleOnly_Bin4->SetMinimum(0); + fOutput4->Add(fhInvMassMultipleOnly_Bin4); + + fhInvMassMultipleOnly_Bin5 = new TH1F("fhInvMassMultipleOnly_Bin5", "D0/D0bar invariant mass Bin5 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2); + fhInvMassMultipleOnly_Bin5->Sumw2(); //Create structure to store sum of squares of weights + fhInvMassMultipleOnly_Bin5->SetMinimum(0); + fOutput5->Add(fhInvMassMultipleOnly_Bin5); + + fScatterP4PID = new TH2F("fScatterP4PID", "Transverse momentum of K vs l-s Pi (D0 + D0bar); Pt of K [GeV/c]; Pt of Pi [GeV/c]",500,0.,5.,500,0.,5.); + fScatterP4PID->SetMinimum(0); + fOutput->Add(fScatterP4PID); + + fPtVsY = new TH2F("fPtVsY", "Pt vs Y PPR Sel. Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.); + fPtVsY->SetMinimum(0); + fOutputC->Add(fPtVsY); + + fPtVsYAll = new TH2F("fPtVsYAll", "Pt vs Y All Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.); + fPtVsYAll->SetMinimum(0); + fOutputC->Add(fPtVsYAll); + + fEventCounter = new TH1F("fEventCounter", "N° of total events; NA; Events",1,0.,1.); + fEventCounter->SetMinimum(0); + fOutputC->Add(fEventCounter); + + fCutDCA = new TH1F("fCutDCA", "DCA of candidate (couple); DCA [cm]; Entries/micron",500,0.,0.05); + fCutDCA->SetMinimum(0); + fOutputC->Add(fCutDCA); + + fCutDCA3 = new TH1F("fCutDCA3", "DCA of candidate (trips); DCA [cm]; Entries/micron",500,0.,0.05); + fCutDCA3->SetMinimum(0); + fOutputC->Add(fCutDCA3); + + fCutDCA2 = new TH1F("fCutDCA2", "DCA of candidate (quads1); DCA [cm]; Entries/micron",500,0.,0.05); + fCutDCA2->SetMinimum(0); + fOutputC->Add(fCutDCA2); + + fCutDCA5 = new TH1F("fCutDCA5", "DCA of candidate (quads2); DCA [cm]; Entries/micron",500,0.,0.05); + fCutDCA5->SetMinimum(0); + fOutputC->Add(fCutDCA5); + + fCutVertexDist2 = new TH1F("fCutVertexDist2", "Distance Vtx doubl.-Primary Vtx; Distance [cm]; Entries/15 micron",500,0.,0.75); + fCutVertexDist2->SetMinimum(0); + fOutputC->Add(fCutVertexDist2); + + fCutVertexDist3 = new TH1F("fCutVertexDist3", "Distance Vtx trips-Primary Vtx; Distance [cm]; Entries/10 micron",500,0.,0.5); + fCutVertexDist3->SetMinimum(0); + fOutputC->Add(fCutVertexDist3); + + fCutVertexDist4 = new TH1F("fCutVertexDist4", "Distance Vtx quads-Primary Vtx; Distance [cm]; Entries/5 micron",500,0.,0.25); + fCutVertexDist4->SetMinimum(0); + fOutputC->Add(fCutVertexDist4); + + fCutCosinePoint = new TH1F("fCutCosinePoint", "Cosine of angle of pointing; Cos(Thetapt.); Entries/10^(-3)",250,0.75,1.); + fCutCosinePoint->SetMinimum(0); + fOutputC->Add(fCutCosinePoint); + + fCutPt = new TH1F("fCutPt", "Pt of candidate D0; Pt [GeV/c]; Entries/5 MeV",3000,0.,15.); + fCutPt->SetMinimum(0); + fOutputC->Add(fCutPt); + + fCutY = new TH1F("fCutY", "Y of candidate D0; Pt [GeV/c]; Entries/5 MeV",900,-9.,9.); + fCutY->SetMinimum(0); + fOutputC->Add(fCutY); + + fPIDSel = new TH1F("fPIDSel", "Ratio of D0 selected by PID for Correction",3,0.,3.); + fPIDSel->SetMinimum(0); + fPIDSel->GetXaxis()->SetBinLabel(1,"D0allhyp All"); + fPIDSel->GetXaxis()->SetBinLabel(2,"D0allhyp PID"); + fPIDSel->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)"); + + fPIDSel_Bin1 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.); + fPIDSel_Bin1->SetMinimum(0); + fPIDSel_Bin1->GetXaxis()->SetBinLabel(1,"D0allhyp All"); + fPIDSel_Bin1->GetXaxis()->SetBinLabel(2,"D0allhyp PID"); + fPIDSel_Bin1->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)"); + + fPIDSel_Bin2 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.); + fPIDSel_Bin2->SetMinimum(0); + fPIDSel_Bin2->GetXaxis()->SetBinLabel(1,"D0allhyp All"); + fPIDSel_Bin2->GetXaxis()->SetBinLabel(2,"D0allhyp PID"); + fPIDSel_Bin2->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)"); + + fPIDSel_Bin3 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.); + fPIDSel_Bin3->SetMinimum(0); + fPIDSel_Bin3->GetXaxis()->SetBinLabel(1,"D0allhyp All"); + fPIDSel_Bin3->GetXaxis()->SetBinLabel(2,"D0allhyp PID"); + fPIDSel_Bin3->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)"); + + fPIDSel_Bin4 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.); + fPIDSel_Bin4->SetMinimum(0); + fPIDSel_Bin4->GetXaxis()->SetBinLabel(1,"D0allhyp All"); + fPIDSel_Bin4->GetXaxis()->SetBinLabel(2,"D0allhyp PID"); + fPIDSel_Bin4->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)"); + + fPIDSel_Bin5 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.); + fPIDSel_Bin5->SetMinimum(0); + fPIDSel_Bin5->GetXaxis()->SetBinLabel(1,"D0allhyp All"); + fPIDSel_Bin5->GetXaxis()->SetBinLabel(2,"D0allhyp PID"); + fPIDSel_Bin5->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)"); + + fMultipleHyps = new TH1F("fMultipleHyps", "N. of hyp. accepted for each candidate (accounted N. times)",8,0.,8.); + fMultipleHyps->SetMinimum(0); + fMultipleHyps->GetXaxis()->SetBinLabel(1,"1 (PPR)"); + fMultipleHyps->GetXaxis()->SetBinLabel(2,"2 (PPR)"); + fMultipleHyps->GetXaxis()->SetBinLabel(3,"3 (PPR)"); + fMultipleHyps->GetXaxis()->SetBinLabel(4,"4 (PPR)"); + fMultipleHyps->GetXaxis()->SetBinLabel(5,"1 (PPR+PID)"); + fMultipleHyps->GetXaxis()->SetBinLabel(6,"2 (PPR+PID)"); + fMultipleHyps->GetXaxis()->SetBinLabel(7,"3 (PPR+PID)"); + fMultipleHyps->GetXaxis()->SetBinLabel(8,"4 (PPR+PID)"); + + fMultipleHypsType = new TH1F("fMultipleHypsType", "Type of hyp. accepted for each candidate",8,0.,8.); + fMultipleHypsType->SetMinimum(0); + fMultipleHypsType->GetXaxis()->SetBinLabel(1,"D0"); + fMultipleHypsType->GetXaxis()->SetBinLabel(2,"D0bar"); + fMultipleHypsType->GetXaxis()->SetBinLabel(3,"2D0"); + fMultipleHypsType->GetXaxis()->SetBinLabel(4,"2D0bar"); + fMultipleHypsType->GetXaxis()->SetBinLabel(5,"D0+D0bar"); + fMultipleHypsType->GetXaxis()->SetBinLabel(6,"2D0+D0bar"); + fMultipleHypsType->GetXaxis()->SetBinLabel(7,"D0+2D0bar"); + fMultipleHypsType->GetXaxis()->SetBinLabel(8,"2D0+2D0bar"); + + fOutputC->Add(fMultipleHyps); + fOutputC->Add(fMultipleHypsType); + + fOutputC->Add(fPIDSel); + fOutput->Add(fPIDSel_Bin1); + fOutput2->Add(fPIDSel_Bin2); + fOutput3->Add(fPIDSel_Bin3); + fOutput4->Add(fPIDSel_Bin4); + fOutput5->Add(fPIDSel_Bin5); + + fPtSel = new TH1F("fPtSel", "Pt of candidates accepted; Pt [GeV/c]; Entries/10 MeV",2000,0.,20.); + fPtSel->SetMinimum(0); + fOutputC->Add(fPtSel); + + return; +} + +//________________________________________________________________________ +void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/) +{ + // Execute analysis for current event: + // heavy flavor candidates selection and histograms + + AliAODEvent *aodIn = dynamic_cast (InputEvent()); + + TClonesArray *inputArrayCharm4Prong = 0; + + if(!aodIn && AODEvent() && IsStandardAOD()) { + // In case there is an AOD handler writing a standard AOD, use the AOD + // event in memory rather than the input (ESD) event. + aodIn = dynamic_cast (AODEvent()); + // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) + // have to taken from the AOD event hold by the AliAODExtension + AliAODHandler* aodHandler = (AliAODHandler*) + ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); + if(aodHandler->GetExtensions()) { + AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); + AliAODEvent *aodFromExt = ext->GetAOD(); + // load D0 candidates + inputArrayCharm4Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong"); + } + } else { + // load D0 candidates + inputArrayCharm4Prong=(TClonesArray*)aodIn->GetList()->FindObject("Charm4Prong"); + } + + if(!inputArrayCharm4Prong) { + printf("AliAnalysisTaskSESelectHF4Prong::UserExec: D0to3Kpi branch not found!\n"); + return; + } + + //print event info +// aodIn->GetHeader()->Print(); + + //Event counter ++ + fEventCounter->Fill(0); + + // fix for temporary bug in ESDfilter + // the AODs with null vertex pointer didn't pass the PhysSel + if(!aodIn->GetPrimaryVertex()) return; + + // primary vertex + AliAODVertex *vtx1 = (AliAODVertex*)aodIn->GetPrimaryVertex(); +// vtx1->Print(); + + // make trkIDtoEntry register (temporary) + Int_t trkIDtoEntry[100000]; + for(Int_t it=0;itGetNumberOfTracks();it++) { + AliAODTrack *track = aodIn->GetTrack(it); + trkIDtoEntry[track->GetID()]=it; + } + + Int_t iOutVerticesHF=0,iOutCharm4Prong=0; + fVerticesHFTClArr->Delete(); + iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast(); + TClonesArray &verticesHFRef = *fVerticesHFTClArr; + fCharm4ProngTClArr->Delete(); + iOutCharm4Prong = fCharm4ProngTClArr->GetEntriesFast(); + TClonesArray &aodCharm4ProngRef = *fCharm4ProngTClArr; + + // loop over D0->K3pi candidates + Int_t nInCharm4Prong = inputArrayCharm4Prong->GetEntriesFast(); + printf("Number of D0->K3pi: %d\n",nInCharm4Prong); + + for (Int_t iCharm4Prong = 0; iCharm4Prong < nInCharm4Prong; iCharm4Prong++) { + AliAODRecoDecayHF4Prong *dIn = (AliAODRecoDecayHF4Prong*)inputArrayCharm4Prong->UncheckedAt(iCharm4Prong); + Bool_t unsetvtx=kFALSE; + + if(!dIn->GetOwnPrimaryVtx()) { + dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables + unsetvtx=kTRUE; + } + + //fill histos of cuts + Double_t dca = dIn->GetDCA(); + Double_t dist2 = dIn->GetDist12toPrim(); + Double_t dist3 = dIn->GetDist3toPrim(); + Double_t dist4 = dIn->GetDist4toPrim(); + Double_t cosine = dIn->CosPointingAngle(); + Double_t ptPart = dIn->Pt(); + Double_t yPart = dIn->YD0(); + Double_t dcatrip = dIn->GetDCA(3); + Double_t dcaquad1 = dIn->GetDCA(2); + Double_t dcaquad2 = dIn->GetDCA(5); + + Double_t ptBinH[6] = {2.5,4.5,6.0,8.0,12.0,25.0}; //bin i has pt between values i and i+1 + + fCutDCA->Fill(dca); + fCutDCA3->Fill(dcatrip); + fCutDCA2->Fill(dcaquad1); + fCutDCA5->Fill(dcaquad2); + fCutVertexDist2->Fill(dist2); + fCutVertexDist3->Fill(dist3); + fCutVertexDist4->Fill(dist4); + fCutCosinePoint->Fill(cosine); + fCutPt->Fill(ptPart); + fCutY->Fill(yPart); + fPtVsYAll->Fill(ptPart,yPart); + + //flags initialization + fSelected = fCuts->IsSelected(dIn,AliRDHFCuts::kCandidate); + Int_t selD01 = fCuts->D01Selected(dIn,AliRDHFCuts::kCandidate); + Int_t selD02 = fCuts->D02Selected(dIn,AliRDHFCuts::kCandidate); + Int_t selD0bar1 = fCuts->D0bar1Selected(dIn,AliRDHFCuts::kCandidate); + Int_t selD0bar2 = fCuts->D0bar2Selected(dIn,AliRDHFCuts::kCandidate); + Int_t flagAccLim = 1; + + //Limited Acceptance + if(ptPart > 5.) { + if (TMath::Abs(yPart) > 0.8) flagAccLim = 0; + } + else { + Double_t maxFiducialY = -0.2/15*ptPart*ptPart+1.9/15*ptPart+0.5; + Double_t minFiducialY = 0.2/15*ptPart*ptPart-1.9/15*ptPart-0.5; + if (yPart < minFiducialY || yPart > maxFiducialY) flagAccLim = 0;; + } + + //number of CANDIDATES (regardless of hypotheses) passing PPR + if(fSelected==1||fSelected==2||fSelected==3) { + fPtVsY->Fill(ptPart,yPart); + fPIDSel->Fill(0); + if (ptPart >= ptBinH[0] && ptPart < ptBinH[1]) fPIDSel_Bin1->Fill(0); + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(0); + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(0); + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(0); + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(0); + } + + + PostData(6,fOutputC); + + //selection + if((fSelected==1||fSelected==2||fSelected==3) && flagAccLim == 1) { + // get daughter AOD tracks + AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0); + AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1); + AliAODTrack *trk2 = (AliAODTrack*)dIn->GetDaughter(2); + AliAODTrack *trk3 = (AliAODTrack*)dIn->GetDaughter(3); + if(!trk0 || !trk1 || !trk2 || !trk3) { + trk0=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(0)]); + trk1=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(1)]); + trk2=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(2)]); + trk3=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(3)]); + } + printf("pt of positive track #1: %f\n",trk0->Pt()); + printf("pt of negative track #1: %f\n",trk1->Pt()); + printf("pt of positive track #2: %f\n",trk2->Pt()); + printf("pt of negative track #2: %f\n",trk3->Pt()); + + dIn->InvMassD0(fmassD0); + dIn->InvMassD0bar(fmassD0bar); + + //fill histos (combining selection form cuts & PID (with rho information)) + Int_t hypD01 = 0, hypD02 = 0, hypD0bar1 = 0, hypD0bar2 = 0; + Int_t pid1 = 0, pid2 = 0, pidbar1 = 0, pidbar2 = 0; + Int_t pidSelection = fCuts->IsSelectedFromPID(dIn, &pid1, &pid2, &pidbar1, &pidbar2); + + //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR may not be the same! + if (pidSelection > 0) { + fPIDSel->Fill(1); + if (ptPart >= ptBinH[0] && ptPart < ptBinH[1]) fPIDSel_Bin1->Fill(1); + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(1); + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(1); + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(1); + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(1); + } + + //number of hypoteses accepted per candidate after PPR + if(selD01+selD02+selD0bar1+selD0bar2 == 1) fMultipleHyps->Fill(0); + if(selD01+selD02+selD0bar1+selD0bar2 == 2) fMultipleHyps->Fill(1); + if(selD01+selD02+selD0bar1+selD0bar2 == 3) fMultipleHyps->Fill(2); + if(selD01+selD02+selD0bar1+selD0bar2 == 4) fMultipleHyps->Fill(3); + + //combine PPD + PID cuts + if (selD01 == 1 && pid1==1) hypD01 = 1; + if (selD02 == 1 && pid2==1) hypD02 = 1; + if (selD0bar1 == 1 && pidbar1==1) hypD0bar1 = 1; + if (selD0bar2 == 1 && pidbar2==1) hypD0bar2 = 1; + + //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR must match (at least one)! + if (hypD01 == 1 || hypD02 == 1 || hypD0bar1 == 1 || hypD0bar2 == 1) { + fPIDSel->Fill(2); + if (ptPart >= ptBinH[0] && ptPart < ptBinH[1]) fPIDSel_Bin1->Fill(2); + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(2); + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(2); + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(2); + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(2); + } + + //number of hypoteses accepted per candidate after PPR and PID + if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 1) fMultipleHyps->Fill(4); + if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 2) fMultipleHyps->Fill(5); + if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 3) fMultipleHyps->Fill(6); + if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 4) fMultipleHyps->Fill(7); + + //type of hypoteses accepted per candidate after PPR and PID + if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill(0); + if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(1); + if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill(2); + if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(3); + if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(4); + if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(5); + if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(6); + if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(7); + + +// +// This section is auxiliary: (multiple hypotheses histos) +// + if (hypD01+hypD02+hypD0bar1+hypD0bar2 == 2 || hypD01+hypD02+hypD0bar1+hypD0bar2 == 3 || hypD01+hypD02+hypD0bar1+hypD0bar2 == 4) + { + + if (ptPart > ptBinH[0]) { + // D01 hyp. + if(hypD01==1) { + if (ptPart < ptBinH[1]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0[0]); + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0[0]); + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0[0]); + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0[0]); + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0[0]); + } + // D02 hyp. + if(hypD02==1) { + if (ptPart < ptBinH[1]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0[1]); + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0[1]); + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0[1]); + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0[1]); + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0[1]); + } + // D0bar1 hyp. + if(hypD0bar1==1) { + if (ptPart < ptBinH[1]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0bar[0]); + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0bar[0]); + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0bar[0]); + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0bar[0]); + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0bar[0]); + } + // D0bar2 hyp. + if(hypD0bar2==1) { + if (ptPart < ptBinH[1]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0bar[1]); + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0bar[1]); + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0bar[1]); + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0bar[1]); + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0bar[1]); + } + } + + } +// +//end of auxiliary section +// + + + //All histos are filled if Pt of candidate is greater than minimum of first bin (in this way: bin1+bin2+...binN = whole) + if (ptPart > ptBinH[0]) { + + // D01 hyp. + if(hypD01==1) { + fPtSel->Fill(ptPart); + fScatterP4PID->Fill(trk1->Pt(),trk3->Pt()); + if (ptPart < ptBinH[1]) {fhInvMassD0Sum_10Mev_Bin1->Fill(fmassD0[0]); + fhInvMassD0Sum_5Mev_Bin1->Fill(fmassD0[0]); + fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0[0]); + fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0[0]);} + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum_10Mev_Bin2->Fill(fmassD0[0]); + fhInvMassD0Sum_5Mev_Bin2->Fill(fmassD0[0]); + fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0[0]); + fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0[0]);} + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum_10Mev_Bin3->Fill(fmassD0[0]); + fhInvMassD0Sum_5Mev_Bin3->Fill(fmassD0[0]); + fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0[0]); + fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0[0]);} + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum_10Mev_Bin4->Fill(fmassD0[0]); + fhInvMassD0Sum_5Mev_Bin4->Fill(fmassD0[0]); + fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0[0]); + fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0[0]);} + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0Sum_10Mev_Bin5->Fill(fmassD0[0]); + fhInvMassD0Sum_5Mev_Bin5->Fill(fmassD0[0]); + fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0[0]); + fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0[0]);} + } + // D02 hyp. + if(hypD02==1) { + fPtSel->Fill(ptPart); + fScatterP4PID->Fill(trk3->Pt(),trk1->Pt()); + if (ptPart < ptBinH[1]) {fhInvMassD0Sum_10Mev_Bin1->Fill(fmassD0[1]); + fhInvMassD0Sum_5Mev_Bin1->Fill(fmassD0[1]); + fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0[1]); + fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0[1]);} + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum_10Mev_Bin2->Fill(fmassD0[1]); + fhInvMassD0Sum_5Mev_Bin2->Fill(fmassD0[1]); + fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0[1]); + fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0[1]);} + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum_10Mev_Bin3->Fill(fmassD0[1]); + fhInvMassD0Sum_5Mev_Bin3->Fill(fmassD0[1]); + fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0[1]); + fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0[1]);} + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum_10Mev_Bin4->Fill(fmassD0[1]); + fhInvMassD0Sum_5Mev_Bin4->Fill(fmassD0[1]); + fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0[1]); + fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0[1]);} + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0Sum_10Mev_Bin5->Fill(fmassD0[1]); + fhInvMassD0Sum_5Mev_Bin5->Fill(fmassD0[1]); + fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0[1]); + fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0[1]);} + } + // D0bar1 hyp. + if(hypD0bar1==1) { + fPtSel->Fill(ptPart); + fScatterP4PID->Fill(trk0->Pt(),trk2->Pt()); + if (ptPart < ptBinH[1]) {fhInvMassD0barSum_10Mev_Bin1->Fill(fmassD0bar[0]); + fhInvMassD0barSum_5Mev_Bin1->Fill(fmassD0bar[0]); + fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0bar[0]); + fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0bar[0]);} + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum_10Mev_Bin2->Fill(fmassD0bar[0]); + fhInvMassD0barSum_5Mev_Bin2->Fill(fmassD0bar[0]); + fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0bar[0]); + fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0bar[0]);} + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum_10Mev_Bin3->Fill(fmassD0bar[0]); + fhInvMassD0barSum_5Mev_Bin3->Fill(fmassD0bar[0]); + fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0bar[0]); + fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0bar[0]);} + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum_10Mev_Bin4->Fill(fmassD0bar[0]); + fhInvMassD0barSum_5Mev_Bin4->Fill(fmassD0bar[0]); + fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0bar[0]); + fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0bar[0]);} + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0barSum_10Mev_Bin5->Fill(fmassD0bar[0]); + fhInvMassD0barSum_5Mev_Bin5->Fill(fmassD0bar[0]); + fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0bar[0]); + fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0bar[0]);} + } + // D0bar2 hyp. + if(hypD0bar2==1) { + fPtSel->Fill(ptPart); + fScatterP4PID->Fill(trk2->Pt(),trk0->Pt()); + if (ptPart < ptBinH[1]) {fhInvMassD0barSum_10Mev_Bin1->Fill(fmassD0bar[1]); + fhInvMassD0barSum_5Mev_Bin1->Fill(fmassD0bar[1]); + fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0bar[1]); + fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0bar[1]);} + else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum_10Mev_Bin2->Fill(fmassD0bar[1]); + fhInvMassD0barSum_5Mev_Bin2->Fill(fmassD0bar[1]); + fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0bar[1]); + fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0bar[1]);} + else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum_10Mev_Bin3->Fill(fmassD0bar[1]); + fhInvMassD0barSum_5Mev_Bin3->Fill(fmassD0bar[1]); + fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0bar[1]); + fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0bar[1]);} + else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum_10Mev_Bin4->Fill(fmassD0bar[1]); + fhInvMassD0barSum_5Mev_Bin4->Fill(fmassD0bar[1]); + fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0bar[1]); + fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0bar[1]);} + else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0barSum_10Mev_Bin5->Fill(fmassD0bar[1]); + fhInvMassD0barSum_5Mev_Bin5->Fill(fmassD0bar[1]); + fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0bar[1]); + fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0bar[1]);} + } + } + + PostData(1,fOutput); + PostData(2,fOutput2); + PostData(3,fOutput3); + PostData(4,fOutput4); + PostData(5,fOutput5); + + // HERE ONE COULD RECALCULATE THE VERTEX USING THE KF PACKAGE + + // clone candidate for output AOD + if(hypD01||hypD02||hypD0bar1||hypD0bar2) { + AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++]) + AliAODVertex(*(dIn->GetSecondaryVtx())); + AliAODRecoDecayHF4Prong *dOut=new(aodCharm4ProngRef[iOutCharm4Prong++]) + AliAODRecoDecayHF4Prong(*dIn); + dOut->SetSecondaryVtx(v); + dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone())); + v->SetParent(dOut); + } + } //end of selection loop (starts with fSelected == 1, 2 or 3) + if(unsetvtx) dIn->UnsetOwnPrimaryVtx(); + } // end loop on D0->K3pi + + printf("Number of selected D0->K3pi: %d\n",iOutCharm4Prong); + + return; +} + +//________________________________________________________________________ +void AliAnalysisTaskSESelectHF4Prong::Terminate(Option_t */*option*/) +{ + // Terminate analysis + // + if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong: Terminate() \n"); + + fOutput = dynamic_cast (GetOutputData(1)); + if (!fOutput) { + printf("ERROR: fOutput not available\n"); + return; + } + fOutput2 = dynamic_cast (GetOutputData(2)); + if (!fOutput2) { + printf("ERROR: fOutput not available\n"); + return; + } + fOutput3 = dynamic_cast (GetOutputData(3)); + if (!fOutput3) { + printf("ERROR: fOutput not available\n"); + return; + } + fOutput4 = dynamic_cast (GetOutputData(4)); + if (!fOutput4) { + printf("ERROR: fOutput not available\n"); + return; + } + if (!fOutput5) { + printf("ERROR: fOutput not available\n"); + return; + } + fOutputC = dynamic_cast (GetOutputData(5)); + if (!fOutputC) { + printf("ERROR: fOutputC not available\n"); + return; + } + + fhInvMassD0Sum_10Mev_Bin1 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin1")); + fhInvMassD0barSum_10Mev_Bin1 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin1")); + fhInvMassSumAll_10Mev_Bin1 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin1")); + fhInvMassD0Sum_5Mev_Bin1 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin1")); + fhInvMassD0barSum_5Mev_Bin1 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin1")); + fhInvMassSumAll_5Mev_Bin1 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin1")); + fhInvMassD0Sum_10Mev_Bin2 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin2")); + fhInvMassD0barSum_10Mev_Bin2 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin2")); + fhInvMassSumAll_10Mev_Bin2 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin2")); + fhInvMassD0Sum_5Mev_Bin2 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin2")); + fhInvMassD0barSum_5Mev_Bin2 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin2")); + fhInvMassSumAll_5Mev_Bin2 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin2")); + fhInvMassD0Sum_10Mev_Bin3 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin3")); + fhInvMassD0barSum_10Mev_Bin3 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin3")); + fhInvMassSumAll_10Mev_Bin3 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin3")); + fhInvMassD0Sum_5Mev_Bin3 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin3")); + fhInvMassD0barSum_5Mev_Bin3 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin3")); + fhInvMassSumAll_5Mev_Bin3 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin3")); + fhInvMassD0Sum_10Mev_Bin4 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin4")); + fhInvMassD0barSum_10Mev_Bin4 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin4")); + fhInvMassSumAll_10Mev_Bin4 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin4")); + fhInvMassD0Sum_5Mev_Bin4 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin4")); + fhInvMassD0barSum_5Mev_Bin4 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin4")); + fhInvMassSumAll_5Mev_Bin4 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin4")); + fhInvMassD0Sum_10Mev_Bin5 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin5")); + fhInvMassD0barSum_10Mev_Bin5 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin5")); + fhInvMassSumAll_10Mev_Bin5 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin5")); + fhInvMassD0Sum_5Mev_Bin5 = dynamic_cast(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin5")); + fhInvMassD0barSum_5Mev_Bin5 = dynamic_cast(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin5")); + fhInvMassSumAll_5Mev_Bin5 = dynamic_cast(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin5")); + + fScatterP4PID = dynamic_cast(fOutput->FindObject("fScatterP4PID")); + fPtVsY = dynamic_cast(fOutputC->FindObject("fPtVsY")); + fPtVsYAll = dynamic_cast(fOutputC->FindObject("fPtVsYAll")); + + fEventCounter = dynamic_cast(fOutputC->FindObject("fEventCounter")); + + fCutDCA = dynamic_cast(fOutputC->FindObject("fCutDCA")); + fCutDCA3 = dynamic_cast(fOutputC->FindObject("fCutDCA3")); + fCutDCA2 = dynamic_cast(fOutputC->FindObject("fCutDCA2")); + fCutDCA5 = dynamic_cast(fOutputC->FindObject("fCutDCA5")); + fCutVertexDist2 = dynamic_cast(fOutputC->FindObject("fCutVertexDist2")); + fCutVertexDist3 = dynamic_cast(fOutputC->FindObject("fCutVertexDist3")); + fCutVertexDist4 = dynamic_cast(fOutputC->FindObject("fCutVertexDist4")); + fCutCosinePoint = dynamic_cast(fOutputC->FindObject("fCutCosinePoint")); + + fCutPt = dynamic_cast(fOutputC->FindObject("fCutPt")); + fCutY = dynamic_cast(fOutputC->FindObject("fCutY")); + fPIDSel = dynamic_cast(fOutputC->FindObject("fPIDSel")); + fPIDSel_Bin1 = dynamic_cast(fOutputC->FindObject("fPIDSel_Bin1")); + fPIDSel_Bin2 = dynamic_cast(fOutputC->FindObject("fPIDSel_Bin2")); + fPIDSel_Bin3 = dynamic_cast(fOutputC->FindObject("fPIDSel_Bin3")); + fPIDSel_Bin4 = dynamic_cast(fOutputC->FindObject("fPIDSel_Bin4")); + fPIDSel_Bin5 = dynamic_cast(fOutputC->FindObject("fPIDSel_Bin5")); +} + diff --git a/PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.h b/PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.h new file mode 100644 index 00000000000..d988cb1a8d7 --- /dev/null +++ b/PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.h @@ -0,0 +1,141 @@ +#ifndef ALIANALYSISTASKSESELECTHF4PRONG_H +#define ALIANALYSISTASKSESELECTHF4PRONG_H + +/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//************************************************************************* +// Class AliAnalysisTaskSESelectHF4Prong +// AliAnalysisTaskSE for the selection of heavy-flavour decay candidates +// and creation of a stand-alone AOD for 4prong D0 decay +// Author: A.Dainese, andrea.dainese@lnl.infn.it +// F.Colamaria, fabio.colamaria@ba.infn.it +//************************************************************************* + +#include +#include +#include +#include +#include +#include +#include + +#include "AliAnalysisTaskSE.h" +#include "AliAODEvent.h" +#include "AliAnalysisManager.h" +#include "AliAnalysisVertexingHF.h" +#include "AliRDHFCuts.h" +#include "AliRDHFCutsD0toKpipipi.h" + + +class AliAnalysisTaskSESelectHF4Prong : public AliAnalysisTaskSE +{ + public: + + AliAnalysisTaskSESelectHF4Prong(); + AliAnalysisTaskSESelectHF4Prong(const char *name,AliRDHFCutsD0toKpipipi* cuts); + virtual ~AliAnalysisTaskSESelectHF4Prong(); + + + // Implementation of interface methods + virtual void UserCreateOutputObjects(); + virtual void Init(); + virtual void LocalInit() {Init();} + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *option); + + private: + + AliAnalysisTaskSESelectHF4Prong(const AliAnalysisTaskSESelectHF4Prong &source); + AliAnalysisTaskSESelectHF4Prong& operator=(const AliAnalysisTaskSESelectHF4Prong& source); + TClonesArray *fVerticesHFTClArr; //! Array of heavy-flavour vertices + TClonesArray *fCharm4ProngTClArr; //! Array of D0->K3pi + + Double_t fmassD0[2]; + Double_t fmassD0bar[2]; + Int_t fSelected; + + TList *fOutput; //! list send on output slot 1 + TList *fOutput2; //! list send on output slot 2 + TList *fOutput3; //! list send on output slot 3 + TList *fOutput4; //! list send on output slot 4 + TList *fOutput5; //! list send on output slot 5 + TList *fOutputC; //! list send on output slot 6 + + //output histograms + TH1F *fhInvMassD0Sum_10Mev_Bin1; //! Invariant mass D01+D02 (good hyp) _10Mev BIN1 + TH1F *fhInvMassD0barSum_10Mev_Bin1; //! Invariant mass D0bar1+D0bar2 (good hyp) _10Mev + TH1F *fhInvMassSumAll_10Mev_Bin1; //! Invariant mass superimpose (good hyp only)_10Mev + TH1F *fhInvMassD0Sum_5Mev_Bin1; //! Invariant mass D01+D02 (good hyp) _5Mev + TH1F *fhInvMassD0barSum_5Mev_Bin1; //! Invariant mass D0bar1+D0bar2 (good hyp) _5Mev + TH1F *fhInvMassSumAll_5Mev_Bin1; //! Invariant mass superimpose (good hyp only)_5Mev + + TH1F *fhInvMassD0Sum_10Mev_Bin2; //! Invariant mass D01+D02 (good hyp) _10Mev BIN2 + TH1F *fhInvMassD0barSum_10Mev_Bin2; //! Invariant mass D0bar1+D0bar2 (good hyp) _10Mev + TH1F *fhInvMassSumAll_10Mev_Bin2; //! Invariant mass superimpose (good hyp only)_10Mev + TH1F *fhInvMassD0Sum_5Mev_Bin2; //! Invariant mass D01+D02 (good hyp) _5Mev + TH1F *fhInvMassD0barSum_5Mev_Bin2; //! Invariant mass D0bar1+D0bar2 (good hyp) _5Mev + TH1F *fhInvMassSumAll_5Mev_Bin2; //! Invariant mass superimpose (good hyp only)_5Mev + + TH1F *fhInvMassD0Sum_10Mev_Bin3; //! Invariant mass D01+D02 (good hyp) _10Mev BIN3 + TH1F *fhInvMassD0barSum_10Mev_Bin3; //! Invariant mass D0bar1+D0bar2 (good hyp) _10Mev + TH1F *fhInvMassSumAll_10Mev_Bin3; //! Invariant mass superimpose (good hyp only)_10Mev + TH1F *fhInvMassD0Sum_5Mev_Bin3; //! Invariant mass D01+D02 (good hyp) _5Mev + TH1F *fhInvMassD0barSum_5Mev_Bin3; //! Invariant mass D0bar1+D0bar2 (good hyp) _5Mev + TH1F *fhInvMassSumAll_5Mev_Bin3; //! Invariant mass superimpose (good hyp only)_5Mev + + TH1F *fhInvMassD0Sum_10Mev_Bin4; //! Invariant mass D01+D02 (good hyp) _10Mev BIN4 + TH1F *fhInvMassD0barSum_10Mev_Bin4; //! Invariant mass D0bar1+D0bar2 (good hyp) _10Mev + TH1F *fhInvMassSumAll_10Mev_Bin4; //! Invariant mass superimpose (good hyp only)_10Mev + TH1F *fhInvMassD0Sum_5Mev_Bin4; //! Invariant mass D01+D02 (good hyp) _5Mev + TH1F *fhInvMassD0barSum_5Mev_Bin4; //! Invariant mass D0bar1+D0bar2 (good hyp) _5Mev + TH1F *fhInvMassSumAll_5Mev_Bin4; //! Invariant mass superimpose (good hyp only)_5Mev + + TH1F *fhInvMassD0Sum_10Mev_Bin5; //! Invariant mass D01+D02 (good hyp) _10Mev BIN5 + TH1F *fhInvMassD0barSum_10Mev_Bin5; //! Invariant mass D0bar1+D0bar2 (good hyp) _10Mev + TH1F *fhInvMassSumAll_10Mev_Bin5; //! Invariant mass superimpose (good hyp only)_10Mev + TH1F *fhInvMassD0Sum_5Mev_Bin5; //! Invariant mass D01+D02 (good hyp) _5Mev + TH1F *fhInvMassD0barSum_5Mev_Bin5; //! Invariant mass D0bar1+D0bar2 (good hyp) _5Mev + TH1F *fhInvMassSumAll_5Mev_Bin5; //! Invariant mass superimpose (good hyp only)_5Mev + + TH1F *fhInvMassMultipleOnly_Bin1; //! Invariant mass superimpose good hyp only for multiple hyps accepted Bin1 + TH1F *fhInvMassMultipleOnly_Bin2; //! Invariant mass superimpose good hyp only for multiple hyps accepted Bin2 + TH1F *fhInvMassMultipleOnly_Bin3; //! Invariant mass superimpose good hyp only for multiple hyps accepted Bin3 + TH1F *fhInvMassMultipleOnly_Bin4; //! Invariant mass superimpose good hyp only for multiple hyps accepted Bin4 + TH1F *fhInvMassMultipleOnly_Bin5; //! Invariant mass superimpose good hyp only for multiple hyps accepted Bin5 + + TH2F *fScatterP4PID; //! K momentum vs like sign Pi momentum after PID + TH2F *fPtVsY; //! Pt vs Y of selected candidates (by PPR cuts) + TH2F *fPtVsYAll; //! Pt vs Y of all candidates + + TH1F *fEventCounter; //! Event Counter + TH1F *fCutDCA; //! DCA histogram doubl. + TH1F *fCutDCA3; //! DCA histogram trips + TH1F *fCutDCA2; //! DCA histogram quads1 + TH1F *fCutDCA5; //! DCA histogram quads2 + TH1F *fCutVertexDist2; //! Vertex doubl. to primary distance + TH1F *fCutVertexDist3; //! Vertex trips to primary distance + TH1F *fCutVertexDist4; //! Vertex quads to primary distance + TH1F *fCutCosinePoint; //! Cosine of pointing angle + TH1F *fCutPt; //! Candidate D0 Pt + TH1F *fCutY; //! Candidate D0 Y + TH1F *fPIDSel; //! PID Selected + TH1F *fPIDSel_Bin1; //! PID Selected Bin1 + TH1F *fPIDSel_Bin2; //! PID Selected Bin2 + TH1F *fPIDSel_Bin3; //! PID Selected Bin3 + TH1F *fPIDSel_Bin4; //! PID Selected Bin4 + TH1F *fPIDSel_Bin5; //! PID Selected Bin5 + TH1F *fMultipleHyps; //! Multiple hypotesis accepted counter + TH1F *fMultipleHypsType; //! Multiple hypotesis accepted counter + + TH1F *fPtSel; //! Pt of selected candidates + + AliRDHFCutsD0toKpipipi *fCuts; //! Cuts container + + AliAnalysisVertexingHF *fVHF; // analysis (used to pass the cuts) + + ClassDef(AliAnalysisTaskSESelectHF4Prong,2); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates +}; + +#endif + diff --git a/PWG3/vertexingHF/AliRDHFCuts.cxx b/PWG3/vertexingHF/AliRDHFCuts.cxx index 49984246b13..2e44a1c62d9 100644 --- a/PWG3/vertexingHF/AliRDHFCuts.cxx +++ b/PWG3/vertexingHF/AliRDHFCuts.cxx @@ -258,11 +258,15 @@ Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) { if(vertex->GetNContributors()GetZ())>fMaxVtxZ) return kFALSE; + if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) { + fWhyRejection=6; + return kFALSE; + } // switch to settings for 1-pad cls in TPC if(fPidHF) { - if(event->GetRunNumber()>121693 && event->GetRunNumber()<136851) + if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) || + event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE); if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE); diff --git a/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.cxx b/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.cxx index 4dc58069af0..7109e5b618d 100644 --- a/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.cxx +++ b/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.cxx @@ -13,13 +13,12 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ - ///////////////////////////////////////////////////////////// // // Class for cuts on AOD reconstructed D0->Kpipipi // -// Author: r.romita@gsi.de, andrea.dainese@pd.infn.it +// Author: r.romita@gsi.de, andrea.dainese@pd.infn.it, +// fabio.colamaria@ba.infn.it ///////////////////////////////////////////////////////////// #include @@ -29,6 +28,7 @@ #include "AliAODRecoDecayHF4Prong.h" #include "AliAODTrack.h" #include "AliESDtrack.h" +#include "AliAODPidHF.h" ClassImp(AliRDHFCutsD0toKpipipi) @@ -189,7 +189,7 @@ Int_t AliRDHFCutsD0toKpipipi::IsSelected(TObject* obj,Int_t selectionLevel) { Double_t ptD=d->Pt(); if(ptDfMaxPtCand) return 0; + if(ptD>fMaxPtCand) return 0; // selection on daughter tracks if(selectionLevel==AliRDHFCuts::kAll || @@ -205,8 +205,8 @@ Int_t AliRDHFCutsD0toKpipipi::IsSelected(TObject* obj,Int_t selectionLevel) { selectionLevel==AliRDHFCuts::kCandidate) { Int_t ptbin=PtBin(d->Pt()); - - Int_t okD0=1,okD0bar=1; + + Int_t okD0=1,okD0bar=1; Double_t mD0[2],mD0bar[2]; Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass(); @@ -218,19 +218,263 @@ Int_t AliRDHFCutsD0toKpipipi::IsSelected(TObject* obj,Int_t selectionLevel) { TMath::Abs(mD0bar[1]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0; if(!okD0 && !okD0bar) return 0; - if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0; + if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)] + || d->GetDCA(3) > fCutsRD[GetGlobalIndex(1,ptbin)] + || d->GetDCA(2) > fCutsRD[GetGlobalIndex(1,ptbin)] + || d->GetDCA(5) > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0; if(d->GetDist12toPrim() < fCutsRD[GetGlobalIndex(2,ptbin)]) return 0; if(d->GetDist3toPrim() < fCutsRD[GetGlobalIndex(3,ptbin)]) return 0; if(d->GetDist4toPrim() < fCutsRD[GetGlobalIndex(4,ptbin)]) return 0; if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(5,ptbin)]) return 0; if(d->Pt() < fCutsRD[GetGlobalIndex(6,ptbin)]) return 0; if(!d->CutRhoMass(mD0,mD0bar,fCutsRD[GetGlobalIndex(0,ptbin)],fCutsRD[GetGlobalIndex(7,ptbin)])) return 0; - + if (okD0) returnvalue=1; //cuts passed as D0 if (okD0bar) returnvalue=2; //cuts passed as D0bar if (okD0 && okD0bar) returnvalue=3; //cuts passed as D0 and D0bar } - + + // selection on PID (from AliAODPidHF) + if(selectionLevel==AliRDHFCuts::kAll || + selectionLevel==AliRDHFCuts::kPID) { + + Int_t selD01 = D01Selected(d,AliRDHFCuts::kCandidate); + Int_t selD02 = D02Selected(d,AliRDHFCuts::kCandidate); + Int_t selD0bar1 = D0bar1Selected(d,AliRDHFCuts::kCandidate); + Int_t selD0bar2 = D0bar2Selected(d,AliRDHFCuts::kCandidate); + + Int_t d01PID = 0, d02PID = 0, d0bar1PID = 0, d0bar2PID = 0; + returnvalue = IsSelectedFromPID(d, &d01PID, &d02PID, &d0bar1PID, &d0bar2PID); //This returnvalue is dummy! Now it's modified as it must be! + +returnvalue = 0; + + if((selD01 == 1 && d01PID == 1)||(selD02 == 1 && d02PID == 1)||(selD0bar1 == 1 && d0bar1PID == 1)||(selD0bar2 == 1 && d0bar2PID == 1)) returnvalue = 1; + } + return returnvalue; +} + +//--------------------------------------------------------------------------- +Int_t AliRDHFCutsD0toKpipipi::IsSelectedFromPID(AliAODRecoDecayHF4Prong *d, Int_t *hyp1, Int_t *hyp2, Int_t *hyp3, Int_t *hyp4) { + // + // Apply selection (using AliAODPidHF methods) + // Mass hypothesis true if each particle is at least compatible with specie of hypothesis + // + + Int_t output=0; + + Int_t matchK[4], matchPi[4]; + Double_t ptlimit[2] = {0.6,0.8}; + AliAODTrack* trk[4]; + trk[0] = (AliAODTrack*)d->GetDaughter(0); + trk[1] = (AliAODTrack*)d->GetDaughter(1); + trk[2] = (AliAODTrack*)d->GetDaughter(2); + trk[3] = (AliAODTrack*)d->GetDaughter(3); +// if(!trk[0] || !trk[1] || !trk[2] || !trk[3]) { //REMOVED (needs also AliAODEvent to be passed, here and in IsSelected method) +// trk[0]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(0)]); +// trk[1]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(1)]); +// trk[2]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(2)]); +// trk[3]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(3)]);} + + AliAODPidHF* PidObj = new AliAODPidHF(); + PidObj->SetAsym(kTRUE); + PidObj->SetPLimit(ptlimit); + PidObj->SetSigma(0,2.); //TPC sigma, in three pT ranges + PidObj->SetSigma(1,1.); + PidObj->SetSigma(2,0.); + PidObj->SetSigma(3,2.); //TOF sigma, whole pT range + PidObj->SetTPC(kTRUE); + PidObj->SetTOF(kTRUE); + + for(Int_t ii=0; ii<4; ii++) { + PidObj->SetSigma(0,2.); + matchK[ii] = PidObj->MatchTPCTOF(trk[ii],1,3,kTRUE); //arguments: track, mode, particle#, compatibility allowed + PidObj->SetSigma(0,2.); + matchPi[ii] = PidObj->MatchTPCTOF(trk[ii],1,2,kTRUE); + } + + //Rho invariant mass under various hypotheses (to be matched with PID infos in order to accet the candidate) + Int_t d01rho03 = 0, d01rho23 = 0, d02rho01 = 0, d02rho12 = 0, d0bar1rho12 = 0, d0bar1rho23 = 0, d0bar2rho01 = 0, d0bar2rho03 = 0; + if(TMath::Abs(0.775 - d->InvMassRho(0,3))<0.1) {d01rho03 = 1; d0bar2rho03 = 1;} + if(TMath::Abs(0.775 - d->InvMassRho(2,3))<0.1) {d01rho23 = 1; d0bar1rho23 = 1;} + if(TMath::Abs(0.775 - d->InvMassRho(0,1))<0.1) {d02rho01 = 1; d0bar2rho01 = 1;} + if(TMath::Abs(0.775 - d->InvMassRho(1,2))<0.1) {d02rho12 = 1; d0bar1rho12 = 1;} + Int_t d01rho = 0, d02rho = 0, d0bar1rho = 0, d0bar2rho = 0; + if(d01rho03==1||d01rho23==1) d01rho = 1; + if(d02rho01==1||d02rho12==1) d02rho = 1; + if(d0bar1rho12==1||d0bar1rho23==1) d0bar1rho = 1; + if(d0bar2rho01==1||d0bar2rho03==1) d0bar2rho = 1; + + //This way because there could be multiple hypotheses accepted + if(d01rho==1 && (matchK[1]>=0 && matchPi[0]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp1 = 1; output = 1;} //d01 hyp + if(d02rho==1 && (matchK[3]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0)) {*hyp2 = 1; output = 1;} //d02 hyp + if(d0bar1rho==1 && (matchK[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp3 = 1; output = 1;} //d0bar1 hyp + if(d0bar2rho==1 && (matchK[2]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[3]>=0)) {*hyp4 = 1; output = 1;} //d0bar2 hyp + + return output; } //--------------------------------------------------------------------------- +Int_t AliRDHFCutsD0toKpipipi::D01Selected(TObject* obj,Int_t selectionLevel) { + // + // Apply selection + // + + if(!fCutsRD){ + cout<<"Cut matrix not inizialized. Exit..."<Pt()); + + Double_t mD0[2]; + Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass(); + + d->InvMassD0(mD0); + if(TMath::Abs(mD0[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1; + } + + return returnvalue; +} +//--------------------------------------------------------------------------- +Int_t AliRDHFCutsD0toKpipipi::D02Selected(TObject* obj,Int_t selectionLevel) { + // + // Apply selection + // + + if(!fCutsRD){ + cout<<"Cut matrix not inizialized. Exit..."<Pt()); + + Double_t mD0[2]; + Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass(); + + d->InvMassD0(mD0); + if(TMath::Abs(mD0[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1; + } + + return returnvalue; +} //--------------------------------------------------------------------------- +Int_t AliRDHFCutsD0toKpipipi::D0bar1Selected(TObject* obj,Int_t selectionLevel) { + // + // Apply selection + // + + if(!fCutsRD){ + cout<<"Cut matrix not inizialized. Exit..."<Pt()); + + Double_t mD0bar[2]; + Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass(); + + d->InvMassD0bar(mD0bar); + if(TMath::Abs(mD0bar[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1; + } + + return returnvalue; +} +//--------------------------------------------------------------------------- +Int_t AliRDHFCutsD0toKpipipi::D0bar2Selected(TObject* obj,Int_t selectionLevel) { + // + // Apply selection + // + + if(!fCutsRD){ + cout<<"Cut matrix not inizialized. Exit..."<Pt()); + + Double_t mD0bar[2]; + Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass(); + + d->InvMassD0bar(mD0bar); + if(TMath::Abs(mD0bar[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1; + } + + return returnvalue; +} +//--------------------------------------------------------------------------- +Bool_t AliRDHFCutsD0toKpipipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const +{ + // + // Checking if D0 is in fiducial acceptance region + // + + if(pt > 5.) { + // applying cut for pt > 5 GeV + AliDebug(4,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt)); + if (TMath::Abs(y) > 0.8){ + return kFALSE; + } + } else { + // appliying smooth cut for pt < 5 GeV + Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; + Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5; + AliDebug(4,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); + if (y < minFiducialY || y > maxFiducialY){ + return kFALSE; + } + } + + return kTRUE; +} + diff --git a/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.h b/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.h index 0b757076166..21229135717 100644 --- a/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.h +++ b/PWG3/vertexingHF/AliRDHFCutsD0toKpipipi.h @@ -3,15 +3,15 @@ /* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -/* $Id$ */ - //*********************************************************** // Class AliRDHFCutsD0toKpipipi // class for cuts on AOD reconstructed D0->Kpipipi -// Author: A.Dainese, andrea.dainese@pd.infn.it +// Author: A.Dainese, andrea.dainese@pd.infn.it +// F.Colamaria, fabio.colamaria@ba.infn.it //*********************************************************** -#include "AliRDHFCuts.h" +#include "AliRDHFCuts.h" +#include "AliAODRecoDecayHF4Prong.h" class AliRDHFCutsD0toKpipipi : public AliRDHFCuts { @@ -27,11 +27,18 @@ class AliRDHFCutsD0toKpipipi : public AliRDHFCuts virtual void GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters); using AliRDHFCuts::IsSelected; - virtual Int_t IsSelected(TObject* obj,Int_t selectionLevel); + virtual Int_t IsSelected(TObject* obj,Int_t selectionLevel); + virtual Int_t IsSelectedFromPID(AliAODRecoDecayHF4Prong *d, Int_t *hyp1, Int_t *hyp2, Int_t *hyp3, Int_t *hyp4); + virtual Int_t D01Selected(TObject* obj,Int_t selectionLevel); + virtual Int_t D02Selected(TObject* obj,Int_t selectionLevel); + virtual Int_t D0bar1Selected(TObject* obj,Int_t selectionLevel); + virtual Int_t D0bar2Selected(TObject* obj,Int_t selectionLevel); Float_t GetMassCut(Int_t iPtBin=0) const { return (GetCuts() ? fCutsRD[GetGlobalIndex(0,iPtBin)] : 1.e6);} - Float_t GetDCACut(Int_t iPtBin=0) const { return (GetCuts() ? fCutsRD[GetGlobalIndex(1,iPtBin)] : 1.e6);} - Bool_t GetUsePID(Int_t iPtBin=0) const { return (GetCuts() ? (Bool_t)(fCutsRD[GetGlobalIndex(8,iPtBin)]) : kFALSE);} + Float_t GetDCACut(Int_t iPtBin=0) const { return (GetCuts() ? fCutsRD[GetGlobalIndex(1,iPtBin)] : 1.e6);} + Bool_t GetUsePID(Int_t iPtBin=0) const { return (GetCuts() ? (Bool_t)(fCutsRD[GetGlobalIndex(8,iPtBin)]) : kFALSE);} + + virtual Bool_t IsInFiducialAcceptance(Double_t pt,Double_t y) const; protected: diff --git a/PWG3/vertexingHF/macros/AddTaskSelectHF4Prong.C b/PWG3/vertexingHF/macros/AddTaskSelectHF4Prong.C new file mode 100644 index 00000000000..0fc33e1a335 --- /dev/null +++ b/PWG3/vertexingHF/macros/AddTaskSelectHF4Prong.C @@ -0,0 +1,63 @@ +AliAnalysisTaskSESelectHF4Prong *AddTaskSelectHF4Prong() +{ + // + // Test macro for the AliAnalysisTaskSE for heavy-flavour selection + // and creation of a stand-alone AOD + // A.Dainese, andrea.dainese@lnl.infn.it + // F.Colamaria, fabio.colamaria@ba.infn.it + // + + // Get the pointer to the existing analysis manager via the static access method. + //============================================================================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTaskSelectHF4Prong", "No analysis manager to connect to."); + return NULL; + } + + // Output + AliAODHandler *aodHandler = new AliAODHandler(); + aodHandler->SetOutputFileName("AliAOD.VertexingHF.sa.root"); + aodHandler->SetCreateNonStandardAOD(); + mgr->SetOutputEventHandler(aodHandler); + + //Cuts loading + TFile* filecuts=new TFile("Charm4ProngCutsDef.root"); + if(!filecuts->IsOpen()){ + cout<<"Input file not found: exit"<Get("Charm4ProngCuts"); + RDHFCharm4Prong->SetName(Form("Charm4ProngCuts%d",1)); + + if(!RDHFCharm4Prong){ + cout<<"Specific AliRDHFCuts not found"<SetDebugLevel(2); + mgr->AddTask(hfTask); + + // Create containers for input/output mgr->ConnectInput(hfTask,0,mgr->GetCommonInputContainer()); + + AliAnalysisDataContainer *contHist = mgr->CreateContainer("histos_bin1",TList::Class(),AliAnalysisManager::kOutputContainer,"HistMassInvAndCuts.root"); + AliAnalysisDataContainer *contHist2 = mgr->CreateContainer("histos_bin2",TList::Class(),AliAnalysisManager::kOutputContainer,"HistMassInvAndCuts.root"); + AliAnalysisDataContainer *contHist3 = mgr->CreateContainer("histos_bin3",TList::Class(),AliAnalysisManager::kOutputContainer,"HistMassInvAndCuts.root"); + AliAnalysisDataContainer *contHist4 = mgr->CreateContainer("histos_bin4",TList::Class(),AliAnalysisManager::kOutputContainer,"HistMassInvAndCuts.root"); + AliAnalysisDataContainer *contHist5 = mgr->CreateContainer("histos_bin5",TList::Class(),AliAnalysisManager::kOutputContainer,"HistMassInvAndCuts.root"); + AliAnalysisDataContainer *contHistCuts = mgr->CreateContainer("histoscuts",TList::Class(),AliAnalysisManager::kOutputContainer,"HistMassInvAndCuts.root"); + + mgr->ConnectOutput(hfTask,0,mgr->GetCommonOutputContainer()); + mgr->ConnectOutput(hfTask,1,contHist); + mgr->ConnectOutput(hfTask,2,contHist2); + mgr->ConnectOutput(hfTask,3,contHist3); + mgr->ConnectOutput(hfTask,4,contHist4); + mgr->ConnectOutput(hfTask,5,contHist5); + mgr->ConnectOutput(hfTask,6,contHistCuts); + + return hfTask; +} diff --git a/PWG3/vertexingHF/macros/MakeCuts4Charm4Prong.C b/PWG3/vertexingHF/macros/MakeCuts4Charm4Prong.C new file mode 100644 index 00000000000..396cb759909 --- /dev/null +++ b/PWG3/vertexingHF/macros/MakeCuts4Charm4Prong.C @@ -0,0 +1,369 @@ +//gSystem->AddIncludePath("-I/$ALICE_ROOT/PWG3/vertexingHF"); + +#include +#include +//#include +#include +#include + +//Author: Chiara Bianchin, cbianchi@pd.infn.it, Fabio Colamaria, fabio.colamaria@ba.infn.it + +//macro to make a .root file which contains an AliRDHFCutsD0toKpipipi for AliAnalysisTaskSESelectHF4Prong task + +void MakeCuts4Charm4Prong(){ + +gSystem->Load("libTree.so"); +gSystem->Load("libGeom.so"); +gSystem->Load("libPhysics.so"); +gSystem->Load("libVMC.so"); +gSystem->Load("libMinuit.so"); +gSystem->Load("libSTEERBase.so"); +gSystem->Load("libESD.so"); +gSystem->Load("libAOD.so"); +gSystem->Load("libANALYSIS.so"); +gSystem->Load("libANALYSISalice.so"); +gSystem->Load("libCORRFW.so"); +gSystem->Load("libPWG3base.so"); +gSystem->Load("libPWG3vertexingHF.so"); +gSystem->Load("libPWG3muon.so"); + + AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi(); + RDHFCharm4Prong->SetName("Charm4ProngCuts"); + RDHFCharm4Prong->SetTitle("Cuts for D0 analysis"); + + AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts(); + esdTrackCuts->SetRequireSigmaToVertex(kFALSE); + //default + esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetRequireITSRefit(kTRUE); + esdTrackCuts->SetMinNClustersITS(4); // default is 5 + //esdTrackCuts->SetMinNClustersTPC(70); + esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, + AliESDtrackCuts::kAny); + // default is kBoth, otherwise kAny + esdTrackCuts->SetMinDCAToVertexXY(0.); + esdTrackCuts->SetPtRange(0.3,1.e10); + + RDHFCharm4Prong->AddTrackCuts(esdTrackCuts); + + + + const Int_t nvars=9; + + const Int_t nptbins=5; + Float_t* ptbins; + ptbins=new Float_t[nptbins+1]; + ptbins[0]=0.; + ptbins[1]=4.5; + ptbins[2]=6.; + ptbins[3]=8.; + ptbins[4]=12.; + ptbins[5]=16.; + + RDHFCharm4Prong->SetPtBins(nptbins+1,ptbins); + + + Float_t** rdcutsvalmine; + rdcutsvalmine=new Float_t*[nvars]; + for(Int_t iv=0;ivSetCuts(nvars,nptbins,rdcutsvalmine); + cout<<"This is the odject I'm going to save:"<PrintAll(); + TFile* fout=new TFile("Charm4ProngCutsDef_16GeV.root","recreate"); //set this!! + fout->cd(); + RDHFCharm4Prong->Write(); + fout->Close(); + +} + +void MakeCuts4Charm4ProngForMaxim(){ + +gSystem->Clear(); +gSystem->Load("libTree.so"); +gSystem->Load("libGeom.so"); +gSystem->Load("libPhysics.so"); +gSystem->Load("libVMC.so"); +gSystem->Load("libMinuit.so"); +gSystem->Load("libSTEERBase.so"); +gSystem->Load("libESD.so"); +gSystem->Load("libAOD.so"); +gSystem->Load("libANALYSIS.so"); +gSystem->Load("libANALYSISalice.so"); +gSystem->Load("libCORRFW.so"); +gSystem->Load("libPWG3base.so"); +gSystem->Load("libPWG3vertexingHF.so"); +gSystem->Load("libPWG3muon.so"); + + AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi(); + RDHFCharm4Prong->SetName("loosercuts"); + RDHFCharm4Prong->SetTitle("Cuts for significance maximization"); + + AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts(); + esdTrackCuts->SetRequireSigmaToVertex(kFALSE); + //default + esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetRequireITSRefit(kTRUE); + //esdTrackCuts->SetMinNClustersITS(4); + //esdTrackCuts->SetMinNClustersTPC(120); + + esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); + esdTrackCuts->SetMinDCAToVertexXY(0.); + esdTrackCuts->SetEtaRange(-0.8,0.8); + esdTrackCuts->SetPtRange(0.3,1.e10); + + RDHFCharm4Prong->AddTrackCuts(esdTrackCuts); + + const Int_t nvars=9; + + const Int_t nptbins=5; //change this when adding pt bins! + Float_t ptbins[nptbins+1]; + ptbins[0]=0.; + ptbins[1]=4.5; + ptbins[2]=6.; + ptbins[3]=8.; + ptbins[4]=12.; + ptbins[5]=16.; //o 25! + + RDHFCharm4Prong->SetPtBins(nptbins+1,ptbins); + + Float_t** rdcutsvalmine; + rdcutsvalmine=new Float_t*[nvars]; + for(Int_t iv=0;ivSetCuts(nvars,nptbins,cutsMatrixTransposeStand); + + Int_t nvarsforopt=RDHFCharm4Prong->GetNVarsForOpt(); + Int_t dim=5; //set this!! + Bool_t *boolforopt; + boolforopt=new Bool_t[nvars]; + if(dim>nvarsforopt){ + cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<GetVarsForOpt(); + }else{ + TString *names; + names=new TString[nvars]; + TString answer=""; + Int_t checktrue=0; + names=RDHFCharm4Prong->GetVarNames(); + for(Int_t i=0;i>answer; + if(answer=="y") { + boolforopt[i]=kTRUE; + checktrue++; + } + else boolforopt[i]=kFALSE; + } + if (checktrue!=dim) { + cout<<"Error! You set "<SetVarsForOpt(dim,boolforopt); + } + } + +// Tight values for variables to be maximized +// (DCA, Dist12, Dist3, Dist4, CosThetaPoint) +// Indexes: variable, ptbin + + Float_t tighterval[5][5]; + + //number of steps for each variable is set in the AddTask arguments (default=8) + + tighterval[0][0]=0.025; + tighterval[1][0]=0.09; + tighterval[2][0]=0.09; + tighterval[3][0]=0.06; + tighterval[4][0]=0.955; + + tighterval[0][1]=0.015; + tighterval[1][1]=0.12; + tighterval[2][1]=0.10; + tighterval[3][1]=0.085; + tighterval[4][1]=0.99; + + tighterval[0][2]=0.025; + tighterval[1][2]=0.0975; + tighterval[2][2]=0.09; + tighterval[3][2]=0.08; + tighterval[4][2]=0.99; + + tighterval[0][3]=0.015; + tighterval[1][3]=0.0975; + tighterval[2][3]=0.07; + tighterval[3][3]=0.075; + tighterval[4][3]=0.9825; + + tighterval[0][4]=0.015; + tighterval[1][4]=0.0975; + tighterval[2][4]=0.07; + tighterval[3][4]=0.075; + tighterval[4][4]=0.9825; + + TString name=""; + Int_t arrdim=dim*nptbins; + cout<<"Will save "<"<",arrdim); + for(Int_t ival=0;ival(name.Data(),tighterval[ival][jpt]); + } + } + + Bool_t flagPID=kFALSE; + RDHFCharm4Prong->SetUsePID(flagPID); + + RDHFCharm4Prong->PrintAll(); + printf("Use PID? %s\n",flagPID ? "yes" : "no"); + +/* + //pid settings + AliAODPidHF* pidObj=new AliAODPidHF(); + //pidObj->SetName("pid4D0"); + Int_t mode=1; + const Int_t nlims=2; + Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c] + Bool_t compat=kTRUE; //effective only for this mode + Bool_t asym=kTRUE; + Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella + pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC + pidObj->SetMatch(mode); + pidObj->SetPLimit(plims,nlims); + pidObj->SetSigma(sigmas); + pidObj->SetCompat(compat); + pidObj->SetTPC(kTRUE); + pidObj->SetTOF(kTRUE); + RDHFCharm4Prong->SetPidHF(pidObj); + + RDHFCharm4Prong->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF +*/ + + //activate pileup rejection (for pp) + RDHFCharm4Prong->SetOptPileup(AliRDHFCuts::kRejectPileupEvent); + + //Do not recalculate the vertex + RDHFCharm4Prong->SetRemoveDaughtersFromPrim(kTRUE); //activate for pp + +/* Per il pp levo la centralità! +//Metto kCentOff! Per il Pb c'era kCentV0M + TString cent=""; + //centrality selection (Pb-Pb) + Float_t minc=20,maxc=80; + RDHFCharm4Prong->SetMinCentrality(minc); + RDHFCharm4Prong->SetMaxCentrality(maxc); + cent=Form("%.0f%.0f",minc,maxc); + RDHFCharm4Prong->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid +*/ + + //temporary +// RDHFCharm4Prong->SetFixRefs(); + + TFile* fout=new TFile("Charm4ProngCutsForMaxim.root","recreate"); + fout->cd(); + RDHFCharm4Prong->Write(); + max.Write(); + fout->Close(); +} -- 2.43.0