/************************************************************************** * 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 #include #include #include #include #include #include #include "TROOT.h" #include #include "TVector3.h" #include #include #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 "AliExternalTrackParam.h" #include "AliESDtrack.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), fCutKFChi2NDF(999999.), fCutKFDeviationFromVtx(999999.), fCutKFDeviationFromVtxV0(0.), fCurrentEvent(-1), fBField(0), fKeepingOnlyPYTHIABkg(kFALSE), fHistoMCLcK0SpGen(0x0), fHistoMCLcK0SpGenAcc(0x0), fHistoMCLcK0SpGenLimAcc(0x0) { // // 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), fCutKFChi2NDF(999999.), fCutKFDeviationFromVtx(999999.), fCutKFDeviationFromVtxV0(0.), fCurrentEvent(-1), fBField(0), fKeepingOnlyPYTHIABkg(kFALSE), fHistoMCLcK0SpGen(0x0), fHistoMCLcK0SpGenAcc(0x0), fHistoMCLcK0SpGenLimAcc(0x0) { // // 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 || fKeepingOnlyPYTHIABkg)) 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 (GetOutputData(1)); if (!fOutput) { AliError("fOutput not available"); return; } //AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries())); //AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries())); //AliDebug(2, Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries())); if(fHistoMCLcK0SpGen) { AliInfo(Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries())); } else { AliInfo("fHistoMCLcK0SpGen not available"); } if(fHistoMCLcK0SpGenAcc) { AliInfo(Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries())); } else { AliInfo("fHistoMCLcK0SpGenAcc not available"); } if(fVariablesTreeSgn) { AliInfo(Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries())); } else { AliInfo("fVariablesTreeSgn not available"); } fOutputKF = dynamic_cast (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 = 86; 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]="TPCProtonProb"; fCandidateVariableNames[28]="TOFProtonProb"; 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"; fCandidateVariableNames[66]="ITSrefitV0pos"; fCandidateVariableNames[67]="ITSrefitV0neg"; fCandidateVariableNames[68]="TPCClV0pos"; fCandidateVariableNames[69]="TPCClV0neg"; fCandidateVariableNames[70]="v0Xcoord"; fCandidateVariableNames[71]="v0Ycoord"; fCandidateVariableNames[72]="v0Zcoord"; fCandidateVariableNames[73]="primVtxX"; fCandidateVariableNames[74]="primVtxY"; fCandidateVariableNames[75]="primVtxZ"; fCandidateVariableNames[76]="ITSclBach"; fCandidateVariableNames[77]="SPDclBach"; fCandidateVariableNames[78]="ITSclV0pos"; fCandidateVariableNames[79]="SPDclV0pos"; fCandidateVariableNames[80]="ITSclV0neg"; fCandidateVariableNames[81]="SPDclV0neg"; fCandidateVariableNames[82]="alphaArmLc"; fCandidateVariableNames[83]="alphaArmLcCharge"; fCandidateVariableNames[84]="ptArmLc"; fCandidateVariableNames[85]="CosThetaStar"; for(Int_t ivar=0; ivarBranch(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", 4, -0.5, 3.5); TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"}; for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){ fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data()); } //fOutput->Add(fVariablesTreeSgn); //fOutput->Add(fVariablesTreeBkg); const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.}; fHistoMCLcK0SpGen = new TH1F("fHistoMCLcK0SpGen", "fHistoMCLcK0SpGen", 14, ptbins); fHistoMCLcK0SpGenAcc = new TH1F("fHistoMCLcK0SpGenAcc", "fHistoMCLcK0SpGenAcc", 14, ptbins); fHistoMCLcK0SpGenLimAcc = new TH1F("fHistoMCLcK0SpGenLimAcc", "fHistoMCLcK0SpGenLimAcc", 14, ptbins); 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); fOutput->Add(fHistoMCLcK0SpGen); fOutput->Add(fHistoMCLcK0SpGenAcc); fOutput->Add(fHistoMCLcK0SpGenLimAcc); 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]); fCounter = new AliNormalizationCounter("NormalizationCounter"); fCounter->Init(); PostData(2, fCounter); // Histograms from KF if (fCallKFVertexing){ 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"); if (fCallKFVertexing){ 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(fHistoMassLcSgnFromAOD); 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); return; } //_________________________________________________ void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *) { // user exec if (!fInputEvent) { AliError("NO EVENT FOUND!"); return; } fCurrentEvent++; AliDebug(2, Form("Processing event = %d", fCurrentEvent)); AliAODEvent* aodEvent = dynamic_cast(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 (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); } Int_t runnumber = aodEvent->GetRunNumber(); if (aodEvent->GetTriggerMask() == 0 && (runnumber >= 195344 && runnumber <= 195677)){ AliDebug(3,"Event rejected because of null trigger mask"); return; } fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo); // mc analysis TClonesArray *mcArray = 0; AliAODMCHeader *mcHeader=0; if (fUseMCInfo) { // MC array need for matching mcArray = dynamic_cast(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; } Double_t zMCVertex = mcHeader->GetVtxZ(); if (TMath::Abs(zMCVertex) > fAnalCuts->GetMaxVtxZ()){ AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fAnalCuts->GetMaxVtxZ(), fAnalCuts->GetMaxVtxZ())); AliInfo(Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fAnalCuts->GetMaxVtxZ(), fAnalCuts->GetMaxVtxZ())); return; } //Printf("Filling MC histo"); FillMCHisto(mcArray); } // 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 fBField = aodEvent->GetMagneticField(); AliKFParticle::SetField(fBField); 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::FillMCHisto(TClonesArray *mcArray){ // method to fill MC histo: how many Lc --> K0S + p are there at MC level for (Int_t iPart=0; iPartGetEntriesFast(); iPart++) { AliAODMCParticle* mcPart = dynamic_cast(mcArray->At(iPart)); if (!mcPart){ AliError("Failed casting particle from MC array!, Skipping particle"); continue; } Int_t pdg = mcPart->GetPdgCode(); if (TMath::Abs(pdg) != 4122){ AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg)); continue; } AliDebug(2, Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg)); Int_t labeldaugh0 = mcPart->GetDaughter(0); Int_t labeldaugh1 = mcPart->GetDaughter(1); if (labeldaugh0 <= 0 || labeldaugh1 <= 0){ AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!")); continue; } else if (labeldaugh1 - labeldaugh0 == 1){ AliDebug(2, Form("Step 1 ok: The MC particle has correct daughters!!")); AliAODMCParticle* daugh0 = dynamic_cast(mcArray->At(labeldaugh0)); AliAODMCParticle* daugh1 = dynamic_cast(mcArray->At(labeldaugh1)); if(!daugh0 || !daugh1){ AliDebug(2,"Particle daughters not properly retrieved!"); return; } Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode()); Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode()); AliAODMCParticle* bachelorMC = daugh0; AliAODMCParticle* v0MC = daugh1; AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1)); if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){ // we are in the case of Lc --> K0 + p; now we have to check if the K0 decays in K0S, and if this goes in pi+pi- /// first, we set the bachelor and the v0: above we assumed first proton and second V0, but we could have to change it: if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) { bachelorMC = daugh1; v0MC = daugh0; } AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters())); if (v0MC->GetNDaughters() != 1) { AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing..."); continue; } else { // So far: Lc --> K0 + p, K0 with 1 daughter AliDebug(2, "Step 2 ok: The K0 does decay in 1 body only! "); Int_t labelK0daugh = v0MC->GetDaughter(0); AliAODMCParticle* partK0S = dynamic_cast(mcArray->At(labelK0daugh)); if(!partK0S){ AliError("Error while casting particle! returning a NULL array"); continue; } else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){ AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies"); continue; } else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies AliDebug(2, "Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies"); Int_t labelK0Sdaugh0 = partK0S->GetDaughter(0); Int_t labelK0Sdaugh1 = partK0S->GetDaughter(1); AliAODMCParticle* daughK0S0 = dynamic_cast(mcArray->At(labelK0Sdaugh0)); AliAODMCParticle* daughK0S1 = dynamic_cast(mcArray->At(labelK0Sdaugh1)); if (!daughK0S0 || ! daughK0S1){ AliDebug(2, "Could not access K0S daughters, continuing..."); continue; } else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access AliDebug(2, "Step 4 ok: Could access K0S daughters, continuing..."); Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode(); Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode(); if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){ AliDebug(2, "The K0S does not decay in pi+pi-, continuing"); //AliInfo("The K0S does not decay in pi+pi-, continuing"); } else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi- if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) { AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt())); if(TMath::Abs(mcPart->Y()) < 0.5) fHistoMCLcK0SpGenLimAcc->Fill(mcPart->Pt()); //AliInfo(Form("\nparticle = %d, Filling MC Gen histo\n", iPart)); fHistoMCLcK0SpGen->Fill(mcPart->Pt()); if(!(TMath::Abs(bachelorMC->Eta()) > 0.9 || bachelorMC->Pt() < 0.1 || TMath::Abs(daughK0S0->Eta()) > 0.9 || daughK0S0->Pt() < 0.1 || TMath::Abs(daughK0S1->Eta()) > 0.9 || daughK0S1->Pt() < 0.1)) { fHistoMCLcK0SpGenAcc->Fill(mcPart->Pt()); } } else { AliDebug(2, "not in fiducial acceptance! Skipping"); continue; } } } } } } } } } // closing loop over mcArray return; } //------------------------------------------------------------------------------- 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(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(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) is a Lc!", iLctopK0s, nCascades)); AliAODMCParticle *partLc = dynamic_cast(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(lcK0spr->Getv0()); AliAODTrack * bachPart = dynamic_cast(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(lcK0spr->Getv0PositiveTrack()); AliAODTrack * v0Neg = dynamic_cast(lcK0spr->Getv0NegativeTrack()); if (!v0Neg || !v0Pos) { 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(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(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; } } else if (fKeepingOnlyPYTHIABkg){ // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event! AliAODTrack *bachelor = (AliAODTrack*)lcK0spr->GetBachelor(); AliAODTrack *v0pos = (AliAODTrack*)lcK0spr->Getv0PositiveTrack(); AliAODTrack *v0neg = (AliAODTrack*)lcK0spr->Getv0NegativeTrack(); if (!bachelor || !v0pos || !v0neg) { AliDebug(2, "Cannot retrieve one of the tracks while checking origin, continuing"); continue; } else { Int_t labelbachelor = TMath::Abs(bachelor->GetLabel()); Int_t labelv0pos = TMath::Abs(v0pos->GetLabel()); Int_t labelv0neg = TMath::Abs(v0neg->GetLabel()); AliAODMCParticle* MCbachelor = (AliAODMCParticle*)mcArray->At(labelbachelor); AliAODMCParticle* MCv0pos = (AliAODMCParticle*)mcArray->At(labelv0pos); AliAODMCParticle* MCv0neg = (AliAODMCParticle*)mcArray->At(labelv0neg); if (!MCbachelor || !MCv0pos || !MCv0neg) { AliDebug(2, "Cannot retrieve MC particle for one of the tracks while checking origin, continuing"); continue; } else { Int_t isBachelorFromPythia = fUtils->CheckOrigin(mcArray, MCbachelor, kTRUE); Int_t isv0posFromPythia = fUtils->CheckOrigin(mcArray, MCv0pos, kTRUE); Int_t isv0negFromPythia = fUtils->CheckOrigin(mcArray, MCv0neg, kTRUE); if (isBachelorFromPythia != 0 && isv0posFromPythia != 0 && isv0negFromPythia != 0){ AliDebug(2, "The candidate is from PYTHIA (i.e. all daughters originate from a quark), keeping it to fill background"); fHistoBackground->Fill(2); } else { AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background"); fHistoBackground->Fill(3); continue; } } } } } } FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray, iLctopK0s); } return; } //________________________________________________________________________ void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part, Int_t isLc, Int_t &nSelectedAnal, AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *mcArray, Int_t iLctopK0s){ // // 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 V0 window cuts!"); if (isLc) Printf("SIGNAL candidate rejected: V0 window cuts"); return; } else AliDebug(2, "Yes: The candidate has passed the mass cuts!"); 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: cascade window cuts"); 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 AliDebug(2, Form("recoAnalysisCuts = %d", cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate) & (AliRDHFCutsLctoV0::kLcToK0Spr))); 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]={-1.}; UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF); AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed)); Double_t probProton = -1.; // Double_t probPion = -1.; // Double_t probKaon = -1.; if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) { AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton])); probProton = probTPCTOF[AliPID::kProton]; // probPion = probTPCTOF[AliPID::kPion]; // probKaon = probTPCTOF[AliPID::kKaon]; } else { // if you don't have both TOF and TPC, try only TPC fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC); AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC"); detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF); AliDebug(2,Form(" detUsed (TPC case) = %d", detUsed)); if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) { probProton = probTPCTOF[AliPID::kProton]; // probPion = probTPCTOF[AliPID::kPion]; // probKaon = probTPCTOF[AliPID::kKaon]; AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton])); } else { AliDebug(2, "Only TPC did not work..."); } // resetting mask to ask for both TPC+TOF fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF); } AliDebug(2, Form("probProton = %f", probProton)); // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) ) Double_t probProtonTPC = -1.; Double_t probProtonTOF = -1.; Double_t pidTPC[AliPID::kSPECIES]={-1.}; Double_t pidTOF[AliPID::kSPECIES]={-1.}; Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC); Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF); if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton]; if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton]; // checking V0 status (on-the-fly vs offline) 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] = probProtonTPC; fCandidateVariables[28] = probProtonTOF; 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 = kBachInvalid; EK0S k0SCode = kK0SInvalid; if (fUseMCInfo) { bachCode = CheckBachelor(part, bachelor, mcArray); 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(); AliDebug(2, Form("v0pos->GetStatus() & AliESDtrack::kITSrefit= %d, v0neg->GetStatus() & AliESDtrack::kITSrefit = %d, v0pos->GetTPCClusterInfo(2, 1)= %f, v0neg->GetTPCClusterInfo(2, 1) = %f", (Int_t)(v0pos->GetStatus() & AliESDtrack::kITSrefit), (Int_t)(v0pos->GetStatus() & AliESDtrack::kITSrefit), v0pos->GetTPCClusterInfo(2, 1), v0neg->GetTPCClusterInfo(2, 1))); fCandidateVariables[66] = v0pos->GetStatus() & AliESDtrack::kITSrefit; fCandidateVariables[67] = v0neg->GetStatus() & AliESDtrack::kITSrefit; fCandidateVariables[68] = v0pos->GetTPCClusterInfo(2, 1); fCandidateVariables[69] = v0neg->GetTPCClusterInfo(2, 1); fCandidateVariables[70] = v0part->Xv(); fCandidateVariables[71] = v0part->Yv(); fCandidateVariables[72] = v0part->Zv(); fCandidateVariables[73] = fVtx1->GetX(); fCandidateVariables[74] = fVtx1->GetY(); fCandidateVariables[75] = fVtx1->GetZ(); fCandidateVariables[76] = bachelor->GetITSNcls(); fCandidateVariables[77] = bachelor->HasPointOnITSLayer(0) + bachelor->HasPointOnITSLayer(1); fCandidateVariables[78] = v0pos->GetITSNcls(); fCandidateVariables[79] = v0pos->HasPointOnITSLayer(0) + v0pos->HasPointOnITSLayer(1); fCandidateVariables[80] = v0neg->GetITSNcls(); fCandidateVariables[81] = v0neg->HasPointOnITSLayer(0) + v0neg->HasPointOnITSLayer(1); TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz()); TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz()); TVector3 momTot(part->Px(), part->Py(), part->Pz()); Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag(); Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag(); Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2); Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) ); Double_t ptArmLc = mom1.Perp(momTot); fCandidateVariables[82] = alphaArmLc; fCandidateVariables[83] = alphaArmLcCharge; fCandidateVariables[84] = ptArmLc; Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // mass K0S PDG Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); // mass Proton PDG Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass(); // mass Lc PDG Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG); Double_t e = part->E(4122); Double_t beta = part->P()/e; Double_t gamma = e/massLcPDG; Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar; fCandidateVariables[85] = cts; if (fUseMCInfo) { if (isLc){ AliDebug(2, Form("Reco particle %d --> Filling Sgn", iLctopK0s)); 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}; Int_t idv0daugh[2] = {-1, -1}; AliExternalTrackParam* esdv0Daugh1 = 0x0; AliExternalTrackParam* esdv0Daugh2 = 0x0; 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; } Double_t xyz[3], pxpypz[3], cv[21]; Short_t sign; aodTrack->GetXYZ(xyz); aodTrack->PxPyPz(pxpypz); aodTrack->GetCovarianceXYZPxPyPz(cv); sign = aodTrack->Charge(); AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign); if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign); else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign); labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel()); idv0daugh[ipr] = aodTrack->GetID(); 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; } } Double_t xn=0., xp=0.;//, dca; AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2)); // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp); AliExternalTrackParam tr1(*esdv0Daugh1); AliExternalTrackParam tr2(*esdv0Daugh2); tr1.PropagateTo(xn, fBField); tr2.PropagateTo(xp, fBField); AliKFParticle daughterKF1(tr1, 211); AliKFParticle daughterKF2(tr2, 211); V0.AddDaughter(positiveV0KF); V0.AddDaughter(negativeV0KF); //V0.AddDaughter(daughterKF1); //V0.AddDaughter(daughterKF2); delete esdv0Daugh1; delete esdv0Daugh2; esdv0Daugh1=0; esdv0Daugh2=0; // 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())) > fCutKFChi2NDF ) { //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 ) < fCutKFDeviationFromVtxV0) { //Printf("Deviation from vertex 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); if (massV0 < 0.4) Printf("\n\n>>>>>>>>>> Found the Funny V0 (mass = %f, sigma = %f, AOD mass = %f): labels of the tracks = %d, %d, id = %d and %d", massV0, sigmaMassV0, v0part->MassK0Short(), labelsv0daugh[0], labelsv0daugh[1], idv0daugh[0], idv0daugh[1]); if (massV0 > 0.55) Printf("\n\n>>>>>>>>>> Found the Funny V0 (mass = %f, , sigma = %f, AOD mass = %f): labels of the tracks = %d, %d, id = %d and %d", massV0, sigmaMassV0, v0part->MassK0Short(), labelsv0daugh[0], labelsv0daugh[1], idv0daugh[0], idv0daugh[1]); Printf("Vertices: KF: x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ()); Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv()); //Printf("Got MC vtx for V0"); if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) { AliAODMCParticle* tmpdaughv01 = dynamic_cast(mcArray->At(labelsv0daugh[0])); AliAODMCParticle* tmpdaughv02 = dynamic_cast(mcArray->At(labelsv0daugh[1])); if (!tmpdaughv01 && labelsv0daugh[0] > 0){ AliDebug(2, "Could not access MC info for first daughter of V0, continuing"); } if (!tmpdaughv02 && labelsv0daugh[1] > 0){ AliDebug(2, "Could not access MC info for second daughter of V0, continuing"); } if(tmpdaughv01){ Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays Double_t yPionMC = tmpdaughv01->Yv(); Double_t zPionMC = tmpdaughv01->Zv(); //Printf("Got MC vtx for Pion"); Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC); } } else { Printf("Not a true V0"); } //massV0=-1;//return -1;// !!!! // 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); TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310); Double_t massPDGK0S = particlePDG->Mass(); V0.SetMassConstraint(massPDGK0S); 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())) > fCutKFChi2NDF) { 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 Lc Chi^2 deviation from primary vertex /* if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) { 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())) > fCutKFChi2NDF) { AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())))); return -1; } if(ftopoConstraint){ V0.SetProductionVertex(Lc); } // After setting the vertex of the V0, getting/filling some info //* 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(mcArray->At(labelsv0daugh[0])); daughv02 = dynamic_cast(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(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)); // Going back to Lc //* 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(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(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(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 && (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(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"); Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC); 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(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(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(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(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(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; }