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)
196 //___________________________________________________________________________
197 AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
198 AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
199 AliAnalysisTaskSE(name),
205 fIsK0sAnalysis(kFALSE),
209 fUseOnTheFlyV0(useOnTheFly),
210 fIsEventSelected(kFALSE),
211 fVariablesTreeSgn(0),
212 fVariablesTreeBkg(0),
213 fCandidateVariables(),
218 fFillOnlySgn(kFALSE),
219 fHistoLcBeforeCuts(0),
220 fHistoFiducialAcceptance(0),
223 fHistoLcpKpiBeforeCuts(0),
226 fHistoDistanceLcToPrimVtx(0),
227 fHistoDistanceV0ToPrimVtx(0),
228 fHistoDistanceV0ToLc(0),
230 fHistoDistanceLcToPrimVtxSgn(0),
231 fHistoDistanceV0ToPrimVtxSgn(0),
232 fHistoDistanceV0ToLcSgn(0),
234 fHistoVtxLcResidualToPrimVtx(0),
235 fHistoVtxV0ResidualToPrimVtx(0),
236 fHistoVtxV0ResidualToLc(0),
239 fHistoDecayLengthV0All(0),
240 fHistoLifeTimeV0All(0),
243 fHistoDecayLengthV0True(0),
244 fHistoLifeTimeV0True(0),
246 fHistoMassV0TrueFromAOD(0),
248 fHistoMassV0TrueK0S(0),
249 fHistoDecayLengthV0TrueK0S(0),
250 fHistoLifeTimeV0TrueK0S(0),
252 fHistoMassV0TrueK0SFromAOD(0),
255 fHistoDecayLengthLcAll(0),
256 fHistoLifeTimeLcAll(0),
259 fHistoDecayLengthLcTrue(0),
260 fHistoLifeTimeLcTrue(0),
262 fHistoMassLcTrueFromAOD(0),
264 fHistoMassV0fromLcAll(0),
265 fHistoDecayLengthV0fromLcAll(0),
266 fHistoLifeTimeV0fromLcAll(0),
268 fHistoMassV0fromLcTrue(0),
269 fHistoDecayLengthV0fromLcTrue(0),
270 fHistoLifeTimeV0fromLcTrue(0),
273 fHistoMassLcSgnFromAOD(0),
274 fHistoDecayLengthLcSgn(0),
275 fHistoLifeTimeLcSgn(0),
277 fHistoMassV0fromLcSgn(0),
278 fHistoDecayLengthV0fromLcSgn(0),
279 fHistoLifeTimeV0fromLcSgn(0),
286 fHistoDecayLengthKFV0(0),
287 fHistoLifeTimeKFV0(0),
290 fHistoDecayLengthKFLc(0),
291 fHistoLifeTimeKFLc(0),
293 fHistoArmenterosPodolanskiV0KF(0),
294 fHistoArmenterosPodolanskiV0KFSgn(0),
295 fHistoArmenterosPodolanskiV0AOD(0),
296 fHistoArmenterosPodolanskiV0AODSgn(0),
300 ftopoConstraint(kTRUE),
301 fCallKFVertexing(kFALSE),
302 fKeepingOnlyHIJINGBkg(kFALSE),
305 fCutKFChi2NDF(999999.),
306 fCutKFDeviationFromVtx(999999.),
307 fCutKFDeviationFromVtxV0(0.),
310 fKeepingOnlyPYTHIABkg(kFALSE)
314 // Constructor. Initialization of Inputs and Outputs
316 Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
318 DefineOutput(1, TList::Class()); // Tree signal + Tree Bkg + histoEvents
319 DefineOutput(2, AliNormalizationCounter::Class()); // normalization counter object
320 DefineOutput(3, TList::Class()); // Cuts
321 DefineOutput(4, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
322 DefineOutput(5, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
323 DefineOutput(6, TList::Class()); // Tree signal + Tree Bkg + histoEvents
327 //___________________________________________________________________________
328 AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
332 Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
362 if(fVariablesTreeSgn){
363 delete fVariablesTreeSgn;
364 fVariablesTreeSgn = 0;
367 if(fVariablesTreeBkg){
368 delete fVariablesTreeBkg;
369 fVariablesTreeBkg = 0;
383 //_________________________________________________
384 void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
389 fIsEventSelected=kFALSE;
391 if (fDebug > 1) AliInfo("Init");
393 fListCuts = new TList();
394 fListCuts->SetOwner();
395 fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
396 PostData(3,fListCuts);
398 if (fUseMCInfo && (fKeepingOnlyHIJINGBkg || fKeepingOnlyPYTHIABkg)) fUtils = new AliVertexingHFUtils();
403 //________________________________________ terminate ___________________________
404 void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
406 // The Terminate() function is the last function to be called during
407 // a query. It always runs on the client, it can be used to present
408 // the results graphically or save the results to file.
410 //AliInfo("Terminate","");
411 AliAnalysisTaskSE::Terminate();
413 fOutput = dynamic_cast<TList*> (GetOutputData(1));
415 AliError("fOutput not available");
418 fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
420 AliError("fOutputKF not available");
427 //___________________________________________________________________________
428 void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
430 AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
434 fOutput = new TList();
436 fOutput->SetName("listTrees");
438 // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
439 const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
440 fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
441 fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
443 fCandidateVariables = new Float_t [nVar];
444 TString * fCandidateVariableNames = new TString[nVar];
445 fCandidateVariableNames[0]="massLc2K0Sp";
446 fCandidateVariableNames[1]="massLc2Lambdapi";
447 fCandidateVariableNames[2]="massK0S";
448 fCandidateVariableNames[3]="massLambda";
449 fCandidateVariableNames[4]="massLambdaBar";
450 fCandidateVariableNames[5]="cosPAK0S";
451 fCandidateVariableNames[6]="dcaV0";
452 fCandidateVariableNames[7]="tImpParBach";
453 fCandidateVariableNames[8]="tImpParV0";
454 fCandidateVariableNames[9]="nSigmaTPCpr";
455 fCandidateVariableNames[10]="nSigmaTPCpi";
456 fCandidateVariableNames[11]="nSigmaTPCka";
457 fCandidateVariableNames[12]="nSigmaTOFpr";
458 fCandidateVariableNames[13]="nSigmaTOFpi";
459 fCandidateVariableNames[14]="nSigmaTOFka";
460 fCandidateVariableNames[15]="bachelorPt";
461 fCandidateVariableNames[16]="V0positivePt";
462 fCandidateVariableNames[17]="V0negativePt";
463 fCandidateVariableNames[18]="dcaV0pos";
464 fCandidateVariableNames[19]="dcaV0neg";
465 fCandidateVariableNames[20]="v0Pt";
466 fCandidateVariableNames[21]="massGamma";
467 fCandidateVariableNames[22]="LcPt";
468 fCandidateVariableNames[23]="combinedProtonProb";
469 fCandidateVariableNames[24]="LcEta";
470 fCandidateVariableNames[25]="V0positiveEta";
471 fCandidateVariableNames[26]="V0negativeEta";
472 fCandidateVariableNames[27]="TPCProtonProb";
473 fCandidateVariableNames[28]="TOFProtonProb";
474 fCandidateVariableNames[29]="bachelorEta";
475 fCandidateVariableNames[30]="LcP";
476 fCandidateVariableNames[31]="bachelorP";
477 fCandidateVariableNames[32]="v0P";
478 fCandidateVariableNames[33]="V0positiveP";
479 fCandidateVariableNames[34]="V0negativeP";
480 fCandidateVariableNames[35]="LcY";
481 fCandidateVariableNames[36]="v0Y";
482 fCandidateVariableNames[37]="bachelorY";
483 fCandidateVariableNames[38]="V0positiveY";
484 fCandidateVariableNames[39]="V0negativeY";
485 fCandidateVariableNames[40]="v0Eta";
486 fCandidateVariableNames[41]="DecayLengthLc";
487 fCandidateVariableNames[42]="DecayLengthK0S";
488 fCandidateVariableNames[43]="CtLc";
489 fCandidateVariableNames[44]="CtK0S";
490 fCandidateVariableNames[45]="bachCode";
491 fCandidateVariableNames[46]="k0SCode";
493 fCandidateVariableNames[47]="V0KFmass";
494 fCandidateVariableNames[48]="V0KFdecayLength";
495 fCandidateVariableNames[49]="V0KFlifeTime";
497 fCandidateVariableNames[50]="V0KFmassErr";
498 fCandidateVariableNames[51]="V0KFdecayTimeErr";
499 fCandidateVariableNames[52]="V0KFlifeTimeErr";
501 fCandidateVariableNames[53]="LcKFmass";
502 fCandidateVariableNames[54]="LcKFdecayLength";
503 fCandidateVariableNames[55]="LcKFlifeTime";
505 fCandidateVariableNames[56]="LcKFmassErr";
506 fCandidateVariableNames[57]="LcKFdecayTimeErr";
507 fCandidateVariableNames[58]="LcKFlifeTimeErr";
509 fCandidateVariableNames[59]="LcKFDistToPrimVtx";
510 fCandidateVariableNames[60]="V0KFDistToPrimVtx";
511 fCandidateVariableNames[61]="V0KFDistToLc";
512 fCandidateVariableNames[62]="alphaArmKF";
513 fCandidateVariableNames[63]="ptArmKF";
514 fCandidateVariableNames[64]="alphaArm";
515 fCandidateVariableNames[65]="ptArm";
517 fCandidateVariableNames[66]="ITSrefitV0pos";
518 fCandidateVariableNames[67]="ITSrefitV0neg";
520 fCandidateVariableNames[68]="TPCClV0pos";
521 fCandidateVariableNames[69]="TPCClV0neg";
523 fCandidateVariableNames[70]="v0Xcoord";
524 fCandidateVariableNames[71]="v0Ycoord";
525 fCandidateVariableNames[72]="v0Zcoord";
526 fCandidateVariableNames[73]="primVtxX";
527 fCandidateVariableNames[74]="primVtxY";
528 fCandidateVariableNames[75]="primVtxZ";
530 fCandidateVariableNames[76]="ITSclBach";
531 fCandidateVariableNames[77]="SPDclBach";
533 fCandidateVariableNames[78]="ITSclV0pos";
534 fCandidateVariableNames[79]="SPDclV0pos";
535 fCandidateVariableNames[80]="ITSclV0neg";
536 fCandidateVariableNames[81]="SPDclV0neg";
538 fCandidateVariableNames[82]="alphaArmLc";
539 fCandidateVariableNames[83]="alphaArmLcCharge";
540 fCandidateVariableNames[84]="ptArmLc";
542 fCandidateVariableNames[85]="CosThetaStar";
545 for(Int_t ivar=0; ivar<nVar; ivar++){
546 fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
547 fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
550 fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
551 TString labelEv[2] = {"NotSelected", "Selected"};
552 for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
553 fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
556 fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
558 fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
559 TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
560 for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
561 fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
564 fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
565 TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
566 for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
567 fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
568 fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
571 fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
572 TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
573 for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
574 fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
577 fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
578 fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
580 TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
581 "kBachDifferentLambdaMother", "kBachCorrectLambdaMother"};
582 TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
583 "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
585 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
586 fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
587 fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
589 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
590 fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
591 fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
594 fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
595 for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
596 fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
599 fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 4, -0.5, 3.5);
600 TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"};
601 for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
602 fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
605 //fOutput->Add(fVariablesTreeSgn);
606 //fOutput->Add(fVariablesTreeBkg);
608 fOutput->Add(fHistoEvents);
609 fOutput->Add(fHistoLc);
610 fOutput->Add(fHistoLcOnTheFly);
611 fOutput->Add(fHistoLcBeforeCuts);
612 fOutput->Add(fHistoFiducialAcceptance);
613 fOutput->Add(fHistoCodesSgn);
614 fOutput->Add(fHistoCodesBkg);
615 fOutput->Add(fHistoLcpKpiBeforeCuts);
616 fOutput->Add(fHistoBackground);
618 PostData(1, fOutput);
619 PostData(4, fVariablesTreeSgn);
620 PostData(5, fVariablesTreeBkg);
621 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
622 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
623 fPIDResponse = inputHandler->GetPIDResponse();
625 if (fAnalCuts->GetIsUsePID()){
627 fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
628 fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
629 fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
630 fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
631 fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
632 fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
634 fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
637 // Setting properties of PID
638 fPIDCombined=new AliPIDCombined;
639 fPIDCombined->SetDefaultTPCPriors();
640 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
641 //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
643 fCounter = new AliNormalizationCounter("NormalizationCounter");
645 PostData(2, fCounter);
647 // Histograms from KF
649 if (fCallKFVertexing){
650 fHistoDistanceLcToPrimVtx = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
651 fHistoDistanceV0ToPrimVtx = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
652 fHistoDistanceV0ToLc = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
654 fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
655 fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
656 fHistoDistanceV0ToLcSgn = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
658 fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
659 fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
660 fHistoVtxV0ResidualToLc = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
662 fHistoMassV0All = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
663 fHistoDecayLengthV0All = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
664 fHistoLifeTimeV0All = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
666 fHistoMassV0True = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
667 fHistoDecayLengthV0True = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
668 fHistoLifeTimeV0True = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
670 fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
672 fHistoMassV0TrueK0S = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
673 fHistoDecayLengthV0TrueK0S = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
674 fHistoLifeTimeV0TrueK0S = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
676 fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
678 fHistoMassLcAll = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
679 fHistoDecayLengthLcAll = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
680 fHistoLifeTimeLcAll = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
682 fHistoMassLcTrue = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
683 fHistoDecayLengthLcTrue = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
684 fHistoLifeTimeLcTrue = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
686 fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
688 fHistoMassV0fromLcAll = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
689 fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
690 fHistoLifeTimeV0fromLcAll = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
692 fHistoMassV0fromLcTrue = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
693 fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
694 fHistoLifeTimeV0fromLcTrue = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
696 fHistoMassLcSgn = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
697 fHistoMassLcSgnFromAOD = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
698 fHistoDecayLengthLcSgn = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
699 fHistoLifeTimeLcSgn = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
701 fHistoMassV0fromLcSgn = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
702 fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
703 fHistoLifeTimeV0fromLcSgn = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
705 fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
706 fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
707 fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
708 TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
709 "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
710 "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
711 "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
713 for (Int_t ibin = 1; ibin <=16; ibin++){
714 fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
715 fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
716 fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
717 fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
720 fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
721 fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
722 fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
724 fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
725 fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
726 fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
728 fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
729 fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
730 fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
731 fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
734 fOutputKF = new TList();
735 fOutputKF->SetOwner();
736 fOutputKF->SetName("listHistoKF");
738 if (fCallKFVertexing){
739 fOutputKF->Add(fHistoDistanceLcToPrimVtx);
740 fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
741 fOutputKF->Add(fHistoDistanceV0ToLc);
743 fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
744 fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
745 fOutputKF->Add(fHistoDistanceV0ToLcSgn);
747 fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
748 fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
749 fOutputKF->Add(fHistoVtxV0ResidualToLc);
751 fOutputKF->Add(fHistoMassV0All);
752 fOutputKF->Add(fHistoDecayLengthV0All);
753 fOutputKF->Add(fHistoLifeTimeV0All);
755 fOutputKF->Add(fHistoMassV0True);
756 fOutputKF->Add(fHistoDecayLengthV0True);
757 fOutputKF->Add(fHistoLifeTimeV0True);
759 fOutputKF->Add(fHistoMassV0TrueFromAOD);
761 fOutputKF->Add(fHistoMassV0TrueK0S);
762 fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
763 fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
765 fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
767 fOutputKF->Add(fHistoMassLcAll);
768 fOutputKF->Add(fHistoDecayLengthLcAll);
769 fOutputKF->Add(fHistoLifeTimeLcAll);
771 fOutputKF->Add(fHistoMassLcTrue);
772 fOutputKF->Add(fHistoDecayLengthLcTrue);
773 fOutputKF->Add(fHistoLifeTimeLcTrue);
775 fOutputKF->Add(fHistoMassLcTrueFromAOD);
777 fOutputKF->Add(fHistoMassV0fromLcAll);
778 fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
779 fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
781 fOutputKF->Add(fHistoMassV0fromLcTrue);
782 fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
783 fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
785 fOutputKF->Add(fHistoMassLcSgn);
786 fOutputKF->Add(fHistoMassLcSgnFromAOD);
787 fOutputKF->Add(fHistoDecayLengthLcSgn);
788 fOutputKF->Add(fHistoLifeTimeLcSgn);
790 fOutputKF->Add(fHistoMassV0fromLcSgn);
791 fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
792 fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
794 fOutputKF->Add(fHistoKF);
795 fOutputKF->Add(fHistoKFV0);
796 fOutputKF->Add(fHistoKFLc);
798 fOutputKF->Add(fHistoMassKFV0);
799 fOutputKF->Add(fHistoDecayLengthKFV0);
800 fOutputKF->Add(fHistoLifeTimeKFV0);
802 fOutputKF->Add(fHistoMassKFLc);
803 fOutputKF->Add(fHistoDecayLengthKFLc);
804 fOutputKF->Add(fHistoLifeTimeKFLc);
806 fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
807 fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
808 fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
809 fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
812 PostData(6, fOutputKF);
817 //_________________________________________________
818 void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
822 AliError("NO EVENT FOUND!");
827 Printf("Processing event = %d", fCurrentEvent);
828 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
829 TClonesArray *arrayLctopKos=0;
831 TClonesArray *array3Prong = 0;
833 if (!aodEvent && AODEvent() && IsStandardAOD()) {
834 // In case there is an AOD handler writing a standard AOD, use the AOD
835 // event in memory rather than the input (ESD) event.
836 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
837 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
838 // have to taken from the AOD event hold by the AliAODExtension
839 AliAODHandler* aodHandler = (AliAODHandler*)
840 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
842 if (aodHandler->GetExtensions()) {
843 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
844 AliAODEvent *aodFromExt = ext->GetAOD();
845 arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
847 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
850 arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
852 array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
855 if ( !fUseMCInfo && fIspA) {
856 fAnalCuts->SetTriggerClass("");
857 fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
860 fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
862 // AOD primary vertex
863 fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
865 if (fVtx1->GetNContributors()<1) return;
867 fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
869 if ( !fIsEventSelected ) {
870 fHistoEvents->Fill(0);
871 return; // don't take into account not selected events
873 fHistoEvents->Fill(1);
875 // Setting magnetic field for KF vertexing
876 fBField = aodEvent->GetMagneticField();
877 AliKFParticle::SetField(fBField);
880 TClonesArray *mcArray = 0;
881 AliAODMCHeader *mcHeader=0;
884 // MC array need for maching
885 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
887 AliError("Could not find Monte-Carlo in AOD");
891 mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
893 AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
898 Int_t nSelectedAnal = 0;
899 if (fIsK0sAnalysis) {
900 MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
901 nSelectedAnal, fAnalCuts,
902 array3Prong, mcHeader);
904 fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
906 PostData(1, fOutput);
907 PostData(2, fCounter);
908 PostData(4, fVariablesTreeSgn);
909 PostData(5, fVariablesTreeBkg);
910 PostData(6, fOutputKF);
914 //-------------------------------------------------------------------------------
915 void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
916 TClonesArray *mcArray,
917 Int_t &nSelectedAnal,
918 AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
919 AliAODMCHeader* aodheader){
920 //Lc prong needed to MatchToMC method
922 Int_t pdgCand = 4122;
923 Int_t pdgDgLctoV0bachelor[2]={2212, 310};
924 Int_t pdgDgV0toDaughters[2]={211, 211};
926 Int_t pdgDgLctopKpi[3]={2212, 321, 211};
928 // loop to search for candidates Lc->p+K+pi
929 Int_t n3Prong = array3Prong->GetEntriesFast();
930 Int_t nCascades= arrayLctopKos->GetEntriesFast();
932 //AliInfo(Form("\n\n\n\n3 prong candidates = %d, ncascades = %d \n\n\n\n\n", n3Prong, nCascades));
933 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
934 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
935 //Filling a control histogram with no cuts
938 // find associated MC particle for Lc -> p+K+pi
939 Int_t mcLabel = d->MatchToMC(4122, mcArray, 3, pdgDgLctopKpi);
940 //Int_t mcLabelTemp = d->MatchToMC(4122, mcArray);
941 //Printf("mcLabel = %d, mcLabelTemp = %d", mcLabel, mcLabelTemp);
944 AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
946 Int_t pdgCode = partLcpKpi->GetPdgCode();
947 AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
948 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
949 fHistoLcpKpiBeforeCuts->Fill(1);
954 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
955 fHistoLcpKpiBeforeCuts->Fill(0);
960 // loop over cascades to search for candidates Lc->p+K0S
962 for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
964 // Lc candidates and K0s from Lc
965 AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
967 AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
971 //Filling a control histogram with no cuts
976 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
977 fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
979 AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
981 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
983 pdgCode = partLc->GetPdgCode();
984 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
985 pdgCode = TMath::Abs(pdgCode);
986 fHistoLcBeforeCuts->Fill(1);
991 fHistoLcBeforeCuts->Fill(0);
996 //if (!lcK0spr->GetSecondaryVtx()) {
997 // AliInfo("No secondary vertex");
1001 if (lcK0spr->GetNDaughters()!=2) {
1002 AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
1006 AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
1007 AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
1008 if (!v0part || !bachPart) {
1009 AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0s));
1014 if (!v0part->GetSecondaryVtx()) {
1015 AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
1019 if (v0part->GetNDaughters()!=2) {
1020 AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
1024 AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
1025 AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
1026 if (!v0Neg || !v0Neg) {
1027 AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
1032 if (v0Pos->Charge() == v0Neg->Charge()) {
1033 AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
1043 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1044 Int_t mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1046 AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1048 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1050 pdgCode = partLc->GetPdgCode();
1051 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1052 pdgCode = TMath::Abs(pdgCode);
1063 AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
1064 AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1066 if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1069 else { // checking if we want to fill the background
1070 if (fKeepingOnlyHIJINGBkg){
1071 // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1072 Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
1073 if (!isCandidateInjected){
1074 AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1075 fHistoBackground->Fill(1);
1078 AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1079 fHistoBackground->Fill(0);
1083 else if (fKeepingOnlyPYTHIABkg){
1084 // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
1085 AliAODTrack *bachelor = (AliAODTrack*)lcK0spr->GetBachelor();
1086 AliAODTrack *v0pos = (AliAODTrack*)lcK0spr->Getv0PositiveTrack();
1087 AliAODTrack *v0neg = (AliAODTrack*)lcK0spr->Getv0NegativeTrack();
1088 if (!bachelor || !v0pos || !v0neg) {
1089 AliDebug(2, "Cannot retrieve one of the tracks while checking origin, continuing");
1093 Int_t labelbachelor = TMath::Abs(bachelor->GetLabel());
1094 Int_t labelv0pos = TMath::Abs(v0pos->GetLabel());
1095 Int_t labelv0neg = TMath::Abs(v0neg->GetLabel());
1096 AliAODMCParticle* MCbachelor = (AliAODMCParticle*)mcArray->At(labelbachelor);
1097 AliAODMCParticle* MCv0pos = (AliAODMCParticle*)mcArray->At(labelv0pos);
1098 AliAODMCParticle* MCv0neg = (AliAODMCParticle*)mcArray->At(labelv0neg);
1099 if (!MCbachelor || !MCv0pos || !MCv0neg) {
1100 AliDebug(2, "Cannot retrieve MC particle for one of the tracks while checking origin, continuing");
1104 Int_t isBachelorFromPythia = fUtils->CheckOrigin(mcArray, MCbachelor, kTRUE);
1105 Int_t isv0posFromPythia = fUtils->CheckOrigin(mcArray, MCv0pos, kTRUE);
1106 Int_t isv0negFromPythia = fUtils->CheckOrigin(mcArray, MCv0neg, kTRUE);
1107 if (isBachelorFromPythia != 0 && isv0posFromPythia != 0 && isv0negFromPythia != 0){
1108 AliDebug(2, "The candidate is from PYTHIA (i.e. all daughters originate from a quark), keeping it to fill background");
1109 fHistoBackground->Fill(2);
1112 AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background");
1113 fHistoBackground->Fill(3);
1122 FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray);
1129 //________________________________________________________________________
1130 void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
1132 Int_t &nSelectedAnal,
1133 AliRDHFCutsLctoV0 *cutsAnal,
1134 TClonesArray *mcArray){
1136 // Fill histos for Lc -> K0S+proton
1140 if (!part->GetOwnPrimaryVtx()) {
1141 //Printf("No primary vertex for Lc found!!");
1142 part->SetOwnPrimaryVtx(fVtx1);
1145 //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1148 Double_t invmassLc = part->InvMassLctoK0sP();
1150 AliAODv0 * v0part = part->Getv0();
1151 Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1152 if (onFlyV0){ // on-the-fly V0
1154 fHistoLcOnTheFly->Fill(2.);
1157 fHistoLcOnTheFly->Fill(0.);
1160 else { // offline V0
1162 fHistoLcOnTheFly->Fill(3.);
1165 fHistoLcOnTheFly->Fill(1.);
1169 Double_t dcaV0 = v0part->GetDCA();
1170 Double_t invmassK0s = v0part->MassK0Short();
1172 if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
1174 fHistoFiducialAcceptance->Fill(3.);
1177 fHistoFiducialAcceptance->Fill(1.);
1182 fHistoFiducialAcceptance->Fill(2.);
1185 fHistoFiducialAcceptance->Fill(0.);
1189 Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1191 if (isInV0window == 0) {
1192 AliDebug(2, "No: The candidate has NOT passed the mass cuts!");
1193 if (isLc) Printf("SIGNAL candidate rejected");
1196 else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1198 Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1200 if (!isInCascadeWindow) {
1201 AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1202 if (isLc) Printf("SIGNAL candidate rejected");
1205 else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1207 Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1208 if (!isCandidateSelectedCuts){
1209 AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1210 if (isLc) Printf("SIGNAL candidate rejected");
1214 AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1217 AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1219 AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1223 //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1224 Double_t probTPCTOF[AliPID::kSPECIES]={-1.};
1226 UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1227 Printf("detUsed (TPCTOF case) = %d", detUsed);
1228 Double_t probProton = -1.;
1229 Double_t probPion = -1.;
1230 Double_t probKaon = -1.;
1231 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1232 Printf("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]);
1233 probProton = probTPCTOF[AliPID::kProton];
1234 probPion = probTPCTOF[AliPID::kPion];
1235 probKaon = probTPCTOF[AliPID::kKaon];
1237 else { // if you don't have both TOF and TPC, try only TPC
1238 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
1239 Printf("We did not find the detector mask for TOF + TPC, let's see only TPC");
1240 detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1241 Printf("detUsed (TPC case) = %d", detUsed);
1242 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1243 probProton = probTPCTOF[AliPID::kProton];
1244 probPion = probTPCTOF[AliPID::kPion];
1245 probKaon = probTPCTOF[AliPID::kKaon];
1246 Printf("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]);
1249 Printf("Only TPC did not work...");
1251 // resetting mask to ask for both TPC+TOF
1252 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1254 Printf("probProton = %f", probProton);
1256 // now we get the TPC and TOF single PID probabilities (only for Proton, or the tree will explode :) )
1257 Double_t probProtonTPC = -1.;
1258 Double_t probProtonTOF = -1.;
1259 Double_t pidTPC[AliPID::kSPECIES]={-1.};
1260 Double_t pidTOF[AliPID::kSPECIES]={-1.};
1261 Int_t respTPC = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTPC, bachelor, AliPID::kSPECIES, pidTPC);
1262 Int_t respTOF = fPIDResponse->ComputePIDProbability(AliPIDResponse::kDetTOF, bachelor, AliPID::kSPECIES, pidTOF);
1263 if (respTPC == AliPIDResponse::kDetPidOk) probProtonTPC = pidTPC[AliPID::kProton];
1264 if (respTOF == AliPIDResponse::kDetPidOk) probProtonTOF = pidTOF[AliPID::kProton];
1266 // checking V0 status (on-the-fly vs offline)
1267 if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1268 AliDebug(2, "On-the-fly discarded");
1272 if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
1276 if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
1278 if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1279 if (isLc) Printf("SIGNAL candidate rejected");
1280 AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1284 AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1287 Int_t pdgCand = 4122;
1288 Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
1289 Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
1290 Int_t isLc2LBarpi=0, isLc2Lpi=0;
1291 Int_t currentLabel = part->GetLabel();
1294 mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1296 if (bachelor->Charge()==-1) isLc2LBarpi=1;
1297 if (bachelor->Charge()==+1) isLc2Lpi=1;
1301 Int_t pdgDg2prong[2] = {211, 211};
1305 labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1306 if (labelK0S>=0) isK0S = 1;
1309 pdgDg2prong[0] = 211;
1310 pdgDg2prong[1] = 2212;
1312 Int_t isLambdaBar = 0;
1313 Int_t lambdaLabel = 0;
1315 lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1316 if (lambdaLabel>=0) {
1317 AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1318 if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1319 else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1323 pdgDg2prong[0] = 11;
1324 pdgDg2prong[1] = 11;
1326 Int_t gammaLabel = 0;
1328 gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1329 if (gammaLabel>=0) {
1330 AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1331 if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1336 if (currentLabel != -1){
1337 AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1338 pdgTemp = tempPart->GetPdgCode();
1340 if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1341 else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1342 else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1343 else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1344 if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1345 else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
1346 else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
1347 else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
1348 //else AliDebug(2, Form("V0 is something else!!"));
1350 Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1351 Double_t invmassLambda = v0part->MassLambda();
1352 Double_t invmassLambdaBar = v0part->MassAntiLambda();
1354 //Double_t nSigmaITSpr=-999.;
1355 Double_t nSigmaTPCpr=-999.;
1356 Double_t nSigmaTOFpr=-999.;
1358 //Double_t nSigmaITSpi=-999.;
1359 Double_t nSigmaTPCpi=-999.;
1360 Double_t nSigmaTOFpi=-999.;
1362 //Double_t nSigmaITSka=-999.;
1363 Double_t nSigmaTPCka=-999.;
1364 Double_t nSigmaTOFka=-999.;
1367 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1368 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1369 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1370 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1371 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1372 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1373 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1374 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1375 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1378 nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
1379 nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
1380 nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
1381 nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
1382 nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
1383 nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
1386 // Fill candidate variable Tree (track selection, V0 invMass selection)
1387 if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
1389 fCandidateVariables[0] = invmassLc;
1390 fCandidateVariables[1] = invmassLc2Lpi;
1391 fCandidateVariables[2] = invmassK0s;
1392 fCandidateVariables[3] = invmassLambda;
1393 fCandidateVariables[4] = invmassLambdaBar;
1394 fCandidateVariables[5] = part->CosV0PointingAngle();
1395 fCandidateVariables[6] = dcaV0;
1396 fCandidateVariables[7] = part->Getd0Prong(0);
1397 fCandidateVariables[8] = part->Getd0Prong(1);
1398 fCandidateVariables[9] = nSigmaTPCpr;
1399 fCandidateVariables[10] = nSigmaTPCpi;
1400 fCandidateVariables[11] = nSigmaTPCka;
1401 fCandidateVariables[12] = nSigmaTOFpr;
1402 fCandidateVariables[13] = nSigmaTOFpi;
1403 fCandidateVariables[14] = nSigmaTOFka;
1404 fCandidateVariables[15] = bachelor->Pt();
1405 AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1406 fCandidateVariables[16] = v0pos->Pt();
1407 AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
1408 fCandidateVariables[17] = v0neg->Pt();
1409 fCandidateVariables[18] = v0part->Getd0Prong(0);
1410 fCandidateVariables[19] = v0part->Getd0Prong(1);
1411 fCandidateVariables[20] = v0part->Pt();
1412 fCandidateVariables[21] = v0part->InvMass2Prongs(0,1,11,11);
1413 fCandidateVariables[22] = part->Pt();
1414 fCandidateVariables[23] = probProton;
1415 fCandidateVariables[24] = part->Eta();
1416 fCandidateVariables[25] = v0pos->Eta();
1417 fCandidateVariables[26] = v0neg->Eta();
1418 fCandidateVariables[27] = probProtonTPC;
1419 fCandidateVariables[28] = probProtonTOF;
1420 fCandidateVariables[29] = bachelor->Eta();
1422 fCandidateVariables[30] = part->P();
1423 fCandidateVariables[31] = bachelor->P();
1424 fCandidateVariables[32] = v0part->P();
1425 fCandidateVariables[33] = v0pos->P();
1426 fCandidateVariables[34] = v0neg->P();
1428 fCandidateVariables[35] = part->Y(4122);
1429 fCandidateVariables[36] = bachelor->Y(2212);
1430 fCandidateVariables[37] = v0part->Y(310);
1431 fCandidateVariables[38] = v0pos->Y(211);
1432 fCandidateVariables[39] = v0neg->Y(211);
1434 fCandidateVariables[40] = v0part->Eta();
1436 fCandidateVariables[41] = part->DecayLength();
1437 fCandidateVariables[42] = part->DecayLengthV0();
1438 fCandidateVariables[43] = part->Ct(4122);
1439 fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
1441 EBachelor bachCode = kBachInvalid;
1442 EK0S k0SCode = kK0SInvalid;
1444 bachCode = CheckBachelor(part, bachelor, mcArray);
1445 k0SCode = CheckK0S(part, v0part, mcArray);
1448 fCandidateVariables[45] = bachCode;
1449 fCandidateVariables[46] = k0SCode;
1451 Double_t V0KF[3] = {-999999, -999999, -999999};
1452 Double_t errV0KF[3] = {-999999, -999999, -999999};
1453 Double_t LcKF[3] = {-999999, -999999, -999999};
1454 Double_t errLcKF[3] = {-999999, -999999, -999999};
1455 Double_t distances[3] = {-999999, -999999, -999999};
1456 Double_t armPolKF[2] = {-999999, -999999};
1458 if (fCallKFVertexing){
1459 Int_t kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
1460 AliDebug(2, Form("Result from KF = %d", kfResult));
1464 for (Int_t i = 0; i< 3; i++){
1465 Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
1469 fCandidateVariables[47] = V0KF[0];
1470 fCandidateVariables[48] = V0KF[1];
1471 fCandidateVariables[49] = V0KF[2];
1473 fCandidateVariables[50] = errV0KF[0];
1474 fCandidateVariables[51] = errV0KF[1];
1475 fCandidateVariables[52] = errV0KF[2];
1477 fCandidateVariables[53] = LcKF[0];
1478 fCandidateVariables[54] = LcKF[1];
1479 fCandidateVariables[55] = LcKF[2];
1481 fCandidateVariables[56] = errLcKF[0];
1482 fCandidateVariables[57] = errLcKF[1];
1483 fCandidateVariables[58] = errLcKF[2];
1485 fCandidateVariables[59] = distances[0];
1486 fCandidateVariables[60] = distances[1];
1487 fCandidateVariables[61] = distances[2];
1488 fCandidateVariables[62] = armPolKF[0];
1489 fCandidateVariables[63] = armPolKF[1];
1490 fCandidateVariables[64] = v0part->AlphaV0();
1491 fCandidateVariables[65] = v0part->PtArmV0();
1493 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)));
1494 fCandidateVariables[66] = v0pos->GetStatus() & AliESDtrack::kITSrefit;
1495 fCandidateVariables[67] = v0neg->GetStatus() & AliESDtrack::kITSrefit;
1496 fCandidateVariables[68] = v0pos->GetTPCClusterInfo(2, 1);
1497 fCandidateVariables[69] = v0neg->GetTPCClusterInfo(2, 1);
1499 fCandidateVariables[70] = v0part->Xv();
1500 fCandidateVariables[71] = v0part->Yv();
1501 fCandidateVariables[72] = v0part->Zv();
1503 fCandidateVariables[73] = fVtx1->GetX();
1504 fCandidateVariables[74] = fVtx1->GetY();
1505 fCandidateVariables[75] = fVtx1->GetZ();
1507 fCandidateVariables[76] = bachelor->GetITSNcls();
1508 fCandidateVariables[77] = bachelor->HasPointOnITSLayer(0) + bachelor->HasPointOnITSLayer(1);
1510 fCandidateVariables[78] = v0pos->GetITSNcls();
1511 fCandidateVariables[79] = v0pos->HasPointOnITSLayer(0) + v0pos->HasPointOnITSLayer(1);
1513 fCandidateVariables[80] = v0neg->GetITSNcls();
1514 fCandidateVariables[81] = v0neg->HasPointOnITSLayer(0) + v0neg->HasPointOnITSLayer(1);
1516 TVector3 mom1(bachelor->Px(), bachelor->Py(), bachelor->Pz());
1517 TVector3 mom2(v0part->Px(), v0part->Py(), v0part->Pz());
1518 TVector3 momTot(part->Px(), part->Py(), part->Pz());
1520 Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
1521 Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
1523 Double_t alphaArmLc = (Ql1 - Ql2)/(Ql1 + Ql2);
1524 Double_t alphaArmLcCharge = ( bachelor->Charge() > 0 ? (Ql1 - Ql2)/(Ql1 + Ql2) : (Ql2 - Ql1)/(Ql1 + Ql2) );
1525 Double_t ptArmLc = mom1.Perp(momTot);
1527 fCandidateVariables[82] = alphaArmLc;
1528 fCandidateVariables[83] = alphaArmLcCharge;
1529 fCandidateVariables[84] = ptArmLc;
1531 Double_t massK0SPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass(); // mass K0S PDG
1532 Double_t massPrPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); // mass Proton PDG
1533 Double_t massLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass(); // mass Lc PDG
1535 Double_t pStar = TMath::Sqrt((massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)*(massLcPDG*massLcPDG-massPrPDG*massPrPDG-massK0SPDG*massK0SPDG)-4.*massPrPDG*massPrPDG*massK0SPDG*massK0SPDG)/(2.*massLcPDG);
1536 Double_t e = part->E(4122);
1537 Double_t beta = part->P()/e;
1538 Double_t gamma = e/massLcPDG;
1540 Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
1542 fCandidateVariables[85] = cts;
1546 AliDebug(2, "Filling Sgn");
1547 fVariablesTreeSgn->Fill();
1548 fHistoCodesSgn->Fill(bachCode, k0SCode);
1551 if (fFillOnlySgn == kFALSE){
1552 AliDebug(2, "Filling Bkg");
1553 fVariablesTreeBkg->Fill();
1554 fHistoCodesBkg->Fill(bachCode, k0SCode);
1559 fVariablesTreeSgn->Fill();
1567 //________________________________________________________________________
1568 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
1569 Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
1570 Double_t* distances, Double_t* armPolKF) {
1573 // method to perform KF vertexing
1574 // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
1577 Int_t codeKFV0 = -1, codeKFLc = -1;
1579 AliKFVertex primVtxCopy;
1580 Int_t nt = 0, ntcheck = 0;
1581 Double_t pos[3] = {0., 0., 0.};
1584 Int_t contr = fVtx1->GetNContributors();
1585 Double_t covmatrix[6] = {0.};
1586 fVtx1->GetCovarianceMatrix(covmatrix);
1587 Double_t chi2 = fVtx1->GetChi2();
1588 AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
1590 // topological constraint
1591 primVtxCopy = AliKFVertex(primaryESDVtxCopy);
1592 nt = primaryESDVtxCopy.GetNContributors();
1595 Int_t pdg[2] = {211, -211};
1596 Int_t pdgLc[2] = {2212, 310};
1598 Int_t pdgDgV0toDaughters[2] = {211, 211};
1600 Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
1602 // the KF vertex for the V0 has to be built with the prongs of the V0!
1603 Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
1604 AliKFParticle V0, positiveV0KF, negativeV0KF;
1605 Int_t labelsv0daugh[2] = {-1, -1};
1606 Int_t idv0daugh[2] = {-1, -1};
1607 AliExternalTrackParam* esdv0Daugh1 = 0x0;
1608 AliExternalTrackParam* esdv0Daugh2 = 0x0;
1609 for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1610 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1612 AliDebug(2, "No V0 daughters available");
1615 Double_t xyz[3], pxpypz[3], cv[21];
1617 aodTrack->GetXYZ(xyz);
1618 aodTrack->PxPyPz(pxpypz);
1619 aodTrack->GetCovarianceXYZPxPyPz(cv);
1620 sign = aodTrack->Charge();
1621 AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
1623 if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1624 else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1625 labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
1626 idv0daugh[ipr] = aodTrack->GetID();
1627 if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
1629 //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
1631 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1632 if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
1633 positiveV0KF = daughterKF;
1636 negativeV0KF = daughterKF;
1640 Double_t xn, xp, dca;
1641 AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
1642 dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
1644 AliExternalTrackParam tr1(*esdv0Daugh1);
1645 AliExternalTrackParam tr2(*esdv0Daugh2);
1646 tr1.PropagateTo(xn, fBField);
1647 tr2.PropagateTo(xp, fBField);
1649 AliKFParticle daughterKF1(tr1, 211);
1650 AliKFParticle daughterKF2(tr2, 211);
1651 V0.AddDaughter(positiveV0KF);
1652 V0.AddDaughter(negativeV0KF);
1653 //V0.AddDaughter(daughterKF1);
1654 //V0.AddDaughter(daughterKF2);
1660 // Checking the quality of the KF V0 vertex
1661 if( V0.GetNDF() < 1 ) {
1662 //Printf("Number of degrees of freedom < 1, continuing");
1665 if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
1666 //Printf("Chi2 per DOF too big, continuing");
1670 if(ftopoConstraint && nt > 0){
1671 for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1672 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1673 //* subtruct daughters from primary vertex
1674 if(!aodTrack->GetUsedForPrimVtxFit()) {
1675 //Printf("Track %d was not used for primary vertex, continuing", i);
1678 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1679 primVtxCopy -= daughterKF;
1684 // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
1686 if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
1687 //Printf("Deviation from vertex too big, continuing");
1692 //* Get V0 invariant mass
1693 Double_t massV0 = 999999, sigmaMassV0 = 999999;
1694 Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
1697 codeKFV0 = 1; // Mass not ok
1698 if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
1700 else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
1702 fHistoMassKFV0->Fill(massV0, sigmaMassV0);
1704 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]);
1705 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]);
1707 Printf("Vertices: KF: x = %f, y = %f, z = %f", V0.GetX(), V0.GetY(), V0.GetZ());
1708 Printf("Vertices: AOD: x = %f, y = %f, z = %f", v0part->Xv(), v0part->Yv(), v0part->Zv());
1710 //Printf("Got MC vtx for V0");
1711 if (fUseMCInfo && TMath::Abs(labelsv0daugh[0] - labelsv0daugh[1]) == 1) {
1712 AliAODMCParticle* tmpdaughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
1713 AliAODMCParticle* tmpdaughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
1714 if (!tmpdaughv01 && labelsv0daugh[0] > 0){
1715 AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
1717 if (!tmpdaughv02 && labelsv0daugh[1] > 0){
1718 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1720 Double_t xPionMC = tmpdaughv01->Xv(); //Production vertex of Pion --> Where K0S decays
1721 Double_t yPionMC = tmpdaughv01->Yv();
1722 Double_t zPionMC = tmpdaughv01->Zv();
1723 //Printf("Got MC vtx for Pion");
1724 Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
1727 Printf("Not a true V0");
1729 //massV0=-1;//return -1;// !!!!
1731 // now use what we just try with the bachelor, to build the Lc
1733 // topological constraint
1734 nt = primVtxCopy.GetNContributors();
1737 Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
1739 Int_t labelsLcdaugh[2] = {-1, -1};
1740 labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
1741 labelsLcdaugh[1] = mcLabelV0;
1743 if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
1744 AliKFParticle daughterKFLc(*bach, pdgLc[0]);
1745 Lc.AddDaughter(daughterKFLc);
1746 TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(310);
1747 Double_t massPDGK0S = particlePDG->Mass();
1748 V0.SetMassConstraint(massPDGK0S);
1750 if( Lc.GetNDF() < 1 ) {
1751 AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
1754 if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1755 AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1759 if(ftopoConstraint && nt > 0){
1760 //* subtruct daughters from primary vertex
1761 if(!bach->GetUsedForPrimVtxFit()) {
1762 AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
1765 primVtxCopy -= daughterKFLc;
1768 /* the V0 was added above, so it is ok to remove it without checking
1769 if(!V0->GetUsedForPrimVtxFit()) {
1770 Printf("Lc: V0 was not used for primary vertex, continuing");
1774 //primVtxCopy -= V0;
1778 // Check Lc Chi^2 deviation from primary vertex
1780 if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
1781 AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
1785 if(ftopoConstraint){
1787 //* Add Lc to primary vertex to improve the primary vertex resolution
1789 Lc.SetProductionVertex(primVtxCopy);
1794 if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) > fCutKFChi2NDF) {
1795 AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1799 if(ftopoConstraint){
1800 V0.SetProductionVertex(Lc);
1803 // After setting the vertex of the V0, getting/filling some info
1805 //* Get V0 decayLength
1806 Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
1807 Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
1809 if (sigmaDecayLengthV0 > 1e19) {
1810 codeKFV0 = 3; // DecayLength not ok
1812 codeKFV0 = 6; // DecayLength and Mass not ok
1813 if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
1815 else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
1818 fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
1820 //* Get V0 life time
1821 Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
1822 Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
1824 if (sigmaLifeTimeV0 > 1e19) {
1825 codeKFV0 = 4; // LifeTime not ok
1826 if (sigmaDecayLengthV0 > 1e19) {
1827 codeKFV0 = 9; // LifeTime and DecayLength not ok
1829 codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
1830 if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
1832 else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
1834 else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
1835 codeKFV0 = 7; // LifeTime and Mass not ok
1836 if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
1838 else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
1841 fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
1843 if (codeKFV0 == -1) codeKFV0 = 0;
1844 fHistoKFV0->Fill(codeKFV0);
1846 AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
1848 fHistoMassV0All->Fill(massV0);
1849 fHistoDecayLengthV0All->Fill(decayLengthV0);
1850 fHistoLifeTimeV0All->Fill(lifeTimeV0);
1852 Double_t qtAlphaV0[2] = {0., 0.};
1853 Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
1854 positiveV0KF.TransportToPoint(vtxV0KF);
1855 negativeV0KF.TransportToPoint(vtxV0KF);
1856 V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
1857 AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
1858 fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
1859 fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
1860 armPolKF[0] = qtAlphaV0[1];
1861 armPolKF[1] = qtAlphaV0[0];
1863 // Checking MC info for V0
1865 AliAODMCParticle *motherV0 = 0x0;
1866 AliAODMCParticle *daughv01 = 0x0;
1867 AliAODMCParticle *daughv02 = 0x0;
1870 daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
1871 daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
1872 if (!daughv01 && labelsv0daugh[0] > 0){
1873 AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
1876 if (!daughv02 && labelsv0daugh[1] > 0){
1877 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1881 if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
1882 AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
1883 motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
1884 if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
1885 if( motherV0->GetNDaughters() == 2 ){
1886 fHistoMassV0True->Fill(massV0);
1887 fHistoDecayLengthV0True->Fill(decayLengthV0);
1888 fHistoLifeTimeV0True->Fill(lifeTimeV0);
1889 fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
1890 if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
1891 fHistoMassV0TrueK0S->Fill(massV0);
1892 fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
1893 fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
1894 fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
1897 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()));
1899 else if (!motherV0){
1900 AliDebug(3, "could not access MC info for mother, continuing");
1903 AliDebug(3, "MC mother is a gluon, continuing");
1907 AliDebug(3, "Background V0!");
1913 AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
1917 //* Get Lc invariant mass
1918 Double_t massLc = 999999, sigmaMassLc= 999999;
1919 Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
1921 AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
1923 codeKFLc = 1; // Mass not ok
1924 if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
1926 else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
1928 fHistoMassKFLc->Fill(massLc, sigmaMassLc);
1930 //* Get Lc Decay length
1931 Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
1932 Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
1934 AliDebug(3, "----> Lc: Could not get decay length, and sigma");
1935 if (sigmaDecayLengthLc > 1e19) {
1936 codeKFLc = 3; // DecayLength not ok
1938 codeKFLc = 6; // DecayLength and Mass not ok
1939 if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
1941 else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
1944 AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
1946 fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
1948 //* Get Lc life time
1949 Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
1950 Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
1952 AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
1953 if (sigmaLifeTimeLc > 1e19) {
1954 codeKFLc = 4; // LifeTime not ok
1955 if (sigmaDecayLengthLc > 1e19) {
1956 codeKFLc = 9; // LifeTime and DecayLength not ok
1958 codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
1959 if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
1961 else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
1963 else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
1964 codeKFLc = 7; // LifeTime and Mass not ok
1965 if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
1967 else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
1971 fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
1973 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));
1975 if (codeKFLc == -1) codeKFLc = 0;
1976 fHistoKFLc->Fill(codeKFLc);
1978 fHistoKF->Fill(codeKFV0, codeKFLc);
1980 // here we fill the histgrams for all the reconstructed KF vertices for the cascade
1981 fHistoMassLcAll->Fill(massLc);
1982 fHistoDecayLengthLcAll->Fill(decayLengthLc);
1983 fHistoLifeTimeLcAll->Fill(lifeTimeLc);
1985 fHistoMassV0fromLcAll->Fill(massV0);
1986 fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
1987 fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
1989 Double_t xV0 = V0.GetX();
1990 Double_t yV0 = V0.GetY();
1991 Double_t zV0 = V0.GetZ();
1993 Double_t xLc = Lc.GetX();
1994 Double_t yLc = Lc.GetY();
1995 Double_t zLc = Lc.GetZ();
1997 Double_t xPrimVtx = primVtxCopy.GetX();
1998 Double_t yPrimVtx = primVtxCopy.GetY();
1999 Double_t zPrimVtx = primVtxCopy.GetZ();
2001 Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2002 (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2003 (zPrimVtx - zLc) * (zPrimVtx - zLc));
2005 Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2006 (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2007 (zPrimVtx - zV0) * (zPrimVtx - zV0));
2009 Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2010 (yLc - yV0)*(yLc - yV0) +
2011 (zLc - zV0)*(zLc - zV0));
2013 //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2015 fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2016 fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2017 fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2019 distances[0] = distanceLcToPrimVtx;
2020 distances[1] = distanceV0ToPrimVtx;
2021 distances[2] = distanceV0ToLc;
2025 AliAODMCParticle *daughv01Lc = 0x0;
2026 AliAODMCParticle *K0S = 0x0;
2027 AliAODMCParticle *daughv02Lc = 0x0;
2029 if (labelsLcdaugh[0] >= 0) {
2030 // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2031 daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2033 AliDebug(3, "Could not access MC info for first daughter of Lc");
2037 AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2038 if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2041 else { // The bachelor is a fake
2045 if (labelsLcdaugh[1] >= 0) {
2046 //Printf("Getting K0S");
2047 K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2049 AliDebug(3, "Could not access MC info for second daughter of Lc");
2053 if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2057 AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2061 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
2062 if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
2063 Int_t iK0 = K0S->GetMother();
2065 Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2067 else { // The K0S has a mother
2068 daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2070 AliDebug(3, "Could not access MC info for second daughter of Lc");
2072 else { // we can access safely the K0S mother in the MC
2073 if( daughv01Lc && (daughv01Lc->GetMother() == daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){ // This is a true cascade! bachelor and V0 come from the same mother
2074 //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
2075 AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
2076 Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
2077 if (motherLc) pdgMum = motherLc->GetPdgCode();
2078 if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
2079 if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
2080 AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
2082 if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
2083 fHistoMassLcTrue->Fill(massLc);
2084 fHistoDecayLengthLcTrue->Fill(decayLengthLc);
2085 fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
2086 fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
2088 fHistoMassV0fromLcTrue->Fill(massV0);
2089 fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2090 fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2092 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...
2093 AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
2095 fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2096 fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2098 fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2099 fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2100 fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2102 fHistoMassLcSgn->Fill(massLc);
2103 fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2105 fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2106 fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2108 fHistoMassV0fromLcSgn->Fill(massV0);
2109 fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2110 fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2113 isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
2116 // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
2117 // First, we compare the vtx of the Lc
2118 Double_t xLcMC = motherLc->Xv();
2119 Double_t yLcMC = motherLc->Yv();
2120 Double_t zLcMC = motherLc->Zv();
2121 //Printf("Got MC vtx for Lc");
2122 Double_t xProtonMC = daughv01Lc->Xv();
2123 Double_t yProtonMC = daughv01Lc->Yv();
2124 Double_t zProtonMC = daughv01Lc->Zv();
2125 //Printf("Got MC vtx for Proton");
2127 Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2128 (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2129 (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2131 // Then, we compare the vtx of the K0s
2132 Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
2133 Double_t yV0MC = motherV0->Yv();
2134 Double_t zV0MC = motherV0->Zv();
2136 //Printf("Got MC vtx for V0");
2137 Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
2138 Double_t yPionMC = daughv01->Yv();
2139 Double_t zPionMC = daughv01->Zv();
2140 //Printf("Got MC vtx for Pion");
2141 Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
2143 Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2144 (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2145 (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2147 // Then, the K0S vertex wrt primary vertex
2149 Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2150 (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2151 (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2153 fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2154 fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2155 fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2157 } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
2158 else if (!motherLc){
2159 AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
2162 AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2164 } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2166 isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2168 } // closing: else { // we can access safely the K0S mother in the MC
2169 } // closing: else { // The K0S has a mother
2170 } // closing isMCLcok
2171 } // closing !isBkgLc
2172 } // closing fUseMCInfo
2174 //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2175 if ( retMV0 == 0 && retMLc == 0){
2177 errV0KF[0] = sigmaMassV0;
2178 V0KF[1] = decayLengthV0;
2179 errV0KF[1] = sigmaDecayLengthV0;
2180 V0KF[2] = lifeTimeV0;
2181 errV0KF[2] = sigmaLifeTimeV0;
2183 errLcKF[0] = sigmaMassLc;
2184 LcKF[1] = decayLengthLc;
2185 errLcKF[1] = sigmaDecayLengthLc;
2186 LcKF[2] = lifeTimeLc;
2187 errLcKF[2] = sigmaLifeTimeLc;
2193 //________________________________________________________________________
2194 AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
2195 AliAODTrack* bachelor,
2196 TClonesArray *mcArray ){
2198 //Printf("In CheckBachelor");
2200 // function to check if the bachelor is a p, if it is a p but not from Lc
2201 // to be filled for background candidates
2203 Int_t label = bachelor->GetLabel();
2208 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
2209 if (!mcpart) return kBachInvalid;
2210 Int_t pdg = mcpart->PdgCode();
2211 if (TMath::Abs(pdg) != 2212) {
2212 AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
2213 return kBachNoProton;
2215 else { // it is a proton
2216 //Int_t labelLc = part->GetLabel();
2217 Int_t labelLc = FindLcLabel(part, mcArray);
2218 //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
2219 Int_t bachelorMotherLabel = mcpart->GetMother();
2220 //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
2221 if (bachelorMotherLabel == -1) {
2222 return kBachPrimary;
2224 AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2225 if (!bachelorMother) return kBachInvalid;
2226 Int_t pdgMother = bachelorMother->PdgCode();
2227 if (TMath::Abs(pdgMother) != 4122) {
2228 AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2229 return kBachNoLambdaMother;
2231 else { // it comes from Lc
2232 if (labelLc != bachelorMotherLabel){
2233 //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));
2234 AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2235 return kBachDifferentLambdaMother;
2237 else { // it comes from the correct Lc
2238 AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2239 return kBachCorrectLambdaMother;
2244 return kBachInvalid;
2248 //________________________________________________________________________
2249 AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
2251 //AliAODTrack* v0part,
2252 TClonesArray *mcArray ){
2254 // function to check if the K0Spart is a p, if it is a p but not from Lc
2255 // to be filled for background candidates
2257 //Printf(" CheckK0S");
2259 //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
2260 //Int_t label = v0part->GetLabel();
2261 Int_t labelFind = FindV0Label(v0part, mcArray);
2262 //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
2263 if (labelFind == -1) {
2267 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
2268 if (!mcpart) return kK0SInvalid;
2269 Int_t pdg = mcpart->PdgCode();
2270 if (TMath::Abs(pdg) != 310) {
2271 AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2272 //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2275 else { // it is a K0S
2276 //Int_t labelLc = part->GetLabel();
2277 Int_t labelLc = FindLcLabel(part, mcArray);
2278 Int_t K0SpartMotherLabel = mcpart->GetMother();
2279 if (K0SpartMotherLabel == -1) {
2280 return kK0SWithoutMother;
2282 AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
2283 if (!K0SpartMother) return kK0SInvalid;
2284 Int_t pdgMotherK0S = K0SpartMother->PdgCode();
2285 if (TMath::Abs(pdgMotherK0S) != 311) {
2286 AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2287 //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2288 return kK0SNotFromK0;
2290 else { // the K0S comes from a K0
2291 Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
2292 if (K0MotherLabel == -1) {
2295 AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
2296 if (!K0Mother) return kK0SInvalid;
2297 Int_t pdgK0Mother = K0Mother->PdgCode();
2298 if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
2299 AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2300 //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2301 return kK0NoLambdaMother;
2303 else { // the K0 comes from Lc
2304 if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
2305 AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2306 //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2307 return kK0DifferentLambdaMother;
2309 else { // it comes from the correct Lc
2310 AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2311 //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2312 return kK0CorrectLambdaMother;
2322 //----------------------------------------------------------------------------
2323 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2326 //Printf(" FindV0Label");
2328 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2330 Int_t labMother[2]={-1, -1};
2331 AliAODMCParticle *part=0;
2332 AliAODMCParticle *mother=0;
2333 Int_t dgLabels = -1;
2335 Int_t ndg = v0part->GetNDaughters();
2337 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2340 for(Int_t i = 0; i < 2; i++) {
2341 AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
2342 dgLabels = trk->GetLabel();
2343 if (dgLabels == -1) {
2344 //printf("daughter with negative label %d\n", dgLabels);
2347 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2349 //printf("no MC particle\n");
2352 labMother[i] = part->GetMother();
2353 if (labMother[i] != -1){
2354 mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2356 //printf("no MC mother particle\n");
2365 if (labMother[0] == labMother[1]) return labMother[0];
2371 //----------------------------------------------------------------------------
2372 Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
2375 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2377 //Printf(" FindLcLabel");
2379 AliAODMCParticle *part=0;
2380 AliAODMCParticle *mother=0;
2381 AliAODMCParticle *grandmother=0;
2382 Int_t dgLabels = -1;
2384 Int_t ndg = cascade->GetNDaughters();
2386 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2389 // daughter 0 --> bachelor, daughter 1 --> V0
2390 AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
2391 dgLabels = trk->GetLabel();
2392 if (dgLabels == -1 ) {
2393 //printf("daughter with negative label %d\n", dgLabels);
2396 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2398 //printf("no MC particle\n");
2401 Int_t labMotherBach = part->GetMother();
2402 if (labMotherBach == -1){
2405 mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2407 //printf("no MC mother particle\n");
2411 AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
2412 dgLabels = FindV0Label(v0, mcArray);
2413 if (dgLabels == -1 ) {
2414 //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
2417 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2419 //printf("no MC particle\n");
2422 Int_t labMotherv0 = part->GetMother();
2423 if (labMotherv0 == -1){
2426 mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2428 //printf("no MC mother particle\n");
2431 Int_t labGrandMotherv0 = mother->GetMother();
2432 if (labGrandMotherv0 == -1){
2435 grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2437 //printf("no MC mother particle\n");
2441 //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2442 if (labGrandMotherv0 == labMotherBach) return labMotherBach;