--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2009, 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 *
+ * appeuear 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. *
+ **************************************************************************/
+
+/* $Id: AliAnalysisTaskSELc2V0bachelorTMVA.cxx 62231 2013-04-29 09:47:26Z fprino $ */
+
+//
+//
+// Base class for Lc2V0 Analysis to be used with TMVA
+//
+//
+
+//------------------------------------------------------------------------------------------
+//
+// Author: C. Zampolli (from AliAnalysisTaskSELc2V0bachelor class by A.De Caro, P. Pagano)
+//
+//------------------------------------------------------------------------------------------
+
+#include <TSystem.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include <TH1F.h>
+#include <TH1I.h>
+#include <TH2F.h>
+#include <TTree.h>
+#include "TROOT.h"
+#include <TDatabasePDG.h>
+#include <AliAnalysisDataSlot.h>
+#include <AliAnalysisDataContainer.h>
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliAODMCHeader.h"
+#include "AliAODHandler.h"
+#include "AliLog.h"
+#include "AliAODVertex.h"
+#include "AliAODRecoDecay.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliAnalysisVertexingHF.h"
+#include "AliESDtrack.h"
+#include "AliAODTrack.h"
+#include "AliAODv0.h"
+#include "AliAODMCParticle.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskSELc2V0bachelorTMVA.h"
+#include "AliNormalizationCounter.h"
+#include "AliAODPidHF.h"
+#include "AliPIDResponse.h"
+#include "AliPIDCombined.h"
+#include "AliTOFPIDResponse.h"
+#include "AliInputEventHandler.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliKFParticle.h"
+#include "AliKFVertex.h"
+#include "AliVertexingHFUtils.h"
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTaskSELc2V0bachelorTMVA)
+
+//__________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA():
+AliAnalysisTaskSE(),
+ fUseMCInfo(kFALSE),
+ fOutput(0),
+ fCEvents(0),
+ fPIDResponse(0),
+ fPIDCombined(0),
+ fIsK0sAnalysis(kFALSE),
+ fCounter(0),
+ fAnalCuts(0),
+ fListCuts(0),
+ fUseOnTheFlyV0(kFALSE),
+ fIsEventSelected(kFALSE),
+ fVariablesTreeSgn(0),
+ fVariablesTreeBkg(0),
+ fCandidateVariables(),
+ fIspA(kFALSE),
+ fHistoEvents(0),
+ fHistoLc(0),
+ fHistoLcOnTheFly(0),
+ fFillOnlySgn(kFALSE),
+fHistoLcBeforeCuts(0),
+fHistoFiducialAcceptance(0),
+fHistoCodesSgn(0),
+fHistoCodesBkg(0),
+ fHistoLcpKpiBeforeCuts(0),
+ fVtx1(0),
+
+ fHistoDistanceLcToPrimVtx(0),
+ fHistoDistanceV0ToPrimVtx(0),
+ fHistoDistanceV0ToLc(0),
+
+ fHistoDistanceLcToPrimVtxSgn(0),
+ fHistoDistanceV0ToPrimVtxSgn(0),
+ fHistoDistanceV0ToLcSgn(0),
+
+ fHistoVtxLcResidualToPrimVtx(0),
+ fHistoVtxV0ResidualToPrimVtx(0),
+ fHistoVtxV0ResidualToLc(0),
+
+ fHistoMassV0All(0),
+ fHistoDecayLengthV0All(0),
+ fHistoLifeTimeV0All(0),
+
+ fHistoMassV0True(0),
+ fHistoDecayLengthV0True(0),
+ fHistoLifeTimeV0True(0),
+
+ fHistoMassV0TrueFromAOD(0),
+
+ fHistoMassV0TrueK0S(0),
+ fHistoDecayLengthV0TrueK0S(0),
+ fHistoLifeTimeV0TrueK0S(0),
+
+ fHistoMassV0TrueK0SFromAOD(0),
+
+ fHistoMassLcAll(0),
+ fHistoDecayLengthLcAll(0),
+ fHistoLifeTimeLcAll(0),
+
+ fHistoMassLcTrue(0),
+ fHistoDecayLengthLcTrue(0),
+ fHistoLifeTimeLcTrue(0),
+
+ fHistoMassLcTrueFromAOD(0),
+
+ fHistoMassV0fromLcAll(0),
+ fHistoDecayLengthV0fromLcAll(0),
+ fHistoLifeTimeV0fromLcAll(0),
+
+ fHistoMassV0fromLcTrue(0),
+ fHistoDecayLengthV0fromLcTrue(0),
+ fHistoLifeTimeV0fromLcTrue(0),
+
+ fHistoMassLcSgn(0),
+ fHistoMassLcSgnFromAOD(0),
+ fHistoDecayLengthLcSgn(0),
+ fHistoLifeTimeLcSgn(0),
+
+ fHistoMassV0fromLcSgn(0),
+ fHistoDecayLengthV0fromLcSgn(0),
+ fHistoLifeTimeV0fromLcSgn(0),
+
+ fHistoKF(0),
+ fHistoKFV0(0),
+ fHistoKFLc(0),
+
+ fHistoMassKFV0(0),
+ fHistoDecayLengthKFV0(0),
+ fHistoLifeTimeKFV0(0),
+
+ fHistoMassKFLc(0),
+ fHistoDecayLengthKFLc(0),
+ fHistoLifeTimeKFLc(0),
+
+ fHistoArmenterosPodolanskiV0KF(0),
+ fHistoArmenterosPodolanskiV0KFSgn(0),
+ fHistoArmenterosPodolanskiV0AOD(0),
+ fHistoArmenterosPodolanskiV0AODSgn(0),
+
+ fOutputKF(0),
+ fmcLabelLc(-1),
+ ftopoConstraint(kTRUE),
+ fCallKFVertexing(kFALSE),
+ fKeepingOnlyHIJINGBkg(kFALSE),
+ fUtils(0),
+ fHistoBackground(0)
+
+{
+ //
+ // Default ctor
+ //
+}
+//___________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
+ AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
+ AliAnalysisTaskSE(name),
+ fUseMCInfo(kFALSE),
+ fOutput(0),
+ fCEvents(0),
+ fPIDResponse(0),
+ fPIDCombined(0),
+ fIsK0sAnalysis(kFALSE),
+ fCounter(0),
+ fAnalCuts(analCuts),
+ fListCuts(0),
+ fUseOnTheFlyV0(useOnTheFly),
+ fIsEventSelected(kFALSE),
+ fVariablesTreeSgn(0),
+ fVariablesTreeBkg(0),
+ fCandidateVariables(),
+ fIspA(kFALSE),
+ fHistoEvents(0),
+ fHistoLc(0),
+ fHistoLcOnTheFly(0),
+ fFillOnlySgn(kFALSE),
+ fHistoLcBeforeCuts(0),
+ fHistoFiducialAcceptance(0),
+ fHistoCodesSgn(0),
+ fHistoCodesBkg(0),
+ fHistoLcpKpiBeforeCuts(0),
+ fVtx1(0),
+
+ fHistoDistanceLcToPrimVtx(0),
+ fHistoDistanceV0ToPrimVtx(0),
+ fHistoDistanceV0ToLc(0),
+
+ fHistoDistanceLcToPrimVtxSgn(0),
+ fHistoDistanceV0ToPrimVtxSgn(0),
+ fHistoDistanceV0ToLcSgn(0),
+
+ fHistoVtxLcResidualToPrimVtx(0),
+ fHistoVtxV0ResidualToPrimVtx(0),
+ fHistoVtxV0ResidualToLc(0),
+
+ fHistoMassV0All(0),
+ fHistoDecayLengthV0All(0),
+ fHistoLifeTimeV0All(0),
+
+ fHistoMassV0True(0),
+ fHistoDecayLengthV0True(0),
+ fHistoLifeTimeV0True(0),
+
+ fHistoMassV0TrueFromAOD(0),
+
+ fHistoMassV0TrueK0S(0),
+ fHistoDecayLengthV0TrueK0S(0),
+ fHistoLifeTimeV0TrueK0S(0),
+
+ fHistoMassV0TrueK0SFromAOD(0),
+
+ fHistoMassLcAll(0),
+ fHistoDecayLengthLcAll(0),
+ fHistoLifeTimeLcAll(0),
+
+ fHistoMassLcTrue(0),
+ fHistoDecayLengthLcTrue(0),
+ fHistoLifeTimeLcTrue(0),
+
+ fHistoMassLcTrueFromAOD(0),
+
+ fHistoMassV0fromLcAll(0),
+ fHistoDecayLengthV0fromLcAll(0),
+ fHistoLifeTimeV0fromLcAll(0),
+
+ fHistoMassV0fromLcTrue(0),
+ fHistoDecayLengthV0fromLcTrue(0),
+ fHistoLifeTimeV0fromLcTrue(0),
+
+ fHistoMassLcSgn(0),
+ fHistoMassLcSgnFromAOD(0),
+ fHistoDecayLengthLcSgn(0),
+ fHistoLifeTimeLcSgn(0),
+
+ fHistoMassV0fromLcSgn(0),
+ fHistoDecayLengthV0fromLcSgn(0),
+ fHistoLifeTimeV0fromLcSgn(0),
+
+ fHistoKF(0),
+ fHistoKFV0(0),
+ fHistoKFLc(0),
+
+ fHistoMassKFV0(0),
+ fHistoDecayLengthKFV0(0),
+ fHistoLifeTimeKFV0(0),
+
+ fHistoMassKFLc(0),
+ fHistoDecayLengthKFLc(0),
+ fHistoLifeTimeKFLc(0),
+
+ fHistoArmenterosPodolanskiV0KF(0),
+ fHistoArmenterosPodolanskiV0KFSgn(0),
+ fHistoArmenterosPodolanskiV0AOD(0),
+ fHistoArmenterosPodolanskiV0AODSgn(0),
+
+ fOutputKF(0),
+ fmcLabelLc(-1),
+ ftopoConstraint(kTRUE),
+ fCallKFVertexing(kFALSE),
+ fKeepingOnlyHIJINGBkg(kFALSE),
+ fUtils(0),
+ fHistoBackground(0)
+
+
+{
+ //
+ // Constructor. Initialization of Inputs and Outputs
+ //
+ Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
+
+ DefineOutput(1, TList::Class()); // Tree signal + Tree Bkg + histoEvents
+ DefineOutput(2, AliNormalizationCounter::Class()); // normalization counter object
+ DefineOutput(3, TList::Class()); // Cuts
+ DefineOutput(4, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
+ DefineOutput(5, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
+ DefineOutput(6, TList::Class()); // Tree signal + Tree Bkg + histoEvents
+
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
+ //
+ // destructor
+ //
+ Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
+
+ if (fOutput) {
+ delete fOutput;
+ fOutput = 0;
+ }
+
+ if (fPIDResponse) {
+ delete fPIDResponse;
+ }
+
+ if (fPIDCombined) {
+ delete fPIDCombined;
+ }
+
+ if (fCounter) {
+ delete fCounter;
+ fCounter = 0;
+ }
+
+ if (fAnalCuts) {
+ delete fAnalCuts;
+ fAnalCuts = 0;
+ }
+
+ if (fListCuts) {
+ delete fListCuts;
+ fListCuts = 0;
+ }
+
+ if(fVariablesTreeSgn){
+ delete fVariablesTreeSgn;
+ fVariablesTreeSgn = 0;
+ }
+
+ if(fVariablesTreeBkg){
+ delete fVariablesTreeBkg;
+ fVariablesTreeBkg = 0;
+ }
+
+ if (fOutputKF) {
+ delete fOutputKF;
+ fOutputKF = 0;
+ }
+
+ if (fUtils) {
+ delete fUtils;
+ fUtils = 0;
+ }
+
+}
+//_________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
+ //
+ // Initialization
+ //
+
+ fIsEventSelected=kFALSE;
+
+ if (fDebug > 1) AliInfo("Init");
+
+ fListCuts = new TList();
+ fListCuts->SetOwner();
+ fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
+ PostData(3,fListCuts);
+
+ if (fUseMCInfo && fKeepingOnlyHIJINGBkg) fUtils = new AliVertexingHFUtils();
+
+ return;
+}
+
+//________________________________________ terminate ___________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
+{
+ // The Terminate() function is the last function to be called during
+ // a query. It always runs on the client, it can be used to present
+ // the results graphically or save the results to file.
+
+ //AliInfo("Terminate","");
+ AliAnalysisTaskSE::Terminate();
+
+ fOutput = dynamic_cast<TList*> (GetOutputData(1));
+ if (!fOutput) {
+ AliError("fOutput not available");
+ return;
+ }
+ fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
+ if (!fOutputKF) {
+ AliError("fOutputKF not available");
+ return;
+ }
+
+ return;
+}
+
+//___________________________________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
+ // output
+ AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
+
+ //slot #1
+ //OpenFile(1);
+ fOutput = new TList();
+ fOutput->SetOwner();
+ fOutput->SetName("listTrees");
+
+ // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
+ const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
+ fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
+ fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
+ Int_t nVar = 66;
+ fCandidateVariables = new Float_t [nVar];
+ TString * fCandidateVariableNames = new TString[nVar];
+ fCandidateVariableNames[0]="massLc2K0Sp";
+ fCandidateVariableNames[1]="massLc2Lambdapi";
+ fCandidateVariableNames[2]="massK0S";
+ fCandidateVariableNames[3]="massLambda";
+ fCandidateVariableNames[4]="massLambdaBar";
+ fCandidateVariableNames[5]="cosPAK0S";
+ fCandidateVariableNames[6]="dcaV0";
+ fCandidateVariableNames[7]="tImpParBach";
+ fCandidateVariableNames[8]="tImpParV0";
+ fCandidateVariableNames[9]="nSigmaTPCpr";
+ fCandidateVariableNames[10]="nSigmaTPCpi";
+ fCandidateVariableNames[11]="nSigmaTPCka";
+ fCandidateVariableNames[12]="nSigmaTOFpr";
+ fCandidateVariableNames[13]="nSigmaTOFpi";
+ fCandidateVariableNames[14]="nSigmaTOFka";
+ fCandidateVariableNames[15]="bachelorPt";
+ fCandidateVariableNames[16]="V0positivePt";
+ fCandidateVariableNames[17]="V0negativePt";
+ fCandidateVariableNames[18]="dcaV0pos";
+ fCandidateVariableNames[19]="dcaV0neg";
+ fCandidateVariableNames[20]="v0Pt";
+ fCandidateVariableNames[21]="massGamma";
+ fCandidateVariableNames[22]="LcPt";
+ fCandidateVariableNames[23]="combinedProtonProb";
+ fCandidateVariableNames[24]="LcEta";
+ fCandidateVariableNames[25]="V0positiveEta";
+ fCandidateVariableNames[26]="V0negativeEta";
+ fCandidateVariableNames[27]="combinedPionProb";
+ fCandidateVariableNames[28]="combinedKaonProb";
+ fCandidateVariableNames[29]="bachelorEta";
+ fCandidateVariableNames[30]="LcP";
+ fCandidateVariableNames[31]="bachelorP";
+ fCandidateVariableNames[32]="v0P";
+ fCandidateVariableNames[33]="V0positiveP";
+ fCandidateVariableNames[34]="V0negativeP";
+ fCandidateVariableNames[35]="LcY";
+ fCandidateVariableNames[36]="v0Y";
+ fCandidateVariableNames[37]="bachelorY";
+ fCandidateVariableNames[38]="V0positiveY";
+ fCandidateVariableNames[39]="V0negativeY";
+ fCandidateVariableNames[40]="v0Eta";
+ fCandidateVariableNames[41]="DecayLengthLc";
+ fCandidateVariableNames[42]="DecayLengthK0S";
+ fCandidateVariableNames[43]="CtLc";
+ fCandidateVariableNames[44]="CtK0S";
+ fCandidateVariableNames[45]="bachCode";
+ fCandidateVariableNames[46]="k0SCode";
+
+ fCandidateVariableNames[47]="V0KFmass";
+ fCandidateVariableNames[48]="V0KFdecayLength";
+ fCandidateVariableNames[49]="V0KFlifeTime";
+
+ fCandidateVariableNames[50]="V0KFmassErr";
+ fCandidateVariableNames[51]="V0KFdecayTimeErr";
+ fCandidateVariableNames[52]="V0KFlifeTimeErr";
+
+ fCandidateVariableNames[53]="LcKFmass";
+ fCandidateVariableNames[54]="LcKFdecayLength";
+ fCandidateVariableNames[55]="LcKFlifeTime";
+
+ fCandidateVariableNames[56]="LcKFmassErr";
+ fCandidateVariableNames[57]="LcKFdecayTimeErr";
+ fCandidateVariableNames[58]="LcKFlifeTimeErr";
+
+ fCandidateVariableNames[59]="LcKFDistToPrimVtx";
+ fCandidateVariableNames[60]="V0KFDistToPrimVtx";
+ fCandidateVariableNames[61]="V0KFDistToLc";
+ fCandidateVariableNames[62]="alphaArmKF";
+ fCandidateVariableNames[63]="ptArmKF";
+ fCandidateVariableNames[64]="alphaArm";
+ fCandidateVariableNames[65]="ptArm";
+
+ for(Int_t ivar=0; ivar<nVar; ivar++){
+ fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
+ fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
+ }
+
+ fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
+ TString labelEv[2] = {"NotSelected", "Selected"};
+ for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
+ fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
+ }
+
+ fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
+
+ fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
+ TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
+ for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
+ fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
+ }
+
+ fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
+ TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
+ for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
+ fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
+ fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
+ }
+
+ fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
+ TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
+ for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
+ fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
+ }
+
+ fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
+ fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
+
+ TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
+ "kBachDifferentLambdaMother", "kBachCorrectLambdaMother"};
+ TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
+ "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
+
+ for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
+ fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
+ fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
+ }
+ for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
+ fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
+ fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
+ }
+
+ fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
+ for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
+ fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
+ }
+
+ fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 2, -0.5, 1.5);
+ TString labelBkg[2] = {"Injected", "Non-injected"};
+ for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
+ fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
+ }
+
+ //fOutput->Add(fVariablesTreeSgn);
+ //fOutput->Add(fVariablesTreeBkg);
+
+ fOutput->Add(fHistoEvents);
+ fOutput->Add(fHistoLc);
+ fOutput->Add(fHistoLcOnTheFly);
+ fOutput->Add(fHistoLcBeforeCuts);
+ fOutput->Add(fHistoFiducialAcceptance);
+ fOutput->Add(fHistoCodesSgn);
+ fOutput->Add(fHistoCodesBkg);
+ fOutput->Add(fHistoLcpKpiBeforeCuts);
+ fOutput->Add(fHistoBackground);
+
+ PostData(1, fOutput);
+ PostData(4, fVariablesTreeSgn);
+ PostData(5, fVariablesTreeBkg);
+ AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+ fPIDResponse = inputHandler->GetPIDResponse();
+
+ if (fAnalCuts->GetIsUsePID()){
+ /*
+ fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
+ fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
+ fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
+ fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
+ fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
+ fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
+ */
+ fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
+ }
+
+ // Setting properties of PID
+ fPIDCombined=new AliPIDCombined;
+ fPIDCombined->SetDefaultTPCPriors();
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
+ //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
+
+ //Printf("Here ok!");
+
+ fCounter = new AliNormalizationCounter("NormalizationCounter");
+ fCounter->Init();
+ PostData(2, fCounter);
+
+ // Histograms from KF
+
+ fHistoDistanceLcToPrimVtx = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
+ fHistoDistanceV0ToPrimVtx = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
+ fHistoDistanceV0ToLc = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
+
+ fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
+ fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
+ fHistoDistanceV0ToLcSgn = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
+
+ fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
+ fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
+ fHistoVtxV0ResidualToLc = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
+
+ fHistoMassV0All = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
+ fHistoDecayLengthV0All = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
+ fHistoLifeTimeV0All = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
+
+ fHistoMassV0True = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
+ fHistoDecayLengthV0True = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
+ fHistoLifeTimeV0True = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
+
+ fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
+
+ fHistoMassV0TrueK0S = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
+ fHistoDecayLengthV0TrueK0S = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
+ fHistoLifeTimeV0TrueK0S = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
+
+ fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
+
+ fHistoMassLcAll = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
+ fHistoDecayLengthLcAll = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
+ fHistoLifeTimeLcAll = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
+
+ fHistoMassLcTrue = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
+ fHistoDecayLengthLcTrue = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
+ fHistoLifeTimeLcTrue = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
+
+ fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
+
+ fHistoMassV0fromLcAll = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
+ fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
+ fHistoLifeTimeV0fromLcAll = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
+
+ fHistoMassV0fromLcTrue = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
+ fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
+ fHistoLifeTimeV0fromLcTrue = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
+
+ fHistoMassLcSgn = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
+ fHistoMassLcSgnFromAOD = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
+ fHistoDecayLengthLcSgn = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
+ fHistoLifeTimeLcSgn = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
+
+ fHistoMassV0fromLcSgn = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
+ fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
+ fHistoLifeTimeV0fromLcSgn = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
+
+ fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
+ fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
+ fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
+ TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
+ "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
+ "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
+ "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
+
+ for (Int_t ibin = 1; ibin <=16; ibin++){
+ fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+ fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+ fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+ fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
+ }
+
+ fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
+ fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
+ fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
+
+ fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
+ fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
+ fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
+
+ fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
+ fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
+ fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
+ fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
+
+ fOutputKF = new TList();
+ fOutputKF->SetOwner();
+ fOutputKF->SetName("listHistoKF");
+
+ fOutputKF->Add(fHistoDistanceLcToPrimVtx);
+ fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
+ fOutputKF->Add(fHistoDistanceV0ToLc);
+
+ fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
+ fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
+ fOutputKF->Add(fHistoDistanceV0ToLcSgn);
+
+ fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
+ fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
+ fOutputKF->Add(fHistoVtxV0ResidualToLc);
+
+ fOutputKF->Add(fHistoMassV0All);
+ fOutputKF->Add(fHistoDecayLengthV0All);
+ fOutputKF->Add(fHistoLifeTimeV0All);
+
+ fOutputKF->Add(fHistoMassV0True);
+ fOutputKF->Add(fHistoDecayLengthV0True);
+ fOutputKF->Add(fHistoLifeTimeV0True);
+
+ fOutputKF->Add(fHistoMassV0TrueFromAOD);
+
+ fOutputKF->Add(fHistoMassV0TrueK0S);
+ fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
+ fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
+
+ fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
+
+ fOutputKF->Add(fHistoMassLcAll);
+ fOutputKF->Add(fHistoDecayLengthLcAll);
+ fOutputKF->Add(fHistoLifeTimeLcAll);
+
+ fOutputKF->Add(fHistoMassLcTrue);
+ fOutputKF->Add(fHistoDecayLengthLcTrue);
+ fOutputKF->Add(fHistoLifeTimeLcTrue);
+
+ fOutputKF->Add(fHistoMassLcTrueFromAOD);
+
+ fOutputKF->Add(fHistoMassV0fromLcAll);
+ fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
+ fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
+
+ fOutputKF->Add(fHistoMassV0fromLcTrue);
+ fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
+ fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
+
+ fOutputKF->Add(fHistoMassLcSgn);
+ fOutputKF->Add(fHistoDecayLengthLcSgn);
+ fOutputKF->Add(fHistoLifeTimeLcSgn);
+
+ fOutputKF->Add(fHistoMassV0fromLcSgn);
+ fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
+ fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
+
+ fOutputKF->Add(fHistoKF);
+ fOutputKF->Add(fHistoKFV0);
+ fOutputKF->Add(fHistoKFLc);
+
+ fOutputKF->Add(fHistoMassKFV0);
+ fOutputKF->Add(fHistoDecayLengthKFV0);
+ fOutputKF->Add(fHistoLifeTimeKFV0);
+
+ fOutputKF->Add(fHistoMassKFLc);
+ fOutputKF->Add(fHistoDecayLengthKFLc);
+ fOutputKF->Add(fHistoLifeTimeKFLc);
+
+ fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
+ fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
+ fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
+ fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
+
+ PostData(6, fOutputKF);
+
+ //Printf("Here ok 1!");
+
+ return;
+}
+
+//_________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
+{
+ // user exec
+ if (!fInputEvent) {
+ AliError("NO EVENT FOUND!");
+ return;
+ }
+
+ AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
+ TClonesArray *arrayLctopKos=0;
+
+ TClonesArray *array3Prong = 0;
+
+ if (!aodEvent && 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.
+ aodEvent = dynamic_cast<AliAODEvent*> (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();
+ arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
+
+ array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
+ }
+ } else {
+ arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
+
+ array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
+ }
+
+ if ( !fUseMCInfo && fIspA) {
+ fAnalCuts->SetTriggerClass("");
+ fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
+ }
+
+ fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
+
+ // AOD primary vertex
+ fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
+ if (!fVtx1) return;
+ if (fVtx1->GetNContributors()<1) return;
+
+ fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
+
+ if ( !fIsEventSelected ) {
+ fHistoEvents->Fill(0);
+ return; // don't take into account not selected events
+ }
+ fHistoEvents->Fill(1);
+
+ // Setting magnetic field for KF vertexing
+ AliKFParticle::SetField(aodEvent->GetMagneticField());
+
+ // mc analysis
+ TClonesArray *mcArray = 0;
+ AliAODMCHeader *mcHeader=0;
+
+ if (fUseMCInfo) {
+ // MC array need for maching
+ mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!mcArray) {
+ AliError("Could not find Monte-Carlo in AOD");
+ return;
+ }
+ // load MC header
+ mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+ if (!mcHeader) {
+ AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
+ return;
+ }
+ }
+
+ Int_t nSelectedAnal = 0;
+ if (fIsK0sAnalysis) {
+ MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
+ nSelectedAnal, fAnalCuts,
+ array3Prong, mcHeader);
+ }
+ fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
+
+ PostData(1, fOutput);
+ PostData(2, fCounter);
+ PostData(4, fVariablesTreeSgn);
+ PostData(5, fVariablesTreeBkg);
+ PostData(6, fOutputKF);
+
+}
+
+//-------------------------------------------------------------------------------
+void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
+ TClonesArray *mcArray,
+ Int_t &nSelectedAnal,
+ AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
+ AliAODMCHeader* aodheader){
+ //Lc prong needed to MatchToMC method
+
+ Int_t pdgCand = 4122;
+ Int_t pdgDgLctoV0bachelor[2]={2212, 310};
+ Int_t pdgDgV0toDaughters[2]={211, 211};
+
+ Int_t pdgDgLctopKpi[3]={2212, 321, 211};
+
+ // loop to search for candidates Lc->p+K+pi
+ Int_t n3Prong = array3Prong->GetEntriesFast();
+ Int_t nCascades= arrayLctopKos->GetEntriesFast();
+
+ //AliInfo(Form("\n\n\n\n3 prong candidates = %d, ncascades = %d \n\n\n\n\n", n3Prong, nCascades));
+ for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
+ AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
+ //Filling a control histogram with no cuts
+ if (fUseMCInfo) {
+
+ // find associated MC particle for Lc -> p+K+pi
+ Int_t mcLabel = d->MatchToMC(4122, mcArray, 3, pdgDgLctopKpi);
+ //Int_t mcLabelTemp = d->MatchToMC(4122, mcArray);
+ //Printf("mcLabel = %d, mcLabelTemp = %d", mcLabel, mcLabelTemp);
+ if (mcLabel >= 0) {
+
+ AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
+ if(partLcpKpi){
+ Int_t pdgCode = partLcpKpi->GetPdgCode();
+ AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
+ //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
+ fHistoLcpKpiBeforeCuts->Fill(1);
+
+ }
+ }
+ else {
+ //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
+ fHistoLcpKpiBeforeCuts->Fill(0);
+ }
+ }
+ }
+
+ // loop over cascades to search for candidates Lc->p+K0S
+
+ for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
+
+ // Lc candidates and K0s from Lc
+ AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
+ if (!lcK0spr) {
+ AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
+ continue;
+ }
+
+ //Filling a control histogram with no cuts
+ if (fUseMCInfo) {
+
+ Int_t pdgCode=-2;
+
+ // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
+ fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
+ if (fmcLabelLc>=0) {
+ AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
+
+ AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
+ if(partLc){
+ pdgCode = partLc->GetPdgCode();
+ if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
+ pdgCode = TMath::Abs(pdgCode);
+ fHistoLcBeforeCuts->Fill(1);
+
+ }
+ }
+ else {
+ fHistoLcBeforeCuts->Fill(0);
+ pdgCode=-1;
+ }
+ }
+
+ //if (!lcK0spr->GetSecondaryVtx()) {
+ // AliInfo("No secondary vertex");
+ //continue;
+ //}
+
+ if (lcK0spr->GetNDaughters()!=2) {
+ AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
+ continue;
+ }
+
+ AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
+ AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
+ if (!v0part || !bachPart) {
+ AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0s));
+ continue;
+ }
+
+
+ if (!v0part->GetSecondaryVtx()) {
+ AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
+ continue;
+ }
+
+ if (v0part->GetNDaughters()!=2) {
+ AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
+ continue;
+ }
+
+ AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
+ AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
+ if (!v0Neg || !v0Neg) {
+ AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
+ continue;
+ }
+
+
+ if (v0Pos->Charge() == v0Neg->Charge()) {
+ AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
+ continue;
+ }
+
+ Int_t isLc = 0;
+
+ if (fUseMCInfo) {
+
+ Int_t pdgCode=-2;
+
+ // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
+ Int_t mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
+ if (mcLabel>=0) {
+ AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
+
+ AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
+ if(partLc){
+ pdgCode = partLc->GetPdgCode();
+ if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
+ pdgCode = TMath::Abs(pdgCode);
+ isLc = 1;
+ fHistoLc->Fill(1);
+
+ }
+ }
+ else {
+ fHistoLc->Fill(0);
+ pdgCode=-1;
+ }
+ }
+ AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
+ AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
+ if (!isLc) {
+ if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
+ continue;
+ }
+ else { // checking if we want to fill the background
+ if (fKeepingOnlyHIJINGBkg){
+ // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
+ //Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(part, aodheader, mcArray);
+ Bool_t isCandidateInjected = CheckInjection(lcK0spr, aodheader, mcArray);
+ if (!isCandidateInjected){
+ AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
+ fHistoBackground->Fill(1);
+ }
+ else {
+ AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
+ fHistoBackground->Fill(0);
+ continue;
+ }
+ }
+ }
+ }
+
+ FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray);
+
+ }
+
+ return;
+
+}
+//________________________________________________________________________
+void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
+ Int_t isLc,
+ Int_t &nSelectedAnal,
+ AliRDHFCutsLctoV0 *cutsAnal,
+ TClonesArray *mcArray){
+ //
+ // Fill histos for Lc -> K0S+proton
+ //
+
+ /*
+ if (!part->GetOwnPrimaryVtx()) {
+ //Printf("No primary vertex for Lc found!!");
+ part->SetOwnPrimaryVtx(fVtx1);
+ }
+ else {
+ //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
+ }
+ */
+ Double_t invmassLc = part->InvMassLctoK0sP();
+
+ AliAODv0 * v0part = part->Getv0();
+ Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
+ if (onFlyV0){ // on-the-fly V0
+ if (isLc) { // Lc
+ fHistoLcOnTheFly->Fill(2.);
+ }
+ else { // not Lc
+ fHistoLcOnTheFly->Fill(0.);
+ }
+ }
+ else { // offline V0
+ if (isLc) { // Lc
+ fHistoLcOnTheFly->Fill(3.);
+ }
+ else { // not Lc
+ fHistoLcOnTheFly->Fill(1.);
+ }
+ }
+
+ Double_t dcaV0 = v0part->GetDCA();
+ Double_t invmassK0s = v0part->MassK0Short();
+
+ if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
+ if (isLc) {
+ fHistoFiducialAcceptance->Fill(3.);
+ }
+ else {
+ fHistoFiducialAcceptance->Fill(1.);
+ }
+ }
+ else {
+ if (isLc) {
+ fHistoFiducialAcceptance->Fill(2.);
+ }
+ else {
+ fHistoFiducialAcceptance->Fill(0.);
+ }
+ }
+
+ Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
+
+ if (isInV0window == 0) {
+ AliDebug(2, "No: The candidate has NOT passed the mass cuts!");
+ if (isLc) Printf("SIGNAL candidate rejected");
+ return;
+ }
+ else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
+
+ // Printf("ciccio!!!!!!!");
+
+ Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
+
+ if (!isInCascadeWindow) {
+ AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
+ if (isLc) Printf("SIGNAL candidate rejected");
+ return;
+ }
+ else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
+
+ Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
+ if (!isCandidateSelectedCuts){
+ AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
+ if (isLc) Printf("SIGNAL candidate rejected");
+ return;
+ }
+ else {
+ AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
+ }
+
+ AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
+ if (!bachelor) {
+ AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
+ return;
+ }
+
+ //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
+ Double_t probTPCTOF[AliPID::kSPECIES]={0.};
+
+ UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
+ Double_t probProton = 0.;
+ Double_t probPion = 0.;
+ Double_t probKaon = 0.;
+ if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
+ probProton = probTPCTOF[AliPID::kProton];
+ probPion = probTPCTOF[AliPID::kPion];
+ probKaon = probTPCTOF[AliPID::kKaon];
+ }
+
+ if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
+ AliDebug(2, "On-the-fly discarded");
+ return;
+ }
+
+ if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
+ nSelectedAnal++;
+ }
+
+ if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
+
+ if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
+ if (isLc) Printf("SIGNAL candidate rejected");
+ AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
+ return;
+ }
+ else {
+ AliDebug(2, "Yes: Analysis cuts kTracks level passed");
+ }
+
+ Int_t pdgCand = 4122;
+ Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
+ Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
+ Int_t isLc2LBarpi=0, isLc2Lpi=0;
+ Int_t currentLabel = part->GetLabel();
+ Int_t mcLabel = 0;
+ if (fUseMCInfo) {
+ mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
+ if (mcLabel>=0) {
+ if (bachelor->Charge()==-1) isLc2LBarpi=1;
+ if (bachelor->Charge()==+1) isLc2Lpi=1;
+ }
+ }
+
+ Int_t pdgDg2prong[2] = {211, 211};
+ Int_t labelK0S = 0;
+ Int_t isK0S = 0;
+ if (fUseMCInfo) {
+ labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
+ if (labelK0S>=0) isK0S = 1;
+ }
+
+ pdgDg2prong[0] = 211;
+ pdgDg2prong[1] = 2212;
+ Int_t isLambda = 0;
+ Int_t isLambdaBar = 0;
+ Int_t lambdaLabel = 0;
+ if (fUseMCInfo) {
+ lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
+ if (lambdaLabel>=0) {
+ AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
+ if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
+ else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
+ }
+ }
+
+ pdgDg2prong[0] = 11;
+ pdgDg2prong[1] = 11;
+ Int_t isGamma = 0;
+ Int_t gammaLabel = 0;
+ if (fUseMCInfo) {
+ gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
+ if (gammaLabel>=0) {
+ AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
+ if (gammaTrack->GetPdgCode()==22) isGamma = 1;
+ }
+ }
+
+ Int_t pdgTemp = -1;
+ if (currentLabel != -1){
+ AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
+ pdgTemp = tempPart->GetPdgCode();
+ }
+ if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
+ else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
+ else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
+ else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
+ if (isK0S) AliDebug(2, Form("V0 is a K0S"));
+ else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
+ else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
+ else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
+ //else AliDebug(2, Form("V0 is something else!!"));
+
+ Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
+ Double_t invmassLambda = v0part->MassLambda();
+ Double_t invmassLambdaBar = v0part->MassAntiLambda();
+
+ //Double_t nSigmaITSpr=-999.;
+ Double_t nSigmaTPCpr=-999.;
+ Double_t nSigmaTOFpr=-999.;
+
+ //Double_t nSigmaITSpi=-999.;
+ Double_t nSigmaTPCpi=-999.;
+ Double_t nSigmaTOFpi=-999.;
+
+ //Double_t nSigmaITSka=-999.;
+ Double_t nSigmaTPCka=-999.;
+ Double_t nSigmaTOFka=-999.;
+
+ /*
+ cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
+ cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
+ cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
+ cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
+ cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
+ cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
+ cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
+ cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
+ cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
+ */
+
+ nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
+ nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
+ nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
+ nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
+ nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
+ nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
+
+
+ // Fill candidate variable Tree (track selection, V0 invMass selection)
+ if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
+
+ fCandidateVariables[0] = invmassLc;
+ fCandidateVariables[1] = invmassLc2Lpi;
+ fCandidateVariables[2] = invmassK0s;
+ fCandidateVariables[3] = invmassLambda;
+ fCandidateVariables[4] = invmassLambdaBar;
+ fCandidateVariables[5] = part->CosV0PointingAngle();
+ fCandidateVariables[6] = dcaV0;
+ fCandidateVariables[7] = part->Getd0Prong(0);
+ fCandidateVariables[8] = part->Getd0Prong(1);
+ fCandidateVariables[9] = nSigmaTPCpr;
+ fCandidateVariables[10] = nSigmaTPCpi;
+ fCandidateVariables[11] = nSigmaTPCka;
+ fCandidateVariables[12] = nSigmaTOFpr;
+ fCandidateVariables[13] = nSigmaTOFpi;
+ fCandidateVariables[14] = nSigmaTOFka;
+ fCandidateVariables[15] = bachelor->Pt();
+ AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
+ fCandidateVariables[16] = v0pos->Pt();
+ AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
+ fCandidateVariables[17] = v0neg->Pt();
+ fCandidateVariables[18] = v0part->Getd0Prong(0);
+ fCandidateVariables[19] = v0part->Getd0Prong(1);
+ fCandidateVariables[20] = v0part->Pt();
+ fCandidateVariables[21] = v0part->InvMass2Prongs(0,1,11,11);
+ fCandidateVariables[22] = part->Pt();
+ fCandidateVariables[23] = probProton;
+ fCandidateVariables[24] = part->Eta();
+ fCandidateVariables[25] = v0pos->Eta();
+ fCandidateVariables[26] = v0neg->Eta();
+ fCandidateVariables[27] = probPion;
+ fCandidateVariables[28] = probKaon;
+ fCandidateVariables[29] = bachelor->Eta();
+
+ fCandidateVariables[30] = part->P();
+ fCandidateVariables[31] = bachelor->P();
+ fCandidateVariables[32] = v0part->P();
+ fCandidateVariables[33] = v0pos->P();
+ fCandidateVariables[34] = v0neg->P();
+
+ fCandidateVariables[35] = part->Y(4122);
+ fCandidateVariables[36] = bachelor->Y(2212);
+ fCandidateVariables[37] = v0part->Y(310);
+ fCandidateVariables[38] = v0pos->Y(211);
+ fCandidateVariables[39] = v0neg->Y(211);
+
+ fCandidateVariables[40] = v0part->Eta();
+
+ fCandidateVariables[41] = part->DecayLength();
+ fCandidateVariables[42] = part->DecayLengthV0();
+ fCandidateVariables[43] = part->Ct(4122);
+ fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
+
+ EBachelor bachCode = CheckBachelor(part, bachelor, mcArray);
+ EK0S k0SCode = CheckK0S(part, v0part, mcArray);
+
+ fCandidateVariables[45] = bachCode;
+ fCandidateVariables[46] = k0SCode;
+
+ Double_t V0KF[3] = {-999999, -999999, -999999};
+ Double_t errV0KF[3] = {-999999, -999999, -999999};
+ Double_t LcKF[3] = {-999999, -999999, -999999};
+ Double_t errLcKF[3] = {-999999, -999999, -999999};
+ Double_t distances[3] = {-999999, -999999, -999999};
+ Double_t armPolKF[2] = {-999999, -999999};
+
+ if (fCallKFVertexing){
+ Int_t kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
+ AliDebug(2, Form("Result from KF = %d", kfResult));
+ }
+
+ /*
+ for (Int_t i = 0; i< 3; i++){
+ Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
+ }
+ */
+
+ fCandidateVariables[47] = V0KF[0];
+ fCandidateVariables[48] = V0KF[1];
+ fCandidateVariables[49] = V0KF[2];
+
+ fCandidateVariables[50] = errV0KF[0];
+ fCandidateVariables[51] = errV0KF[1];
+ fCandidateVariables[52] = errV0KF[2];
+
+ fCandidateVariables[53] = LcKF[0];
+ fCandidateVariables[54] = LcKF[1];
+ fCandidateVariables[55] = LcKF[2];
+
+ fCandidateVariables[56] = errLcKF[0];
+ fCandidateVariables[57] = errLcKF[1];
+ fCandidateVariables[58] = errLcKF[2];
+
+ fCandidateVariables[59] = distances[0];
+ fCandidateVariables[60] = distances[1];
+ fCandidateVariables[61] = distances[2];
+ fCandidateVariables[62] = armPolKF[0];
+ fCandidateVariables[63] = armPolKF[1];
+ fCandidateVariables[64] = v0part->AlphaV0();
+ fCandidateVariables[65] = v0part->PtArmV0();
+
+ if (fUseMCInfo) {
+ if (isLc){
+ AliDebug(2, "Filling Sgn");
+ fVariablesTreeSgn->Fill();
+ fHistoCodesSgn->Fill(bachCode, k0SCode);
+ }
+ else {
+ if (fFillOnlySgn == kFALSE){
+ AliDebug(2, "Filling Bkg");
+ fVariablesTreeBkg->Fill();
+ fHistoCodesBkg->Fill(bachCode, k0SCode);
+ }
+ }
+ }
+ else {
+ fVariablesTreeSgn->Fill();
+ }
+ }
+
+ return;
+
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
+ Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
+ Double_t* distances, Double_t* armPolKF) {
+
+ //
+ // method to perform KF vertexing
+ // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
+ //
+
+ Int_t codeKFV0 = -1, codeKFLc = -1;
+
+ AliKFVertex primVtxCopy;
+ Int_t nt = 0, ntcheck = 0;
+ Double_t pos[3] = {0., 0., 0.};
+
+ fVtx1->GetXYZ(pos);
+ Int_t contr = fVtx1->GetNContributors();
+ Double_t covmatrix[6] = {0.};
+ fVtx1->GetCovarianceMatrix(covmatrix);
+ Double_t chi2 = fVtx1->GetChi2();
+ AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
+
+ // topological constraint
+ primVtxCopy = AliKFVertex(primaryESDVtxCopy);
+ nt = primaryESDVtxCopy.GetNContributors();
+ ntcheck = nt;
+
+ Int_t pdg[2] = {211, -211};
+ Int_t pdgLc[2] = {2212, 310};
+
+ Int_t pdgDgV0toDaughters[2] = {211, 211};
+
+ Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
+
+ // the KF vertex for the V0 has to be built with the prongs of the V0!
+ Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
+ AliKFParticle V0, positiveV0KF, negativeV0KF;;
+ Int_t labelsv0daugh[2] = {-1, -1};
+ for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
+ AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
+ if(!aodTrack) {
+ AliDebug(2, "No V0 daughters available");
+ return -1;
+ }
+ labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
+ if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
+
+ //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
+ AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
+ if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
+ positiveV0KF = daughterKF;
+ }
+ else {
+ negativeV0KF = daughterKF;
+ }
+ V0.AddDaughter(daughterKF);
+ }
+
+ // Checking the quality of the KF V0 vertex
+ if( V0.GetNDF() < 1 ) {
+ //Printf("Number of degrees of freedom < 1, continuing");
+ return -1;
+ }
+ if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > 3. ) {
+ //Printf("Chi2 per DOF too big, continuing");
+ return -1;
+ }
+
+ if(ftopoConstraint && nt > 0){
+ for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
+ AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
+ //* subtruct daughters from primary vertex
+ if(!aodTrack->GetUsedForPrimVtxFit()) {
+ //Printf("Track %d was not used for primary vertex, continuing", i);
+ continue;
+ }
+ AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
+ primVtxCopy -= daughterKF;
+ ntcheck--;
+ }
+ }
+
+ /*
+ // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
+ if( V0.GetDeviationFromVertex( primVtxCopy ) >3. ) {
+ //Printf("Deviation from vertex too big, continuing");
+ return -1;
+ }
+ */
+
+ if(ftopoConstraint){
+ if(ntcheck>0) { // number of constraints left with which the vertex was built, after removing V0 daughters
+ //* Add V0 to primary vertex to improve the primary vertex resolution
+ primVtxCopy += V0;
+ V0.SetProductionVertex(primVtxCopy); // FIXME: is this correct here?
+ }
+ }
+
+ //* Check chi^2
+ if( TMath::Sqrt( TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) {
+ AliDebug(2, "Final Chi2 per DOF too big, continuing");
+ return -1;
+ }
+
+ //* Get V0 invariant mass
+ Double_t massV0 = 999999, sigmaMassV0 = 999999;
+ Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
+ if( retMV0 ) {
+ if (massV0 < 0) {
+ codeKFV0 = 1; // Mass not ok
+ if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
+ }
+ else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
+ }
+ fHistoMassKFV0->Fill(massV0, sigmaMassV0);
+
+ //* Get V0 decayLength
+ Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
+ Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
+ if( retDLV0 ) {
+ if (sigmaDecayLengthV0 > 1e19) {
+ codeKFV0 = 3; // DecayLength not ok
+ if (massV0 < 0) {
+ codeKFV0 = 6; // DecayLength and Mass not ok
+ if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
+ }
+ else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
+ }
+ }
+ fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
+
+ //* Get V0 life time
+ Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
+ Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
+ if( retTLV0 ) {
+ if (sigmaLifeTimeV0 > 1e19) {
+ codeKFV0 = 4; // LifeTime not ok
+ if (sigmaDecayLengthV0 > 1e19) {
+ codeKFV0 = 9; // LifeTime and DecayLength not ok
+ if (massV0 < 0) {
+ codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
+ if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
+ }
+ else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
+ }
+ else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
+ codeKFV0 = 7; // LifeTime and Mass not ok
+ if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
+ }
+ else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
+ }
+ }
+ fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
+
+ if (codeKFV0 == -1) codeKFV0 = 0;
+ fHistoKFV0->Fill(codeKFV0);
+
+ AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
+
+ fHistoMassV0All->Fill(massV0);
+ fHistoDecayLengthV0All->Fill(decayLengthV0);
+ fHistoLifeTimeV0All->Fill(lifeTimeV0);
+
+ Double_t qtAlphaV0[2] = {0., 0.};
+ Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
+ positiveV0KF.TransportToPoint(vtxV0KF);
+ negativeV0KF.TransportToPoint(vtxV0KF);
+ V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
+ AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
+ fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
+ fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
+ armPolKF[0] = qtAlphaV0[1];
+ armPolKF[1] = qtAlphaV0[0];
+
+ // Checking MC info for V0
+
+ AliAODMCParticle *motherV0 = 0x0;
+ AliAODMCParticle *daughv01 = 0x0;
+ AliAODMCParticle *daughv02 = 0x0;
+
+ if (fUseMCInfo) {
+ daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
+ daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
+ if (!daughv01 && labelsv0daugh[0] > 0){
+ AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
+ isMCokV0 = kFALSE;
+ }
+ if (!daughv02 && labelsv0daugh[1] > 0){
+ AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
+ isMCokV0 = kFALSE;
+ }
+ if (isMCokV0){
+ if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
+ AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
+ motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
+ if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
+ if( motherV0->GetNDaughters() == 2 ){
+ fHistoMassV0True->Fill(massV0);
+ fHistoDecayLengthV0True->Fill(decayLengthV0);
+ fHistoLifeTimeV0True->Fill(lifeTimeV0);
+ fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
+ if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
+ fHistoMassV0TrueK0S->Fill(massV0);
+ fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
+ fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
+ fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
+ }
+ }
+ AliDebug(2, Form("PDG V0 = %d, pi = %d, pj = %d, ndaughters = %d, mc mass = %f, reco mass = %f, v0 mass = %f", motherV0->GetPdgCode(), daughv01->GetPdgCode(), daughv02->GetPdgCode(), motherV0->GetNDaughters(), motherV0->GetCalcMass(), massV0, v0part->MassK0Short()));
+ }
+ else if (!motherV0){
+ AliDebug(3, "could not access MC info for mother, continuing");
+ }
+ else {
+ AliDebug(3, "MC mother is a gluon, continuing");
+ }
+ }
+ else {
+ AliDebug(3, "Background V0!");
+ isBkgV0 = kTRUE;
+ }
+ }
+ }
+
+ AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
+
+ // now use what we just try with the bachelor, to build the Lc
+
+ // topological constraint
+ nt = primVtxCopy.GetNContributors();
+ ntcheck = nt;
+
+ Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
+ AliKFParticle Lc;
+ Int_t labelsLcdaugh[2] = {-1, -1};
+ labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
+ labelsLcdaugh[1] = mcLabelV0;
+
+ if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
+ AliKFParticle daughterKFLc(*bach, pdgLc[0]);
+ Lc.AddDaughter(daughterKFLc);
+ Lc.AddDaughter(V0);
+ if( Lc.GetNDF() < 1 ) {
+ AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
+ return -1;
+ }
+ if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) >3. ) {
+ AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
+ return -1;
+ }
+
+ if(ftopoConstraint && nt > 0){
+ //* subtruct daughters from primary vertex
+ if(!bach->GetUsedForPrimVtxFit()) {
+ AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
+ }
+ else{
+ primVtxCopy -= daughterKFLc;
+ ntcheck--;
+ }
+ /* the V0 was added above, so it is ok to remove it without checking
+ if(!V0->GetUsedForPrimVtxFit()) {
+ Printf("Lc: V0 was not used for primary vertex, continuing");
+ continue;
+ }
+ */
+ primVtxCopy -= V0;
+ ntcheck--;
+ }
+
+ //* Check V0 Chi^2 deviation from primary vertex
+ if( Lc.GetDeviationFromVertex( primVtxCopy ) >3. ) {
+ AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
+ return -1;
+ }
+
+ if(ftopoConstraint){
+ if(ntcheck>0) {
+ //* Add Lc to primary vertex to improve the primary vertex resolution
+ primVtxCopy += Lc;
+ Lc.SetProductionVertex(primVtxCopy);
+ }
+ }
+
+ //* Check chi^2
+ if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) >3. ) {
+ AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
+ return -1;
+ }
+
+ //* Get Lc invariant mass
+ Double_t massLc = 999999, sigmaMassLc= 999999;
+ Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
+ if( retMLc ) {
+ AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
+ if (massLc < 0) {
+ codeKFLc = 1; // Mass not ok
+ if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
+ }
+ else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
+ }
+ fHistoMassKFLc->Fill(massLc, sigmaMassLc);
+
+ //* Get Lc Decay length
+ Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
+ Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
+ if( retDLLc ) {
+ AliDebug(3, "----> Lc: Could not get decay length, and sigma");
+ if (sigmaDecayLengthLc > 1e19) {
+ codeKFLc = 3; // DecayLength not ok
+ if (massLc < 0) {
+ codeKFLc = 6; // DecayLength and Mass not ok
+ if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
+ }
+ else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
+ }
+ }
+ AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
+
+ fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
+
+ //* Get Lc life time
+ Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
+ Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
+ if( retTLLc ) {
+ AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
+ if (sigmaLifeTimeLc > 1e19) {
+ codeKFLc = 4; // LifeTime not ok
+ if (sigmaDecayLengthLc > 1e19) {
+ codeKFLc = 9; // LifeTime and DecayLength not ok
+ if (massLc < 0) {
+ codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
+ if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
+ }
+ else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
+ }
+ else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
+ codeKFLc = 7; // LifeTime and Mass not ok
+ if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
+ }
+ else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
+ }
+ }
+
+ fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
+
+ AliDebug(2, Form("Lc: mass = %f (error = %e), decay length = %f (error = %e), life time = %f (error = %e) --> codeKFLc = %d", massLc, sigmaMassLc, decayLengthLc, sigmaDecayLengthLc, lifeTimeLc, sigmaLifeTimeLc, codeKFLc));
+
+ if (codeKFLc == -1) codeKFLc = 0;
+ fHistoKFLc->Fill(codeKFLc);
+
+ fHistoKF->Fill(codeKFV0, codeKFLc);
+
+ // here we fill the histgrams for all the reconstructed KF vertices for the cascade
+ fHistoMassLcAll->Fill(massLc);
+ fHistoDecayLengthLcAll->Fill(decayLengthLc);
+ fHistoLifeTimeLcAll->Fill(lifeTimeLc);
+
+ fHistoMassV0fromLcAll->Fill(massV0);
+ fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
+ fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
+
+ Double_t xV0 = V0.GetX();
+ Double_t yV0 = V0.GetY();
+ Double_t zV0 = V0.GetZ();
+
+ Double_t xLc = Lc.GetX();
+ Double_t yLc = Lc.GetY();
+ Double_t zLc = Lc.GetZ();
+
+ Double_t xPrimVtx = primVtxCopy.GetX();
+ Double_t yPrimVtx = primVtxCopy.GetY();
+ Double_t zPrimVtx = primVtxCopy.GetZ();
+
+ Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
+ (yPrimVtx - yLc) * (yPrimVtx - yLc) +
+ (zPrimVtx - zLc) * (zPrimVtx - zLc));
+
+ Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
+ (yPrimVtx - yV0) * (yPrimVtx - yV0) +
+ (zPrimVtx - zV0) * (zPrimVtx - zV0));
+
+ Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
+ (yLc - yV0)*(yLc - yV0) +
+ (zLc - zV0)*(zLc - zV0));
+
+ //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
+
+ fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
+ fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
+ fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
+
+ distances[0] = distanceLcToPrimVtx;
+ distances[1] = distanceV0ToPrimVtx;
+ distances[2] = distanceV0ToLc;
+
+ if (fUseMCInfo) {
+
+ AliAODMCParticle *daughv01Lc = 0x0;
+ AliAODMCParticle *K0S = 0x0;
+ AliAODMCParticle *daughv02Lc = 0x0;
+
+ if (labelsLcdaugh[0] >= 0) {
+ // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
+ daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
+ if (!daughv01Lc){
+ AliDebug(3, "Could not access MC info for first daughter of Lc");
+ isMCokLc = kFALSE;
+ }
+ else {
+ AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
+ if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
+ }
+ }
+ else { // The bachelor is a fake
+ isBkgLc = kTRUE;
+ }
+
+ if (labelsLcdaugh[1] >= 0) {
+ //Printf("Getting K0S");
+ K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
+ if (!K0S) {
+ AliDebug(3, "Could not access MC info for second daughter of Lc");
+ isMCokLc = kFALSE;
+ }
+ else{
+ if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
+ }
+ }
+ else{
+ AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
+ isBkgLc = kTRUE;
+ }
+
+ if (!isBkgLc){ // so far, we only checked that the V0 and the bachelor are not fake, and in particular, we know that the V0 is a K0S since we used the MatchToMC method
+ if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
+ Int_t iK0 = K0S->GetMother();
+ if (iK0 < 0) {
+ Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
+ }
+ else { // The K0S has a mother
+ daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
+ if (!daughv02Lc){
+ AliDebug(3, "Could not access MC info for second daughter of Lc");
+ }
+ else { // we can access safely the K0S mother in the MC
+ if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 ){ // This is a true cascade! bachelor and V0 come from the same mother
+ //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
+ AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
+ Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
+ if (motherLc) pdgMum = motherLc->GetPdgCode();
+ if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
+ if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
+ AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
+
+ if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
+ fHistoMassLcTrue->Fill(massLc);
+ fHistoDecayLengthLcTrue->Fill(decayLengthLc);
+ fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
+ fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
+
+ fHistoMassV0fromLcTrue->Fill(massV0);
+ fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
+ fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
+
+ if (TMath::Abs(motherLc->GetPdgCode()) == 4122 && TMath::Abs(motherV0->GetPdgCode()) == 310 && TMath::Abs(daughv01Lc->GetPdgCode()) == 2212){ // This is Lc --> K0S + p (the check on the PDG code of the V0 is useless, since we used MathcToMC with it, but fine...
+ AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
+
+ fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
+ fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
+
+ fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
+ fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
+ fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
+
+ fHistoMassLcSgn->Fill(massLc);
+ fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
+
+ fHistoDecayLengthLcSgn->Fill(decayLengthLc);
+ fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
+
+ fHistoMassV0fromLcSgn->Fill(massV0);
+ fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
+ fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
+ }
+ else {
+ isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
+ }
+
+ // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
+ // First, we compare the vtx of the Lc
+ Double_t xLcMC = motherLc->Xv();
+ Double_t yLcMC = motherLc->Yv();
+ Double_t zLcMC = motherLc->Zv();
+ //Printf("Got MC vtx for Lc");
+ Double_t xProtonMC = daughv01Lc->Xv();
+ Double_t yProtonMC = daughv01Lc->Yv();
+ Double_t zProtonMC = daughv01Lc->Zv();
+ //Printf("Got MC vtx for Proton");
+
+ Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
+ (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
+ (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
+
+ // Then, we compare the vtx of the K0s
+ Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
+ Double_t yV0MC = motherV0->Yv();
+ Double_t zV0MC = motherV0->Zv();
+ //Printf("Got MC vtx for V0");
+ Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
+ Double_t yPionMC = daughv01->Yv();
+ Double_t zPionMC = daughv01->Zv();
+ //Printf("Got MC vtx for Pion");
+
+ Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
+ (yV0MC - yPionMC) * (yV0MC - yPionMC) +
+ (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
+
+ // Then, the K0S vertex wrt primary vertex
+
+ Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
+ (yPionMC - yLcMC) * (yPionMC - yLcMC) +
+ (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
+
+ fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
+ fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
+ fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
+
+ } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
+ else if (!motherLc){
+ AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
+ }
+ else {
+ AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
+ }
+ } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
+ else {
+ isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
+ }
+ } // closing: else { // we can access safely the K0S mother in the MC
+ } // closing: else { // The K0S has a mother
+ } // closing isMCLcok
+ } // closing !isBkgLc
+ } // closing fUseMCInfo
+
+ //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
+ if ( retMV0 == 0 && retMLc == 0){
+ V0KF[0] = massV0;
+ errV0KF[0] = sigmaMassV0;
+ V0KF[1] = decayLengthV0;
+ errV0KF[1] = sigmaDecayLengthV0;
+ V0KF[2] = lifeTimeV0;
+ errV0KF[2] = sigmaLifeTimeV0;
+ LcKF[0] = massLc;
+ errLcKF[0] = sigmaMassLc;
+ LcKF[1] = decayLengthLc;
+ errLcKF[1] = sigmaDecayLengthLc;
+ LcKF[2] = lifeTimeLc;
+ errLcKF[2] = sigmaLifeTimeLc;
+ }
+
+ return 1;
+
+}
+//________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
+ AliAODTrack* bachelor,
+ TClonesArray *mcArray ){
+
+ //Printf("In CheckBachelor");
+
+ // function to check if the bachelor is a p, if it is a p but not from Lc
+ // to be filled for background candidates
+
+ Int_t label = bachelor->GetLabel();
+ if (label == -1) {
+ return kBachFake;
+ }
+
+ AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
+ if (!mcpart) return kBachInvalid;
+ Int_t pdg = mcpart->PdgCode();
+ if (TMath::Abs(pdg) != 2212) {
+ AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
+ return kBachNoProton;
+ }
+ else { // it is a proton
+ //Int_t labelLc = part->GetLabel();
+ Int_t labelLc = FindLcLabel(part, mcArray);
+ //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
+ Int_t bachelorMotherLabel = mcpart->GetMother();
+ //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
+ if (bachelorMotherLabel == -1) {
+ return kBachPrimary;
+ }
+ AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
+ if (!bachelorMother) return kBachInvalid;
+ Int_t pdgMother = bachelorMother->PdgCode();
+ if (TMath::Abs(pdgMother) != 4122) {
+ AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
+ return kBachNoLambdaMother;
+ }
+ else { // it comes from Lc
+ if (labelLc != bachelorMotherLabel){
+ //AliInfo(Form("The proton comes from a Lc, but it is not the candidate we are analyzing (label Lc = %d, label p mother = %d", labelLc, bachelorMotherLabel));
+ AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
+ return kBachDifferentLambdaMother;
+ }
+ else { // it comes from the correct Lc
+ AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
+ return kBachCorrectLambdaMother;
+ }
+ }
+ }
+
+ return kBachInvalid;
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
+ AliAODv0* v0part,
+ //AliAODTrack* v0part,
+ TClonesArray *mcArray ){
+
+ // function to check if the K0Spart is a p, if it is a p but not from Lc
+ // to be filled for background candidates
+
+ //Printf(" CheckK0S");
+
+ //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
+ //Int_t label = v0part->GetLabel();
+ Int_t labelFind = FindV0Label(v0part, mcArray);
+ //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
+ if (labelFind == -1) {
+ return kK0SFake;
+ }
+
+ AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
+ if (!mcpart) return kK0SInvalid;
+ Int_t pdg = mcpart->PdgCode();
+ if (TMath::Abs(pdg) != 310) {
+ AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
+ //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
+ return kK0SNoK0S;
+ }
+ else { // it is a K0S
+ //Int_t labelLc = part->GetLabel();
+ Int_t labelLc = FindLcLabel(part, mcArray);
+ Int_t K0SpartMotherLabel = mcpart->GetMother();
+ if (K0SpartMotherLabel == -1) {
+ return kK0SWithoutMother;
+ }
+ AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
+ if (!K0SpartMother) return kK0SInvalid;
+ Int_t pdgMotherK0S = K0SpartMother->PdgCode();
+ if (TMath::Abs(pdgMotherK0S) != 311) {
+ AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
+ //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
+ return kK0SNotFromK0;
+ }
+ else { // the K0S comes from a K0
+ Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
+ if (K0MotherLabel == -1) {
+ return kK0Primary;
+ }
+ AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
+ if (!K0Mother) return kK0SInvalid;
+ Int_t pdgK0Mother = K0Mother->PdgCode();
+ if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
+ AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
+ //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
+ return kK0NoLambdaMother;
+ }
+ else { // the K0 comes from Lc
+ if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
+ AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
+ //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
+ return kK0DifferentLambdaMother;
+ }
+ else { // it comes from the correct Lc
+ AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
+ //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
+ return kK0CorrectLambdaMother;
+ }
+ }
+ }
+ }
+
+ return kK0SInvalid;
+
+}
+
+//----------------------------------------------------------------------------
+Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
+{
+
+ //Printf(" FindV0Label");
+
+ // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
+
+ Int_t labMother[2]={-1, -1};
+ AliAODMCParticle *part=0;
+ AliAODMCParticle *mother=0;
+ Int_t dgLabels = -1;
+
+ Int_t ndg = v0part->GetNDaughters();
+ if (ndg != 2) {
+ AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
+ }
+
+ for(Int_t i = 0; i < 2; i++) {
+ AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
+ dgLabels = trk->GetLabel();
+ if (dgLabels == -1) {
+ //printf("daughter with negative label %d\n", dgLabels);
+ return -1;
+ }
+ part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
+ if (!part) {
+ //printf("no MC particle\n");
+ return -1;
+ }
+ labMother[i] = part->GetMother();
+ if (labMother[i] != -1){
+ mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
+ if(!mother) {
+ //printf("no MC mother particle\n");
+ return -1;
+ }
+ }
+ else {
+ return -1;
+ }
+ }
+
+ if (labMother[0] == labMother[1]) return labMother[0];
+
+ else return -1;
+
+}
+
+//----------------------------------------------------------------------------
+Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
+{
+
+ // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
+
+ //Printf(" FindLcLabel");
+
+ AliAODMCParticle *part=0;
+ AliAODMCParticle *mother=0;
+ AliAODMCParticle *grandmother=0;
+ Int_t dgLabels = -1;
+
+ Int_t ndg = cascade->GetNDaughters();
+ if (ndg != 2) {
+ AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
+ }
+
+ // daughter 0 --> bachelor, daughter 1 --> V0
+ AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
+ dgLabels = trk->GetLabel();
+ if (dgLabels == -1 ) {
+ //printf("daughter with negative label %d\n", dgLabels);
+ return -1;
+ }
+ part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
+ if (!part) {
+ //printf("no MC particle\n");
+ return -1;
+ }
+ Int_t labMotherBach = part->GetMother();
+ if (labMotherBach == -1){
+ return -1;
+ }
+ mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
+ if(!mother) {
+ //printf("no MC mother particle\n");
+ return -1;
+ }
+
+ AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
+ dgLabels = FindV0Label(v0, mcArray);
+ if (dgLabels == -1 ) {
+ //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
+ return -1;
+ }
+ part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
+ if (!part) {
+ //printf("no MC particle\n");
+ return -1;
+ }
+ Int_t labMotherv0 = part->GetMother();
+ if (labMotherv0 == -1){
+ return -1;
+ }
+ mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
+ if(!mother) {
+ //printf("no MC mother particle\n");
+ return -1;
+ }
+ Int_t labGrandMotherv0 = mother->GetMother();
+ if (labGrandMotherv0 == -1){
+ return -1;
+ }
+ grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
+ if(!grandmother) {
+ //printf("no MC mother particle\n");
+ return -1;
+ }
+
+ //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
+ if (labGrandMotherv0 == labMotherBach) return labMotherBach;
+
+ else return -1;
+
+}
+
+//----------------------------------------------------------------------------
+Bool_t AliAnalysisTaskSELc2V0bachelorTMVA::CheckInjection(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
+
+ // temporary method until the same is in AliVertexingHFUtils
+
+ AliAODTrack* bach = cand->GetBachelor();
+ if(fUtils->IsTrackInjected(bach, header, arrayMC)) {
+ AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
+ return kTRUE;
+ }
+ AliAODv0* v0 = cand->Getv0();
+ Int_t nprongs = v0->GetNProngs();
+ for(Int_t i = 0; i < nprongs; i++){
+ AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
+ if(fUtils->IsTrackInjected(daugh,header,arrayMC)) {
+ AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
+ return kTRUE;
+ }
+ }
+ return kFALSE;
+}
+
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff