]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskSELc2V0bachelorTMVA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix on the V0 daughter label for TMVA analysis (C. Zampolli)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSELc2V0bachelorTMVA.cxx
CommitLineData
0abe5247 1/**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/* $Id: AliAnalysisTaskSELc2V0bachelorTMVA.cxx 62231 2013-04-29 09:47:26Z fprino $ */
17
18//
19//
20// Base class for Lc2V0 Analysis to be used with TMVA
21//
22//
23
24//------------------------------------------------------------------------------------------
25//
26// Author: C. Zampolli (from AliAnalysisTaskSELc2V0bachelor class by A.De Caro, P. Pagano)
27//
28//------------------------------------------------------------------------------------------
29
30#include <TSystem.h>
31#include <TParticle.h>
32#include <TParticlePDG.h>
33#include <TH1F.h>
34#include <TH1I.h>
35#include <TH2F.h>
36#include <TTree.h>
37#include "TROOT.h"
38#include <TDatabasePDG.h>
ef87f0d6 39#include "TVector3.h"
40
0abe5247 41#include <AliAnalysisDataSlot.h>
42#include <AliAnalysisDataContainer.h>
43#include "AliStack.h"
44#include "AliMCEvent.h"
45#include "AliAnalysisManager.h"
46#include "AliAODMCHeader.h"
47#include "AliAODHandler.h"
48#include "AliLog.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"
56#include "AliAODv0.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"
ef87f0d6 69#include "AliExternalTrackParam.h"
70#include "AliESDtrack.h"
0abe5247 71
72using std::cout;
73using std::endl;
74
75ClassImp(AliAnalysisTaskSELc2V0bachelorTMVA)
76
77//__________________________________________________________________________
78AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA():
79AliAnalysisTaskSE(),
80 fUseMCInfo(kFALSE),
81 fOutput(0),
82 fCEvents(0),
83 fPIDResponse(0),
84 fPIDCombined(0),
85 fIsK0sAnalysis(kFALSE),
86 fCounter(0),
87 fAnalCuts(0),
88 fListCuts(0),
89 fUseOnTheFlyV0(kFALSE),
90 fIsEventSelected(kFALSE),
91 fVariablesTreeSgn(0),
92 fVariablesTreeBkg(0),
93 fCandidateVariables(),
94 fIspA(kFALSE),
95 fHistoEvents(0),
96 fHistoLc(0),
97 fHistoLcOnTheFly(0),
98 fFillOnlySgn(kFALSE),
ef87f0d6 99 fHistoLcBeforeCuts(0),
100 fHistoFiducialAcceptance(0),
101 fHistoCodesSgn(0),
102 fHistoCodesBkg(0),
0abe5247 103 fHistoLcpKpiBeforeCuts(0),
104 fVtx1(0),
105
106 fHistoDistanceLcToPrimVtx(0),
107 fHistoDistanceV0ToPrimVtx(0),
108 fHistoDistanceV0ToLc(0),
109
110 fHistoDistanceLcToPrimVtxSgn(0),
111 fHistoDistanceV0ToPrimVtxSgn(0),
112 fHistoDistanceV0ToLcSgn(0),
113
114 fHistoVtxLcResidualToPrimVtx(0),
115 fHistoVtxV0ResidualToPrimVtx(0),
116 fHistoVtxV0ResidualToLc(0),
117
118 fHistoMassV0All(0),
119 fHistoDecayLengthV0All(0),
120 fHistoLifeTimeV0All(0),
121
122 fHistoMassV0True(0),
123 fHistoDecayLengthV0True(0),
124 fHistoLifeTimeV0True(0),
125
126 fHistoMassV0TrueFromAOD(0),
127
128 fHistoMassV0TrueK0S(0),
129 fHistoDecayLengthV0TrueK0S(0),
130 fHistoLifeTimeV0TrueK0S(0),
131
132 fHistoMassV0TrueK0SFromAOD(0),
133
134 fHistoMassLcAll(0),
135 fHistoDecayLengthLcAll(0),
136 fHistoLifeTimeLcAll(0),
137
138 fHistoMassLcTrue(0),
139 fHistoDecayLengthLcTrue(0),
140 fHistoLifeTimeLcTrue(0),
141
142 fHistoMassLcTrueFromAOD(0),
143
144 fHistoMassV0fromLcAll(0),
145 fHistoDecayLengthV0fromLcAll(0),
146 fHistoLifeTimeV0fromLcAll(0),
147
148 fHistoMassV0fromLcTrue(0),
149 fHistoDecayLengthV0fromLcTrue(0),
150 fHistoLifeTimeV0fromLcTrue(0),
151
152 fHistoMassLcSgn(0),
153 fHistoMassLcSgnFromAOD(0),
154 fHistoDecayLengthLcSgn(0),
155 fHistoLifeTimeLcSgn(0),
156
157 fHistoMassV0fromLcSgn(0),
158 fHistoDecayLengthV0fromLcSgn(0),
159 fHistoLifeTimeV0fromLcSgn(0),
160
161 fHistoKF(0),
162 fHistoKFV0(0),
163 fHistoKFLc(0),
164
165 fHistoMassKFV0(0),
166 fHistoDecayLengthKFV0(0),
167 fHistoLifeTimeKFV0(0),
168
169 fHistoMassKFLc(0),
170 fHistoDecayLengthKFLc(0),
171 fHistoLifeTimeKFLc(0),
172
173 fHistoArmenterosPodolanskiV0KF(0),
174 fHistoArmenterosPodolanskiV0KFSgn(0),
175 fHistoArmenterosPodolanskiV0AOD(0),
176 fHistoArmenterosPodolanskiV0AODSgn(0),
177
178 fOutputKF(0),
179 fmcLabelLc(-1),
180 ftopoConstraint(kTRUE),
181 fCallKFVertexing(kFALSE),
182 fKeepingOnlyHIJINGBkg(kFALSE),
183 fUtils(0),
ef87f0d6 184 fHistoBackground(0),
185 fCutKFChi2NDF(999999.),
186 fCutKFDeviationFromVtx(999999.),
187 fCutKFDeviationFromVtxV0(0.),
188 fCurrentEvent(-1),
ae48be26 189 fBField(0),
4539d15b 190 fKeepingOnlyPYTHIABkg(kFALSE),
191 fHistoMCLcK0Sp(0x0)
0abe5247 192{
193 //
194 // Default ctor
195 //
196}
197//___________________________________________________________________________
198AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
199 AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
200 AliAnalysisTaskSE(name),
201 fUseMCInfo(kFALSE),
202 fOutput(0),
203 fCEvents(0),
204 fPIDResponse(0),
205 fPIDCombined(0),
206 fIsK0sAnalysis(kFALSE),
207 fCounter(0),
208 fAnalCuts(analCuts),
209 fListCuts(0),
210 fUseOnTheFlyV0(useOnTheFly),
211 fIsEventSelected(kFALSE),
212 fVariablesTreeSgn(0),
213 fVariablesTreeBkg(0),
214 fCandidateVariables(),
215 fIspA(kFALSE),
216 fHistoEvents(0),
217 fHistoLc(0),
218 fHistoLcOnTheFly(0),
219 fFillOnlySgn(kFALSE),
220 fHistoLcBeforeCuts(0),
221 fHistoFiducialAcceptance(0),
222 fHistoCodesSgn(0),
223 fHistoCodesBkg(0),
224 fHistoLcpKpiBeforeCuts(0),
225 fVtx1(0),
226
227 fHistoDistanceLcToPrimVtx(0),
228 fHistoDistanceV0ToPrimVtx(0),
229 fHistoDistanceV0ToLc(0),
230
231 fHistoDistanceLcToPrimVtxSgn(0),
232 fHistoDistanceV0ToPrimVtxSgn(0),
233 fHistoDistanceV0ToLcSgn(0),
234
235 fHistoVtxLcResidualToPrimVtx(0),
236 fHistoVtxV0ResidualToPrimVtx(0),
237 fHistoVtxV0ResidualToLc(0),
238
239 fHistoMassV0All(0),
240 fHistoDecayLengthV0All(0),
241 fHistoLifeTimeV0All(0),
242
243 fHistoMassV0True(0),
244 fHistoDecayLengthV0True(0),
245 fHistoLifeTimeV0True(0),
246
247 fHistoMassV0TrueFromAOD(0),
248
249 fHistoMassV0TrueK0S(0),
250 fHistoDecayLengthV0TrueK0S(0),
251 fHistoLifeTimeV0TrueK0S(0),
252
253 fHistoMassV0TrueK0SFromAOD(0),
254
255 fHistoMassLcAll(0),
256 fHistoDecayLengthLcAll(0),
257 fHistoLifeTimeLcAll(0),
258
259 fHistoMassLcTrue(0),
260 fHistoDecayLengthLcTrue(0),
261 fHistoLifeTimeLcTrue(0),
262
263 fHistoMassLcTrueFromAOD(0),
264
265 fHistoMassV0fromLcAll(0),
266 fHistoDecayLengthV0fromLcAll(0),
267 fHistoLifeTimeV0fromLcAll(0),
268
269 fHistoMassV0fromLcTrue(0),
270 fHistoDecayLengthV0fromLcTrue(0),
271 fHistoLifeTimeV0fromLcTrue(0),
272
273 fHistoMassLcSgn(0),
274 fHistoMassLcSgnFromAOD(0),
275 fHistoDecayLengthLcSgn(0),
276 fHistoLifeTimeLcSgn(0),
277
278 fHistoMassV0fromLcSgn(0),
279 fHistoDecayLengthV0fromLcSgn(0),
280 fHistoLifeTimeV0fromLcSgn(0),
281
282 fHistoKF(0),
283 fHistoKFV0(0),
284 fHistoKFLc(0),
285
286 fHistoMassKFV0(0),
287 fHistoDecayLengthKFV0(0),
288 fHistoLifeTimeKFV0(0),
289
290 fHistoMassKFLc(0),
291 fHistoDecayLengthKFLc(0),
292 fHistoLifeTimeKFLc(0),
293
294 fHistoArmenterosPodolanskiV0KF(0),
295 fHistoArmenterosPodolanskiV0KFSgn(0),
296 fHistoArmenterosPodolanskiV0AOD(0),
297 fHistoArmenterosPodolanskiV0AODSgn(0),
298
299 fOutputKF(0),
300 fmcLabelLc(-1),
301 ftopoConstraint(kTRUE),
302 fCallKFVertexing(kFALSE),
303 fKeepingOnlyHIJINGBkg(kFALSE),
304 fUtils(0),
ef87f0d6 305 fHistoBackground(0),
306 fCutKFChi2NDF(999999.),
307 fCutKFDeviationFromVtx(999999.),
308 fCutKFDeviationFromVtxV0(0.),
309 fCurrentEvent(-1),
ae48be26 310 fBField(0),
4539d15b 311 fKeepingOnlyPYTHIABkg(kFALSE),
312 fHistoMCLcK0Sp(0x0)
0abe5247 313
314{
315 //
316 // Constructor. Initialization of Inputs and Outputs
317 //
318 Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
319
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
326
327}
328
329//___________________________________________________________________________
330AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
331 //
332 // destructor
333 //
334 Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
335
336 if (fOutput) {
337 delete fOutput;
338 fOutput = 0;
339 }
340
341 if (fPIDResponse) {
342 delete fPIDResponse;
343 }
344
345 if (fPIDCombined) {
346 delete fPIDCombined;
347 }
348
349 if (fCounter) {
350 delete fCounter;
351 fCounter = 0;
352 }
353
354 if (fAnalCuts) {
355 delete fAnalCuts;
356 fAnalCuts = 0;
357 }
358
359 if (fListCuts) {
360 delete fListCuts;
361 fListCuts = 0;
362 }
363
364 if(fVariablesTreeSgn){
365 delete fVariablesTreeSgn;
366 fVariablesTreeSgn = 0;
367 }
368
369 if(fVariablesTreeBkg){
370 delete fVariablesTreeBkg;
371 fVariablesTreeBkg = 0;
372 }
373
374 if (fOutputKF) {
375 delete fOutputKF;
376 fOutputKF = 0;
377 }
378
379 if (fUtils) {
380 delete fUtils;
381 fUtils = 0;
382 }
383
384}
385//_________________________________________________
386void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
387 //
388 // Initialization
389 //
390
391 fIsEventSelected=kFALSE;
392
393 if (fDebug > 1) AliInfo("Init");
394
395 fListCuts = new TList();
396 fListCuts->SetOwner();
397 fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
398 PostData(3,fListCuts);
399
ae48be26 400 if (fUseMCInfo && (fKeepingOnlyHIJINGBkg || fKeepingOnlyPYTHIABkg)) fUtils = new AliVertexingHFUtils();
0abe5247 401
402 return;
403}
404
405//________________________________________ terminate ___________________________
406void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
407{
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.
411
412 //AliInfo("Terminate","");
413 AliAnalysisTaskSE::Terminate();
414
415 fOutput = dynamic_cast<TList*> (GetOutputData(1));
416 if (!fOutput) {
417 AliError("fOutput not available");
418 return;
419 }
4539d15b 420
421
422 AliDebug(2, Form("At MC level, %f Lc --> K0S + p were found", fHistoMCLcK0Sp->GetEntries()));
423
0abe5247 424 fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
425 if (!fOutputKF) {
426 AliError("fOutputKF not available");
427 return;
428 }
429
430 return;
431}
432
433//___________________________________________________________________________
434void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
435 // output
436 AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
437
438 //slot #1
439 //OpenFile(1);
440 fOutput = new TList();
441 fOutput->SetOwner();
442 fOutput->SetName("listTrees");
443
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");
ef87f0d6 448 Int_t nVar = 86;
0abe5247 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";
ae48be26 478 fCandidateVariableNames[27]="TPCProtonProb";
479 fCandidateVariableNames[28]="TOFProtonProb";
0abe5247 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";
498
499 fCandidateVariableNames[47]="V0KFmass";
500 fCandidateVariableNames[48]="V0KFdecayLength";
501 fCandidateVariableNames[49]="V0KFlifeTime";
502
503 fCandidateVariableNames[50]="V0KFmassErr";
504 fCandidateVariableNames[51]="V0KFdecayTimeErr";
505 fCandidateVariableNames[52]="V0KFlifeTimeErr";
506
507 fCandidateVariableNames[53]="LcKFmass";
508 fCandidateVariableNames[54]="LcKFdecayLength";
509 fCandidateVariableNames[55]="LcKFlifeTime";
510
511 fCandidateVariableNames[56]="LcKFmassErr";
512 fCandidateVariableNames[57]="LcKFdecayTimeErr";
513 fCandidateVariableNames[58]="LcKFlifeTimeErr";
514
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";
522
ef87f0d6 523 fCandidateVariableNames[66]="ITSrefitV0pos";
524 fCandidateVariableNames[67]="ITSrefitV0neg";
525
526 fCandidateVariableNames[68]="TPCClV0pos";
527 fCandidateVariableNames[69]="TPCClV0neg";
528
529 fCandidateVariableNames[70]="v0Xcoord";
530 fCandidateVariableNames[71]="v0Ycoord";
531 fCandidateVariableNames[72]="v0Zcoord";
532 fCandidateVariableNames[73]="primVtxX";
533 fCandidateVariableNames[74]="primVtxY";
534 fCandidateVariableNames[75]="primVtxZ";
535
536 fCandidateVariableNames[76]="ITSclBach";
537 fCandidateVariableNames[77]="SPDclBach";
538
539 fCandidateVariableNames[78]="ITSclV0pos";
540 fCandidateVariableNames[79]="SPDclV0pos";
541 fCandidateVariableNames[80]="ITSclV0neg";
542 fCandidateVariableNames[81]="SPDclV0neg";
543
544 fCandidateVariableNames[82]="alphaArmLc";
545 fCandidateVariableNames[83]="alphaArmLcCharge";
546 fCandidateVariableNames[84]="ptArmLc";
547
548 fCandidateVariableNames[85]="CosThetaStar";
549
550
0abe5247 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()));
554 }
555
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());
560 }
561
562 fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
563
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());
568 }
569
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());
575 }
576
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());
581 }
582
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);
585
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"};
590
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());
594 }
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());
598 }
599
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());
603 }
604
ae48be26 605 fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 4, -0.5, 3.5);
606 TString labelBkg[4] = {"Injected", "Non-injected", "Non-PYTHIA", "PYTHIA"};
0abe5247 607 for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
608 fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
609 }
610
611 //fOutput->Add(fVariablesTreeSgn);
612 //fOutput->Add(fVariablesTreeBkg);
613
4539d15b 614 const Float_t ptbins[15] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 17., 25., 35.};
615
616 fHistoMCLcK0Sp = new TH1F("fHistoMCLcK0Sp", "fHistoMCLcK0Sp", 14, ptbins);
617
0abe5247 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);
4539d15b 627 fOutput->Add(fHistoMCLcK0Sp);
0abe5247 628
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();
635
636 if (fAnalCuts->GetIsUsePID()){
637 /*
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);
644 */
645 fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
646 }
647
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]);
653
0abe5247 654 fCounter = new AliNormalizationCounter("NormalizationCounter");
655 fCounter->Init();
656 PostData(2, fCounter);
657
658 // Histograms from KF
659
ef87f0d6 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.);
664
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.);
668
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.);
672
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);
676
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);
680
681 fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
682
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);
686
687 fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
688
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);
692
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);
696
697 fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
698
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);
702
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);
706
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);
711
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);
715
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"};
723
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());
729 }
730
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);
734
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);
738
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);
0abe5247 743 }
744
0abe5247 745 fOutputKF = new TList();
746 fOutputKF->SetOwner();
747 fOutputKF->SetName("listHistoKF");
748
ef87f0d6 749 if (fCallKFVertexing){
750 fOutputKF->Add(fHistoDistanceLcToPrimVtx);
751 fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
752 fOutputKF->Add(fHistoDistanceV0ToLc);
753
754 fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
755 fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
756 fOutputKF->Add(fHistoDistanceV0ToLcSgn);
757
758 fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
759 fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
760 fOutputKF->Add(fHistoVtxV0ResidualToLc);
761
762 fOutputKF->Add(fHistoMassV0All);
763 fOutputKF->Add(fHistoDecayLengthV0All);
764 fOutputKF->Add(fHistoLifeTimeV0All);
765
766 fOutputKF->Add(fHistoMassV0True);
767 fOutputKF->Add(fHistoDecayLengthV0True);
768 fOutputKF->Add(fHistoLifeTimeV0True);
769
770 fOutputKF->Add(fHistoMassV0TrueFromAOD);
771
772 fOutputKF->Add(fHistoMassV0TrueK0S);
773 fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
774 fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
775
776 fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
777
778 fOutputKF->Add(fHistoMassLcAll);
779 fOutputKF->Add(fHistoDecayLengthLcAll);
780 fOutputKF->Add(fHistoLifeTimeLcAll);
781
782 fOutputKF->Add(fHistoMassLcTrue);
783 fOutputKF->Add(fHistoDecayLengthLcTrue);
784 fOutputKF->Add(fHistoLifeTimeLcTrue);
785
786 fOutputKF->Add(fHistoMassLcTrueFromAOD);
787
788 fOutputKF->Add(fHistoMassV0fromLcAll);
789 fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
790 fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
791
792 fOutputKF->Add(fHistoMassV0fromLcTrue);
793 fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
794 fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
795
796 fOutputKF->Add(fHistoMassLcSgn);
797 fOutputKF->Add(fHistoMassLcSgnFromAOD);
798 fOutputKF->Add(fHistoDecayLengthLcSgn);
799 fOutputKF->Add(fHistoLifeTimeLcSgn);
800
801 fOutputKF->Add(fHistoMassV0fromLcSgn);
802 fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
803 fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
804
805 fOutputKF->Add(fHistoKF);
806 fOutputKF->Add(fHistoKFV0);
807 fOutputKF->Add(fHistoKFLc);
808
809 fOutputKF->Add(fHistoMassKFV0);
810 fOutputKF->Add(fHistoDecayLengthKFV0);
811 fOutputKF->Add(fHistoLifeTimeKFV0);
812
813 fOutputKF->Add(fHistoMassKFLc);
814 fOutputKF->Add(fHistoDecayLengthKFLc);
815 fOutputKF->Add(fHistoLifeTimeKFLc);
816
817 fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
818 fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
819 fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
820 fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
821 }
0abe5247 822
823 PostData(6, fOutputKF);
824
0abe5247 825 return;
826}
827
828//_________________________________________________
829void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
830{
831 // user exec
832 if (!fInputEvent) {
833 AliError("NO EVENT FOUND!");
834 return;
835 }
836
ef87f0d6 837 fCurrentEvent++;
4539d15b 838 AliDebug(2, Form("Processing event = %d", fCurrentEvent));
0abe5247 839 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
840 TClonesArray *arrayLctopKos=0;
841
842 TClonesArray *array3Prong = 0;
843
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());
852
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");
857
858 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
859 }
860 } else {
861 arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
862
863 array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
864 }
865
866 if ( !fUseMCInfo && fIspA) {
867 fAnalCuts->SetTriggerClass("");
868 fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
869 }
870
871 fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
872
0abe5247 873 // mc analysis
874 TClonesArray *mcArray = 0;
875 AliAODMCHeader *mcHeader=0;
876
877 if (fUseMCInfo) {
878 // MC array need for maching
879 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
880 if (!mcArray) {
881 AliError("Could not find Monte-Carlo in AOD");
882 return;
883 }
884 // load MC header
885 mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
886 if (!mcHeader) {
887 AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
888 return;
889 }
4539d15b 890 //Printf("Filling MC histo");
891 FillMCHisto(mcArray);
892 }
893
894 // AOD primary vertex
895 fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
896 if (!fVtx1) return;
897 if (fVtx1->GetNContributors()<1) return;
898
899 fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
900
901 if ( !fIsEventSelected ) {
902 fHistoEvents->Fill(0);
903 return; // don't take into account not selected events
0abe5247 904 }
4539d15b 905 fHistoEvents->Fill(1);
906
907 // Setting magnetic field for KF vertexing
908 fBField = aodEvent->GetMagneticField();
909 AliKFParticle::SetField(fBField);
0abe5247 910
911 Int_t nSelectedAnal = 0;
912 if (fIsK0sAnalysis) {
913 MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
914 nSelectedAnal, fAnalCuts,
915 array3Prong, mcHeader);
916 }
917 fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
918
919 PostData(1, fOutput);
920 PostData(2, fCounter);
921 PostData(4, fVariablesTreeSgn);
922 PostData(5, fVariablesTreeBkg);
923 PostData(6, fOutputKF);
924
4539d15b 925}
926//-------------------------------------------------------------------------------
927void AliAnalysisTaskSELc2V0bachelorTMVA::FillMCHisto(TClonesArray *mcArray){
928
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));
932 if (!mcPart){
933 AliError("Failed casting particle from MC array!, Skipping particle");
934 AliInfo("Failed casting particle from MC array!, Skipping particle");
935 continue;
936 }
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));
940 continue;
941 }
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!!"));
947 continue;
948 }
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) {
962 bachelorMC = daugh1;
963 v0MC = daugh0;
964 }
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...");
968 continue;
969 }
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));
974 if(!partK0S){
975 AliError("Error while casting particle! returning a NULL array");
976 AliInfo("Error while casting particle! returning a NULL array");
977 continue;
978 }
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");
982 continue;
983 }
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...");
992 continue;
993 }
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");
1001 }
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());
1006 }
1007 else {
1008 AliDebug(2, "not in fiducial acceptance! Skipping");
1009 continue;
1010 }
1011 }
1012 }
1013 }
1014 }
1015 }
1016 }
1017 }
1018 } // closing loop over mcArray
1019
1020 return;
1021
0abe5247 1022}
1023
1024//-------------------------------------------------------------------------------
1025void 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
1031
1032 Int_t pdgCand = 4122;
1033 Int_t pdgDgLctoV0bachelor[2]={2212, 310};
1034 Int_t pdgDgV0toDaughters[2]={211, 211};
1035
1036 Int_t pdgDgLctopKpi[3]={2212, 321, 211};
1037
1038 // loop to search for candidates Lc->p+K+pi
1039 Int_t n3Prong = array3Prong->GetEntriesFast();
1040 Int_t nCascades= arrayLctopKos->GetEntriesFast();
1041
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
1046 if (fUseMCInfo) {
1047
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);
1052 if (mcLabel >= 0) {
1053
1054 AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1055 if(partLcpKpi){
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);
1060
1061 }
1062 }
1063 else {
1064 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
1065 fHistoLcpKpiBeforeCuts->Fill(0);
1066 }
1067 }
1068 }
1069
1070 // loop over cascades to search for candidates Lc->p+K0S
1071
1072 for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
1073
1074 // Lc candidates and K0s from Lc
1075 AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
1076 if (!lcK0spr) {
1077 AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
1078 continue;
1079 }
1080
1081 //Filling a control histogram with no cuts
1082 if (fUseMCInfo) {
1083
1084 Int_t pdgCode=-2;
1085
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) {
4539d15b 1089 AliDebug(2, Form("----> cascade number %d (total cascade number = %d) is a Lc!", iLctopK0s, nCascades));
0abe5247 1090
1091 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
1092 if(partLc){
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);
1097
1098 }
1099 }
1100 else {
1101 fHistoLcBeforeCuts->Fill(0);
1102 pdgCode=-1;
1103 }
1104 }
1105
1106 //if (!lcK0spr->GetSecondaryVtx()) {
1107 // AliInfo("No secondary vertex");
1108 //continue;
1109 //}
1110
1111 if (lcK0spr->GetNDaughters()!=2) {
1112 AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
1113 continue;
1114 }
1115
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));
1120 continue;
1121 }
1122
1123
1124 if (!v0part->GetSecondaryVtx()) {
1125 AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
1126 continue;
1127 }
1128
1129 if (v0part->GetNDaughters()!=2) {
1130 AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
1131 continue;
1132 }
1133
1134 AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
1135 AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
ff22c824 1136 if (!v0Neg || !v0Pos) {
0abe5247 1137 AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
1138 continue;
1139 }
1140
1141
1142 if (v0Pos->Charge() == v0Neg->Charge()) {
1143 AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
1144 continue;
1145 }
1146
1147 Int_t isLc = 0;
1148
1149 if (fUseMCInfo) {
1150
1151 Int_t pdgCode=-2;
1152
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);
1155 if (mcLabel>=0) {
1156 AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1157
1158 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1159 if(partLc){
1160 pdgCode = partLc->GetPdgCode();
1161 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1162 pdgCode = TMath::Abs(pdgCode);
1163 isLc = 1;
1164 fHistoLc->Fill(1);
1165
1166 }
1167 }
1168 else {
1169 fHistoLc->Fill(0);
1170 pdgCode=-1;
1171 }
1172 }
1173 AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
1174 AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1175 if (!isLc) {
1176 if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1177 continue;
1178 }
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!
122f753d 1182 Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
0abe5247 1183 if (!isCandidateInjected){
1184 AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1185 fHistoBackground->Fill(1);
1186 }
1187 else {
1188 AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1189 fHistoBackground->Fill(0);
1190 continue;
1191 }
1192 }
ae48be26 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");
1200 continue;
1201 }
1202 else {
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");
1211 continue;
1212 }
1213 else {
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);
1220 }
1221 else {
1222 AliDebug(2, "The candidate is NOT from PYTHIA, we skip it when filling background");
1223 fHistoBackground->Fill(3);
1224 continue;
1225 }
1226 }
1227 }
1228 }
0abe5247 1229 }
1230 }
1231
1232 FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray);
1233
1234 }
1235
1236 return;
1237
1238}
1239//________________________________________________________________________
1240void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
1241 Int_t isLc,
1242 Int_t &nSelectedAnal,
1243 AliRDHFCutsLctoV0 *cutsAnal,
1244 TClonesArray *mcArray){
1245 //
1246 // Fill histos for Lc -> K0S+proton
1247 //
1248
1249 /*
1250 if (!part->GetOwnPrimaryVtx()) {
1251 //Printf("No primary vertex for Lc found!!");
1252 part->SetOwnPrimaryVtx(fVtx1);
1253 }
1254 else {
1255 //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1256 }
1257 */
1258 Double_t invmassLc = part->InvMassLctoK0sP();
1259
1260 AliAODv0 * v0part = part->Getv0();
1261 Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1262 if (onFlyV0){ // on-the-fly V0
1263 if (isLc) { // Lc
1264 fHistoLcOnTheFly->Fill(2.);
1265 }
1266 else { // not Lc
1267 fHistoLcOnTheFly->Fill(0.);
1268 }
1269 }
1270 else { // offline V0
1271 if (isLc) { // Lc
1272 fHistoLcOnTheFly->Fill(3.);
1273 }
1274 else { // not Lc
1275 fHistoLcOnTheFly->Fill(1.);
1276 }
1277 }
1278
1279 Double_t dcaV0 = v0part->GetDCA();
1280 Double_t invmassK0s = v0part->MassK0Short();
1281
1282 if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
1283 if (isLc) {
1284 fHistoFiducialAcceptance->Fill(3.);
1285 }
1286 else {
1287 fHistoFiducialAcceptance->Fill(1.);
1288 }
1289 }
1290 else {
1291 if (isLc) {
1292 fHistoFiducialAcceptance->Fill(2.);
1293 }
1294 else {
1295 fHistoFiducialAcceptance->Fill(0.);
1296 }
1297 }
1298
1299 Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1300
1301 if (isInV0window == 0) {
1302 AliDebug(2, "No: The candidate has NOT passed the mass cuts!");
1303 if (isLc) Printf("SIGNAL candidate rejected");
1304 return;
1305 }
1306 else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1307
0abe5247 1308 Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1309
1310 if (!isInCascadeWindow) {
1311 AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1312 if (isLc) Printf("SIGNAL candidate rejected");
1313 return;
1314 }
1315 else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1316
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");
1321 return;
1322 }
1323 else {
1324 AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1325 }
1326
1327 AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1328 if (!bachelor) {
1329 AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1330 return;
1331 }
1332
1333 //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
ae48be26 1334 Double_t probTPCTOF[AliPID::kSPECIES]={-1.};
0abe5247 1335
1336 UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
2b0b36b0 1337 AliDebug(2, Form("detUsed (TPCTOF case) = %d", detUsed));
ae48be26 1338 Double_t probProton = -1.;
6be55058 1339 // Double_t probPion = -1.;
1340 // Double_t probKaon = -1.;
0abe5247 1341 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
2b0b36b0 1342 AliDebug(2, Form("We have found the detector mask for TOF + TPC: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
0abe5247 1343 probProton = probTPCTOF[AliPID::kProton];
6be55058 1344 // probPion = probTPCTOF[AliPID::kPion];
1345 // probKaon = probTPCTOF[AliPID::kKaon];
0abe5247 1346 }
ae48be26 1347 else { // if you don't have both TOF and TPC, try only TPC
1348 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
2b0b36b0 1349 AliDebug(2, "We did not find the detector mask for TOF + TPC, let's see only TPC");
ae48be26 1350 detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
2b0b36b0 1351 AliDebug(2,Form(" detUsed (TPC case) = %d", detUsed));
ae48be26 1352 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask()) {
1353 probProton = probTPCTOF[AliPID::kProton];
6be55058 1354 // probPion = probTPCTOF[AliPID::kPion];
1355 // probKaon = probTPCTOF[AliPID::kKaon];
2b0b36b0 1356 AliDebug(2, Form("TPC only worked: probProton will be set to %f", probTPCTOF[AliPID::kProton]));
ae48be26 1357 }
1358 else {
2b0b36b0 1359 AliDebug(2, "Only TPC did not work...");
ae48be26 1360 }
1361 // resetting mask to ask for both TPC+TOF
1362 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
1363 }
2b0b36b0 1364 AliDebug(2, Form("probProton = %f", probProton));
ae48be26 1365
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];
1375
1376 // checking V0 status (on-the-fly vs offline)
0abe5247 1377 if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1378 AliDebug(2, "On-the-fly discarded");
1379 return;
1380 }
1381
1382 if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
1383 nSelectedAnal++;
1384 }
1385
1386 if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
1387
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");
1391 return;
1392 }
1393 else {
1394 AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1395 }
1396
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();
1402 Int_t mcLabel = 0;
1403 if (fUseMCInfo) {
1404 mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1405 if (mcLabel>=0) {
1406 if (bachelor->Charge()==-1) isLc2LBarpi=1;
1407 if (bachelor->Charge()==+1) isLc2Lpi=1;
1408 }
1409 }
1410
1411 Int_t pdgDg2prong[2] = {211, 211};
1412 Int_t labelK0S = 0;
1413 Int_t isK0S = 0;
1414 if (fUseMCInfo) {
1415 labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1416 if (labelK0S>=0) isK0S = 1;
1417 }
1418
1419 pdgDg2prong[0] = 211;
1420 pdgDg2prong[1] = 2212;
1421 Int_t isLambda = 0;
1422 Int_t isLambdaBar = 0;
1423 Int_t lambdaLabel = 0;
1424 if (fUseMCInfo) {
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;
1430 }
1431 }
1432
1433 pdgDg2prong[0] = 11;
1434 pdgDg2prong[1] = 11;
1435 Int_t isGamma = 0;
1436 Int_t gammaLabel = 0;
1437 if (fUseMCInfo) {
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;
1442 }
1443 }
1444
1445 Int_t pdgTemp = -1;
1446 if (currentLabel != -1){
1447 AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1448 pdgTemp = tempPart->GetPdgCode();
1449 }
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!!"));
1459
1460 Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1461 Double_t invmassLambda = v0part->MassLambda();
1462 Double_t invmassLambdaBar = v0part->MassAntiLambda();
1463
1464 //Double_t nSigmaITSpr=-999.;
1465 Double_t nSigmaTPCpr=-999.;
1466 Double_t nSigmaTOFpr=-999.;
1467
1468 //Double_t nSigmaITSpi=-999.;
1469 Double_t nSigmaTPCpi=-999.;
1470 Double_t nSigmaTOFpi=-999.;
1471
1472 //Double_t nSigmaITSka=-999.;
1473 Double_t nSigmaTPCka=-999.;
1474 Double_t nSigmaTOFka=-999.;
1475
1476 /*
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);
1486 */
1487
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));
1494
1495
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) {
1498
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();
ae48be26 1528 fCandidateVariables[27] = probProtonTPC;
1529 fCandidateVariables[28] = probProtonTOF;
0abe5247 1530 fCandidateVariables[29] = bachelor->Eta();
1531
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();
1537
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);
1543
1544 fCandidateVariables[40] = v0part->Eta();
1545
1546 fCandidateVariables[41] = part->DecayLength();
1547 fCandidateVariables[42] = part->DecayLengthV0();
1548 fCandidateVariables[43] = part->Ct(4122);
1549 fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
1550
ef87f0d6 1551 EBachelor bachCode = kBachInvalid;
1552 EK0S k0SCode = kK0SInvalid;
1553 if (fUseMCInfo) {
1554 bachCode = CheckBachelor(part, bachelor, mcArray);
1555 k0SCode = CheckK0S(part, v0part, mcArray);
1556 }
0abe5247 1557
1558 fCandidateVariables[45] = bachCode;
1559 fCandidateVariables[46] = k0SCode;
1560
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};
1567
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));
1571 }
1572
1573 /*
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]);
1576 }
1577 */
1578
1579 fCandidateVariables[47] = V0KF[0];
1580 fCandidateVariables[48] = V0KF[1];
1581 fCandidateVariables[49] = V0KF[2];
1582
1583 fCandidateVariables[50] = errV0KF[0];
1584 fCandidateVariables[51] = errV0KF[1];
1585 fCandidateVariables[52] = errV0KF[2];
1586
1587 fCandidateVariables[53] = LcKF[0];
1588 fCandidateVariables[54] = LcKF[1];
1589 fCandidateVariables[55] = LcKF[2];
1590
1591 fCandidateVariables[56] = errLcKF[0];
1592 fCandidateVariables[57] = errLcKF[1];
1593 fCandidateVariables[58] = errLcKF[2];
1594
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();
1602
ef87f0d6 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);
1608
1609 fCandidateVariables[70] = v0part->Xv();
1610 fCandidateVariables[71] = v0part->Yv();
1611 fCandidateVariables[72] = v0part->Zv();
1612
1613 fCandidateVariables[73] = fVtx1->GetX();
1614 fCandidateVariables[74] = fVtx1->GetY();
1615 fCandidateVariables[75] = fVtx1->GetZ();
1616
1617 fCandidateVariables[76] = bachelor->GetITSNcls();
1618 fCandidateVariables[77] = bachelor->HasPointOnITSLayer(0) + bachelor->HasPointOnITSLayer(1);
1619
1620 fCandidateVariables[78] = v0pos->GetITSNcls();
1621 fCandidateVariables[79] = v0pos->HasPointOnITSLayer(0) + v0pos->HasPointOnITSLayer(1);
1622
1623 fCandidateVariables[80] = v0neg->GetITSNcls();
1624 fCandidateVariables[81] = v0neg->HasPointOnITSLayer(0) + v0neg->HasPointOnITSLayer(1);
1625
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());
1629
1630 Double_t Ql1 = mom1.Dot(momTot)/momTot.Mag();
1631 Double_t Ql2 = mom2.Dot(momTot)/momTot.Mag();
1632
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);
1636
1637 fCandidateVariables[82] = alphaArmLc;
1638 fCandidateVariables[83] = alphaArmLcCharge;
1639 fCandidateVariables[84] = ptArmLc;
1640
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
1644
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;
1649
1650 Double_t cts = (Ql1/gamma-beta*TMath::Sqrt(pStar*pStar+massPrPDG*massPrPDG))/pStar;
1651
1652 fCandidateVariables[85] = cts;
1653
0abe5247 1654 if (fUseMCInfo) {
1655 if (isLc){
1656 AliDebug(2, "Filling Sgn");
1657 fVariablesTreeSgn->Fill();
1658 fHistoCodesSgn->Fill(bachCode, k0SCode);
1659 }
1660 else {
1661 if (fFillOnlySgn == kFALSE){
1662 AliDebug(2, "Filling Bkg");
1663 fVariablesTreeBkg->Fill();
1664 fHistoCodesBkg->Fill(bachCode, k0SCode);
1665 }
1666 }
1667 }
1668 else {
1669 fVariablesTreeSgn->Fill();
1670 }
1671 }
1672
1673 return;
1674
1675}
1676
1677//________________________________________________________________________
1678Int_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) {
1681
1682 //
1683 // method to perform KF vertexing
1684 // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
1685 //
1686
1687 Int_t codeKFV0 = -1, codeKFLc = -1;
1688
1689 AliKFVertex primVtxCopy;
1690 Int_t nt = 0, ntcheck = 0;
1691 Double_t pos[3] = {0., 0., 0.};
1692
1693 fVtx1->GetXYZ(pos);
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");
1699
1700 // topological constraint
1701 primVtxCopy = AliKFVertex(primaryESDVtxCopy);
1702 nt = primaryESDVtxCopy.GetNContributors();
1703 ntcheck = nt;
1704
1705 Int_t pdg[2] = {211, -211};
1706 Int_t pdgLc[2] = {2212, 310};
1707
1708 Int_t pdgDgV0toDaughters[2] = {211, 211};
1709
1710 Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
1711
1712 // the KF vertex for the V0 has to be built with the prongs of the V0!
1713 Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
ef87f0d6 1714 AliKFParticle V0, positiveV0KF, negativeV0KF;
0abe5247 1715 Int_t labelsv0daugh[2] = {-1, -1};
ef87f0d6 1716 Int_t idv0daugh[2] = {-1, -1};
1717 AliExternalTrackParam* esdv0Daugh1 = 0x0;
1718 AliExternalTrackParam* esdv0Daugh2 = 0x0;
0abe5247 1719 for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1720 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1721 if(!aodTrack) {
1722 AliDebug(2, "No V0 daughters available");
1723 return -1;
1724 }
ef87f0d6 1725 Double_t xyz[3], pxpypz[3], cv[21];
1726 Short_t sign;
1727 aodTrack->GetXYZ(xyz);
1728 aodTrack->PxPyPz(pxpypz);
1729 aodTrack->GetCovarianceXYZPxPyPz(cv);
1730 sign = aodTrack->Charge();
1731 AliExternalTrackParam tmp1( xyz, pxpypz, cv, sign);
1732
1733 if (ipr == 0) esdv0Daugh1 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
1734 else esdv0Daugh2 = new AliExternalTrackParam( xyz, pxpypz, cv, sign);
0abe5247 1735 labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
ef87f0d6 1736 idv0daugh[ipr] = aodTrack->GetID();
0abe5247 1737 if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
1738
1739 //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
ef87f0d6 1740
0abe5247 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;
1744 }
1745 else {
1746 negativeV0KF = daughterKF;
1747 }
0abe5247 1748 }
ef87f0d6 1749
6be55058 1750 Double_t xn=0., xp=0.;//, dca;
ef87f0d6 1751 AliDebug(2, Form("bField = %f, esdv0Daugh1 = %p, esdv0Daugh2 = %p", fBField, esdv0Daugh1, esdv0Daugh2));
6be55058 1752 // dca = esdv0Daugh1->GetDCA(esdv0Daugh2, fBField, xn, xp);
ef87f0d6 1753
1754 AliExternalTrackParam tr1(*esdv0Daugh1);
1755 AliExternalTrackParam tr2(*esdv0Daugh2);
1756 tr1.PropagateTo(xn, fBField);
1757 tr2.PropagateTo(xp, fBField);
1758
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);
1765
1766 delete esdv0Daugh1;
1767 delete esdv0Daugh2;
1768 esdv0Daugh1=0;
1769 esdv0Daugh2=0;
0abe5247 1770 // Checking the quality of the KF V0 vertex
1771 if( V0.GetNDF() < 1 ) {
1772 //Printf("Number of degrees of freedom < 1, continuing");
1773 return -1;
1774 }
ef87f0d6 1775 if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > fCutKFChi2NDF ) {
0abe5247 1776 //Printf("Chi2 per DOF too big, continuing");
1777 return -1;
1778 }
1779
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);
1786 continue;
1787 }
1788 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1789 primVtxCopy -= daughterKF;
1790 ntcheck--;
1791 }
1792 }
1793
0abe5247 1794 // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
ef87f0d6 1795 /*
1796 if( V0.GetDeviationFromVertex( primVtxCopy ) < fCutKFDeviationFromVtxV0) {
0abe5247 1797 //Printf("Deviation from vertex too big, continuing");
0abe5247 1798 return -1;
1799 }
ef87f0d6 1800 */
0abe5247 1801
1802 //* Get V0 invariant mass
1803 Double_t massV0 = 999999, sigmaMassV0 = 999999;
1804 Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
1805 if( retMV0 ) {
1806 if (massV0 < 0) {
1807 codeKFV0 = 1; // Mass not ok
1808 if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
1809 }
1810 else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
ef87f0d6 1811 }
0abe5247 1812 fHistoMassKFV0->Fill(massV0, sigmaMassV0);
1813
ef87f0d6 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]);
1816
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());
1819
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");
1826 }
1827 if (!tmpdaughv02 && labelsv0daugh[1] > 0){
1828 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1829 }
ff22c824 1830 if(tmpdaughv01){
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);
1836 }
ef87f0d6 1837 }
1838 else {
1839 Printf("Not a true V0");
1840 }
1841 //massV0=-1;//return -1;// !!!!
1842
1843 // now use what we just try with the bachelor, to build the Lc
1844
1845 // topological constraint
1846 nt = primVtxCopy.GetNContributors();
1847 ntcheck = nt;
1848
1849 Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
1850 AliKFParticle Lc;
1851 Int_t labelsLcdaugh[2] = {-1, -1};
1852 labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
1853 labelsLcdaugh[1] = mcLabelV0;
1854
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);
1861 Lc.AddDaughter(V0);
1862 if( Lc.GetNDF() < 1 ) {
1863 AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
1864 return -1;
1865 }
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()))));
1868 return -1;
1869 }
1870
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");
1875 }
1876 else{
1877 primVtxCopy -= daughterKFLc;
1878 ntcheck--;
1879 }
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");
1883 continue;
1884 }
1885 */
1886 //primVtxCopy -= V0;
1887 //ntcheck--;
1888 }
1889
1890 // Check Lc Chi^2 deviation from primary vertex
1891 /*
1892 if( Lc.GetDeviationFromVertex( primVtxCopy ) > fCutKFDeviationFromVtx) {
1893 AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
1894 return -1;
1895 }
1896
1897 if(ftopoConstraint){
1898 if(ntcheck>0) {
a3c87206 1899 // Add Lc to primary vertex to improve the primary vertex resolution
ef87f0d6 1900 primVtxCopy += Lc;
1901 Lc.SetProductionVertex(primVtxCopy);
1902 }
1903 }
1904 */
1905 //* Check chi^2
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()))));
1908 return -1;
1909 }
1910
1911 if(ftopoConstraint){
1912 V0.SetProductionVertex(Lc);
1913 }
1914
1915 // After setting the vertex of the V0, getting/filling some info
1916
0abe5247 1917 //* Get V0 decayLength
1918 Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
1919 Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
1920 if( retDLV0 ) {
1921 if (sigmaDecayLengthV0 > 1e19) {
1922 codeKFV0 = 3; // DecayLength not ok
1923 if (massV0 < 0) {
1924 codeKFV0 = 6; // DecayLength and Mass not ok
1925 if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
1926 }
1927 else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
1928 }
1929 }
1930 fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
1931
1932 //* Get V0 life time
1933 Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
1934 Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
1935 if( retTLV0 ) {
1936 if (sigmaLifeTimeV0 > 1e19) {
1937 codeKFV0 = 4; // LifeTime not ok
1938 if (sigmaDecayLengthV0 > 1e19) {
1939 codeKFV0 = 9; // LifeTime and DecayLength not ok
1940 if (massV0 < 0) {
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
1943 }
1944 else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
1945 }
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
1949 }
1950 else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
1951 }
1952 }
1953 fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
1954
1955 if (codeKFV0 == -1) codeKFV0 = 0;
1956 fHistoKFV0->Fill(codeKFV0);
1957
1958 AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
1959
1960 fHistoMassV0All->Fill(massV0);
1961 fHistoDecayLengthV0All->Fill(decayLengthV0);
1962 fHistoLifeTimeV0All->Fill(lifeTimeV0);
1963
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];
1974
1975 // Checking MC info for V0
1976
1977 AliAODMCParticle *motherV0 = 0x0;
1978 AliAODMCParticle *daughv01 = 0x0;
1979 AliAODMCParticle *daughv02 = 0x0;
1980
1981 if (fUseMCInfo) {
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");
1986 isMCokV0 = kFALSE;
1987 }
1988 if (!daughv02 && labelsv0daugh[1] > 0){
1989 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1990 isMCokV0 = kFALSE;
1991 }
1992 if (isMCokV0){
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());
2007 }
2008 }
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()));
2010 }
2011 else if (!motherV0){
2012 AliDebug(3, "could not access MC info for mother, continuing");
2013 }
2014 else {
2015 AliDebug(3, "MC mother is a gluon, continuing");
2016 }
2017 }
2018 else {
2019 AliDebug(3, "Background V0!");
2020 isBkgV0 = kTRUE;
2021 }
2022 }
2023 }
2024
2025 AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
2026
ef87f0d6 2027 // Going back to Lc
0abe5247 2028
2029 //* Get Lc invariant mass
2030 Double_t massLc = 999999, sigmaMassLc= 999999;
2031 Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
2032 if( retMLc ) {
2033 AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
2034 if (massLc < 0) {
2035 codeKFLc = 1; // Mass not ok
2036 if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
2037 }
2038 else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
2039 }
2040 fHistoMassKFLc->Fill(massLc, sigmaMassLc);
2041
2042 //* Get Lc Decay length
2043 Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
2044 Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
2045 if( retDLLc ) {
2046 AliDebug(3, "----> Lc: Could not get decay length, and sigma");
2047 if (sigmaDecayLengthLc > 1e19) {
2048 codeKFLc = 3; // DecayLength not ok
2049 if (massLc < 0) {
2050 codeKFLc = 6; // DecayLength and Mass not ok
2051 if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
2052 }
2053 else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
2054 }
2055 }
2056 AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
2057
2058 fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
2059
2060 //* Get Lc life time
2061 Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
2062 Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
2063 if( retTLLc ) {
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
2069 if (massLc < 0) {
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
2072 }
2073 else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
2074 }
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
2078 }
2079 else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
2080 }
2081 }
2082
2083 fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
2084
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));
2086
2087 if (codeKFLc == -1) codeKFLc = 0;
2088 fHistoKFLc->Fill(codeKFLc);
2089
2090 fHistoKF->Fill(codeKFV0, codeKFLc);
2091
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);
2096
2097 fHistoMassV0fromLcAll->Fill(massV0);
2098 fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
2099 fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
2100
2101 Double_t xV0 = V0.GetX();
2102 Double_t yV0 = V0.GetY();
2103 Double_t zV0 = V0.GetZ();
2104
2105 Double_t xLc = Lc.GetX();
2106 Double_t yLc = Lc.GetY();
2107 Double_t zLc = Lc.GetZ();
2108
2109 Double_t xPrimVtx = primVtxCopy.GetX();
2110 Double_t yPrimVtx = primVtxCopy.GetY();
2111 Double_t zPrimVtx = primVtxCopy.GetZ();
2112
2113 Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
2114 (yPrimVtx - yLc) * (yPrimVtx - yLc) +
2115 (zPrimVtx - zLc) * (zPrimVtx - zLc));
2116
2117 Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
2118 (yPrimVtx - yV0) * (yPrimVtx - yV0) +
2119 (zPrimVtx - zV0) * (zPrimVtx - zV0));
2120
2121 Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
2122 (yLc - yV0)*(yLc - yV0) +
2123 (zLc - zV0)*(zLc - zV0));
2124
2125 //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
2126
2127 fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
2128 fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
2129 fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
2130
2131 distances[0] = distanceLcToPrimVtx;
2132 distances[1] = distanceV0ToPrimVtx;
2133 distances[2] = distanceV0ToLc;
2134
2135 if (fUseMCInfo) {
2136
2137 AliAODMCParticle *daughv01Lc = 0x0;
2138 AliAODMCParticle *K0S = 0x0;
2139 AliAODMCParticle *daughv02Lc = 0x0;
2140
2141 if (labelsLcdaugh[0] >= 0) {
2142 // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
2143 daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
2144 if (!daughv01Lc){
2145 AliDebug(3, "Could not access MC info for first daughter of Lc");
2146 isMCokLc = kFALSE;
2147 }
2148 else {
2149 AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
2150 if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
2151 }
2152 }
2153 else { // The bachelor is a fake
2154 isBkgLc = kTRUE;
2155 }
2156
2157 if (labelsLcdaugh[1] >= 0) {
2158 //Printf("Getting K0S");
2159 K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
2160 if (!K0S) {
2161 AliDebug(3, "Could not access MC info for second daughter of Lc");
2162 isMCokLc = kFALSE;
2163 }
2164 else{
2165 if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
2166 }
2167 }
2168 else{
2169 AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
2170 isBkgLc = kTRUE;
2171 }
2172
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();
2176 if (iK0 < 0) {
2177 Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
2178 }
2179 else { // The K0S has a mother
2180 daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
2181 if (!daughv02Lc){
2182 AliDebug(3, "Could not access MC info for second daughter of Lc");
2183 }
2184 else { // we can access safely the K0S mother in the MC
ae48be26 2185 if( daughv01Lc && (daughv01Lc->GetMother() == daughv02Lc->GetMother()) && (daughv01Lc->GetMother()>=0) ){ // This is a true cascade! bachelor and V0 come from the same mother
0abe5247 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));
2193
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());
2199
2200 fHistoMassV0fromLcTrue->Fill(massV0);
2201 fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
2202 fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
2203
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()));
2206
2207 fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
2208 fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
2209
2210 fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
2211 fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
2212 fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
2213
2214 fHistoMassLcSgn->Fill(massLc);
2215 fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
2216
2217 fHistoDecayLengthLcSgn->Fill(decayLengthLc);
2218 fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
2219
2220 fHistoMassV0fromLcSgn->Fill(massV0);
2221 fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
2222 fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
2223 }
2224 else {
2225 isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
2226 }
2227
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");
2238
2239 Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
2240 (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
2241 (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
2242
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();
ef87f0d6 2247
0abe5247 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");
ef87f0d6 2253 Printf("Vertices: MC: x = %f, y = %f, z = %f", xPionMC, yPionMC, zPionMC);
0abe5247 2254
2255 Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
2256 (yV0MC - yPionMC) * (yV0MC - yPionMC) +
2257 (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
2258
2259 // Then, the K0S vertex wrt primary vertex
2260
2261 Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
2262 (yPionMC - yLcMC) * (yPionMC - yLcMC) +
2263 (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
2264
2265 fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
2266 fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
2267 fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
2268
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");
2272 }
2273 else {
2274 AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
2275 }
2276 } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
2277 else {
2278 isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
2279 }
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
2285
2286 //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
2287 if ( retMV0 == 0 && retMLc == 0){
2288 V0KF[0] = massV0;
2289 errV0KF[0] = sigmaMassV0;
2290 V0KF[1] = decayLengthV0;
2291 errV0KF[1] = sigmaDecayLengthV0;
2292 V0KF[2] = lifeTimeV0;
2293 errV0KF[2] = sigmaLifeTimeV0;
2294 LcKF[0] = massLc;
2295 errLcKF[0] = sigmaMassLc;
2296 LcKF[1] = decayLengthLc;
2297 errLcKF[1] = sigmaDecayLengthLc;
2298 LcKF[2] = lifeTimeLc;
2299 errLcKF[2] = sigmaLifeTimeLc;
2300 }
2301
2302 return 1;
2303
2304}
2305//________________________________________________________________________
2306AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
2307 AliAODTrack* bachelor,
2308 TClonesArray *mcArray ){
2309
2310 //Printf("In CheckBachelor");
2311
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
2314
2315 Int_t label = bachelor->GetLabel();
2316 if (label == -1) {
2317 return kBachFake;
2318 }
2319
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;
2326 }
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;
2335 }
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;
2342 }
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;
2348 }
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;
2352 }
2353 }
2354 }
2355
2356 return kBachInvalid;
2357
2358}
2359
2360//________________________________________________________________________
2361AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
2362 AliAODv0* v0part,
2363 //AliAODTrack* v0part,
2364 TClonesArray *mcArray ){
2365
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
2368
2369 //Printf(" CheckK0S");
2370
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) {
2376 return kK0SFake;
2377 }
2378
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));
2385 return kK0SNoK0S;
2386 }
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;
2393 }
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;
2401 }
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) {
2405 return kK0Primary;
2406 }
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;
2414 }
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;
2420 }
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;
2425 }
2426 }
2427 }
2428 }
2429
2430 return kK0SInvalid;
2431
2432}
2433
2434//----------------------------------------------------------------------------
2435Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2436{
2437
2438 //Printf(" FindV0Label");
2439
2440 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2441
2442 Int_t labMother[2]={-1, -1};
2443 AliAODMCParticle *part=0;
2444 AliAODMCParticle *mother=0;
2445 Int_t dgLabels = -1;
2446
2447 Int_t ndg = v0part->GetNDaughters();
2448 if (ndg != 2) {
2449 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2450 }
2451
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);
2457 return -1;
2458 }
2459 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2460 if (!part) {
2461 //printf("no MC particle\n");
2462 return -1;
2463 }
2464 labMother[i] = part->GetMother();
2465 if (labMother[i] != -1){
2466 mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2467 if(!mother) {
2468 //printf("no MC mother particle\n");
2469 return -1;
2470 }
2471 }
2472 else {
2473 return -1;
2474 }
2475 }
2476
2477 if (labMother[0] == labMother[1]) return labMother[0];
2478
2479 else return -1;
2480
2481}
2482
2483//----------------------------------------------------------------------------
2484Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
2485{
2486
2487 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2488
2489 //Printf(" FindLcLabel");
2490
2491 AliAODMCParticle *part=0;
2492 AliAODMCParticle *mother=0;
2493 AliAODMCParticle *grandmother=0;
2494 Int_t dgLabels = -1;
2495
2496 Int_t ndg = cascade->GetNDaughters();
2497 if (ndg != 2) {
2498 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2499 }
2500
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);
2506 return -1;
2507 }
2508 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2509 if (!part) {
2510 //printf("no MC particle\n");
2511 return -1;
2512 }
2513 Int_t labMotherBach = part->GetMother();
2514 if (labMotherBach == -1){
2515 return -1;
2516 }
2517 mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2518 if(!mother) {
2519 //printf("no MC mother particle\n");
2520 return -1;
2521 }
2522
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);
2527 return -1;
2528 }
2529 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2530 if (!part) {
2531 //printf("no MC particle\n");
2532 return -1;
2533 }
2534 Int_t labMotherv0 = part->GetMother();
2535 if (labMotherv0 == -1){
2536 return -1;
2537 }
2538 mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2539 if(!mother) {
2540 //printf("no MC mother particle\n");
2541 return -1;
2542 }
2543 Int_t labGrandMotherv0 = mother->GetMother();
2544 if (labGrandMotherv0 == -1){
2545 return -1;
2546 }
2547 grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2548 if(!grandmother) {
2549 //printf("no MC mother particle\n");
2550 return -1;
2551 }
2552
2553 //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2554 if (labGrandMotherv0 == labMotherBach) return labMotherBach;
2555
2556 else return -1;
2557
2558}
0abe5247 2559
2560
2561
2562
2563
2564