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),
197 //___________________________________________________________________________
198 AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
199 AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
200 AliAnalysisTaskSE(name),
206 fIsK0sAnalysis(kFALSE),
210 fUseOnTheFlyV0(useOnTheFly),
211 fIsEventSelected(kFALSE),
212 fVariablesTreeSgn(0),
213 fVariablesTreeBkg(0),
214 fCandidateVariables(),
219 fFillOnlySgn(kFALSE),
220 fHistoLcBeforeCuts(0),
221 fHistoFiducialAcceptance(0),
224 fHistoLcpKpiBeforeCuts(0),
227 fHistoDistanceLcToPrimVtx(0),
228 fHistoDistanceV0ToPrimVtx(0),
229 fHistoDistanceV0ToLc(0),
231 fHistoDistanceLcToPrimVtxSgn(0),
232 fHistoDistanceV0ToPrimVtxSgn(0),
233 fHistoDistanceV0ToLcSgn(0),
235 fHistoVtxLcResidualToPrimVtx(0),
236 fHistoVtxV0ResidualToPrimVtx(0),
237 fHistoVtxV0ResidualToLc(0),
240 fHistoDecayLengthV0All(0),
241 fHistoLifeTimeV0All(0),
244 fHistoDecayLengthV0True(0),
245 fHistoLifeTimeV0True(0),
247 fHistoMassV0TrueFromAOD(0),
249 fHistoMassV0TrueK0S(0),
250 fHistoDecayLengthV0TrueK0S(0),
251 fHistoLifeTimeV0TrueK0S(0),
253 fHistoMassV0TrueK0SFromAOD(0),
256 fHistoDecayLengthLcAll(0),
257 fHistoLifeTimeLcAll(0),
260 fHistoDecayLengthLcTrue(0),
261 fHistoLifeTimeLcTrue(0),
263 fHistoMassLcTrueFromAOD(0),
265 fHistoMassV0fromLcAll(0),
266 fHistoDecayLengthV0fromLcAll(0),
267 fHistoLifeTimeV0fromLcAll(0),
269 fHistoMassV0fromLcTrue(0),
270 fHistoDecayLengthV0fromLcTrue(0),
271 fHistoLifeTimeV0fromLcTrue(0),
274 fHistoMassLcSgnFromAOD(0),
275 fHistoDecayLengthLcSgn(0),
276 fHistoLifeTimeLcSgn(0),
278 fHistoMassV0fromLcSgn(0),
279 fHistoDecayLengthV0fromLcSgn(0),
280 fHistoLifeTimeV0fromLcSgn(0),
287 fHistoDecayLengthKFV0(0),
288 fHistoLifeTimeKFV0(0),
291 fHistoDecayLengthKFLc(0),
292 fHistoLifeTimeKFLc(0),
294 fHistoArmenterosPodolanskiV0KF(0),
295 fHistoArmenterosPodolanskiV0KFSgn(0),
296 fHistoArmenterosPodolanskiV0AOD(0),
297 fHistoArmenterosPodolanskiV0AODSgn(0),
301 ftopoConstraint(kTRUE),
302 fCallKFVertexing(kFALSE),
303 fKeepingOnlyHIJINGBkg(kFALSE),
306 fCutKFChi2NDF(999999.),
307 fCutKFDeviationFromVtx(999999.),
308 fCutKFDeviationFromVtxV0(0.),
311 fKeepingOnlyPYTHIABkg(kFALSE),
316 // Constructor. Initialization of Inputs and Outputs
318 Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
320 DefineOutput(1, TList::Class()); // Tree signal + Tree Bkg + histoEvents
321 DefineOutput(2, AliNormalizationCounter::Class()); // normalization counter object
322 DefineOutput(3, TList::Class()); // Cuts
323 DefineOutput(4, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
324 DefineOutput(5, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
325 DefineOutput(6, TList::Class()); // Tree signal + Tree Bkg + histoEvents
329 //___________________________________________________________________________
330 AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
334 Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
364 if(fVariablesTreeSgn){
365 delete fVariablesTreeSgn;
366 fVariablesTreeSgn = 0;
369 if(fVariablesTreeBkg){
370 delete fVariablesTreeBkg;
371 fVariablesTreeBkg = 0;
385 //_________________________________________________
386 void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
391 fIsEventSelected=kFALSE;
393 if (fDebug > 1) AliInfo("Init");
395 fListCuts = new TList();
396 fListCuts->SetOwner();
397 fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
398 PostData(3,fListCuts);
400 if (fUseMCInfo && (fKeepingOnlyHIJINGBkg || fKeepingOnlyPYTHIABkg)) fUtils = new AliVertexingHFUtils();
405 //________________________________________ terminate ___________________________
406 void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
408 // The Terminate() function is the last function to be called during
409 // a query. It always runs on the client, it can be used to present
410 // the results graphically or save the results to file.
412 //AliInfo("Terminate","");
413 AliAnalysisTaskSE::Terminate();
415 fOutput = dynamic_cast<TList*> (GetOutputData(1));
417 AliError("fOutput not available");
422 AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0Sp->GetEntries()));
424 fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
426 AliError("fOutputKF not available");
433 //___________________________________________________________________________
434 void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
436 AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
440 fOutput = new TList();
442 fOutput->SetName("listTrees");
444 // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
445 const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
446 fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
447 fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
449 fCandidateVariables = new Float_t [nVar];
450 TString * fCandidateVariableNames = new TString[nVar];
451 fCandidateVariableNames[0]="massLc2K0Sp";
452 fCandidateVariableNames[1]="massLc2Lambdapi";
453 fCandidateVariableNames[2]="massK0S";
454 fCandidateVariableNames[3]="massLambda";
455 fCandidateVariableNames[4]="massLambdaBar";
456 fCandidateVariableNames[5]="cosPAK0S";
457 fCandidateVariableNames[6]="dcaV0";
458 fCandidateVariableNames[7]="tImpParBach";
459 fCandidateVariableNames[8]="tImpParV0";
460 fCandidateVariableNames[9]="nSigmaTPCpr";
461 fCandidateVariableNames[10]="nSigmaTPCpi";
462 fCandidateVariableNames[11]="nSigmaTPCka";
463 fCandidateVariableNames[12]="nSigmaTOFpr";
464 fCandidateVariableNames[13]="nSigmaTOFpi";
465 fCandidateVariableNames[14]="nSigmaTOFka";
466 fCandidateVariableNames[15]="bachelorPt";
467 fCandidateVariableNames[16]="V0positivePt";
468 fCandidateVariableNames[17]="V0negativePt";
469 fCandidateVariableNames[18]="dcaV0pos";
470 fCandidateVariableNames[19]="dcaV0neg";
471 fCandidateVariableNames[20]="v0Pt";
472 fCandidateVariableNames[21]="massGamma";
473 fCandidateVariableNames[22]="LcPt";
474 fCandidateVariableNames[23]="combinedProtonProb";
475 fCandidateVariableNames[24]="LcEta";
476 fCandidateVariableNames[25]="V0positiveEta";
477 fCandidateVariableNames[26]="V0negativeEta";
478 fCandidateVariableNames[27]="TPCProtonProb";
479 fCandidateVariableNames[28]="TOFProtonProb";
480 fCandidateVariableNames[29]="bachelorEta";
481 fCandidateVariableNames[30]="LcP";
482 fCandidateVariableNames[31]="bachelorP";
483 fCandidateVariableNames[32]="v0P";
484 fCandidateVariableNames[33]="V0positiveP";
485 fCandidateVariableNames[34]="V0negativeP";
486 fCandidateVariableNames[35]="LcY";
487 fCandidateVariableNames[36]="v0Y";
488 fCandidateVariableNames[37]="bachelorY";
489 fCandidateVariableNames[38]="V0positiveY";
490 fCandidateVariableNames[39]="V0negativeY";
491 fCandidateVariableNames[40]="v0Eta";
492 fCandidateVariableNames[41]="DecayLengthLc";
493 fCandidateVariableNames[42]="DecayLengthK0S";
494 fCandidateVariableNames[43]="CtLc";
495 fCandidateVariableNames[44]="CtK0S";
496 fCandidateVariableNames[45]="bachCode";
497 fCandidateVariableNames[46]="k0SCode";
499 fCandidateVariableNames[47]="V0KFmass";
500 fCandidateVariableNames[48]="V0KFdecayLength";
501 fCandidateVariableNames[49]="V0KFlifeTime";
503 fCandidateVariableNames[50]="V0KFmassErr";
504 fCandidateVariableNames[51]="V0KFdecayTimeErr";
505 fCandidateVariableNames[52]="V0KFlifeTimeErr";
507 fCandidateVariableNames[53]="LcKFmass";
508 fCandidateVariableNames[54]="LcKFdecayLength";
509 fCandidateVariableNames[55]="LcKFlifeTime";
511 fCandidateVariableNames[56]="LcKFmassErr";
512 fCandidateVariableNames[57]="LcKFdecayTimeErr";
513 fCandidateVariableNames[58]="LcKFlifeTimeErr";
515 fCandidateVariableNames[59]="LcKFDistToPrimVtx";
516 fCandidateVariableNames[60]="V0KFDistToPrimVtx";
517 fCandidateVariableNames[61]="V0KFDistToLc";
518 fCandidateVariableNames[62]="alphaArmKF";
519 fCandidateVariableNames[63]="ptArmKF";
520 fCandidateVariableNames[64]="alphaArm";
521 fCandidateVariableNames[65]="ptArm";
523 fCandidateVariableNames[66]="ITSrefitV0pos";
524 fCandidateVariableNames[67]="ITSrefitV0neg";
526 fCandidateVariableNames[68]="TPCClV0pos";
527 fCandidateVariableNames[69]="TPCClV0neg";
529 fCandidateVariableNames[70]="v0Xcoord";
530 fCandidateVariableNames[71]="v0Ycoord";
531 fCandidateVariableNames[72]="v0Zcoord";
532 fCandidateVariableNames[73]="primVtxX";
533 fCandidateVariableNames[74]="primVtxY";
534 fCandidateVariableNames[75]="primVtxZ";
536 fCandidateVariableNames[76]="ITSclBach";
537 fCandidateVariableNames[77]="SPDclBach";
539 fCandidateVariableNames[78]="ITSclV0pos";
540 fCandidateVariableNames[79]="SPDclV0pos";
541 fCandidateVariableNames[80]="ITSclV0neg";
542 fCandidateVariableNames[81]="SPDclV0neg";
544 fCandidateVariableNames[82]="alphaArmLc";
545 fCandidateVariableNames[83]="alphaArmLcCharge";
546 fCandidateVariableNames[84]="ptArmLc";
548 fCandidateVariableNames[85]="CosThetaStar";
551 for(Int_t ivar=0; ivar<nVar; ivar++){
552 fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
553 fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
556 fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
557 TString labelEv[2] = {"NotSelected", "Selected"};
558 for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
559 fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
562 fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
564 fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
565 TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
566 for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
567 fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
570 fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
571 TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
572 for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
573 fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
574 fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
577 fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
578 TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
579 for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
580 fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
583 fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
584 fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
586 TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
587 "kBachDifferentLambdaMother", "kBachCorrectLambdaMother"};
588 TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
589 "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
591 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
592 fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
593 fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
595 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
596 fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
597 fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
600 fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
601 for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
602 fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
605 fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 4, -0.5, 3.5);
606 TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"};
607 for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
608 fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
611 //fOutput->Add(fVariablesTreeSgn);
612 //fOutput->Add(fVariablesTreeBkg);
614 const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.};
616 fHistoMCLcK0Sp = new TH1F("fHistoMCLcK0Sp", "fHistoMCLcK0Sp", 14, ptbins);
618 fOutput->Add(fHistoEvents);
619 fOutput->Add(fHistoLc);
620 fOutput->Add(fHistoLcOnTheFly);
621 fOutput->Add(fHistoLcBeforeCuts);
622 fOutput->Add(fHistoFiducialAcceptance);
623 fOutput->Add(fHistoCodesSgn);
624 fOutput->Add(fHistoCodesBkg);
625 fOutput->Add(fHistoLcpKpiBeforeCuts);
626 fOutput->Add(fHistoBackground);
627 fOutput->Add(fHistoMCLcK0Sp);
629 PostData(1, fOutput);
630 PostData(4, fVariablesTreeSgn);
631 PostData(5, fVariablesTreeBkg);
632 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
633 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
634 fPIDResponse = inputHandler->GetPIDResponse();
636 if (fAnalCuts->GetIsUsePID()){
638 fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
639 fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
640 fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
641 fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
642 fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
643 fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
645 fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
648 // Setting properties of PID
649 fPIDCombined=new AliPIDCombined;
650 fPIDCombined->SetDefaultTPCPriors();
651 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
652 //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
654 fCounter = new AliNormalizationCounter("NormalizationCounter");
656 PostData(2, fCounter);
658 // Histograms from KF
660 if (fCallKFVertexing){
661 fHistoDistanceLcToPrimVtx = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
662 fHistoDistanceV0ToPrimVtx = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
663 fHistoDistanceV0ToLc = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
665 fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
666 fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
667 fHistoDistanceV0ToLcSgn = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
669 fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
670 fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
671 fHistoVtxV0ResidualToLc = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
673 fHistoMassV0All = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
674 fHistoDecayLengthV0All = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
675 fHistoLifeTimeV0All = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
677 fHistoMassV0True = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
678 fHistoDecayLengthV0True = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
679 fHistoLifeTimeV0True = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
681 fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
683 fHistoMassV0TrueK0S = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
684 fHistoDecayLengthV0TrueK0S = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
685 fHistoLifeTimeV0TrueK0S = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
687 fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
689 fHistoMassLcAll = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
690 fHistoDecayLengthLcAll = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
691 fHistoLifeTimeLcAll = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
693 fHistoMassLcTrue = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
694 fHistoDecayLengthLcTrue = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
695 fHistoLifeTimeLcTrue = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
697 fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
699 fHistoMassV0fromLcAll = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
700 fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
701 fHistoLifeTimeV0fromLcAll = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
703 fHistoMassV0fromLcTrue = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
704 fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
705 fHistoLifeTimeV0fromLcTrue = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
707 fHistoMassLcSgn = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
708 fHistoMassLcSgnFromAOD = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
709 fHistoDecayLengthLcSgn = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
710 fHistoLifeTimeLcSgn = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
712 fHistoMassV0fromLcSgn = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
713 fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
714 fHistoLifeTimeV0fromLcSgn = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
716 fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
717 fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
718 fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
719 TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
720 "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
721 "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
722 "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
724 for (Int_t ibin = 1; ibin <=16; ibin++){
725 fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
726 fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
727 fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
728 fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
731 fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
732 fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
733 fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
735 fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
736 fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
737 fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
739 fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
740 fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
741 fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
742 fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
745 fOutputKF = new TList();
746 fOutputKF->SetOwner();
747 fOutputKF->SetName("listHistoKF");
749 if (fCallKFVertexing){
750 fOutputKF->Add(fHistoDistanceLcToPrimVtx);
751 fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
752 fOutputKF->Add(fHistoDistanceV0ToLc);
754 fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
755 fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
756 fOutputKF->Add(fHistoDistanceV0ToLcSgn);
758 fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
759 fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
760 fOutputKF->Add(fHistoVtxV0ResidualToLc);
762 fOutputKF->Add(fHistoMassV0All);
763 fOutputKF->Add(fHistoDecayLengthV0All);
764 fOutputKF->Add(fHistoLifeTimeV0All);
766 fOutputKF->Add(fHistoMassV0True);
767 fOutputKF->Add(fHistoDecayLengthV0True);
768 fOutputKF->Add(fHistoLifeTimeV0True);
770 fOutputKF->Add(fHistoMassV0TrueFromAOD);
772 fOutputKF->Add(fHistoMassV0TrueK0S);
773 fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
774 fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
776 fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
778 fOutputKF->Add(fHistoMassLcAll);
779 fOutputKF->Add(fHistoDecayLengthLcAll);
780 fOutputKF->Add(fHistoLifeTimeLcAll);
782 fOutputKF->Add(fHistoMassLcTrue);
783 fOutputKF->Add(fHistoDecayLengthLcTrue);
784 fOutputKF->Add(fHistoLifeTimeLcTrue);
786 fOutputKF->Add(fHistoMassLcTrueFromAOD);
788 fOutputKF->Add(fHistoMassV0fromLcAll);
789 fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
790 fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
792 fOutputKF->Add(fHistoMassV0fromLcTrue);
793 fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
794 fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
796 fOutputKF->Add(fHistoMassLcSgn);
797 fOutputKF->Add(fHistoMassLcSgnFromAOD);
798 fOutputKF->Add(fHistoDecayLengthLcSgn);
799 fOutputKF->Add(fHistoLifeTimeLcSgn);
801 fOutputKF->Add(fHistoMassV0fromLcSgn);
802 fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
803 fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
805 fOutputKF->Add(fHistoKF);
806 fOutputKF->Add(fHistoKFV0);
807 fOutputKF->Add(fHistoKFLc);
809 fOutputKF->Add(fHistoMassKFV0);
810 fOutputKF->Add(fHistoDecayLengthKFV0);
811 fOutputKF->Add(fHistoLifeTimeKFV0);
813 fOutputKF->Add(fHistoMassKFLc);
814 fOutputKF->Add(fHistoDecayLengthKFLc);
815 fOutputKF->Add(fHistoLifeTimeKFLc);
817 fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
818 fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
819 fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
820 fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
823 PostData(6, fOutputKF);
828 //_________________________________________________
829 void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
833 AliError("NO EVENT FOUND!");
838 AliDebug(2, Form("Processing event = %d", fCurrentEvent));
839 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
840 TClonesArray *arrayLctopKos=0;
842 TClonesArray *array3Prong = 0;
844 if (!aodEvent && AODEvent() && IsStandardAOD()) {
845 // In case there is an AOD handler writing a standard AOD, use the AOD
846 // event in memory rather than the input (ESD) event.
847 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
848 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
849 // have to taken from the AOD event hold by the AliAODExtension
850 AliAODHandler* aodHandler = (AliAODHandler*)
851 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
853 if (aodHandler->GetExtensions()) {
854 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
855 AliAODEvent *aodFromExt = ext->GetAOD();
856 arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
858 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
861 arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
863 array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
866 if ( !fUseMCInfo && fIspA) {
867 fAnalCuts->SetTriggerClass("");
868 fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
871 fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
874 TClonesArray *mcArray = 0;
875 AliAODMCHeader *mcHeader=0;
878 // MC array need for maching
879 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
881 AliError("Could not find Monte-Carlo in AOD");
885 mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
887 AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
890 //Printf("Filling MC histo");
891 FillMCHisto(mcArray);
894 // AOD primary vertex
895 fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
897 if (fVtx1->GetNContributors()<1) return;
899 fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
901 if ( !fIsEventSelected ) {
902 fHistoEvents->Fill(0);
903 return; // don't take into account not selected events
905 fHistoEvents->Fill(1);
907 // Setting magnetic field for KF vertexing
908 fBField = aodEvent->GetMagneticField();
909 AliKFParticle::SetField(fBField);
911 Int_t nSelectedAnal = 0;
912 if (fIsK0sAnalysis) {
913 MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
914 nSelectedAnal, fAnalCuts,
915 array3Prong, mcHeader);
917 fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
919 PostData(1, fOutput);
920 PostData(2, fCounter);
921 PostData(4, fVariablesTreeSgn);
922 PostData(5, fVariablesTreeBkg);
923 PostData(6, fOutputKF);
926 //-------------------------------------------------------------------------------
927 void AliAnalysisTaskSELc2V0bachelorTMVA::FillMCHisto(TClonesArray *mcArray){
929 // method to fill MC histo: how many Lc --> K0S + p are there at MC level
930 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
931 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
933 AliError("Failed casting particle from MC array!, Skipping particle");
934 AliInfo("Failed casting particle from MC array!, Skipping particle");
937 Int_t pdg = mcPart->GetPdgCode();
938 if (TMath::Abs(pdg) != 4122){
939 AliDebug(2, Form("MC particle %d is not a Lc: its pdg code is %d", iPart, pdg));
942 AliInfo(Form("Step 0 ok: MC particle %d is a Lc: its pdg code is %d", iPart, pdg));
943 Int_t labeldaugh0 = mcPart->GetDaughter(0);
944 Int_t labeldaugh1 = mcPart->GetDaughter(1);
945 if (labeldaugh0 <= 0 || labeldaugh1 <= 0){
946 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
949 else if (labeldaugh1 - labeldaugh0 == 1){
950 AliInfo(Form("Step 1 ok: The MC particle has correct daughters!!"));
951 AliAODMCParticle* daugh0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh0));
952 AliAODMCParticle* daugh1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labeldaugh1));
953 Int_t pdgCodeDaugh0 = TMath::Abs(daugh0->GetPdgCode());
954 Int_t pdgCodeDaugh1 = TMath::Abs(daugh1->GetPdgCode());
955 AliAODMCParticle* bachelorMC = daugh0;
956 AliAODMCParticle* v0MC = daugh1;
957 AliDebug(2, Form("pdgCodeDaugh0 = %d, pdgCodeDaugh1 = %d", pdgCodeDaugh0, pdgCodeDaugh1));
958 if ((pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) || (pdgCodeDaugh0 == 2212 && pdgCodeDaugh1 == 311)){
959 // 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-
960 /// first, we set the bachelor and the v0: above we assumed first proton and second V0, but we could have to change it:
961 if (pdgCodeDaugh0 == 311 && pdgCodeDaugh1 == 2212) {
965 AliDebug(2, Form("Number of Daughters of v0 = %d", v0MC->GetNDaughters()));
966 if (v0MC->GetNDaughters() != 1) {
967 AliDebug(2, "The K0 does not decay in 1 body only! Impossible... Continuing...");
970 else { // So far: Lc --> K0 + p, K0 with 1 daughter
971 AliInfo("Step 2 ok: The K0 does decay in 1 body only! ");
972 Int_t labelK0daugh = v0MC->GetDaughter(0);
973 AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0daugh));
975 AliError("Error while casting particle! returning a NULL array");
976 AliInfo("Error while casting particle! returning a NULL array");
979 else { // So far: Lc --> K0 + p, K0 with 1 daughter that we can access
980 if (partK0S->GetNDaughters() != 2 || TMath::Abs(partK0S->GetPdgCode() != 310)){
981 AliDebug(2, "The K0 daughter is not a K0S or does not decay in 2 bodies");
984 else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies
985 AliInfo("Step 3 ok: The K0 daughter is a K0S and does decay in 2 bodies");
986 Int_t labelK0Sdaugh0 = partK0S->GetDaughter(0);
987 Int_t labelK0Sdaugh1 = partK0S->GetDaughter(1);
988 AliAODMCParticle* daughK0S0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh0));
989 AliAODMCParticle* daughK0S1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelK0Sdaugh1));
990 if (!daughK0S0 || ! daughK0S1){
991 AliDebug(2, "Could not access K0S daughters, continuing...");
994 else { // So far: Lc --> K0 + p, K0 --> K0S, K0S in 2 bodies that we can access
995 AliInfo("Step 4 ok: Could access K0S daughters, continuing...");
996 Int_t pdgK0Sdaugh0 = daughK0S0->GetPdgCode();
997 Int_t pdgK0Sdaugh1 = daughK0S1->GetPdgCode();
998 if (TMath::Abs(pdgK0Sdaugh0) != 211 || TMath::Abs(pdgK0Sdaugh1) != 211){
999 AliDebug(2, "The K0S does not decay in pi+pi-, continuing");
1000 AliInfo("The K0S does not decay in pi+pi-, continuing");
1002 else { // Full chain: Lc --> K0 + p, K0 --> K0S, K0S --> pi+pi-
1003 if (fAnalCuts->IsInFiducialAcceptance(mcPart->Pt(), mcPart->Y())) {
1004 AliDebug(2, Form("----> Filling histo with pt = %f", mcPart->Pt()));
1005 fHistoMCLcK0Sp->Fill(mcPart->Pt());
1008 AliDebug(2, "not in fiducial acceptance! Skipping");
1018 } // closing loop over mcArray
1024 //-------------------------------------------------------------------------------
1025 void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
1026 TClonesArray *mcArray,
1027 Int_t &nSelectedAnal,
1028 AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
1029 AliAODMCHeader* aodheader){
1030 //Lc prong needed to MatchToMC method
1032 Int_t pdgCand = 4122;
1033 Int_t pdgDgLctoV0bachelor[2]={2212, 310};
1034 Int_t pdgDgV0toDaughters[2]={211, 211};
1036 Int_t pdgDgLctopKpi[3]={2212, 321, 211};
1038 // loop to search for candidates Lc->p+K+pi
1039 Int_t n3Prong = array3Prong->GetEntriesFast();
1040 Int_t nCascades= arrayLctopKos->GetEntriesFast();
1042 //AliInfo(Form("\n\n\n\n3 prong candidates = %d, ncascades = %d \n\n\n\n\n", n3Prong, nCascades));
1043 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
1044 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
1045 //Filling a control histogram with no cuts
1048 // find associated MC particle for Lc -> p+K+pi
1049 Int_t mcLabel = d->MatchToMC(4122, mcArray, 3, pdgDgLctopKpi);
1050 //Int_t mcLabelTemp = d->MatchToMC(4122, mcArray);
1051 //Printf("mcLabel = %d, mcLabelTemp = %d", mcLabel, mcLabelTemp);
1054 AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1056 Int_t pdgCode = partLcpKpi->GetPdgCode();
1057 AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1058 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1059 fHistoLcpKpiBeforeCuts->Fill(1);
1064 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
1065 fHistoLcpKpiBeforeCuts->Fill(0);
1070 // loop over cascades to search for candidates Lc->p+K0S
1072 for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
1074 // Lc candidates and K0s from Lc
1075 AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
1077 AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
1081 //Filling a control histogram with no cuts
1086 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1087 fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1088 if (fmcLabelLc>=0) {
1089 AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades));
1091 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
1093 pdgCode = partLc->GetPdgCode();
1094 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
1095 pdgCode = TMath::Abs(pdgCode);
1096 fHistoLcBeforeCuts->Fill(1);
1101 fHistoLcBeforeCuts->Fill(0);
1106 //if (!lcK0spr->GetSecondaryVtx()) {
1107 // AliInfo("No secondary vertex");
1111 if (lcK0spr->GetNDaughters()!=2) {
1112 AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
1116 AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
1117 AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
1118 if (!v0part || !bachPart) {
1119 AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0s));
1124 if (!v0part->GetSecondaryVtx()) {
1125 AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
1129 if (v0part->GetNDaughters()!=2) {
1130 AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
1134 AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
1135 AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
1136 if (!v0Neg || !v0Pos) {
1137 AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
1142 if (v0Pos->Charge() == v0Neg->Charge()) {
1143 AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
1153 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1154 Int_t mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1156 AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1158 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1160 pdgCode = partLc->GetPdgCode();
1161 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1162 pdgCode = TMath::Abs(pdgCode);
1173 AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
1174 AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1176 if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1179 else { // checking if we want to fill the background
1180 if (fKeepingOnlyHIJINGBkg){
1181 // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1182 Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
1183 if (!isCandidateInjected){
1184 AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1185 fHistoBackground->Fill(1);
1188 AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1189 fHistoBackground->Fill(0);
1193 else if (fKeepingOnlyPYTHIABkg){
1194 // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1195 AliAODTrack *bachelor = (AliAODTrack*)lcK0spr->GetBachelor();
1196 AliAODTrack *v0pos = (AliAODTrack*)lcK0spr->Getv0PositiveTrack();
1197 AliAODTrack *v0neg = (AliAODTrack*)lcK0spr->Getv0NegativeTrack();
1198 if (!bachelor || !v0pos || !v0neg) {
1199 AliDebug(2, "Cannot retrieve one of the tracks while checking origin, continuing");
1203 Int_t labelbachelor = TMath::Abs(bachelor->GetLabel());
1204 Int_t labelv0pos = TMath::Abs(v0pos->GetLabel());
1205 Int_t labelv0neg = TMath::Abs(v0neg->GetLabel());
1206 AliAODMCParticle* MCbachelor = (AliAODMCParticle*)mcArray->At(labelbachelor);
1207 AliAODMCParticle* MCv0pos = (AliAODMCParticle*)mcArray->At(labelv0pos);
1208 AliAODMCParticle* MCv0neg = (AliAODMCParticle*)mcArray->At(labelv0neg);
1209 if (!MCbachelor || !MCv0pos || !MCv0neg) {
1210 AliDebug(2, "Cannot retrieve MC particle for one of the tracks while checking origin, continuing");
1214 Int_t isBachelorFromPythia = fUtils->CheckOrigin(mcArray, MCbachelor, kTRUE);
1215 Int_t isv0posFromPythia = fUtils->CheckOrigin(mcArray, MCv0pos, kTRUE);
1216 Int_t isv0negFromPythia = fUtils->CheckOrigin(mcArray, MCv0neg, kTRUE);
1217 if (isBachelorFromPythia != 0 && isv0posFromPythia != 0 && isv0negFromPythia != 0){
1218 AliDebug(2, "The candidate is from PYTHIA (i.e. all daughters originate from a quark), keeping it to fill background");
1219 fHistoBackground->Fill(2);
1222 AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background");
1223 fHistoBackground->Fill(3);
1232 FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray);
1239 //________________________________________________________________________
1240 void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
1242 Int_t &nSelectedAnal,
1243 AliRDHFCutsLctoV0 *cutsAnal,
1244 TClonesArray *mcArray){
1246 // Fill histos for Lc -> K0S+proton
1250 if (!part->GetOwnPrimaryVtx()) {
1251 //Printf("No primary vertex for Lc found!!");
1252 part->SetOwnPrimaryVtx(fVtx1);
1255 //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1258 Double_t invmassLc = part->InvMassLctoK0sP();
1260 AliAODv0 * v0part = part->Getv0();
1261 Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1262 if (onFlyV0){ // on-the-fly V0
1264 fHistoLcOnTheFly->Fill(2.);
1267 fHistoLcOnTheFly->Fill(0.);
1270 else { // offline V0
1272 fHistoLcOnTheFly->Fill(3.);
1275 fHistoLcOnTheFly->Fill(1.);
1279 Double_t dcaV0 = v0part->GetDCA();
1280 Double_t invmassK0s = v0part->MassK0Short();
1282 if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
1284 fHistoFiducialAcceptance->Fill(3.);
1287 fHistoFiducialAcceptance->Fill(1.);
1292 fHistoFiducialAcceptance->Fill(2.);
1295 fHistoFiducialAcceptance->Fill(0.);
1299 Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1301 if (isInV0window == 0) {
1302 AliDebug(2, "No: The candidate has NOT passed the mass cuts!");
1303 if (isLc) Printf("SIGNAL candidate rejected");
1306 else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1308 Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1310 if (!isInCascadeWindow) {
1311 AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1312 if (isLc) Printf("SIGNAL candidate rejected");
1315 else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1317 Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1318 if (!isCandidateSelectedCuts){
1319 AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1320 if (isLc) Printf("SIGNAL candidate rejected");
1324 AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1327 AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1329 AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1333 //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1334 Double_t probTPCTOF[AliPID::kSPECIES]={-1.};
1336 UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1337 AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed));
1338 Double_t probProton = -1.;
1339 // Double_t probPion = -1.;
1340 // Double_t probKaon = -1.;
1341 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1342 AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1343 probProton = probTPCTOF[AliPID::kProton];
1344 // probPion = probTPCTOF[AliPID::kPion];
1345 // probKaon = probTPCTOF[AliPID::kKaon];
1347 else { // if you don't have both TOF and TPC, try only TPC
1348 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1349 AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC");
1350 detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1351 AliDebug(2,Form(" detUsed (TPC case) = %d", detUsed));
1352 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1353 probProton = probTPCTOF[AliPID::kProton];
1354 // probPion = probTPCTOF[AliPID::kPion];
1355 // probKaon = probTPCTOF[AliPID::kKaon];
1356 AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
1359 AliDebug(2, "Only TPC did not work...");
1361 // resetting mask to ask for both TPC+TOF
1362 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1364 AliDebug(2, Form("probProton = %f", probProton));
1366 // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) )
1367 Double_t probProtonTPC = -1.;
1368 Double_t probProtonTOF = -1.;
1369 Double_t pidTPC[AliPID::kSPECIES]={-1.};
1370 Double_t pidTOF[AliPID::kSPECIES]={-1.};
1371 Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC);
1372 Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF);
1373 if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton];
1374 if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton];
1376 // checking V0 status (on-the-fly vs offline)
1377 if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1378 AliDebug(2, "On-the-fly discarded");
1382 if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
1386 if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
1388 if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1389 if (isLc) Printf("SIGNAL candidate rejected");
1390 AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1394 AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1397 Int_t pdgCand = 4122;
1398 Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
1399 Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
1400 Int_t isLc2LBarpi=0, isLc2Lpi=0;
1401 Int_t currentLabel = part->GetLabel();
1404 mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1406 if (bachelor->Charge()==-1) isLc2LBarpi=1;
1407 if (bachelor->Charge()==+1) isLc2Lpi=1;
1411 Int_t pdgDg2prong[2] = {211, 211};
1415 labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1416 if (labelK0S>=0) isK0S = 1;
1419 pdgDg2prong[0] = 211;
1420 pdgDg2prong[1] = 2212;
1422 Int_t isLambdaBar = 0;
1423 Int_t lambdaLabel = 0;
1425 lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1426 if (lambdaLabel>=0) {
1427 AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1428 if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1429 else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1433 pdgDg2prong[0] = 11;
1434 pdgDg2prong[1] = 11;
1436 Int_t gammaLabel = 0;
1438 gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1439 if (gammaLabel>=0) {
1440 AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1441 if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1446 if (currentLabel != -1){
1447 AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1448 pdgTemp = tempPart->GetPdgCode();
1450 if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1451 else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1452 else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1453 else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1454 if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1455 else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
1456 else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
1457 else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
1458 //else AliDebug(2, Form("V0 is something else!!"));
1460 Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1461 Double_t invmassLambda = v0part->MassLambda();
1462 Double_t invmassLambdaBar = v0part->MassAntiLambda();
1464 //Double_t nSigmaITSpr=-999.;
1465 Double_t nSigmaTPCpr=-999.;
1466 Double_t nSigmaTOFpr=-999.;
1468 //Double_t nSigmaITSpi=-999.;
1469 Double_t nSigmaTPCpi=-999.;
1470 Double_t nSigmaTOFpi=-999.;
1472 //Double_t nSigmaITSka=-999.;
1473 Double_t nSigmaTPCka=-999.;
1474 Double_t nSigmaTOFka=-999.;
1477 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1478 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1479 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1480 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1481 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1482 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1483 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1484 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1485 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1488 nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
1489 nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
1490 nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
1491 nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
1492 nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
1493 nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
1496 // Fill candidate variable Tree (track selection, V0 invMass selection)
1497 if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
1499 fCandidateVariables[0] = invmassLc;
1500 fCandidateVariables[1] = invmassLc2Lpi;
1501 fCandidateVariables[2] = invmassK0s;
1502 fCandidateVariables[3] = invmassLambda;
1503 fCandidateVariables[4] = invmassLambdaBar;
1504 fCandidateVariables[5] = part->CosV0PointingAngle();
1505 fCandidateVariables[6] = dcaV0;
1506 fCandidateVariables[7] = part->Getd0Prong(0);
1507 fCandidateVariables[8] = part->Getd0Prong(1);
1508 fCandidateVariables[9] = nSigmaTPCpr;
1509 fCandidateVariables[10] = nSigmaTPCpi;
1510 fCandidateVariables[11] = nSigmaTPCka;
1511 fCandidateVariables[12] = nSigmaTOFpr;
1512 fCandidateVariables[13] = nSigmaTOFpi;
1513 fCandidateVariables[14] = nSigmaTOFka;
1514 fCandidateVariables[15] = bachelor->Pt();
1515 AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1516 fCandidateVariables[16] = v0pos->Pt();
1517 AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
1518 fCandidateVariables[17] = v0neg->Pt();
1519 fCandidateVariables[18] = v0part->Getd0Prong(0);
1520 fCandidateVariables[19] = v0part->Getd0Prong(1);
1521 fCandidateVariables[20] = v0part->Pt();
1522 fCandidateVariables[21] = v0part->InvMass2Prongs(0,1,11,11);
1523 fCandidateVariables[22] = part->Pt();
1524 fCandidateVariables[23] = probProton;
1525 fCandidateVariables[24] = part->Eta();
1526 fCandidateVariables[25] = v0pos->Eta();
1527 fCandidateVariables[26] = v0neg->Eta();
1528 fCandidateVariables[27] = probProtonTPC;
1529 fCandidateVariables[28] = probProtonTOF;
1530 fCandidateVariables[29] = bachelor->Eta();
1532 fCandidateVariables[30] = part->P();
1533 fCandidateVariables[31] = bachelor->P();
1534 fCandidateVariables[32] = v0part->P();
1535 fCandidateVariables[33] = v0pos->P();
1536 fCandidateVariables[34] = v0neg->P();
1538 fCandidateVariables[35] = part->Y(4122);
1539 fCandidateVariables[36] = bachelor->Y(2212);
1540 fCandidateVariables[37] = v0part->Y(310);
1541 fCandidateVariables[38] = v0pos->Y(211);
1542 fCandidateVariables[39] = v0neg->Y(211);
1544 fCandidateVariables[40] = v0part->Eta();
1546 fCandidateVariables[41] = part->DecayLength();
1547 fCandidateVariables[42] = part->DecayLengthV0();
1548 fCandidateVariables[43] = part->Ct(4122);
1549 fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
1551 EBachelor bachCode = kBachInvalid;
1552 EK0S k0SCode = kK0SInvalid;
1554 bachCode = CheckBachelor(part, bachelor, mcArray);
1555 k0SCode = CheckK0S(part, v0part, mcArray);
1558 fCandidateVariables[45] = bachCode;
1559 fCandidateVariables[46] = k0SCode;
1561 Double_t V0KF[3] = {-999999, -999999, -999999};
1562 Double_t errV0KF[3] = {-999999, -999999, -999999};
1563 Double_t LcKF[3] = {-999999, -999999, -999999};
1564 Double_t errLcKF[3] = {-999999, -999999, -999999};
1565 Double_t distances[3] = {-999999, -999999, -999999};
1566 Double_t armPolKF[2] = {-999999, -999999};
1568 if (fCallKFVertexing){
1569 Int_t kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
1570 AliDebug(2, Form("Result from KF = %d", kfResult));
1574 for (Int_t i = 0; i< 3; i++){
1575 Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
1579 fCandidateVariables[47] = V0KF[0];
1580 fCandidateVariables[48] = V0KF[1];
1581 fCandidateVariables[49] = V0KF[2];
1583 fCandidateVariables[50] = errV0KF[0];
1584 fCandidateVariables[51] = errV0KF[1];
1585 fCandidateVariables[52] = errV0KF[2];
1587 fCandidateVariables[53] = LcKF[0];
1588 fCandidateVariables[54] = LcKF[1];
1589 fCandidateVariables[55] = LcKF[2];
1591 fCandidateVariables[56] = errLcKF[0];
1592 fCandidateVariables[57] = errLcKF[1];
1593 fCandidateVariables[58] = errLcKF[2];
1595 fCandidateVariables[59] = distances[0];
1596 fCandidateVariables[60] = distances[1];
1597 fCandidateVariables[61] = distances[2];
1598 fCandidateVariables[62] = armPolKF[0];
1599 fCandidateVariables[63] = armPolKF[1];
1600 fCandidateVariables[64] = v0part->AlphaV0();
1601 fCandidateVariables[65] = v0part->PtArmV0();
1603 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)));
1604 fCandidateVariables[66] = v0pos->GetStatus() & AliESDtrack::kITSrefit;
1605 fCandidateVariables[67] = v0neg->GetStatus() & AliESDtrack::kITSrefit;
1606 fCandidateVariables[68] = v0pos->GetTPCClusterInfo(2, 1);
1607 fCandidateVariables[69] = v0neg->GetTPCClusterInfo(2, 1);
1609 fCandidateVariables[70] = v0part->Xv();
1610 fCandidateVariables[71] = v0part->Yv();
1611 fCandidateVariables[72] = v0part->Zv();
1613 fCandidateVariables[73] = fVtx1->GetX();
1614 fCandidateVariables[74] = fVtx1->GetY();
1615 fCandidateVariables[75] = fVtx1->GetZ();
1617 fCandidateVariables[76] = bachelor->GetITSNcls();
1618 fCandidateVariables[77] = bachelor->HasPointOnITSLayer(0) + bachelor->HasPointOnITSLayer(1);
1620 fCandidateVariables[78] = v0pos->GetITSNcls();
1621 fCandidateVariables[79] = v0pos->HasPointOnITSLayer(0) + v0pos->HasPointOnITSLayer(1);
1623 fCandidateVariables[80] = v0neg->GetITSNcls();
1624 fCandidateVariables[81] = v0neg->HasPointOnITSLayer(0) + v0neg->HasPointOnITSLayer(1);
1626 TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz());
1627 TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz());
1628 TVector3 momTot(part->Px(), part->Py(), part->Pz());
1630 Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
1631 Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
1633 Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2);
1634 Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) );
1635 Double_t ptArmLc = mom1.Perp(momTot);
1637 fCandidateVariables[82] = alphaArmLc;
1638 fCandidateVariables[83] = alphaArmLcCharge;
1639 fCandidateVariables[84] = ptArmLc;
1641 Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // mass K0S PDG
1642 Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); // mass Proton PDG
1643 Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass(); // mass Lc PDG
1645 Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG);
1646 Double_t e = part->E(4122);
1647 Double_t beta = part->P()/e;
1648 Double_t gamma = e/massLcPDG;
1650 Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
1652 fCandidateVariables[85] = cts;
1656 AliDebug(2, "Filling Sgn");
1657 fVariablesTreeSgn->Fill();
1658 fHistoCodesSgn->Fill(bachCode, k0SCode);
1661 if (fFillOnlySgn == kFALSE){
1662 AliDebug(2, "Filling Bkg");
1663 fVariablesTreeBkg->Fill();
1664 fHistoCodesBkg->Fill(bachCode, k0SCode);
1669 fVariablesTreeSgn->Fill();
1677 //________________________________________________________________________
1678 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
1679 Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
1680 Double_t* distances, Double_t* armPolKF) {
1683 // method to perform KF vertexing
1684 // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
1687 Int_t codeKFV0 = -1, codeKFLc = -1;
1689 AliKFVertex primVtxCopy;
1690 Int_t nt = 0, ntcheck = 0;
1691 Double_t pos[3] = {0., 0., 0.};
1694 Int_t contr = fVtx1->GetNContributors();
1695 Double_t covmatrix[6] = {0.};
1696 fVtx1->GetCovarianceMatrix(covmatrix);
1697 Double_t chi2 = fVtx1->GetChi2();
1698 AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
1700 // topological constraint
1701 primVtxCopy = AliKFVertex(primaryESDVtxCopy);
1702 nt = primaryESDVtxCopy.GetNContributors();
1705 Int_t pdg[2] = {211, -211};
1706 Int_t pdgLc[2] = {2212, 310};
1708 Int_t pdgDgV0toDaughters[2] = {211, 211};
1710 Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
1712 // the KF vertex for the V0 has to be built with the prongs of the V0!
1713 Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
1714 AliKFParticle V0, positiveV0KF, negativeV0KF;
1715 Int_t labelsv0daugh[2] = {-1, -1};
1716 Int_t idv0daugh[2] = {-1, -1};
1717 AliExternalTrackParam* esdv0Daugh1 = 0x0;
1718 AliExternalTrackParam* esdv0Daugh2 = 0x0;
1719 for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1720 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1722 AliDebug(2, "No V0 daughters available");
1725 Double_t xyz[3], pxpypz[3], cv[21];
1727 aodTrack->GetXYZ(xyz);
1728 aodTrack->PxPyPz(pxpypz);
1729 aodTrack->GetCovarianceXYZPxPyPz(cv);
1730 sign = aodTrack->Charge();
1731 AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
1733 if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1734 else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1735 labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
1736 idv0daugh[ipr] = aodTrack->GetID();
1737 if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
1739 //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
1741 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1742 if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
1743 positiveV0KF = daughterKF;
1746 negativeV0KF = daughterKF;
1750 Double_t xn=0., xp=0.;//, dca;
1751 AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
1752 // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
1754 AliExternalTrackParam tr1(*esdv0Daugh1);
1755 AliExternalTrackParam tr2(*esdv0Daugh2);
1756 tr1.PropagateTo(xn, fBField);
1757 tr2.PropagateTo(xp, fBField);
1759 AliKFParticle daughterKF1(tr1, 211);
1760 AliKFParticle daughterKF2(tr2, 211);
1761 V0.AddDaughter(positiveV0KF);
1762 V0.AddDaughter(negativeV0KF);
1763 //V0.AddDaughter(daughterKF1);
1764 //V0.AddDaughter(daughterKF2);
1770 // Checking the quality of the KF V0 vertex
1771 if( V0.GetNDF() < 1 ) {
1772 //Printf("Number of degrees of freedom < 1, continuing");
1775 if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
1776 //Printf("Chi2 per DOF too big, continuing");
1780 if(ftopoConstraint && nt > 0){
1781 for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1782 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1783 //* subtruct daughters from primary vertex
1784 if(!aodTrack->GetUsedForPrimVtxFit()) {
1785 //Printf("Track %d was not used for primary vertex, continuing", i);
1788 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1789 primVtxCopy -= daughterKF;
1794 // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
1796 if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
1797 //Printf("Deviation from vertex too big, continuing");
1802 //* Get V0 invariant mass
1803 Double_t massV0 = 999999, sigmaMassV0 = 999999;
1804 Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
1807 codeKFV0 = 1; // Mass not ok
1808 if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
1810 else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
1812 fHistoMassKFV0->Fill(massV0, sigmaMassV0);
1814 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]);
1815 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]);
1817 Printf("Vertices: KF: x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ());
1818 Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv());
1820 //Printf("Got MC vtx for V0");
1821 if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) {
1822 AliAODMCParticle* tmpdaughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
1823 AliAODMCParticle* tmpdaughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
1824 if (!tmpdaughv01 && labelsv0daugh[0] > 0){
1825 AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
1827 if (!tmpdaughv02 && labelsv0daugh[1] > 0){
1828 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1831 Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays
1832 Double_t yPionMC = tmpdaughv01->Yv();
1833 Double_t zPionMC = tmpdaughv01->Zv();
1834 //Printf("Got MC vtx for Pion");
1835 Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
1839 Printf("Not a true V0");
1841 //massV0=-1;//return -1;// !!!!
1843 // now use what we just try with the bachelor, to build the Lc
1845 // topological constraint
1846 nt = primVtxCopy.GetNContributors();
1849 Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
1851 Int_t labelsLcdaugh[2] = {-1, -1};
1852 labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
1853 labelsLcdaugh[1] = mcLabelV0;
1855 if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
1856 AliKFParticle daughterKFLc(*bach, pdgLc[0]);
1857 Lc.AddDaughter(daughterKFLc);
1858 TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310);
1859 Double_t massPDGK0S = particlePDG->Mass();
1860 V0.SetMassConstraint(massPDGK0S);
1862 if( Lc.GetNDF() < 1 ) {
1863 AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
1866 if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1867 AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1871 if(ftopoConstraint && nt > 0){
1872 //* subtruct daughters from primary vertex
1873 if(!bach->GetUsedForPrimVtxFit()) {
1874 AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
1877 primVtxCopy -= daughterKFLc;
1880 /* the V0 was added above, so it is ok to remove it without checking
1881 if(!V0->GetUsedForPrimVtxFit()) {
1882 Printf("Lc: V0 was not used for primary vertex, continuing");
1886 //primVtxCopy -= V0;
1890 // Check Lc Chi^2 deviation from primary vertex
1892 if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
1893 AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
1897 if(ftopoConstraint){
1899 // Add Lc to primary vertex to improve the primary vertex resolution
1901 Lc.SetProductionVertex(primVtxCopy);
1906 if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1907 AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1911 if(ftopoConstraint){
1912 V0.SetProductionVertex(Lc);
1915 // After setting the vertex of the V0, getting/filling some info
1917 //* Get V0 decayLength
1918 Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
1919 Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
1921 if (sigmaDecayLengthV0 > 1e19) {
1922 codeKFV0 = 3; // DecayLength not ok
1924 codeKFV0 = 6; // DecayLength and Mass not ok
1925 if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
1927 else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
1930 fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
1932 //* Get V0 life time
1933 Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
1934 Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
1936 if (sigmaLifeTimeV0 > 1e19) {
1937 codeKFV0 = 4; // LifeTime not ok
1938 if (sigmaDecayLengthV0 > 1e19) {
1939 codeKFV0 = 9; // LifeTime and DecayLength not ok
1941 codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
1942 if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
1944 else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
1946 else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
1947 codeKFV0 = 7; // LifeTime and Mass not ok
1948 if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
1950 else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
1953 fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
1955 if (codeKFV0 == -1) codeKFV0 = 0;
1956 fHistoKFV0->Fill(codeKFV0);
1958 AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
1960 fHistoMassV0All->Fill(massV0);
1961 fHistoDecayLengthV0All->Fill(decayLengthV0);
1962 fHistoLifeTimeV0All->Fill(lifeTimeV0);
1964 Double_t qtAlphaV0[2] = {0., 0.};
1965 Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
1966 positiveV0KF.TransportToPoint(vtxV0KF);
1967 negativeV0KF.TransportToPoint(vtxV0KF);
1968 V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
1969 AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
1970 fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
1971 fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
1972 armPolKF[0] = qtAlphaV0[1];
1973 armPolKF[1] = qtAlphaV0[0];
1975 // Checking MC info for V0
1977 AliAODMCParticle *motherV0 = 0x0;
1978 AliAODMCParticle *daughv01 = 0x0;
1979 AliAODMCParticle *daughv02 = 0x0;
1982 daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
1983 daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
1984 if (!daughv01 && labelsv0daugh[0] > 0){
1985 AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
1988 if (!daughv02 && labelsv0daugh[1] > 0){
1989 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1993 if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
1994 AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
1995 motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
1996 if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
1997 if( motherV0->GetNDaughters() == 2 ){
1998 fHistoMassV0True->Fill(massV0);
1999 fHistoDecayLengthV0True->Fill(decayLengthV0);
2000 fHistoLifeTimeV0True->Fill(lifeTimeV0);
2001 fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
2002 if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
2003 fHistoMassV0TrueK0S->Fill(massV0);
2004 fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
2005 fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
2006 fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
2009 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()));
2011 else if (!motherV0){
2012 AliDebug(3, "could not access MC info for mother, continuing");
2015 AliDebug(3, "MC mother is a gluon, continuing");
2019 AliDebug(3, "Background V0!");
2025 AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
2029 //* Get Lc invariant mass
2030 Double_t massLc = 999999, sigmaMassLc= 999999;
2031 Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
2033 AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
2035 codeKFLc = 1; // Mass not ok
2036 if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
2038 else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
2040 fHistoMassKFLc->Fill(massLc, sigmaMassLc);
2042 //* Get Lc Decay length
2043 Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
2044 Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
2046 AliDebug(3, "----> Lc: Could not get decay length, and sigma");
2047 if (sigmaDecayLengthLc > 1e19) {
2048 codeKFLc = 3; // DecayLength not ok
2050 codeKFLc = 6; // DecayLength and Mass not ok
2051 if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
2053 else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
2056 AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
2058 fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
2060 //* Get Lc life time
2061 Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
2062 Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
2064 AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
2065 if (sigmaLifeTimeLc > 1e19) {
2066 codeKFLc = 4; // LifeTime not ok
2067 if (sigmaDecayLengthLc > 1e19) {
2068 codeKFLc = 9; // LifeTime and DecayLength not ok
2070 codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
2071 if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
2073 else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
2075 else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
2076 codeKFLc = 7; // LifeTime and Mass not ok
2077 if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
2079 else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
2083 fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
2085 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));
2087 if (codeKFLc == -1) codeKFLc = 0;
2088 fHistoKFLc->Fill(codeKFLc);
2090 fHistoKF->Fill(codeKFV0, codeKFLc);
2092 // here we fill the histgrams for all the reconstructed KF vertices for the cascade
2093 fHistoMassLcAll->Fill(massLc);
2094 fHistoDecayLengthLcAll->Fill(decayLengthLc);
2095 fHistoLifeTimeLcAll->Fill(lifeTimeLc);
2097 fHistoMassV0fromLcAll->Fill(massV0);
2098 fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
2099 fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
2101 Double_t xV0 = V0.GetX();
2102 Double_t yV0 = V0.GetY();
2103 Double_t zV0 = V0.GetZ();
2105 Double_t xLc = Lc.GetX();
2106 Double_t yLc = Lc.GetY();
2107 Double_t zLc = Lc.GetZ();
2109 Double_t xPrimVtx = primVtxCopy.GetX();
2110 Double_t yPrimVtx = primVtxCopy.GetY();
2111 Double_t zPrimVtx = primVtxCopy.GetZ();
2113 Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2114 (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2115 (zPrimVtx - zLc) * (zPrimVtx - zLc));
2117 Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2118 (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2119 (zPrimVtx - zV0) * (zPrimVtx - zV0));
2121 Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2122 (yLc - yV0)*(yLc - yV0) +
2123 (zLc - zV0)*(zLc - zV0));
2125 //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2127 fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2128 fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2129 fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2131 distances[0] = distanceLcToPrimVtx;
2132 distances[1] = distanceV0ToPrimVtx;
2133 distances[2] = distanceV0ToLc;
2137 AliAODMCParticle *daughv01Lc = 0x0;
2138 AliAODMCParticle *K0S = 0x0;
2139 AliAODMCParticle *daughv02Lc = 0x0;
2141 if (labelsLcdaugh[0] >= 0) {
2142 // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2143 daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2145 AliDebug(3, "Could not access MC info for first daughter of Lc");
2149 AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2150 if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2153 else { // The bachelor is a fake
2157 if (labelsLcdaugh[1] >= 0) {
2158 //Printf("Getting K0S");
2159 K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2161 AliDebug(3, "Could not access MC info for second daughter of Lc");
2165 if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2169 AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2173 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
2174 if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
2175 Int_t iK0 = K0S->GetMother();
2177 Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2179 else { // The K0S has a mother
2180 daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2182 AliDebug(3, "Could not access MC info for second daughter of Lc");
2184 else { // we can access safely the K0S mother in the MC
2185 if( daughv01Lc && (daughv01Lc->GetMother() == daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){ // This is a true cascade! bachelor and V0 come from the same mother
2186 //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
2187 AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
2188 Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
2189 if (motherLc) pdgMum = motherLc->GetPdgCode();
2190 if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
2191 if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
2192 AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
2194 if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
2195 fHistoMassLcTrue->Fill(massLc);
2196 fHistoDecayLengthLcTrue->Fill(decayLengthLc);
2197 fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
2198 fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
2200 fHistoMassV0fromLcTrue->Fill(massV0);
2201 fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2202 fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2204 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...
2205 AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
2207 fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2208 fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2210 fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2211 fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2212 fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2214 fHistoMassLcSgn->Fill(massLc);
2215 fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2217 fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2218 fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2220 fHistoMassV0fromLcSgn->Fill(massV0);
2221 fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2222 fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2225 isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
2228 // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
2229 // First, we compare the vtx of the Lc
2230 Double_t xLcMC = motherLc->Xv();
2231 Double_t yLcMC = motherLc->Yv();
2232 Double_t zLcMC = motherLc->Zv();
2233 //Printf("Got MC vtx for Lc");
2234 Double_t xProtonMC = daughv01Lc->Xv();
2235 Double_t yProtonMC = daughv01Lc->Yv();
2236 Double_t zProtonMC = daughv01Lc->Zv();
2237 //Printf("Got MC vtx for Proton");
2239 Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2240 (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2241 (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2243 // Then, we compare the vtx of the K0s
2244 Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
2245 Double_t yV0MC = motherV0->Yv();
2246 Double_t zV0MC = motherV0->Zv();
2248 //Printf("Got MC vtx for V0");
2249 Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2250 Double_t yPionMC = daughv01->Yv();
2251 Double_t zPionMC = daughv01->Zv();
2252 //Printf("Got MC vtx for Pion");
2253 Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2255 Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2256 (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2257 (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2259 // Then, the K0S vertex wrt primary vertex
2261 Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2262 (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2263 (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2265 fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2266 fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2267 fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2269 } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
2270 else if (!motherLc){
2271 AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
2274 AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2276 } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2278 isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2280 } // closing: else { // we can access safely the K0S mother in the MC
2281 } // closing: else { // The K0S has a mother
2282 } // closing isMCLcok
2283 } // closing !isBkgLc
2284 } // closing fUseMCInfo
2286 //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2287 if ( retMV0 == 0 && retMLc == 0){
2289 errV0KF[0] = sigmaMassV0;
2290 V0KF[1] = decayLengthV0;
2291 errV0KF[1] = sigmaDecayLengthV0;
2292 V0KF[2] = lifeTimeV0;
2293 errV0KF[2] = sigmaLifeTimeV0;
2295 errLcKF[0] = sigmaMassLc;
2296 LcKF[1] = decayLengthLc;
2297 errLcKF[1] = sigmaDecayLengthLc;
2298 LcKF[2] = lifeTimeLc;
2299 errLcKF[2] = sigmaLifeTimeLc;
2305 //________________________________________________________________________
2306 AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
2307 AliAODTrack* bachelor,
2308 TClonesArray *mcArray ){
2310 //Printf("In CheckBachelor");
2312 // function to check if the bachelor is a p, if it is a p but not from Lc
2313 // to be filled for background candidates
2315 Int_t label = bachelor->GetLabel();
2320 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
2321 if (!mcpart) return kBachInvalid;
2322 Int_t pdg = mcpart->PdgCode();
2323 if (TMath::Abs(pdg) != 2212) {
2324 AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
2325 return kBachNoProton;
2327 else { // it is a proton
2328 //Int_t labelLc = part->GetLabel();
2329 Int_t labelLc = FindLcLabel(part, mcArray);
2330 //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
2331 Int_t bachelorMotherLabel = mcpart->GetMother();
2332 //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
2333 if (bachelorMotherLabel == -1) {
2334 return kBachPrimary;
2336 AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2337 if (!bachelorMother) return kBachInvalid;
2338 Int_t pdgMother = bachelorMother->PdgCode();
2339 if (TMath::Abs(pdgMother) != 4122) {
2340 AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2341 return kBachNoLambdaMother;
2343 else { // it comes from Lc
2344 if (labelLc != bachelorMotherLabel){
2345 //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));
2346 AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2347 return kBachDifferentLambdaMother;
2349 else { // it comes from the correct Lc
2350 AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2351 return kBachCorrectLambdaMother;
2356 return kBachInvalid;
2360 //________________________________________________________________________
2361 AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
2363 //AliAODTrack* v0part,
2364 TClonesArray *mcArray ){
2366 // function to check if the K0Spart is a p, if it is a p but not from Lc
2367 // to be filled for background candidates
2369 //Printf(" CheckK0S");
2371 //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
2372 //Int_t label = v0part->GetLabel();
2373 Int_t labelFind = FindV0Label(v0part, mcArray);
2374 //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
2375 if (labelFind == -1) {
2379 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
2380 if (!mcpart) return kK0SInvalid;
2381 Int_t pdg = mcpart->PdgCode();
2382 if (TMath::Abs(pdg) != 310) {
2383 AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2384 //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2387 else { // it is a K0S
2388 //Int_t labelLc = part->GetLabel();
2389 Int_t labelLc = FindLcLabel(part, mcArray);
2390 Int_t K0SpartMotherLabel = mcpart->GetMother();
2391 if (K0SpartMotherLabel == -1) {
2392 return kK0SWithoutMother;
2394 AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
2395 if (!K0SpartMother) return kK0SInvalid;
2396 Int_t pdgMotherK0S = K0SpartMother->PdgCode();
2397 if (TMath::Abs(pdgMotherK0S) != 311) {
2398 AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2399 //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2400 return kK0SNotFromK0;
2402 else { // the K0S comes from a K0
2403 Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
2404 if (K0MotherLabel == -1) {
2407 AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
2408 if (!K0Mother) return kK0SInvalid;
2409 Int_t pdgK0Mother = K0Mother->PdgCode();
2410 if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
2411 AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2412 //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2413 return kK0NoLambdaMother;
2415 else { // the K0 comes from Lc
2416 if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
2417 AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2418 //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2419 return kK0DifferentLambdaMother;
2421 else { // it comes from the correct Lc
2422 AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2423 //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2424 return kK0CorrectLambdaMother;
2434 //----------------------------------------------------------------------------
2435 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2438 //Printf(" FindV0Label");
2440 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2442 Int_t labMother[2]={-1, -1};
2443 AliAODMCParticle *part=0;
2444 AliAODMCParticle *mother=0;
2445 Int_t dgLabels = -1;
2447 Int_t ndg = v0part->GetNDaughters();
2449 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2452 for(Int_t i = 0; i < 2; i++) {
2453 AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
2454 dgLabels = trk->GetLabel();
2455 if (dgLabels == -1) {
2456 //printf("daughter with negative label %d\n", dgLabels);
2459 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2461 //printf("no MC particle\n");
2464 labMother[i] = part->GetMother();
2465 if (labMother[i] != -1){
2466 mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2468 //printf("no MC mother particle\n");
2477 if (labMother[0] == labMother[1]) return labMother[0];
2483 //----------------------------------------------------------------------------
2484 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
2487 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2489 //Printf(" FindLcLabel");
2491 AliAODMCParticle *part=0;
2492 AliAODMCParticle *mother=0;
2493 AliAODMCParticle *grandmother=0;
2494 Int_t dgLabels = -1;
2496 Int_t ndg = cascade->GetNDaughters();
2498 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2501 // daughter 0 --> bachelor, daughter 1 --> V0
2502 AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
2503 dgLabels = trk->GetLabel();
2504 if (dgLabels == -1 ) {
2505 //printf("daughter with negative label %d\n", dgLabels);
2508 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2510 //printf("no MC particle\n");
2513 Int_t labMotherBach = part->GetMother();
2514 if (labMotherBach == -1){
2517 mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2519 //printf("no MC mother particle\n");
2523 AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
2524 dgLabels = FindV0Label(v0, mcArray);
2525 if (dgLabels == -1 ) {
2526 //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
2529 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2531 //printf("no MC particle\n");
2534 Int_t labMotherv0 = part->GetMother();
2535 if (labMotherv0 == -1){
2538 mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2540 //printf("no MC mother particle\n");
2543 Int_t labGrandMotherv0 = mother->GetMother();
2544 if (labGrandMotherv0 == -1){
2547 grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2549 //printf("no MC mother particle\n");
2553 //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2554 if (labGrandMotherv0 == labMotherBach) return labMotherBach;