1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appeuear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 /* $Id: AliAnalysisTaskSELc2V0bachelorTMVA.cxx 62231 2013-04-29 09:47:26Z fprino $ */
20 // Base class for Lc2V0 Analysis to be used with TMVA
24 //------------------------------------------------------------------------------------------
26 // Author: C. Zampolli (from AliAnalysisTaskSELc2V0bachelor class by A.De Caro, P. Pagano)
28 //------------------------------------------------------------------------------------------
31 #include <TParticle.h>
32 #include <TParticlePDG.h>
38 #include <TDatabasePDG.h>
41 #include <AliAnalysisDataSlot.h>
42 #include <AliAnalysisDataContainer.h>
44 #include "AliMCEvent.h"
45 #include "AliAnalysisManager.h"
46 #include "AliAODMCHeader.h"
47 #include "AliAODHandler.h"
49 #include "AliAODVertex.h"
50 #include "AliAODRecoDecay.h"
51 #include "AliAODRecoDecayHF.h"
52 #include "AliAODRecoCascadeHF.h"
53 #include "AliAnalysisVertexingHF.h"
54 #include "AliESDtrack.h"
55 #include "AliAODTrack.h"
57 #include "AliAODMCParticle.h"
58 #include "AliAnalysisTaskSE.h"
59 #include "AliAnalysisTaskSELc2V0bachelorTMVA.h"
60 #include "AliNormalizationCounter.h"
61 #include "AliAODPidHF.h"
62 #include "AliPIDResponse.h"
63 #include "AliPIDCombined.h"
64 #include "AliTOFPIDResponse.h"
65 #include "AliInputEventHandler.h"
66 #include "AliAODRecoDecayHF3Prong.h"
67 #include "AliKFParticle.h"
68 #include "AliKFVertex.h"
69 #include "AliExternalTrackParam.h"
70 #include "AliESDtrack.h"
75 ClassImp(AliAnalysisTaskSELc2V0bachelorTMVA)
77 //__________________________________________________________________________
78 AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA():
85 fIsK0sAnalysis(kFALSE),
89 fUseOnTheFlyV0(kFALSE),
90 fIsEventSelected(kFALSE),
93 fCandidateVariables(),
99 fHistoLcBeforeCuts(0),
100 fHistoFiducialAcceptance(0),
103 fHistoLcpKpiBeforeCuts(0),
106 fHistoDistanceLcToPrimVtx(0),
107 fHistoDistanceV0ToPrimVtx(0),
108 fHistoDistanceV0ToLc(0),
110 fHistoDistanceLcToPrimVtxSgn(0),
111 fHistoDistanceV0ToPrimVtxSgn(0),
112 fHistoDistanceV0ToLcSgn(0),
114 fHistoVtxLcResidualToPrimVtx(0),
115 fHistoVtxV0ResidualToPrimVtx(0),
116 fHistoVtxV0ResidualToLc(0),
119 fHistoDecayLengthV0All(0),
120 fHistoLifeTimeV0All(0),
123 fHistoDecayLengthV0True(0),
124 fHistoLifeTimeV0True(0),
126 fHistoMassV0TrueFromAOD(0),
128 fHistoMassV0TrueK0S(0),
129 fHistoDecayLengthV0TrueK0S(0),
130 fHistoLifeTimeV0TrueK0S(0),
132 fHistoMassV0TrueK0SFromAOD(0),
135 fHistoDecayLengthLcAll(0),
136 fHistoLifeTimeLcAll(0),
139 fHistoDecayLengthLcTrue(0),
140 fHistoLifeTimeLcTrue(0),
142 fHistoMassLcTrueFromAOD(0),
144 fHistoMassV0fromLcAll(0),
145 fHistoDecayLengthV0fromLcAll(0),
146 fHistoLifeTimeV0fromLcAll(0),
148 fHistoMassV0fromLcTrue(0),
149 fHistoDecayLengthV0fromLcTrue(0),
150 fHistoLifeTimeV0fromLcTrue(0),
153 fHistoMassLcSgnFromAOD(0),
154 fHistoDecayLengthLcSgn(0),
155 fHistoLifeTimeLcSgn(0),
157 fHistoMassV0fromLcSgn(0),
158 fHistoDecayLengthV0fromLcSgn(0),
159 fHistoLifeTimeV0fromLcSgn(0),
166 fHistoDecayLengthKFV0(0),
167 fHistoLifeTimeKFV0(0),
170 fHistoDecayLengthKFLc(0),
171 fHistoLifeTimeKFLc(0),
173 fHistoArmenterosPodolanskiV0KF(0),
174 fHistoArmenterosPodolanskiV0KFSgn(0),
175 fHistoArmenterosPodolanskiV0AOD(0),
176 fHistoArmenterosPodolanskiV0AODSgn(0),
180 ftopoConstraint(kTRUE),
181 fCallKFVertexing(kFALSE),
182 fKeepingOnlyHIJINGBkg(kFALSE),
185 fCutKFChi2NDF(999999.),
186 fCutKFDeviationFromVtx(999999.),
187 fCutKFDeviationFromVtxV0(0.),
190 fKeepingOnlyPYTHIABkg(kFALSE),
191 fHistoMCLcK0SpGen(0x0),
192 fHistoMCLcK0SpGenAcc(0x0),
193 fHistoMCLcK0SpGenLimAcc(0x0)
199 //___________________________________________________________________________
200 AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
201 AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
202 AliAnalysisTaskSE(name),
208 fIsK0sAnalysis(kFALSE),
212 fUseOnTheFlyV0(useOnTheFly),
213 fIsEventSelected(kFALSE),
214 fVariablesTreeSgn(0),
215 fVariablesTreeBkg(0),
216 fCandidateVariables(),
221 fFillOnlySgn(kFALSE),
222 fHistoLcBeforeCuts(0),
223 fHistoFiducialAcceptance(0),
226 fHistoLcpKpiBeforeCuts(0),
229 fHistoDistanceLcToPrimVtx(0),
230 fHistoDistanceV0ToPrimVtx(0),
231 fHistoDistanceV0ToLc(0),
233 fHistoDistanceLcToPrimVtxSgn(0),
234 fHistoDistanceV0ToPrimVtxSgn(0),
235 fHistoDistanceV0ToLcSgn(0),
237 fHistoVtxLcResidualToPrimVtx(0),
238 fHistoVtxV0ResidualToPrimVtx(0),
239 fHistoVtxV0ResidualToLc(0),
242 fHistoDecayLengthV0All(0),
243 fHistoLifeTimeV0All(0),
246 fHistoDecayLengthV0True(0),
247 fHistoLifeTimeV0True(0),
249 fHistoMassV0TrueFromAOD(0),
251 fHistoMassV0TrueK0S(0),
252 fHistoDecayLengthV0TrueK0S(0),
253 fHistoLifeTimeV0TrueK0S(0),
255 fHistoMassV0TrueK0SFromAOD(0),
258 fHistoDecayLengthLcAll(0),
259 fHistoLifeTimeLcAll(0),
262 fHistoDecayLengthLcTrue(0),
263 fHistoLifeTimeLcTrue(0),
265 fHistoMassLcTrueFromAOD(0),
267 fHistoMassV0fromLcAll(0),
268 fHistoDecayLengthV0fromLcAll(0),
269 fHistoLifeTimeV0fromLcAll(0),
271 fHistoMassV0fromLcTrue(0),
272 fHistoDecayLengthV0fromLcTrue(0),
273 fHistoLifeTimeV0fromLcTrue(0),
276 fHistoMassLcSgnFromAOD(0),
277 fHistoDecayLengthLcSgn(0),
278 fHistoLifeTimeLcSgn(0),
280 fHistoMassV0fromLcSgn(0),
281 fHistoDecayLengthV0fromLcSgn(0),
282 fHistoLifeTimeV0fromLcSgn(0),
289 fHistoDecayLengthKFV0(0),
290 fHistoLifeTimeKFV0(0),
293 fHistoDecayLengthKFLc(0),
294 fHistoLifeTimeKFLc(0),
296 fHistoArmenterosPodolanskiV0KF(0),
297 fHistoArmenterosPodolanskiV0KFSgn(0),
298 fHistoArmenterosPodolanskiV0AOD(0),
299 fHistoArmenterosPodolanskiV0AODSgn(0),
303 ftopoConstraint(kTRUE),
304 fCallKFVertexing(kFALSE),
305 fKeepingOnlyHIJINGBkg(kFALSE),
308 fCutKFChi2NDF(999999.),
309 fCutKFDeviationFromVtx(999999.),
310 fCutKFDeviationFromVtxV0(0.),
313 fKeepingOnlyPYTHIABkg(kFALSE),
314 fHistoMCLcK0SpGen(0x0),
315 fHistoMCLcK0SpGenAcc(0x0),
316 fHistoMCLcK0SpGenLimAcc(0x0)
320 // Constructor. Initialization of Inputs and Outputs
322 Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
324 DefineOutput(1, TList::Class()); // Tree signal + Tree Bkg + histoEvents
325 DefineOutput(2, AliNormalizationCounter::Class()); // normalization counter object
326 DefineOutput(3, TList::Class()); // Cuts
327 DefineOutput(4, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
328 DefineOutput(5, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
329 DefineOutput(6, TList::Class()); // Tree signal + Tree Bkg + histoEvents
333 //___________________________________________________________________________
334 AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
338 Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
368 if(fVariablesTreeSgn){
369 delete fVariablesTreeSgn;
370 fVariablesTreeSgn = 0;
373 if(fVariablesTreeBkg){
374 delete fVariablesTreeBkg;
375 fVariablesTreeBkg = 0;
389 //_________________________________________________
390 void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
395 fIsEventSelected=kFALSE;
397 if (fDebug > 1) AliInfo("Init");
399 fListCuts = new TList();
400 fListCuts->SetOwner();
401 fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
402 PostData(3,fListCuts);
404 if (fUseMCInfo && (fKeepingOnlyHIJINGBkg || fKeepingOnlyPYTHIABkg)) fUtils = new AliVertexingHFUtils();
409 //________________________________________ terminate ___________________________
410 void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
412 // The Terminate() function is the last function to be called during
413 // a query. It always runs on the client, it can be used to present
414 // the results graphically or save the results to file.
416 AliInfo("Terminate");
417 AliAnalysisTaskSE::Terminate();
419 fOutput = dynamic_cast<TList*> (GetOutputData(1));
421 AliError("fOutput not available");
426 //AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries()));
427 //AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries()));
428 //AliDebug(2, Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries()));
429 if(fHistoMCLcK0SpGen) {
430 AliInfo(Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0SpGen->GetEntries()));
432 AliInfo("fHistoMCLcK0SpGen not available");
434 if(fHistoMCLcK0SpGenAcc) {
435 AliInfo(Form("At MC level, %f Lc --> K0S + p were found in the acceptance", fHistoMCLcK0SpGenAcc->GetEntries()));
437 AliInfo("fHistoMCLcK0SpGenAcc not available");
439 if(fVariablesTreeSgn) {
440 AliInfo(Form("At Reco level, %lld Lc --> K0S + p were found", fVariablesTreeSgn->GetEntries()));
442 AliInfo("fVariablesTreeSgn not available");
445 fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
447 AliError("fOutputKF not available");
454 //___________________________________________________________________________
455 void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
457 AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
461 fOutput = new TList();
463 fOutput->SetName("listTrees");
465 // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
466 const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
467 fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
468 fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
470 fCandidateVariables = new Float_t [nVar];
471 TString * fCandidateVariableNames = new TString[nVar];
472 fCandidateVariableNames[0]="massLc2K0Sp";
473 fCandidateVariableNames[1]="massLc2Lambdapi";
474 fCandidateVariableNames[2]="massK0S";
475 fCandidateVariableNames[3]="massLambda";
476 fCandidateVariableNames[4]="massLambdaBar";
477 fCandidateVariableNames[5]="cosPAK0S";
478 fCandidateVariableNames[6]="dcaV0";
479 fCandidateVariableNames[7]="tImpParBach";
480 fCandidateVariableNames[8]="tImpParV0";
481 fCandidateVariableNames[9]="nSigmaTPCpr";
482 fCandidateVariableNames[10]="nSigmaTPCpi";
483 fCandidateVariableNames[11]="nSigmaTPCka";
484 fCandidateVariableNames[12]="nSigmaTOFpr";
485 fCandidateVariableNames[13]="nSigmaTOFpi";
486 fCandidateVariableNames[14]="nSigmaTOFka";
487 fCandidateVariableNames[15]="bachelorPt";
488 fCandidateVariableNames[16]="V0positivePt";
489 fCandidateVariableNames[17]="V0negativePt";
490 fCandidateVariableNames[18]="dcaV0pos";
491 fCandidateVariableNames[19]="dcaV0neg";
492 fCandidateVariableNames[20]="v0Pt";
493 fCandidateVariableNames[21]="massGamma";
494 fCandidateVariableNames[22]="LcPt";
495 fCandidateVariableNames[23]="combinedProtonProb";
496 fCandidateVariableNames[24]="LcEta";
497 fCandidateVariableNames[25]="V0positiveEta";
498 fCandidateVariableNames[26]="V0negativeEta";
499 fCandidateVariableNames[27]="TPCProtonProb";
500 fCandidateVariableNames[28]="TOFProtonProb";
501 fCandidateVariableNames[29]="bachelorEta";
502 fCandidateVariableNames[30]="LcP";
503 fCandidateVariableNames[31]="bachelorP";
504 fCandidateVariableNames[32]="v0P";
505 fCandidateVariableNames[33]="V0positiveP";
506 fCandidateVariableNames[34]="V0negativeP";
507 fCandidateVariableNames[35]="LcY";
508 fCandidateVariableNames[36]="v0Y";
509 fCandidateVariableNames[37]="bachelorY";
510 fCandidateVariableNames[38]="V0positiveY";
511 fCandidateVariableNames[39]="V0negativeY";
512 fCandidateVariableNames[40]="v0Eta";
513 fCandidateVariableNames[41]="DecayLengthLc";
514 fCandidateVariableNames[42]="DecayLengthK0S";
515 fCandidateVariableNames[43]="CtLc";
516 fCandidateVariableNames[44]="CtK0S";
517 fCandidateVariableNames[45]="bachCode";
518 fCandidateVariableNames[46]="k0SCode";
520 fCandidateVariableNames[47]="V0KFmass";
521 fCandidateVariableNames[48]="V0KFdecayLength";
522 fCandidateVariableNames[49]="V0KFlifeTime";
524 fCandidateVariableNames[50]="V0KFmassErr";
525 fCandidateVariableNames[51]="V0KFdecayTimeErr";
526 fCandidateVariableNames[52]="V0KFlifeTimeErr";
528 fCandidateVariableNames[53]="LcKFmass";
529 fCandidateVariableNames[54]="LcKFdecayLength";
530 fCandidateVariableNames[55]="LcKFlifeTime";
532 fCandidateVariableNames[56]="LcKFmassErr";
533 fCandidateVariableNames[57]="LcKFdecayTimeErr";
534 fCandidateVariableNames[58]="LcKFlifeTimeErr";
536 fCandidateVariableNames[59]="LcKFDistToPrimVtx";
537 fCandidateVariableNames[60]="V0KFDistToPrimVtx";
538 fCandidateVariableNames[61]="V0KFDistToLc";
539 fCandidateVariableNames[62]="alphaArmKF";
540 fCandidateVariableNames[63]="ptArmKF";
541 fCandidateVariableNames[64]="alphaArm";
542 fCandidateVariableNames[65]="ptArm";
544 fCandidateVariableNames[66]="ITSrefitV0pos";
545 fCandidateVariableNames[67]="ITSrefitV0neg";
547 fCandidateVariableNames[68]="TPCClV0pos";
548 fCandidateVariableNames[69]="TPCClV0neg";
550 fCandidateVariableNames[70]="v0Xcoord";
551 fCandidateVariableNames[71]="v0Ycoord";
552 fCandidateVariableNames[72]="v0Zcoord";
553 fCandidateVariableNames[73]="primVtxX";
554 fCandidateVariableNames[74]="primVtxY";
555 fCandidateVariableNames[75]="primVtxZ";
557 fCandidateVariableNames[76]="ITSclBach";
558 fCandidateVariableNames[77]="SPDclBach";
560 fCandidateVariableNames[78]="ITSclV0pos";
561 fCandidateVariableNames[79]="SPDclV0pos";
562 fCandidateVariableNames[80]="ITSclV0neg";
563 fCandidateVariableNames[81]="SPDclV0neg";
565 fCandidateVariableNames[82]="alphaArmLc";
566 fCandidateVariableNames[83]="alphaArmLcCharge";
567 fCandidateVariableNames[84]="ptArmLc";
569 fCandidateVariableNames[85]="CosThetaStar";
572 for(Int_t ivar=0; ivar<nVar; ivar++){
573 fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
574 fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
577 fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
578 TString labelEv[2] = {"NotSelected", "Selected"};
579 for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
580 fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
583 fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
585 fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
586 TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
587 for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
588 fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
591 fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
592 TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
593 for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
594 fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
595 fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
598 fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
599 TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
600 for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
601 fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
604 fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
605 fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
607 TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
608 "kBachDifferentLambdaMother", "kBachCorrectLambdaMother"};
609 TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
610 "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
612 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
613 fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
614 fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
616 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
617 fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
618 fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
621 fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
622 for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
623 fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
626 fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 4, -0.5, 3.5);
627 TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"};
628 for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
629 fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
632 //fOutput->Add(fVariablesTreeSgn);
633 //fOutput->Add(fVariablesTreeBkg);
635 const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.};
637 fHistoMCLcK0SpGen = new TH1F("fHistoMCLcK0SpGen", "fHistoMCLcK0SpGen", 14, ptbins);
638 fHistoMCLcK0SpGenAcc = new TH1F("fHistoMCLcK0SpGenAcc", "fHistoMCLcK0SpGenAcc", 14, ptbins);
639 fHistoMCLcK0SpGenLimAcc = new TH1F("fHistoMCLcK0SpGenLimAcc", "fHistoMCLcK0SpGenLimAcc", 14, ptbins);
641 fOutput->Add(fHistoEvents);
642 fOutput->Add(fHistoLc);
643 fOutput->Add(fHistoLcOnTheFly);
644 fOutput->Add(fHistoLcBeforeCuts);
645 fOutput->Add(fHistoFiducialAcceptance);
646 fOutput->Add(fHistoCodesSgn);
647 fOutput->Add(fHistoCodesBkg);
648 fOutput->Add(fHistoLcpKpiBeforeCuts);
649 fOutput->Add(fHistoBackground);
650 fOutput->Add(fHistoMCLcK0SpGen);
651 fOutput->Add(fHistoMCLcK0SpGenAcc);
652 fOutput->Add(fHistoMCLcK0SpGenLimAcc);
654 PostData(1, fOutput);
655 PostData(4, fVariablesTreeSgn);
656 PostData(5, fVariablesTreeBkg);
657 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
658 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
659 fPIDResponse = inputHandler->GetPIDResponse();
661 if (fAnalCuts->GetIsUsePID()){
663 fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
664 fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
665 fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
666 fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
667 fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
668 fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
670 fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
673 // Setting properties of PID
674 fPIDCombined=new AliPIDCombined;
675 fPIDCombined->SetDefaultTPCPriors();
676 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
677 //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
679 fCounter = new AliNormalizationCounter("NormalizationCounter");
681 PostData(2, fCounter);
683 // Histograms from KF
685 if (fCallKFVertexing){
686 fHistoDistanceLcToPrimVtx = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
687 fHistoDistanceV0ToPrimVtx = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
688 fHistoDistanceV0ToLc = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
690 fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
691 fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
692 fHistoDistanceV0ToLcSgn = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
694 fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
695 fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
696 fHistoVtxV0ResidualToLc = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
698 fHistoMassV0All = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
699 fHistoDecayLengthV0All = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
700 fHistoLifeTimeV0All = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
702 fHistoMassV0True = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
703 fHistoDecayLengthV0True = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
704 fHistoLifeTimeV0True = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
706 fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
708 fHistoMassV0TrueK0S = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
709 fHistoDecayLengthV0TrueK0S = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
710 fHistoLifeTimeV0TrueK0S = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
712 fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
714 fHistoMassLcAll = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
715 fHistoDecayLengthLcAll = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
716 fHistoLifeTimeLcAll = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
718 fHistoMassLcTrue = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
719 fHistoDecayLengthLcTrue = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
720 fHistoLifeTimeLcTrue = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
722 fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
724 fHistoMassV0fromLcAll = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
725 fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
726 fHistoLifeTimeV0fromLcAll = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
728 fHistoMassV0fromLcTrue = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
729 fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
730 fHistoLifeTimeV0fromLcTrue = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
732 fHistoMassLcSgn = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
733 fHistoMassLcSgnFromAOD = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
734 fHistoDecayLengthLcSgn = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
735 fHistoLifeTimeLcSgn = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
737 fHistoMassV0fromLcSgn = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
738 fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
739 fHistoLifeTimeV0fromLcSgn = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
741 fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
742 fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
743 fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
744 TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
745 "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
746 "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
747 "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
749 for (Int_t ibin = 1; ibin <=16; ibin++){
750 fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
751 fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
752 fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
753 fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
756 fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
757 fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
758 fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
760 fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
761 fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
762 fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
764 fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
765 fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
766 fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
767 fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
770 fOutputKF = new TList();
771 fOutputKF->SetOwner();
772 fOutputKF->SetName("listHistoKF");
774 if (fCallKFVertexing){
775 fOutputKF->Add(fHistoDistanceLcToPrimVtx);
776 fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
777 fOutputKF->Add(fHistoDistanceV0ToLc);
779 fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
780 fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
781 fOutputKF->Add(fHistoDistanceV0ToLcSgn);
783 fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
784 fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
785 fOutputKF->Add(fHistoVtxV0ResidualToLc);
787 fOutputKF->Add(fHistoMassV0All);
788 fOutputKF->Add(fHistoDecayLengthV0All);
789 fOutputKF->Add(fHistoLifeTimeV0All);
791 fOutputKF->Add(fHistoMassV0True);
792 fOutputKF->Add(fHistoDecayLengthV0True);
793 fOutputKF->Add(fHistoLifeTimeV0True);
795 fOutputKF->Add(fHistoMassV0TrueFromAOD);
797 fOutputKF->Add(fHistoMassV0TrueK0S);
798 fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
799 fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
801 fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
803 fOutputKF->Add(fHistoMassLcAll);
804 fOutputKF->Add(fHistoDecayLengthLcAll);
805 fOutputKF->Add(fHistoLifeTimeLcAll);
807 fOutputKF->Add(fHistoMassLcTrue);
808 fOutputKF->Add(fHistoDecayLengthLcTrue);
809 fOutputKF->Add(fHistoLifeTimeLcTrue);
811 fOutputKF->Add(fHistoMassLcTrueFromAOD);
813 fOutputKF->Add(fHistoMassV0fromLcAll);
814 fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
815 fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
817 fOutputKF->Add(fHistoMassV0fromLcTrue);
818 fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
819 fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
821 fOutputKF->Add(fHistoMassLcSgn);
822 fOutputKF->Add(fHistoMassLcSgnFromAOD);
823 fOutputKF->Add(fHistoDecayLengthLcSgn);
824 fOutputKF->Add(fHistoLifeTimeLcSgn);
826 fOutputKF->Add(fHistoMassV0fromLcSgn);
827 fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
828 fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
830 fOutputKF->Add(fHistoKF);
831 fOutputKF->Add(fHistoKFV0);
832 fOutputKF->Add(fHistoKFLc);
834 fOutputKF->Add(fHistoMassKFV0);
835 fOutputKF->Add(fHistoDecayLengthKFV0);
836 fOutputKF->Add(fHistoLifeTimeKFV0);
838 fOutputKF->Add(fHistoMassKFLc);
839 fOutputKF->Add(fHistoDecayLengthKFLc);
840 fOutputKF->Add(fHistoLifeTimeKFLc);
842 fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
843 fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
844 fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
845 fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
848 PostData(6, fOutputKF);
853 //_________________________________________________
854 void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
858 AliError("NO EVENT FOUND!");
863 AliDebug(2, Form("Processing event = %d", fCurrentEvent));
864 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
865 TClonesArray *arrayLctopKos=0;
867 TClonesArray *array3Prong = 0;
869 if (!aodEvent && AODEvent() && IsStandardAOD()) {
870 // In case there is an AOD handler writing a standard AOD, use the AOD
871 // event in memory rather than the input (ESD) event.
872 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
873 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
874 // have to taken from the AOD event hold by the AliAODExtension
875 AliAODHandler* aodHandler = (AliAODHandler*)
876 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
878 if (aodHandler->GetExtensions()) {
879 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
880 AliAODEvent *aodFromExt = ext->GetAOD();
881 arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
883 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
886 arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
888 array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
891 if ( !fUseMCInfo && fIspA) {
892 fAnalCuts->SetTriggerClass("");
893 fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
896 Int_t runnumber = aodEvent->GetRunNumber();
897 if (aodEvent->GetTriggerMask() == 0 && (runnumber >= 195344 && runnumber <= 195677)){
898 AliDebug(3,"Event rejected because of null trigger mask");
902 fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
905 TClonesArray *mcArray = 0;
906 AliAODMCHeader *mcHeader=0;
909 // MC array need for matching
910 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
912 AliError("Could not find Monte-Carlo in AOD");
916 mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
918 AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
921 //Printf("Filling MC histo");
922 FillMCHisto(mcArray);
925 // AOD primary vertex
926 fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
928 if (fVtx1->GetNContributors()<1) return;
930 fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
932 if ( !fIsEventSelected ) {
933 fHistoEvents->Fill(0);
934 return; // don't take into account not selected events
936 fHistoEvents->Fill(1);
938 // Setting magnetic field for KF vertexing
939 fBField = aodEvent->GetMagneticField();
940 AliKFParticle::SetField(fBField);
942 Int_t nSelectedAnal = 0;
943 if (fIsK0sAnalysis) {
944 MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
945 nSelectedAnal, fAnalCuts,
946 array3Prong, mcHeader);
948 fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
950 PostData(1, fOutput);
951 PostData(2, fCounter);
952 PostData(4, fVariablesTreeSgn);
953 PostData(5, fVariablesTreeBkg);
954 PostData(6, fOutputKF);
957 //-------------------------------------------------------------------------------
958 void AliAnalysisTaskSELc2V0bachelorTMVA::FillMCHisto(TClonesArray *mcArray){
960 // method to fill MC histo: how many Lc --> K0S + p are there at MC level
961 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
962 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
964 AliError("Failed casting particle from MC array!, Skipping particle");
967 Int_t pdg = mcPart->GetPdgCode();
968 if (TMath::Abs(pdg) != 4122){
969 AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg));
972 AliDebug(2, Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg));
973 Int_t labeldaugh0 = mcPart->GetDaughter(0);
974 Int_t labeldaugh1 = mcPart->GetDaughter(1);
975 if (labeldaugh0 <= 0 || labeldaugh1 <= 0){
976 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
979 else if (labeldaugh1 - labeldaugh0 == 1){
980 AliDebug(2, Form("Step 1 ok: The MC particle has correct daughters!!"));
981 AliAODMCParticle* daugh0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh0));
982 AliAODMCParticle* daugh1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh1));
983 if(!daugh0 || !daugh1){
984 AliDebug(2,"Particle daughters not properly retrieved!");
987 Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode());
988 Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode());
989 AliAODMCParticle* bachelorMC = daugh0;
990 AliAODMCParticle* v0MC = daugh1;
991 AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1));
992 if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){
993 // 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-
994 /// first, we set the bachelor and the v0: above we assumed first proton and second V0, but we could have to change it:
995 if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) {
999 AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters()));
1000 if (v0MC->GetNDaughters() != 1) {
1001 AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing...");
1004 else { // So far: Lc --> K0 + p, K0 with 1 daughter
1005 AliDebug(2, "Step 2 ok: The K0 does decay in 1 body only! ");
1006 Int_t labelK0daugh = v0MC->GetDaughter(0);
1007 AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0daugh));
1009 AliError("Error while casting particle! returning a NULL array");
1012 else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access
1013 if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){
1014 AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies");
1017 else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies
1018 AliDebug(2, "Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies");
1019 Int_t labelK0Sdaugh0 = partK0S->GetDaughter(0);
1020 Int_t labelK0Sdaugh1 = partK0S->GetDaughter(1);
1021 AliAODMCParticle* daughK0S0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh0));
1022 AliAODMCParticle* daughK0S1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh1));
1023 if (!daughK0S0 || ! daughK0S1){
1024 AliDebug(2, "Could not access K0S daughters, continuing...");
1027 else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access
1028 AliDebug(2, "Step 4 ok: Could access K0S daughters, continuing...");
1029 Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode();
1030 Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode();
1031 if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){
1032 AliDebug(2, "The K0S does not decay in pi+pi-, continuing");
1033 //AliInfo("The K0S does not decay in pi+pi-, continuing");
1035 else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi-
1036 if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) {
1037 AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt()));
1038 if(TMath::Abs(mcPart->Y()) < 0.5) fHistoMCLcK0SpGenLimAcc->Fill(mcPart->Pt());
1039 //AliInfo(Form("\nparticle = %d, Filling MC Gen histo\n", iPart));
1040 fHistoMCLcK0SpGen->Fill(mcPart->Pt());
1041 if(!(TMath::Abs(bachelorMC->Eta()) > 0.9 || bachelorMC->Pt() < 0.1 ||
1042 TMath::Abs(daughK0S0->Eta()) > 0.9 || daughK0S0->Pt() < 0.1 ||
1043 TMath::Abs(daughK0S1->Eta()) > 0.9 || daughK0S1->Pt() < 0.1)) {
1044 fHistoMCLcK0SpGenAcc->Fill(mcPart->Pt());
1048 AliDebug(2, "not in fiducial acceptance! Skipping");
1058 } // closing loop over mcArray
1064 //-------------------------------------------------------------------------------
1065 void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
1066 TClonesArray *mcArray,
1067 Int_t &nSelectedAnal,
1068 AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
1069 AliAODMCHeader* aodheader){
1070 //Lc prong needed to MatchToMC method
1072 Int_t pdgCand = 4122;
1073 Int_t pdgDgLctoV0bachelor[2]={2212, 310};
1074 Int_t pdgDgV0toDaughters[2]={211, 211};
1076 Int_t pdgDgLctopKpi[3]={2212, 321, 211};
1078 // loop to search for candidates Lc->p+K+pi
1079 Int_t n3Prong = array3Prong->GetEntriesFast();
1080 Int_t nCascades= arrayLctopKos->GetEntriesFast();
1082 //AliInfo(Form("\n\n\n\n3 prong candidates = %d, ncascades = %d \n\n\n\n\n", n3Prong, nCascades));
1083 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
1084 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
1085 //Filling a control histogram with no cuts
1088 // find associated MC particle for Lc -> p+K+pi
1089 Int_t mcLabel = d->MatchToMC(4122, mcArray, 3, pdgDgLctopKpi);
1090 //Int_t mcLabelTemp = d->MatchToMC(4122, mcArray);
1091 //Printf("mcLabel = %d, mcLabelTemp = %d", mcLabel, mcLabelTemp);
1094 AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1096 Int_t pdgCode = partLcpKpi->GetPdgCode();
1097 AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1098 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1099 fHistoLcpKpiBeforeCuts->Fill(1);
1104 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
1105 fHistoLcpKpiBeforeCuts->Fill(0);
1110 // loop over cascades to search for candidates Lc->p+K0S
1112 for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
1114 // Lc candidates and K0s from Lc
1115 AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
1117 AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
1121 //Filling a control histogram with no cuts
1126 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1127 fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1128 if (fmcLabelLc>=0) {
1129 AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades));
1131 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
1133 pdgCode = partLc->GetPdgCode();
1134 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
1135 pdgCode = TMath::Abs(pdgCode);
1136 fHistoLcBeforeCuts->Fill(1);
1141 fHistoLcBeforeCuts->Fill(0);
1146 //if (!lcK0spr->GetSecondaryVtx()) {
1147 // AliInfo("No secondary vertex");
1151 if (lcK0spr->GetNDaughters()!=2) {
1152 AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
1156 AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
1157 AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
1158 if (!v0part || !bachPart) {
1159 AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0s));
1164 if (!v0part->GetSecondaryVtx()) {
1165 AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
1169 if (v0part->GetNDaughters()!=2) {
1170 AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
1174 AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
1175 AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
1176 if (!v0Neg || !v0Pos) {
1177 AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
1182 if (v0Pos->Charge() == v0Neg->Charge()) {
1183 AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
1193 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1194 Int_t mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1196 AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1198 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1200 pdgCode = partLc->GetPdgCode();
1201 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1202 pdgCode = TMath::Abs(pdgCode);
1213 AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
1214 AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1216 if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1219 else { // checking if we want to fill the background
1220 if (fKeepingOnlyHIJINGBkg){
1221 // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1222 Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
1223 if (!isCandidateInjected){
1224 AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1225 fHistoBackground->Fill(1);
1228 AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1229 fHistoBackground->Fill(0);
1233 else if (fKeepingOnlyPYTHIABkg){
1234 // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1235 AliAODTrack *bachelor = (AliAODTrack*)lcK0spr->GetBachelor();
1236 AliAODTrack *v0pos = (AliAODTrack*)lcK0spr->Getv0PositiveTrack();
1237 AliAODTrack *v0neg = (AliAODTrack*)lcK0spr->Getv0NegativeTrack();
1238 if (!bachelor || !v0pos || !v0neg) {
1239 AliDebug(2, "Cannot retrieve one of the tracks while checking origin, continuing");
1243 Int_t labelbachelor = TMath::Abs(bachelor->GetLabel());
1244 Int_t labelv0pos = TMath::Abs(v0pos->GetLabel());
1245 Int_t labelv0neg = TMath::Abs(v0neg->GetLabel());
1246 AliAODMCParticle* MCbachelor = (AliAODMCParticle*)mcArray->At(labelbachelor);
1247 AliAODMCParticle* MCv0pos = (AliAODMCParticle*)mcArray->At(labelv0pos);
1248 AliAODMCParticle* MCv0neg = (AliAODMCParticle*)mcArray->At(labelv0neg);
1249 if (!MCbachelor || !MCv0pos || !MCv0neg) {
1250 AliDebug(2, "Cannot retrieve MC particle for one of the tracks while checking origin, continuing");
1254 Int_t isBachelorFromPythia = fUtils->CheckOrigin(mcArray, MCbachelor, kTRUE);
1255 Int_t isv0posFromPythia = fUtils->CheckOrigin(mcArray, MCv0pos, kTRUE);
1256 Int_t isv0negFromPythia = fUtils->CheckOrigin(mcArray, MCv0neg, kTRUE);
1257 if (isBachelorFromPythia != 0 && isv0posFromPythia != 0 && isv0negFromPythia != 0){
1258 AliDebug(2, "The candidate is from PYTHIA (i.e. all daughters originate from a quark), keeping it to fill background");
1259 fHistoBackground->Fill(2);
1262 AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background");
1263 fHistoBackground->Fill(3);
1272 FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray, iLctopK0s);
1279 //________________________________________________________________________
1280 void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
1282 Int_t &nSelectedAnal,
1283 AliRDHFCutsLctoV0 *cutsAnal,
1284 TClonesArray *mcArray, Int_t iLctopK0s){
1286 // Fill histos for Lc -> K0S+proton
1290 if (!part->GetOwnPrimaryVtx()) {
1291 //Printf("No primary vertex for Lc found!!");
1292 part->SetOwnPrimaryVtx(fVtx1);
1295 //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1298 Double_t invmassLc = part->InvMassLctoK0sP();
1300 AliAODv0 * v0part = part->Getv0();
1301 Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1302 if (onFlyV0){ // on-the-fly V0
1304 fHistoLcOnTheFly->Fill(2.);
1307 fHistoLcOnTheFly->Fill(0.);
1310 else { // offline V0
1312 fHistoLcOnTheFly->Fill(3.);
1315 fHistoLcOnTheFly->Fill(1.);
1319 Double_t dcaV0 = v0part->GetDCA();
1320 Double_t invmassK0s = v0part->MassK0Short();
1322 if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
1324 fHistoFiducialAcceptance->Fill(3.);
1327 fHistoFiducialAcceptance->Fill(1.);
1332 fHistoFiducialAcceptance->Fill(2.);
1335 fHistoFiducialAcceptance->Fill(0.);
1339 Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1341 if (isInV0window == 0) {
1342 AliDebug(2, "No: The candidate has NOT passed the V0 window cuts!");
1343 if (isLc) Printf("SIGNAL candidate rejected: V0 window cuts");
1346 else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1348 Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1350 if (!isInCascadeWindow) {
1351 AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1352 if (isLc) Printf("SIGNAL candidate rejected: cascade window cuts");
1355 else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1357 Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1358 AliDebug(2, Form("recoAnalysisCuts = %d", cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate) & (AliRDHFCutsLctoV0::kLcToK0Spr)));
1359 if (!isCandidateSelectedCuts){
1360 AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1361 if (isLc) Printf("SIGNAL candidate rejected");
1365 AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1368 AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1370 AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1374 //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1375 Double_t probTPCTOF[AliPID::kSPECIES]={-1.};
1377 UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1378 AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed));
1379 Double_t probProton = -1.;
1380 // Double_t probPion = -1.;
1381 // Double_t probKaon = -1.;
1382 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1383 AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1384 probProton = probTPCTOF[AliPID::kProton];
1385 // probPion = probTPCTOF[AliPID::kPion];
1386 // probKaon = probTPCTOF[AliPID::kKaon];
1388 else { // if you don't have both TOF and TPC, try only TPC
1389 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1390 AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC");
1391 detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1392 AliDebug(2,Form(" detUsed (TPC case) = %d", detUsed));
1393 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1394 probProton = probTPCTOF[AliPID::kProton];
1395 // probPion = probTPCTOF[AliPID::kPion];
1396 // probKaon = probTPCTOF[AliPID::kKaon];
1397 AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1400 AliDebug(2, "Only TPC did not work...");
1402 // resetting mask to ask for both TPC+TOF
1403 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1405 AliDebug(2, Form("probProton = %f", probProton));
1407 // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) )
1408 Double_t probProtonTPC = -1.;
1409 Double_t probProtonTOF = -1.;
1410 Double_t pidTPC[AliPID::kSPECIES]={-1.};
1411 Double_t pidTOF[AliPID::kSPECIES]={-1.};
1412 Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC);
1413 Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF);
1414 if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton];
1415 if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton];
1417 // checking V0 status (on-the-fly vs offline)
1418 if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1419 AliDebug(2, "On-the-fly discarded");
1423 if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
1427 if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
1429 if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1430 if (isLc) Printf("SIGNAL candidate rejected");
1431 AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1435 AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1438 Int_t pdgCand = 4122;
1439 Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
1440 Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
1441 Int_t isLc2LBarpi=0, isLc2Lpi=0;
1442 Int_t currentLabel = part->GetLabel();
1445 mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1447 if (bachelor->Charge()==-1) isLc2LBarpi=1;
1448 if (bachelor->Charge()==+1) isLc2Lpi=1;
1452 Int_t pdgDg2prong[2] = {211, 211};
1456 labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1457 if (labelK0S>=0) isK0S = 1;
1460 pdgDg2prong[0] = 211;
1461 pdgDg2prong[1] = 2212;
1463 Int_t isLambdaBar = 0;
1464 Int_t lambdaLabel = 0;
1466 lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1467 if (lambdaLabel>=0) {
1468 AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1469 if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1470 else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1474 pdgDg2prong[0] = 11;
1475 pdgDg2prong[1] = 11;
1477 Int_t gammaLabel = 0;
1479 gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1480 if (gammaLabel>=0) {
1481 AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1482 if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1487 if (currentLabel != -1){
1488 AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1489 pdgTemp = tempPart->GetPdgCode();
1491 if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1492 else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1493 else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1494 else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1495 if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1496 else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
1497 else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
1498 else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
1499 //else AliDebug(2, Form("V0 is something else!!"));
1501 Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1502 Double_t invmassLambda = v0part->MassLambda();
1503 Double_t invmassLambdaBar = v0part->MassAntiLambda();
1505 //Double_t nSigmaITSpr=-999.;
1506 Double_t nSigmaTPCpr=-999.;
1507 Double_t nSigmaTOFpr=-999.;
1509 //Double_t nSigmaITSpi=-999.;
1510 Double_t nSigmaTPCpi=-999.;
1511 Double_t nSigmaTOFpi=-999.;
1513 //Double_t nSigmaITSka=-999.;
1514 Double_t nSigmaTPCka=-999.;
1515 Double_t nSigmaTOFka=-999.;
1518 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1519 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1520 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1521 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1522 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1523 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1524 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1525 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1526 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1529 nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
1530 nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
1531 nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
1532 nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
1533 nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
1534 nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
1537 // Fill candidate variable Tree (track selection, V0 invMass selection)
1538 if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
1540 fCandidateVariables[0] = invmassLc;
1541 fCandidateVariables[1] = invmassLc2Lpi;
1542 fCandidateVariables[2] = invmassK0s;
1543 fCandidateVariables[3] = invmassLambda;
1544 fCandidateVariables[4] = invmassLambdaBar;
1545 fCandidateVariables[5] = part->CosV0PointingAngle();
1546 fCandidateVariables[6] = dcaV0;
1547 fCandidateVariables[7] = part->Getd0Prong(0);
1548 fCandidateVariables[8] = part->Getd0Prong(1);
1549 fCandidateVariables[9] = nSigmaTPCpr;
1550 fCandidateVariables[10] = nSigmaTPCpi;
1551 fCandidateVariables[11] = nSigmaTPCka;
1552 fCandidateVariables[12] = nSigmaTOFpr;
1553 fCandidateVariables[13] = nSigmaTOFpi;
1554 fCandidateVariables[14] = nSigmaTOFka;
1555 fCandidateVariables[15] = bachelor->Pt();
1556 AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1557 fCandidateVariables[16] = v0pos->Pt();
1558 AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
1559 fCandidateVariables[17] = v0neg->Pt();
1560 fCandidateVariables[18] = v0part->Getd0Prong(0);
1561 fCandidateVariables[19] = v0part->Getd0Prong(1);
1562 fCandidateVariables[20] = v0part->Pt();
1563 fCandidateVariables[21] = v0part->InvMass2Prongs(0,1,11,11);
1564 fCandidateVariables[22] = part->Pt();
1565 fCandidateVariables[23] = probProton;
1566 fCandidateVariables[24] = part->Eta();
1567 fCandidateVariables[25] = v0pos->Eta();
1568 fCandidateVariables[26] = v0neg->Eta();
1569 fCandidateVariables[27] = probProtonTPC;
1570 fCandidateVariables[28] = probProtonTOF;
1571 fCandidateVariables[29] = bachelor->Eta();
1573 fCandidateVariables[30] = part->P();
1574 fCandidateVariables[31] = bachelor->P();
1575 fCandidateVariables[32] = v0part->P();
1576 fCandidateVariables[33] = v0pos->P();
1577 fCandidateVariables[34] = v0neg->P();
1579 fCandidateVariables[35] = part->Y(4122);
1580 fCandidateVariables[36] = bachelor->Y(2212);
1581 fCandidateVariables[37] = v0part->Y(310);
1582 fCandidateVariables[38] = v0pos->Y(211);
1583 fCandidateVariables[39] = v0neg->Y(211);
1585 fCandidateVariables[40] = v0part->Eta();
1587 fCandidateVariables[41] = part->DecayLength();
1588 fCandidateVariables[42] = part->DecayLengthV0();
1589 fCandidateVariables[43] = part->Ct(4122);
1590 fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
1592 EBachelor bachCode = kBachInvalid;
1593 EK0S k0SCode = kK0SInvalid;
1595 bachCode = CheckBachelor(part, bachelor, mcArray);
1596 k0SCode = CheckK0S(part, v0part, mcArray);
1599 fCandidateVariables[45] = bachCode;
1600 fCandidateVariables[46] = k0SCode;
1602 Double_t V0KF[3] = {-999999, -999999, -999999};
1603 Double_t errV0KF[3] = {-999999, -999999, -999999};
1604 Double_t LcKF[3] = {-999999, -999999, -999999};
1605 Double_t errLcKF[3] = {-999999, -999999, -999999};
1606 Double_t distances[3] = {-999999, -999999, -999999};
1607 Double_t armPolKF[2] = {-999999, -999999};
1609 if (fCallKFVertexing){
1610 Int_t kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
1611 AliDebug(2, Form("Result from KF = %d", kfResult));
1615 for (Int_t i = 0; i< 3; i++){
1616 Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
1620 fCandidateVariables[47] = V0KF[0];
1621 fCandidateVariables[48] = V0KF[1];
1622 fCandidateVariables[49] = V0KF[2];
1624 fCandidateVariables[50] = errV0KF[0];
1625 fCandidateVariables[51] = errV0KF[1];
1626 fCandidateVariables[52] = errV0KF[2];
1628 fCandidateVariables[53] = LcKF[0];
1629 fCandidateVariables[54] = LcKF[1];
1630 fCandidateVariables[55] = LcKF[2];
1632 fCandidateVariables[56] = errLcKF[0];
1633 fCandidateVariables[57] = errLcKF[1];
1634 fCandidateVariables[58] = errLcKF[2];
1636 fCandidateVariables[59] = distances[0];
1637 fCandidateVariables[60] = distances[1];
1638 fCandidateVariables[61] = distances[2];
1639 fCandidateVariables[62] = armPolKF[0];
1640 fCandidateVariables[63] = armPolKF[1];
1641 fCandidateVariables[64] = v0part->AlphaV0();
1642 fCandidateVariables[65] = v0part->PtArmV0();
1644 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)));
1645 fCandidateVariables[66] = v0pos->GetStatus() & AliESDtrack::kITSrefit;
1646 fCandidateVariables[67] = v0neg->GetStatus() & AliESDtrack::kITSrefit;
1647 fCandidateVariables[68] = v0pos->GetTPCClusterInfo(2, 1);
1648 fCandidateVariables[69] = v0neg->GetTPCClusterInfo(2, 1);
1650 fCandidateVariables[70] = v0part->Xv();
1651 fCandidateVariables[71] = v0part->Yv();
1652 fCandidateVariables[72] = v0part->Zv();
1654 fCandidateVariables[73] = fVtx1->GetX();
1655 fCandidateVariables[74] = fVtx1->GetY();
1656 fCandidateVariables[75] = fVtx1->GetZ();
1658 fCandidateVariables[76] = bachelor->GetITSNcls();
1659 fCandidateVariables[77] = bachelor->HasPointOnITSLayer(0) + bachelor->HasPointOnITSLayer(1);
1661 fCandidateVariables[78] = v0pos->GetITSNcls();
1662 fCandidateVariables[79] = v0pos->HasPointOnITSLayer(0) + v0pos->HasPointOnITSLayer(1);
1664 fCandidateVariables[80] = v0neg->GetITSNcls();
1665 fCandidateVariables[81] = v0neg->HasPointOnITSLayer(0) + v0neg->HasPointOnITSLayer(1);
1667 TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz());
1668 TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz());
1669 TVector3 momTot(part->Px(), part->Py(), part->Pz());
1671 Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
1672 Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
1674 Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2);
1675 Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) );
1676 Double_t ptArmLc = mom1.Perp(momTot);
1678 fCandidateVariables[82] = alphaArmLc;
1679 fCandidateVariables[83] = alphaArmLcCharge;
1680 fCandidateVariables[84] = ptArmLc;
1682 Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // mass K0S PDG
1683 Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); // mass Proton PDG
1684 Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass(); // mass Lc PDG
1686 Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG);
1687 Double_t e = part->E(4122);
1688 Double_t beta = part->P()/e;
1689 Double_t gamma = e/massLcPDG;
1691 Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
1693 fCandidateVariables[85] = cts;
1697 AliDebug(2, Form("Reco particle %d --> Filling Sgn", iLctopK0s));
1698 fVariablesTreeSgn->Fill();
1699 fHistoCodesSgn->Fill(bachCode, k0SCode);
1702 if (fFillOnlySgn == kFALSE){
1703 AliDebug(2, "Filling Bkg");
1704 fVariablesTreeBkg->Fill();
1705 fHistoCodesBkg->Fill(bachCode, k0SCode);
1710 fVariablesTreeSgn->Fill();
1718 //________________________________________________________________________
1719 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
1720 Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
1721 Double_t* distances, Double_t* armPolKF) {
1724 // method to perform KF vertexing
1725 // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
1728 Int_t codeKFV0 = -1, codeKFLc = -1;
1730 AliKFVertex primVtxCopy;
1731 Int_t nt = 0, ntcheck = 0;
1732 Double_t pos[3] = {0., 0., 0.};
1735 Int_t contr = fVtx1->GetNContributors();
1736 Double_t covmatrix[6] = {0.};
1737 fVtx1->GetCovarianceMatrix(covmatrix);
1738 Double_t chi2 = fVtx1->GetChi2();
1739 AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
1741 // topological constraint
1742 primVtxCopy = AliKFVertex(primaryESDVtxCopy);
1743 nt = primaryESDVtxCopy.GetNContributors();
1746 Int_t pdg[2] = {211, -211};
1747 Int_t pdgLc[2] = {2212, 310};
1749 Int_t pdgDgV0toDaughters[2] = {211, 211};
1751 Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
1753 // the KF vertex for the V0 has to be built with the prongs of the V0!
1754 Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
1755 AliKFParticle V0, positiveV0KF, negativeV0KF;
1756 Int_t labelsv0daugh[2] = {-1, -1};
1757 Int_t idv0daugh[2] = {-1, -1};
1758 AliExternalTrackParam* esdv0Daugh1 = 0x0;
1759 AliExternalTrackParam* esdv0Daugh2 = 0x0;
1760 for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1761 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1763 AliDebug(2, "No V0 daughters available");
1766 Double_t xyz[3], pxpypz[3], cv[21];
1768 aodTrack->GetXYZ(xyz);
1769 aodTrack->PxPyPz(pxpypz);
1770 aodTrack->GetCovarianceXYZPxPyPz(cv);
1771 sign = aodTrack->Charge();
1772 AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
1774 if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1775 else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1776 labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
1777 idv0daugh[ipr] = aodTrack->GetID();
1778 if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
1780 //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
1782 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1783 if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
1784 positiveV0KF = daughterKF;
1787 negativeV0KF = daughterKF;
1791 Double_t xn=0., xp=0.;//, dca;
1792 AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
1793 // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
1795 AliExternalTrackParam tr1(*esdv0Daugh1);
1796 AliExternalTrackParam tr2(*esdv0Daugh2);
1797 tr1.PropagateTo(xn, fBField);
1798 tr2.PropagateTo(xp, fBField);
1800 AliKFParticle daughterKF1(tr1, 211);
1801 AliKFParticle daughterKF2(tr2, 211);
1802 V0.AddDaughter(positiveV0KF);
1803 V0.AddDaughter(negativeV0KF);
1804 //V0.AddDaughter(daughterKF1);
1805 //V0.AddDaughter(daughterKF2);
1811 // Checking the quality of the KF V0 vertex
1812 if( V0.GetNDF() < 1 ) {
1813 //Printf("Number of degrees of freedom < 1, continuing");
1816 if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
1817 //Printf("Chi2 per DOF too big, continuing");
1821 if(ftopoConstraint && nt > 0){
1822 for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1823 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1824 //* subtruct daughters from primary vertex
1825 if(!aodTrack->GetUsedForPrimVtxFit()) {
1826 //Printf("Track %d was not used for primary vertex, continuing", i);
1829 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1830 primVtxCopy -= daughterKF;
1835 // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
1837 if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
1838 //Printf("Deviation from vertex too big, continuing");
1843 //* Get V0 invariant mass
1844 Double_t massV0 = 999999, sigmaMassV0 = 999999;
1845 Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
1848 codeKFV0 = 1; // Mass not ok
1849 if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
1851 else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
1853 fHistoMassKFV0->Fill(massV0, sigmaMassV0);
1855 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]);
1856 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]);
1858 Printf("Vertices: KF: x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ());
1859 Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv());
1861 //Printf("Got MC vtx for V0");
1862 if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) {
1863 AliAODMCParticle* tmpdaughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
1864 AliAODMCParticle* tmpdaughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
1865 if (!tmpdaughv01 && labelsv0daugh[0] > 0){
1866 AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
1868 if (!tmpdaughv02 && labelsv0daugh[1] > 0){
1869 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1872 Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays
1873 Double_t yPionMC = tmpdaughv01->Yv();
1874 Double_t zPionMC = tmpdaughv01->Zv();
1875 //Printf("Got MC vtx for Pion");
1876 Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
1880 Printf("Not a true V0");
1882 //massV0=-1;//return -1;// !!!!
1884 // now use what we just try with the bachelor, to build the Lc
1886 // topological constraint
1887 nt = primVtxCopy.GetNContributors();
1890 Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
1892 Int_t labelsLcdaugh[2] = {-1, -1};
1893 labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
1894 labelsLcdaugh[1] = mcLabelV0;
1896 if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
1897 AliKFParticle daughterKFLc(*bach, pdgLc[0]);
1898 Lc.AddDaughter(daughterKFLc);
1899 TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310);
1900 Double_t massPDGK0S = particlePDG->Mass();
1901 V0.SetMassConstraint(massPDGK0S);
1903 if( Lc.GetNDF() < 1 ) {
1904 AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
1907 if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1908 AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1912 if(ftopoConstraint && nt > 0){
1913 //* subtruct daughters from primary vertex
1914 if(!bach->GetUsedForPrimVtxFit()) {
1915 AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
1918 primVtxCopy -= daughterKFLc;
1921 /* the V0 was added above, so it is ok to remove it without checking
1922 if(!V0->GetUsedForPrimVtxFit()) {
1923 Printf("Lc: V0 was not used for primary vertex, continuing");
1927 //primVtxCopy -= V0;
1931 // Check Lc Chi^2 deviation from primary vertex
1933 if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
1934 AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
1938 if(ftopoConstraint){
1940 // Add Lc to primary vertex to improve the primary vertex resolution
1942 Lc.SetProductionVertex(primVtxCopy);
1947 if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1948 AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1952 if(ftopoConstraint){
1953 V0.SetProductionVertex(Lc);
1956 // After setting the vertex of the V0, getting/filling some info
1958 //* Get V0 decayLength
1959 Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
1960 Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
1962 if (sigmaDecayLengthV0 > 1e19) {
1963 codeKFV0 = 3; // DecayLength not ok
1965 codeKFV0 = 6; // DecayLength and Mass not ok
1966 if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
1968 else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
1971 fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
1973 //* Get V0 life time
1974 Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
1975 Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
1977 if (sigmaLifeTimeV0 > 1e19) {
1978 codeKFV0 = 4; // LifeTime not ok
1979 if (sigmaDecayLengthV0 > 1e19) {
1980 codeKFV0 = 9; // LifeTime and DecayLength not ok
1982 codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
1983 if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
1985 else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
1987 else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
1988 codeKFV0 = 7; // LifeTime and Mass not ok
1989 if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
1991 else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
1994 fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
1996 if (codeKFV0 == -1) codeKFV0 = 0;
1997 fHistoKFV0->Fill(codeKFV0);
1999 AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
2001 fHistoMassV0All->Fill(massV0);
2002 fHistoDecayLengthV0All->Fill(decayLengthV0);
2003 fHistoLifeTimeV0All->Fill(lifeTimeV0);
2005 Double_t qtAlphaV0[2] = {0., 0.};
2006 Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
2007 positiveV0KF.TransportToPoint(vtxV0KF);
2008 negativeV0KF.TransportToPoint(vtxV0KF);
2009 V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
2010 AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
2011 fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2012 fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2013 armPolKF[0] = qtAlphaV0[1];
2014 armPolKF[1] = qtAlphaV0[0];
2016 // Checking MC info for V0
2018 AliAODMCParticle *motherV0 = 0x0;
2019 AliAODMCParticle *daughv01 = 0x0;
2020 AliAODMCParticle *daughv02 = 0x0;
2023 daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
2024 daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
2025 if (!daughv01 && labelsv0daugh[0] > 0){
2026 AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
2029 if (!daughv02 && labelsv0daugh[1] > 0){
2030 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
2034 if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
2035 AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
2036 motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
2037 if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
2038 if( motherV0->GetNDaughters() == 2 ){
2039 fHistoMassV0True->Fill(massV0);
2040 fHistoDecayLengthV0True->Fill(decayLengthV0);
2041 fHistoLifeTimeV0True->Fill(lifeTimeV0);
2042 fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
2043 if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
2044 fHistoMassV0TrueK0S->Fill(massV0);
2045 fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
2046 fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
2047 fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
2050 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()));
2052 else if (!motherV0){
2053 AliDebug(3, "could not access MC info for mother, continuing");
2056 AliDebug(3, "MC mother is a gluon, continuing");
2060 AliDebug(3, "Background V0!");
2066 AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
2070 //* Get Lc invariant mass
2071 Double_t massLc = 999999, sigmaMassLc= 999999;
2072 Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
2074 AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
2076 codeKFLc = 1; // Mass not ok
2077 if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
2079 else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
2081 fHistoMassKFLc->Fill(massLc, sigmaMassLc);
2083 //* Get Lc Decay length
2084 Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
2085 Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
2087 AliDebug(3, "----> Lc: Could not get decay length, and sigma");
2088 if (sigmaDecayLengthLc > 1e19) {
2089 codeKFLc = 3; // DecayLength not ok
2091 codeKFLc = 6; // DecayLength and Mass not ok
2092 if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
2094 else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
2097 AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
2099 fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
2101 //* Get Lc life time
2102 Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
2103 Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
2105 AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
2106 if (sigmaLifeTimeLc > 1e19) {
2107 codeKFLc = 4; // LifeTime not ok
2108 if (sigmaDecayLengthLc > 1e19) {
2109 codeKFLc = 9; // LifeTime and DecayLength not ok
2111 codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
2112 if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2114 else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
2116 else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
2117 codeKFLc = 7; // LifeTime and Mass not ok
2118 if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
2120 else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
2124 fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
2126 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));
2128 if (codeKFLc == -1) codeKFLc = 0;
2129 fHistoKFLc->Fill(codeKFLc);
2131 fHistoKF->Fill(codeKFV0, codeKFLc);
2133 // here we fill the histgrams for all the reconstructed KF vertices for the cascade
2134 fHistoMassLcAll->Fill(massLc);
2135 fHistoDecayLengthLcAll->Fill(decayLengthLc);
2136 fHistoLifeTimeLcAll->Fill(lifeTimeLc);
2138 fHistoMassV0fromLcAll->Fill(massV0);
2139 fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
2140 fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
2142 Double_t xV0 = V0.GetX();
2143 Double_t yV0 = V0.GetY();
2144 Double_t zV0 = V0.GetZ();
2146 Double_t xLc = Lc.GetX();
2147 Double_t yLc = Lc.GetY();
2148 Double_t zLc = Lc.GetZ();
2150 Double_t xPrimVtx = primVtxCopy.GetX();
2151 Double_t yPrimVtx = primVtxCopy.GetY();
2152 Double_t zPrimVtx = primVtxCopy.GetZ();
2154 Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2155 (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2156 (zPrimVtx - zLc) * (zPrimVtx - zLc));
2158 Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2159 (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2160 (zPrimVtx - zV0) * (zPrimVtx - zV0));
2162 Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2163 (yLc - yV0)*(yLc - yV0) +
2164 (zLc - zV0)*(zLc - zV0));
2166 //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2168 fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2169 fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2170 fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2172 distances[0] = distanceLcToPrimVtx;
2173 distances[1] = distanceV0ToPrimVtx;
2174 distances[2] = distanceV0ToLc;
2178 AliAODMCParticle *daughv01Lc = 0x0;
2179 AliAODMCParticle *K0S = 0x0;
2180 AliAODMCParticle *daughv02Lc = 0x0;
2182 if (labelsLcdaugh[0] >= 0) {
2183 // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2184 daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2186 AliDebug(3, "Could not access MC info for first daughter of Lc");
2190 AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2191 if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2194 else { // The bachelor is a fake
2198 if (labelsLcdaugh[1] >= 0) {
2199 //Printf("Getting K0S");
2200 K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2202 AliDebug(3, "Could not access MC info for second daughter of Lc");
2206 if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2210 AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2214 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
2215 if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
2216 Int_t iK0 = K0S->GetMother();
2218 Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2220 else { // The K0S has a mother
2221 daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2223 AliDebug(3, "Could not access MC info for second daughter of Lc");
2225 else { // we can access safely the K0S mother in the MC
2226 if( daughv01Lc && (daughv01Lc->GetMother() == daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){ // This is a true cascade! bachelor and V0 come from the same mother
2227 //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
2228 AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
2229 Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
2230 if (motherLc) pdgMum = motherLc->GetPdgCode();
2231 if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
2232 if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
2233 AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
2235 if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
2236 fHistoMassLcTrue->Fill(massLc);
2237 fHistoDecayLengthLcTrue->Fill(decayLengthLc);
2238 fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
2239 fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
2241 fHistoMassV0fromLcTrue->Fill(massV0);
2242 fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2243 fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2245 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...
2246 AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
2248 fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2249 fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2251 fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2252 fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2253 fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2255 fHistoMassLcSgn->Fill(massLc);
2256 fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2258 fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2259 fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2261 fHistoMassV0fromLcSgn->Fill(massV0);
2262 fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2263 fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2266 isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
2269 // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
2270 // First, we compare the vtx of the Lc
2271 Double_t xLcMC = motherLc->Xv();
2272 Double_t yLcMC = motherLc->Yv();
2273 Double_t zLcMC = motherLc->Zv();
2274 //Printf("Got MC vtx for Lc");
2275 Double_t xProtonMC = daughv01Lc->Xv();
2276 Double_t yProtonMC = daughv01Lc->Yv();
2277 Double_t zProtonMC = daughv01Lc->Zv();
2278 //Printf("Got MC vtx for Proton");
2280 Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2281 (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2282 (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2284 // Then, we compare the vtx of the K0s
2285 Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
2286 Double_t yV0MC = motherV0->Yv();
2287 Double_t zV0MC = motherV0->Zv();
2289 //Printf("Got MC vtx for V0");
2290 Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2291 Double_t yPionMC = daughv01->Yv();
2292 Double_t zPionMC = daughv01->Zv();
2293 //Printf("Got MC vtx for Pion");
2294 Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2296 Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2297 (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2298 (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2300 // Then, the K0S vertex wrt primary vertex
2302 Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2303 (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2304 (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2306 fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2307 fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2308 fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2310 } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
2311 else if (!motherLc){
2312 AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
2315 AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2317 } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2319 isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2321 } // closing: else { // we can access safely the K0S mother in the MC
2322 } // closing: else { // The K0S has a mother
2323 } // closing isMCLcok
2324 } // closing !isBkgLc
2325 } // closing fUseMCInfo
2327 //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2328 if ( retMV0 == 0 && retMLc == 0){
2330 errV0KF[0] = sigmaMassV0;
2331 V0KF[1] = decayLengthV0;
2332 errV0KF[1] = sigmaDecayLengthV0;
2333 V0KF[2] = lifeTimeV0;
2334 errV0KF[2] = sigmaLifeTimeV0;
2336 errLcKF[0] = sigmaMassLc;
2337 LcKF[1] = decayLengthLc;
2338 errLcKF[1] = sigmaDecayLengthLc;
2339 LcKF[2] = lifeTimeLc;
2340 errLcKF[2] = sigmaLifeTimeLc;
2346 //________________________________________________________________________
2347 AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
2348 AliAODTrack* bachelor,
2349 TClonesArray *mcArray ){
2351 //Printf("In CheckBachelor");
2353 // function to check if the bachelor is a p, if it is a p but not from Lc
2354 // to be filled for background candidates
2356 Int_t label = bachelor->GetLabel();
2361 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
2362 if (!mcpart) return kBachInvalid;
2363 Int_t pdg = mcpart->PdgCode();
2364 if (TMath::Abs(pdg) != 2212) {
2365 AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
2366 return kBachNoProton;
2368 else { // it is a proton
2369 //Int_t labelLc = part->GetLabel();
2370 Int_t labelLc = FindLcLabel(part, mcArray);
2371 //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
2372 Int_t bachelorMotherLabel = mcpart->GetMother();
2373 //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
2374 if (bachelorMotherLabel == -1) {
2375 return kBachPrimary;
2377 AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2378 if (!bachelorMother) return kBachInvalid;
2379 Int_t pdgMother = bachelorMother->PdgCode();
2380 if (TMath::Abs(pdgMother) != 4122) {
2381 AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2382 return kBachNoLambdaMother;
2384 else { // it comes from Lc
2385 if (labelLc != bachelorMotherLabel){
2386 //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));
2387 AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2388 return kBachDifferentLambdaMother;
2390 else { // it comes from the correct Lc
2391 AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2392 return kBachCorrectLambdaMother;
2397 return kBachInvalid;
2401 //________________________________________________________________________
2402 AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
2404 //AliAODTrack* v0part,
2405 TClonesArray *mcArray ){
2407 // function to check if the K0Spart is a p, if it is a p but not from Lc
2408 // to be filled for background candidates
2410 //Printf(" CheckK0S");
2412 //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
2413 //Int_t label = v0part->GetLabel();
2414 Int_t labelFind = FindV0Label(v0part, mcArray);
2415 //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
2416 if (labelFind == -1) {
2420 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
2421 if (!mcpart) return kK0SInvalid;
2422 Int_t pdg = mcpart->PdgCode();
2423 if (TMath::Abs(pdg) != 310) {
2424 AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2425 //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2428 else { // it is a K0S
2429 //Int_t labelLc = part->GetLabel();
2430 Int_t labelLc = FindLcLabel(part, mcArray);
2431 Int_t K0SpartMotherLabel = mcpart->GetMother();
2432 if (K0SpartMotherLabel == -1) {
2433 return kK0SWithoutMother;
2435 AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
2436 if (!K0SpartMother) return kK0SInvalid;
2437 Int_t pdgMotherK0S = K0SpartMother->PdgCode();
2438 if (TMath::Abs(pdgMotherK0S) != 311) {
2439 AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2440 //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2441 return kK0SNotFromK0;
2443 else { // the K0S comes from a K0
2444 Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
2445 if (K0MotherLabel == -1) {
2448 AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
2449 if (!K0Mother) return kK0SInvalid;
2450 Int_t pdgK0Mother = K0Mother->PdgCode();
2451 if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
2452 AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2453 //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2454 return kK0NoLambdaMother;
2456 else { // the K0 comes from Lc
2457 if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
2458 AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2459 //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2460 return kK0DifferentLambdaMother;
2462 else { // it comes from the correct Lc
2463 AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2464 //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2465 return kK0CorrectLambdaMother;
2475 //----------------------------------------------------------------------------
2476 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2479 //Printf(" FindV0Label");
2481 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2483 Int_t labMother[2]={-1, -1};
2484 AliAODMCParticle *part=0;
2485 AliAODMCParticle *mother=0;
2486 Int_t dgLabels = -1;
2488 Int_t ndg = v0part->GetNDaughters();
2490 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2493 for(Int_t i = 0; i < 2; i++) {
2494 AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
2495 dgLabels = trk->GetLabel();
2496 if (dgLabels == -1) {
2497 //printf("daughter with negative label %d\n", dgLabels);
2500 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2502 //printf("no MC particle\n");
2505 labMother[i] = part->GetMother();
2506 if (labMother[i] != -1){
2507 mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2509 //printf("no MC mother particle\n");
2518 if (labMother[0] == labMother[1]) return labMother[0];
2524 //----------------------------------------------------------------------------
2525 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
2528 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2530 //Printf(" FindLcLabel");
2532 AliAODMCParticle *part=0;
2533 AliAODMCParticle *mother=0;
2534 AliAODMCParticle *grandmother=0;
2535 Int_t dgLabels = -1;
2537 Int_t ndg = cascade->GetNDaughters();
2539 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2542 // daughter 0 --> bachelor, daughter 1 --> V0
2543 AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
2544 dgLabels = trk->GetLabel();
2545 if (dgLabels == -1 ) {
2546 //printf("daughter with negative label %d\n", dgLabels);
2549 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2551 //printf("no MC particle\n");
2554 Int_t labMotherBach = part->GetMother();
2555 if (labMotherBach == -1){
2558 mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2560 //printf("no MC mother particle\n");
2564 AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
2565 dgLabels = FindV0Label(v0, mcArray);
2566 if (dgLabels == -1 ) {
2567 //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
2570 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2572 //printf("no MC particle\n");
2575 Int_t labMotherv0 = part->GetMother();
2576 if (labMotherv0 == -1){
2579 mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2581 //printf("no MC mother particle\n");
2584 Int_t labGrandMotherv0 = mother->GetMother();
2585 if (labGrandMotherv0 == -1){
2588 grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2590 //printf("no MC mother particle\n");
2594 //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2595 if (labGrandMotherv0 == labMotherBach) return labMotherBach;