Update in Lc task (ChiaraZ)
[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>
39#include <AliAnalysisDataSlot.h>
40#include <AliAnalysisDataContainer.h>
41#include "AliStack.h"
42#include "AliMCEvent.h"
43#include "AliAnalysisManager.h"
44#include "AliAODMCHeader.h"
45#include "AliAODHandler.h"
46#include "AliLog.h"
47#include "AliAODVertex.h"
48#include "AliAODRecoDecay.h"
49#include "AliAODRecoDecayHF.h"
50#include "AliAODRecoCascadeHF.h"
51#include "AliAnalysisVertexingHF.h"
52#include "AliESDtrack.h"
53#include "AliAODTrack.h"
54#include "AliAODv0.h"
55#include "AliAODMCParticle.h"
56#include "AliAnalysisTaskSE.h"
57#include "AliAnalysisTaskSELc2V0bachelorTMVA.h"
58#include "AliNormalizationCounter.h"
59#include "AliAODPidHF.h"
60#include "AliPIDResponse.h"
61#include "AliPIDCombined.h"
62#include "AliTOFPIDResponse.h"
63#include "AliInputEventHandler.h"
64#include "AliAODRecoDecayHF3Prong.h"
65#include "AliKFParticle.h"
66#include "AliKFVertex.h"
67#include "AliVertexingHFUtils.h"
68
69using std::cout;
70using std::endl;
71
72ClassImp(AliAnalysisTaskSELc2V0bachelorTMVA)
73
74//__________________________________________________________________________
75AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA():
76AliAnalysisTaskSE(),
77 fUseMCInfo(kFALSE),
78 fOutput(0),
79 fCEvents(0),
80 fPIDResponse(0),
81 fPIDCombined(0),
82 fIsK0sAnalysis(kFALSE),
83 fCounter(0),
84 fAnalCuts(0),
85 fListCuts(0),
86 fUseOnTheFlyV0(kFALSE),
87 fIsEventSelected(kFALSE),
88 fVariablesTreeSgn(0),
89 fVariablesTreeBkg(0),
90 fCandidateVariables(),
91 fIspA(kFALSE),
92 fHistoEvents(0),
93 fHistoLc(0),
94 fHistoLcOnTheFly(0),
95 fFillOnlySgn(kFALSE),
96fHistoLcBeforeCuts(0),
97fHistoFiducialAcceptance(0),
98fHistoCodesSgn(0),
99fHistoCodesBkg(0),
100 fHistoLcpKpiBeforeCuts(0),
101 fVtx1(0),
102
103 fHistoDistanceLcToPrimVtx(0),
104 fHistoDistanceV0ToPrimVtx(0),
105 fHistoDistanceV0ToLc(0),
106
107 fHistoDistanceLcToPrimVtxSgn(0),
108 fHistoDistanceV0ToPrimVtxSgn(0),
109 fHistoDistanceV0ToLcSgn(0),
110
111 fHistoVtxLcResidualToPrimVtx(0),
112 fHistoVtxV0ResidualToPrimVtx(0),
113 fHistoVtxV0ResidualToLc(0),
114
115 fHistoMassV0All(0),
116 fHistoDecayLengthV0All(0),
117 fHistoLifeTimeV0All(0),
118
119 fHistoMassV0True(0),
120 fHistoDecayLengthV0True(0),
121 fHistoLifeTimeV0True(0),
122
123 fHistoMassV0TrueFromAOD(0),
124
125 fHistoMassV0TrueK0S(0),
126 fHistoDecayLengthV0TrueK0S(0),
127 fHistoLifeTimeV0TrueK0S(0),
128
129 fHistoMassV0TrueK0SFromAOD(0),
130
131 fHistoMassLcAll(0),
132 fHistoDecayLengthLcAll(0),
133 fHistoLifeTimeLcAll(0),
134
135 fHistoMassLcTrue(0),
136 fHistoDecayLengthLcTrue(0),
137 fHistoLifeTimeLcTrue(0),
138
139 fHistoMassLcTrueFromAOD(0),
140
141 fHistoMassV0fromLcAll(0),
142 fHistoDecayLengthV0fromLcAll(0),
143 fHistoLifeTimeV0fromLcAll(0),
144
145 fHistoMassV0fromLcTrue(0),
146 fHistoDecayLengthV0fromLcTrue(0),
147 fHistoLifeTimeV0fromLcTrue(0),
148
149 fHistoMassLcSgn(0),
150 fHistoMassLcSgnFromAOD(0),
151 fHistoDecayLengthLcSgn(0),
152 fHistoLifeTimeLcSgn(0),
153
154 fHistoMassV0fromLcSgn(0),
155 fHistoDecayLengthV0fromLcSgn(0),
156 fHistoLifeTimeV0fromLcSgn(0),
157
158 fHistoKF(0),
159 fHistoKFV0(0),
160 fHistoKFLc(0),
161
162 fHistoMassKFV0(0),
163 fHistoDecayLengthKFV0(0),
164 fHistoLifeTimeKFV0(0),
165
166 fHistoMassKFLc(0),
167 fHistoDecayLengthKFLc(0),
168 fHistoLifeTimeKFLc(0),
169
170 fHistoArmenterosPodolanskiV0KF(0),
171 fHistoArmenterosPodolanskiV0KFSgn(0),
172 fHistoArmenterosPodolanskiV0AOD(0),
173 fHistoArmenterosPodolanskiV0AODSgn(0),
174
175 fOutputKF(0),
176 fmcLabelLc(-1),
177 ftopoConstraint(kTRUE),
178 fCallKFVertexing(kFALSE),
179 fKeepingOnlyHIJINGBkg(kFALSE),
180 fUtils(0),
181 fHistoBackground(0)
182
183{
184 //
185 // Default ctor
186 //
187}
188//___________________________________________________________________________
189AliAnalysisTaskSELc2V0bachelorTMVA::AliAnalysisTaskSELc2V0bachelorTMVA(const Char_t* name,
190 AliRDHFCutsLctoV0* analCuts, Bool_t useOnTheFly) :
191 AliAnalysisTaskSE(name),
192 fUseMCInfo(kFALSE),
193 fOutput(0),
194 fCEvents(0),
195 fPIDResponse(0),
196 fPIDCombined(0),
197 fIsK0sAnalysis(kFALSE),
198 fCounter(0),
199 fAnalCuts(analCuts),
200 fListCuts(0),
201 fUseOnTheFlyV0(useOnTheFly),
202 fIsEventSelected(kFALSE),
203 fVariablesTreeSgn(0),
204 fVariablesTreeBkg(0),
205 fCandidateVariables(),
206 fIspA(kFALSE),
207 fHistoEvents(0),
208 fHistoLc(0),
209 fHistoLcOnTheFly(0),
210 fFillOnlySgn(kFALSE),
211 fHistoLcBeforeCuts(0),
212 fHistoFiducialAcceptance(0),
213 fHistoCodesSgn(0),
214 fHistoCodesBkg(0),
215 fHistoLcpKpiBeforeCuts(0),
216 fVtx1(0),
217
218 fHistoDistanceLcToPrimVtx(0),
219 fHistoDistanceV0ToPrimVtx(0),
220 fHistoDistanceV0ToLc(0),
221
222 fHistoDistanceLcToPrimVtxSgn(0),
223 fHistoDistanceV0ToPrimVtxSgn(0),
224 fHistoDistanceV0ToLcSgn(0),
225
226 fHistoVtxLcResidualToPrimVtx(0),
227 fHistoVtxV0ResidualToPrimVtx(0),
228 fHistoVtxV0ResidualToLc(0),
229
230 fHistoMassV0All(0),
231 fHistoDecayLengthV0All(0),
232 fHistoLifeTimeV0All(0),
233
234 fHistoMassV0True(0),
235 fHistoDecayLengthV0True(0),
236 fHistoLifeTimeV0True(0),
237
238 fHistoMassV0TrueFromAOD(0),
239
240 fHistoMassV0TrueK0S(0),
241 fHistoDecayLengthV0TrueK0S(0),
242 fHistoLifeTimeV0TrueK0S(0),
243
244 fHistoMassV0TrueK0SFromAOD(0),
245
246 fHistoMassLcAll(0),
247 fHistoDecayLengthLcAll(0),
248 fHistoLifeTimeLcAll(0),
249
250 fHistoMassLcTrue(0),
251 fHistoDecayLengthLcTrue(0),
252 fHistoLifeTimeLcTrue(0),
253
254 fHistoMassLcTrueFromAOD(0),
255
256 fHistoMassV0fromLcAll(0),
257 fHistoDecayLengthV0fromLcAll(0),
258 fHistoLifeTimeV0fromLcAll(0),
259
260 fHistoMassV0fromLcTrue(0),
261 fHistoDecayLengthV0fromLcTrue(0),
262 fHistoLifeTimeV0fromLcTrue(0),
263
264 fHistoMassLcSgn(0),
265 fHistoMassLcSgnFromAOD(0),
266 fHistoDecayLengthLcSgn(0),
267 fHistoLifeTimeLcSgn(0),
268
269 fHistoMassV0fromLcSgn(0),
270 fHistoDecayLengthV0fromLcSgn(0),
271 fHistoLifeTimeV0fromLcSgn(0),
272
273 fHistoKF(0),
274 fHistoKFV0(0),
275 fHistoKFLc(0),
276
277 fHistoMassKFV0(0),
278 fHistoDecayLengthKFV0(0),
279 fHistoLifeTimeKFV0(0),
280
281 fHistoMassKFLc(0),
282 fHistoDecayLengthKFLc(0),
283 fHistoLifeTimeKFLc(0),
284
285 fHistoArmenterosPodolanskiV0KF(0),
286 fHistoArmenterosPodolanskiV0KFSgn(0),
287 fHistoArmenterosPodolanskiV0AOD(0),
288 fHistoArmenterosPodolanskiV0AODSgn(0),
289
290 fOutputKF(0),
291 fmcLabelLc(-1),
292 ftopoConstraint(kTRUE),
293 fCallKFVertexing(kFALSE),
294 fKeepingOnlyHIJINGBkg(kFALSE),
295 fUtils(0),
296 fHistoBackground(0)
297
298
299{
300 //
301 // Constructor. Initialization of Inputs and Outputs
302 //
303 Info("AliAnalysisTaskSELc2V0bachelorTMVA","Calling Constructor");
304
305 DefineOutput(1, TList::Class()); // Tree signal + Tree Bkg + histoEvents
306 DefineOutput(2, AliNormalizationCounter::Class()); // normalization counter object
307 DefineOutput(3, TList::Class()); // Cuts
308 DefineOutput(4, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
309 DefineOutput(5, TTree::Class()); // Tree signal + Tree Bkg + histoEvents
310 DefineOutput(6, TList::Class()); // Tree signal + Tree Bkg + histoEvents
311
312}
313
314//___________________________________________________________________________
315AliAnalysisTaskSELc2V0bachelorTMVA::~AliAnalysisTaskSELc2V0bachelorTMVA() {
316 //
317 // destructor
318 //
319 Info("~AliAnalysisTaskSELc2V0bachelorTMVA","Calling Destructor");
320
321 if (fOutput) {
322 delete fOutput;
323 fOutput = 0;
324 }
325
326 if (fPIDResponse) {
327 delete fPIDResponse;
328 }
329
330 if (fPIDCombined) {
331 delete fPIDCombined;
332 }
333
334 if (fCounter) {
335 delete fCounter;
336 fCounter = 0;
337 }
338
339 if (fAnalCuts) {
340 delete fAnalCuts;
341 fAnalCuts = 0;
342 }
343
344 if (fListCuts) {
345 delete fListCuts;
346 fListCuts = 0;
347 }
348
349 if(fVariablesTreeSgn){
350 delete fVariablesTreeSgn;
351 fVariablesTreeSgn = 0;
352 }
353
354 if(fVariablesTreeBkg){
355 delete fVariablesTreeBkg;
356 fVariablesTreeBkg = 0;
357 }
358
359 if (fOutputKF) {
360 delete fOutputKF;
361 fOutputKF = 0;
362 }
363
364 if (fUtils) {
365 delete fUtils;
366 fUtils = 0;
367 }
368
369}
370//_________________________________________________
371void AliAnalysisTaskSELc2V0bachelorTMVA::Init() {
372 //
373 // Initialization
374 //
375
376 fIsEventSelected=kFALSE;
377
378 if (fDebug > 1) AliInfo("Init");
379
380 fListCuts = new TList();
381 fListCuts->SetOwner();
382 fListCuts->Add(new AliRDHFCutsLctoV0(*fAnalCuts));
383 PostData(3,fListCuts);
384
385 if (fUseMCInfo && fKeepingOnlyHIJINGBkg) fUtils = new AliVertexingHFUtils();
386
387 return;
388}
389
390//________________________________________ terminate ___________________________
391void AliAnalysisTaskSELc2V0bachelorTMVA::Terminate(Option_t*)
392{
393 // The Terminate() function is the last function to be called during
394 // a query. It always runs on the client, it can be used to present
395 // the results graphically or save the results to file.
396
397 //AliInfo("Terminate","");
398 AliAnalysisTaskSE::Terminate();
399
400 fOutput = dynamic_cast<TList*> (GetOutputData(1));
401 if (!fOutput) {
402 AliError("fOutput not available");
403 return;
404 }
405 fOutputKF = dynamic_cast<TList*> (GetOutputData(6));
406 if (!fOutputKF) {
407 AliError("fOutputKF not available");
408 return;
409 }
410
411 return;
412}
413
414//___________________________________________________________________________
415void AliAnalysisTaskSELc2V0bachelorTMVA::UserCreateOutputObjects() {
416 // output
417 AliInfo(Form("CreateOutputObjects of task %s\n", GetName()));
418
419 //slot #1
420 //OpenFile(1);
421 fOutput = new TList();
422 fOutput->SetOwner();
423 fOutput->SetName("listTrees");
424
425 // Output slot 1: list of 2 trees (Sgn + Bkg) of the candidate variables
426 const char* nameoutput = GetOutputSlot(1)->GetContainer()->GetName();
427 fVariablesTreeSgn = new TTree(Form("%s_Sgn", nameoutput), "Candidates variables tree, Signal");
428 fVariablesTreeBkg = new TTree(Form("%s_Bkg", nameoutput), "Candidates variables tree, Background");
429 Int_t nVar = 66;
430 fCandidateVariables = new Float_t [nVar];
431 TString * fCandidateVariableNames = new TString[nVar];
432 fCandidateVariableNames[0]="massLc2K0Sp";
433 fCandidateVariableNames[1]="massLc2Lambdapi";
434 fCandidateVariableNames[2]="massK0S";
435 fCandidateVariableNames[3]="massLambda";
436 fCandidateVariableNames[4]="massLambdaBar";
437 fCandidateVariableNames[5]="cosPAK0S";
438 fCandidateVariableNames[6]="dcaV0";
439 fCandidateVariableNames[7]="tImpParBach";
440 fCandidateVariableNames[8]="tImpParV0";
441 fCandidateVariableNames[9]="nSigmaTPCpr";
442 fCandidateVariableNames[10]="nSigmaTPCpi";
443 fCandidateVariableNames[11]="nSigmaTPCka";
444 fCandidateVariableNames[12]="nSigmaTOFpr";
445 fCandidateVariableNames[13]="nSigmaTOFpi";
446 fCandidateVariableNames[14]="nSigmaTOFka";
447 fCandidateVariableNames[15]="bachelorPt";
448 fCandidateVariableNames[16]="V0positivePt";
449 fCandidateVariableNames[17]="V0negativePt";
450 fCandidateVariableNames[18]="dcaV0pos";
451 fCandidateVariableNames[19]="dcaV0neg";
452 fCandidateVariableNames[20]="v0Pt";
453 fCandidateVariableNames[21]="massGamma";
454 fCandidateVariableNames[22]="LcPt";
455 fCandidateVariableNames[23]="combinedProtonProb";
456 fCandidateVariableNames[24]="LcEta";
457 fCandidateVariableNames[25]="V0positiveEta";
458 fCandidateVariableNames[26]="V0negativeEta";
459 fCandidateVariableNames[27]="combinedPionProb";
460 fCandidateVariableNames[28]="combinedKaonProb";
461 fCandidateVariableNames[29]="bachelorEta";
462 fCandidateVariableNames[30]="LcP";
463 fCandidateVariableNames[31]="bachelorP";
464 fCandidateVariableNames[32]="v0P";
465 fCandidateVariableNames[33]="V0positiveP";
466 fCandidateVariableNames[34]="V0negativeP";
467 fCandidateVariableNames[35]="LcY";
468 fCandidateVariableNames[36]="v0Y";
469 fCandidateVariableNames[37]="bachelorY";
470 fCandidateVariableNames[38]="V0positiveY";
471 fCandidateVariableNames[39]="V0negativeY";
472 fCandidateVariableNames[40]="v0Eta";
473 fCandidateVariableNames[41]="DecayLengthLc";
474 fCandidateVariableNames[42]="DecayLengthK0S";
475 fCandidateVariableNames[43]="CtLc";
476 fCandidateVariableNames[44]="CtK0S";
477 fCandidateVariableNames[45]="bachCode";
478 fCandidateVariableNames[46]="k0SCode";
479
480 fCandidateVariableNames[47]="V0KFmass";
481 fCandidateVariableNames[48]="V0KFdecayLength";
482 fCandidateVariableNames[49]="V0KFlifeTime";
483
484 fCandidateVariableNames[50]="V0KFmassErr";
485 fCandidateVariableNames[51]="V0KFdecayTimeErr";
486 fCandidateVariableNames[52]="V0KFlifeTimeErr";
487
488 fCandidateVariableNames[53]="LcKFmass";
489 fCandidateVariableNames[54]="LcKFdecayLength";
490 fCandidateVariableNames[55]="LcKFlifeTime";
491
492 fCandidateVariableNames[56]="LcKFmassErr";
493 fCandidateVariableNames[57]="LcKFdecayTimeErr";
494 fCandidateVariableNames[58]="LcKFlifeTimeErr";
495
496 fCandidateVariableNames[59]="LcKFDistToPrimVtx";
497 fCandidateVariableNames[60]="V0KFDistToPrimVtx";
498 fCandidateVariableNames[61]="V0KFDistToLc";
499 fCandidateVariableNames[62]="alphaArmKF";
500 fCandidateVariableNames[63]="ptArmKF";
501 fCandidateVariableNames[64]="alphaArm";
502 fCandidateVariableNames[65]="ptArm";
503
504 for(Int_t ivar=0; ivar<nVar; ivar++){
505 fVariablesTreeSgn->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
506 fVariablesTreeBkg->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/f",fCandidateVariableNames[ivar].Data()));
507 }
508
509 fHistoEvents = new TH1F("fHistoEvents", "fHistoEvents", 2, -0.5, 1.5);
510 TString labelEv[2] = {"NotSelected", "Selected"};
511 for (Int_t ibin = 1; ibin <= fHistoEvents->GetNbinsX(); ibin++){
512 fHistoEvents->GetXaxis()->SetBinLabel(ibin, labelEv[ibin-1].Data());
513 }
514
515 fHistoLc = new TH1F("fHistoLc", "fHistoLc", 2, -0.5, 1.5);
516
517 fHistoLcOnTheFly = new TH1F("fHistoLcOnTheFly", "fHistoLcOnTheFly", 4, -0.5, 3.5);
518 TString labelOnTheFly[4] = {"OnTheFly-Bkg", "OfflineBkg", "OnTheFly-Sgn", "OfflineSgn"};
519 for (Int_t ibin = 1; ibin <= fHistoLcOnTheFly->GetNbinsX(); ibin++){
520 fHistoLcOnTheFly->GetXaxis()->SetBinLabel(ibin, labelOnTheFly[ibin-1].Data());
521 }
522
523 fHistoLcBeforeCuts = new TH1F("fHistoLcBeforeCuts", "fHistoLcBeforeCuts", 2, -0.5, 1.5);
524 TString labelBeforeCuts[2] = {"Bkg", "Sgn"};
525 for (Int_t ibin = 1; ibin <= fHistoLcBeforeCuts->GetNbinsX(); ibin++){
526 fHistoLcBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
527 fHistoLc->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
528 }
529
530 fHistoFiducialAcceptance = new TH1F("fHistoFiducialAcceptance", "fHistoFiducialAcceptance", 4, -0.5, 3.5);
531 TString labelAcc[4] = {"NotOk-Bkg", "Ok-Bkg", "NotOk-Sgn", "Ok-Sgn"};
532 for (Int_t ibin = 1; ibin <= fHistoFiducialAcceptance->GetNbinsX(); ibin++){
533 fHistoFiducialAcceptance->GetXaxis()->SetBinLabel(ibin, labelAcc[ibin-1].Data());
534 }
535
536 fHistoCodesSgn = new TH2F("fHistoCodesSgn", "fHistoCodes for Signal; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
537 fHistoCodesBkg = new TH2F("fHistoCodesBkg", "fHistoCodes for Background; bachelor; K0S", 7, -1.5, 5.5, 9, -1.5, 7.5);
538
539 TString labelx[7] = { "kBachInvalid", "kBachFake", "kBachNoProton", "kBachPrimary", "kBachNoLambdaMother",
540 "kBachDifferentLambdaMother", "kBachCorrectLambdaMother"};
541 TString labely[9] = { "kK0SInvalid", "kK0SFake", "kK0SNoK0S", "kK0SWithoutMother", "kK0SNotFromK0",
542 "kK0Primary", "kK0NoLambdaMother", "kK0DifferentLambdaMother", "kK0CorrectLambdaMother"};
543
544 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsX(); ibin++){
545 fHistoCodesSgn->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
546 fHistoCodesBkg->GetXaxis()->SetBinLabel(ibin, labelx[ibin-1].Data());
547 }
548 for (Int_t ibin = 1; ibin <= fHistoCodesSgn->GetNbinsY(); ibin++){
549 fHistoCodesSgn->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
550 fHistoCodesBkg->GetYaxis()->SetBinLabel(ibin, labely[ibin-1].Data());
551 }
552
553 fHistoLcpKpiBeforeCuts = new TH1F("fHistoLcpKpiBeforeCuts", "fHistoLcpKpiBeforeCuts", 2, -0.5, 1.5);
554 for (Int_t ibin = 1; ibin <= fHistoLcpKpiBeforeCuts->GetNbinsX(); ibin++){
555 fHistoLcpKpiBeforeCuts->GetXaxis()->SetBinLabel(ibin, labelBeforeCuts[ibin-1].Data());
556 }
557
558 fHistoBackground = new TH1F("fHistoBackground", "fHistoBackground", 2, -0.5, 1.5);
559 TString labelBkg[2] = {"Injected", "Non-injected"};
560 for (Int_t ibin = 1; ibin <= fHistoBackground->GetNbinsX(); ibin++){
561 fHistoBackground->GetXaxis()->SetBinLabel(ibin, labelBkg[ibin-1].Data());
562 }
563
564 //fOutput->Add(fVariablesTreeSgn);
565 //fOutput->Add(fVariablesTreeBkg);
566
567 fOutput->Add(fHistoEvents);
568 fOutput->Add(fHistoLc);
569 fOutput->Add(fHistoLcOnTheFly);
570 fOutput->Add(fHistoLcBeforeCuts);
571 fOutput->Add(fHistoFiducialAcceptance);
572 fOutput->Add(fHistoCodesSgn);
573 fOutput->Add(fHistoCodesBkg);
574 fOutput->Add(fHistoLcpKpiBeforeCuts);
575 fOutput->Add(fHistoBackground);
576
577 PostData(1, fOutput);
578 PostData(4, fVariablesTreeSgn);
579 PostData(5, fVariablesTreeBkg);
580 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
581 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
582 fPIDResponse = inputHandler->GetPIDResponse();
583
584 if (fAnalCuts->GetIsUsePID()){
585 /*
586 fAnalCuts->GetPidHF()->SetPidResponse(fPIDResponse);
587 fAnalCuts->GetPidV0pos()->SetPidResponse(fPIDResponse);
588 fAnalCuts->GetPidV0neg()->SetPidResponse(fPIDResponse);
589 fAnalCuts->GetPidHF()->SetOldPid(kFALSE);
590 fAnalCuts->GetPidV0pos()->SetOldPid(kFALSE);
591 fAnalCuts->GetPidV0neg()->SetOldPid(kFALSE);
592 */
593 fAnalCuts->SetUsePID(kFALSE); // I don't want to use the PID through the cut object, but I will use the PID response directly!!!
594 }
595
596 // Setting properties of PID
597 fPIDCombined=new AliPIDCombined;
598 fPIDCombined->SetDefaultTPCPriors();
599 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
600 //fPIDCombined->SetPriorDistribution((AliPID::EParticleType)ispec,fPriors[ispec]);
601
602 //Printf("Here ok!");
603
604 fCounter = new AliNormalizationCounter("NormalizationCounter");
605 fCounter->Init();
606 PostData(2, fCounter);
607
608 // Histograms from KF
609
610 fHistoDistanceLcToPrimVtx = new TH1D("fHistoDistanceLcToPrimVtx", "Lc distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
611 fHistoDistanceV0ToPrimVtx = new TH1D("fHistoDistanceV0ToPrimVtx", "V0 distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
612 fHistoDistanceV0ToLc = new TH1D("fHistoDistanceV0ToLc", "V0 distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
613
614 fHistoDistanceLcToPrimVtxSgn = new TH1D("fHistoDistanceLcToPrimVtxSgn", "Lc Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 1);
615 fHistoDistanceV0ToPrimVtxSgn = new TH1D("fHistoDistanceV0ToPrimVtxSgn", "V0 Sgn distance to Prim Vertex from KF; distance [cm]", 1000, 0., 100.);
616 fHistoDistanceV0ToLcSgn = new TH1D("fHistoDistanceV0ToLcSgn", "V0 Sgn distance to Lc Vertex from KF; distance [cm]", 1000, 0., 100.);
617
618 fHistoVtxLcResidualToPrimVtx = new TH1D("fHistoVtxLcResidualToPrimVtx", "Residual between MC and KF (MC - KF): Lc to Prim Vtx; distance [cm]", 1000, -5., 5.);
619 fHistoVtxV0ResidualToPrimVtx = new TH1D("fHistoVtxV0ResidualToPrimVtx", "Residual between MC and KF (MC - KF): V0 to Prim Vtx; distance [cm]", 1000, -5., 5.);
620 fHistoVtxV0ResidualToLc = new TH1D("fHistoVtxV0ResidualToLc", "Residual between MC and KF: V0 to Lc (MC - KF); distance [cm]", 1000, -5., 5.);
621
622 fHistoMassV0All = new TH1D("fHistoMassV0All", "V0 Mass; mass", 500, 0.4, 0.6);
623 fHistoDecayLengthV0All = new TH1D("fHistoDecayLengthV0All", "V0 Decay Length; decayLength", 500, -10, 10.0);
624 fHistoLifeTimeV0All = new TH1D("fHistoLifeTimeV0All", "V0 Life Time; lifeTime", 500, -10.0, 10.0);
625
626 fHistoMassV0True = new TH1D("fHistoMassV0True", "True V0 Mass; mass", 500, 0.4, 0.6);
627 fHistoDecayLengthV0True = new TH1D("fHistoDecayLengthV0True", "True V0 Decay Length; decayLength", 500, -10, 10.0);
628 fHistoLifeTimeV0True = new TH1D("fHistoLifeTimeV0True", "True V0 Life Time; lifeTime", 500, -10.0, 10.0);
629
630 fHistoMassV0TrueFromAOD = new TH1D("fHistoMassV0TrueFormAOD", "True V0 Mass (AOD); mass", 500, 0.4, 0.6);
631
632 fHistoMassV0TrueK0S = new TH1D("fHistoMassV0TrueK0S", "True V0-K0S Mass; mass", 500, 0.4, 0.6);
633 fHistoDecayLengthV0TrueK0S = new TH1D("fHistoDecayLengthV0TrueK0S", "True V0-K0S Decay Length; decayLength", 500, -10, 10.0);
634 fHistoLifeTimeV0TrueK0S = new TH1D("fHistoLifeTimeV0TrueK0S", "True V0-K0S Life Time; lifeTime", 500, -10.0, 10.0);
635
636 fHistoMassV0TrueK0SFromAOD = new TH1D("fHistoMassV0TrueK0SFormAOD", "True V0-K0S Mass (AOD); mass", 500, 0.4, 0.6);
637
638 fHistoMassLcAll = new TH1D("fHistoMassLcAll", "Lc Mass; mass", 500, 2.0, 3.0);
639 fHistoDecayLengthLcAll = new TH1D("fHistoDecayLengthLcAll", "Lc Decay Length; decayLenght", 100000, -0.1, 0.1);
640 fHistoLifeTimeLcAll = new TH1D("fHistoLifeTimeLcAll", "Lc Life Time; lifeTime", 100000, -0.1, 0.1);
641
642 fHistoMassLcTrue = new TH1D("fHistoMassLcTrue", "True Lc Mass; mass", 500, 2.0, 3.0);
643 fHistoDecayLengthLcTrue = new TH1D("fHistoDecayLengthLcTrue", "True Lc Decay Length; decayLength", 100000, -0.1, 0.1);
644 fHistoLifeTimeLcTrue = new TH1D("fHistoLifeTimeLcTrue", "True Lc Life Time; lifeTime", 100000, -0.1, 0.1);
645
646 fHistoMassLcTrueFromAOD = new TH1D("fHistoMassLcTrueFromAOD", "True Lc Mass (AOD); mass", 500, 2.0, 3.0);
647
648 fHistoMassV0fromLcAll = new TH1D("fHistoMassV0fromLcAll", "V0 mass from Lc built in KF; mass", 500, 0.4, 0.6);
649 fHistoDecayLengthV0fromLcAll = new TH1D("fHistoDecayLengthV0fromLcAll", "V0 Decay Length from Lc built in KF; decayLength", 500, 0, 10.0);
650 fHistoLifeTimeV0fromLcAll = new TH1D("fHistoLifeTimeV0fromLcAll", "V0 Life Time from Lc built in KF; lifeTime", 500, 0.0, 3.0);
651
652 fHistoMassV0fromLcTrue = new TH1D("fHistoMassV0fromLcTrue", "V0 mass from true Lc built in KF; mass", 500, 0.4, 0.6);
653 fHistoDecayLengthV0fromLcTrue= new TH1D("fHistoDecayLengthV0fromLcTrue", "V0 Decay Length from true Lc built in KF; decayLength", 500, 0, 10.0);
654 fHistoLifeTimeV0fromLcTrue = new TH1D("fHistoLifeTimeV0fromLcTrue", "V0 Life Time from true Lc built in KF; lifeTime", 500, 0.0, 3.0);
655
656 fHistoMassLcSgn = new TH1D("fHistoMassLcSgn", "True Lc Signal Mass; mass", 500, 2.0, 3.0);
657 fHistoMassLcSgnFromAOD = new TH1D("fHistoMassLcSgnFromAOD", "True Lc Signal Mass (AOD); mass", 500, 2.0, 3.0);
658 fHistoDecayLengthLcSgn = new TH1D("fHistoDecayLengthLcSgn", "True Lc Signal Decay Length; decayLength", 100000, -0.1, 0.1);
659 fHistoLifeTimeLcSgn = new TH1D("fHistoLifeTimeLcSgn", "True Lc Signal Life Time; lifeTime", 100000, -0.1, 0.1);
660
661 fHistoMassV0fromLcSgn = new TH1D("fHistoMassV0fromLcSgn", "V0 from True Lc Signal Mass; mass", 500, 0.4, 0.6);
662 fHistoDecayLengthV0fromLcSgn = new TH1D("fHistoDecayLengthV0fromLcSgn", "V0 True Lc Signal Decay Length; decayLength", 500, 0, 10.0);
663 fHistoLifeTimeV0fromLcSgn = new TH1D("fHistoLifeTimeV0fromLcSgn", "V0 True Lc Signal Life Time; lifeTime", 500, 0.0, 3.0);
664
665 fHistoKF = new TH2D("fHistoKF", "Summary from KF; V0 KF; Lc KF", 16, -0.5, 15.5, 16, -0.5, 15.5);
666 fHistoKFV0 = new TH1D("fHistoKFV0", "Summary from KF; V0 KF", 16, -0.5, 15.5);
667 fHistoKFLc = new TH1D("fHistoKFLc", "Summary from KF; V0 KF", 16, -0.5, 15.5);
668 TString axisLabel[16] = {"AllOk", "M_NotOk", "Sm_NotOk", "Dl_NotOk",
669 "Lt_NotOk", "M_Sm_NotOk", "M_Dl_NotOk", "M_Lt_NotOk",
670 "Dl_Sm_NotOk", "Dl_Lt_NotOk", "Sm_Lt_NotOk", "M_Sm_Dl_NotOk",
671 "M_Sm_Lt_NotOk", "Sm_Dl_Lt_NotOk", "M_Dl_Lt_NotOk", "All_NotOk"};
672
673 for (Int_t ibin = 1; ibin <=16; ibin++){
674 fHistoKF->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
675 fHistoKF->GetYaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
676 fHistoKFV0->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
677 fHistoKFLc->GetXaxis()->SetBinLabel(ibin, axisLabel[ibin-1].Data());
678 }
679
680 fHistoMassKFV0 = new TH2D("fHistoMassKFV0", "mass vs sigmaMass for V0; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
681 fHistoDecayLengthKFV0 = new TH2D("fHistoDecayLengthKFV0", "decayLength vs sigmaDecayLength for V0; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
682 fHistoLifeTimeKFV0 = new TH2D("fHistoLifeTimeKFV0", "lifeTime vs sigmalifeTime for V0; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
683
684 fHistoMassKFLc = new TH2D("fHistoMassKFLc", "mass vs sigmaMass for Lc; mass; sigmaMass", 500, 0.4, 0.6, 500, 0., 10);
685 fHistoDecayLengthKFLc = new TH2D("fHistoDecayLengthKFLc", "decayLength vs sigmaDecayLength for Lc; decayLength; sigmaDecayLength", 500, -10, 10, 500, 0., 10);
686 fHistoLifeTimeKFLc = new TH2D("fHistoLifeTimeKFLc", "lifeTime vs sigmalifeTime for Lc; lifeTime; sigmaLifeTime", 500, -10, 10, 500, 0., 10);
687
688 fHistoArmenterosPodolanskiV0KF = new TH2D("fHistoArmenterosPodolanskiV0KF", "V0 ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
689 fHistoArmenterosPodolanskiV0KFSgn = new TH2D("fHistoArmenterosPodolanskiV0KFSgn", "V0 (signal) ArmenterosPodolanski from KF; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
690 fHistoArmenterosPodolanskiV0AOD = new TH2D("fHistoArmenterosPodolanskiV0AOD", "V0 ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
691 fHistoArmenterosPodolanskiV0AODSgn = new TH2D("fHistoArmenterosPodolanskiV0AODSgn", "V0 (signal) ArmenterosPodolanski from AOD; #alpha; Qt", 1000, -1, 1, 1000, 0, 1);
692
693 fOutputKF = new TList();
694 fOutputKF->SetOwner();
695 fOutputKF->SetName("listHistoKF");
696
697 fOutputKF->Add(fHistoDistanceLcToPrimVtx);
698 fOutputKF->Add(fHistoDistanceV0ToPrimVtx);
699 fOutputKF->Add(fHistoDistanceV0ToLc);
700
701 fOutputKF->Add(fHistoDistanceLcToPrimVtxSgn);
702 fOutputKF->Add(fHistoDistanceV0ToPrimVtxSgn);
703 fOutputKF->Add(fHistoDistanceV0ToLcSgn);
704
705 fOutputKF->Add(fHistoVtxLcResidualToPrimVtx);
706 fOutputKF->Add(fHistoVtxV0ResidualToPrimVtx);
707 fOutputKF->Add(fHistoVtxV0ResidualToLc);
708
709 fOutputKF->Add(fHistoMassV0All);
710 fOutputKF->Add(fHistoDecayLengthV0All);
711 fOutputKF->Add(fHistoLifeTimeV0All);
712
713 fOutputKF->Add(fHistoMassV0True);
714 fOutputKF->Add(fHistoDecayLengthV0True);
715 fOutputKF->Add(fHistoLifeTimeV0True);
716
717 fOutputKF->Add(fHistoMassV0TrueFromAOD);
718
719 fOutputKF->Add(fHistoMassV0TrueK0S);
720 fOutputKF->Add(fHistoDecayLengthV0TrueK0S);
721 fOutputKF->Add(fHistoLifeTimeV0TrueK0S);
722
723 fOutputKF->Add(fHistoMassV0TrueK0SFromAOD);
724
725 fOutputKF->Add(fHistoMassLcAll);
726 fOutputKF->Add(fHistoDecayLengthLcAll);
727 fOutputKF->Add(fHistoLifeTimeLcAll);
728
729 fOutputKF->Add(fHistoMassLcTrue);
730 fOutputKF->Add(fHistoDecayLengthLcTrue);
731 fOutputKF->Add(fHistoLifeTimeLcTrue);
732
733 fOutputKF->Add(fHistoMassLcTrueFromAOD);
734
735 fOutputKF->Add(fHistoMassV0fromLcAll);
736 fOutputKF->Add(fHistoDecayLengthV0fromLcAll);
737 fOutputKF->Add(fHistoLifeTimeV0fromLcAll);
738
739 fOutputKF->Add(fHistoMassV0fromLcTrue);
740 fOutputKF->Add(fHistoDecayLengthV0fromLcTrue);
741 fOutputKF->Add(fHistoLifeTimeV0fromLcTrue);
742
743 fOutputKF->Add(fHistoMassLcSgn);
122f753d 744 fOutputKF->Add(fHistoMassLcSgnFromAOD);
0abe5247 745 fOutputKF->Add(fHistoDecayLengthLcSgn);
746 fOutputKF->Add(fHistoLifeTimeLcSgn);
747
748 fOutputKF->Add(fHistoMassV0fromLcSgn);
749 fOutputKF->Add(fHistoDecayLengthV0fromLcSgn);
750 fOutputKF->Add(fHistoLifeTimeV0fromLcSgn);
751
752 fOutputKF->Add(fHistoKF);
753 fOutputKF->Add(fHistoKFV0);
754 fOutputKF->Add(fHistoKFLc);
755
756 fOutputKF->Add(fHistoMassKFV0);
757 fOutputKF->Add(fHistoDecayLengthKFV0);
758 fOutputKF->Add(fHistoLifeTimeKFV0);
759
760 fOutputKF->Add(fHistoMassKFLc);
761 fOutputKF->Add(fHistoDecayLengthKFLc);
762 fOutputKF->Add(fHistoLifeTimeKFLc);
763
764 fOutputKF->Add(fHistoArmenterosPodolanskiV0KF);
765 fOutputKF->Add(fHistoArmenterosPodolanskiV0KFSgn);
766 fOutputKF->Add(fHistoArmenterosPodolanskiV0AOD);
767 fOutputKF->Add(fHistoArmenterosPodolanskiV0AODSgn);
768
769 PostData(6, fOutputKF);
770
771 //Printf("Here ok 1!");
772
773 return;
774}
775
776//_________________________________________________
777void AliAnalysisTaskSELc2V0bachelorTMVA::UserExec(Option_t *)
778{
779 // user exec
780 if (!fInputEvent) {
781 AliError("NO EVENT FOUND!");
782 return;
783 }
784
785 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
786 TClonesArray *arrayLctopKos=0;
787
788 TClonesArray *array3Prong = 0;
789
790 if (!aodEvent && AODEvent() && IsStandardAOD()) {
791 // In case there is an AOD handler writing a standard AOD, use the AOD
792 // event in memory rather than the input (ESD) event.
793 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
794 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
795 // have to taken from the AOD event hold by the AliAODExtension
796 AliAODHandler* aodHandler = (AliAODHandler*)
797 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
798
799 if (aodHandler->GetExtensions()) {
800 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
801 AliAODEvent *aodFromExt = ext->GetAOD();
802 arrayLctopKos=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
803
804 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
805 }
806 } else {
807 arrayLctopKos=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
808
809 array3Prong=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
810 }
811
812 if ( !fUseMCInfo && fIspA) {
813 fAnalCuts->SetTriggerClass("");
814 fAnalCuts->SetTriggerMask(AliVEvent::kINT7);
815 }
816
817 fCounter->StoreEvent(aodEvent,fAnalCuts,fUseMCInfo);
818
819 // AOD primary vertex
820 fVtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
821 if (!fVtx1) return;
822 if (fVtx1->GetNContributors()<1) return;
823
824 fIsEventSelected = fAnalCuts->IsEventSelected(aodEvent);
825
826 if ( !fIsEventSelected ) {
827 fHistoEvents->Fill(0);
828 return; // don't take into account not selected events
829 }
830 fHistoEvents->Fill(1);
831
832 // Setting magnetic field for KF vertexing
833 AliKFParticle::SetField(aodEvent->GetMagneticField());
834
835 // mc analysis
836 TClonesArray *mcArray = 0;
837 AliAODMCHeader *mcHeader=0;
838
839 if (fUseMCInfo) {
840 // MC array need for maching
841 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
842 if (!mcArray) {
843 AliError("Could not find Monte-Carlo in AOD");
844 return;
845 }
846 // load MC header
847 mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
848 if (!mcHeader) {
849 AliError("AliAnalysisTaskSELc2V0bachelorTMVA::UserExec: MC header branch not found!\n");
850 return;
851 }
852 }
853
854 Int_t nSelectedAnal = 0;
855 if (fIsK0sAnalysis) {
856 MakeAnalysisForLc2prK0S(arrayLctopKos, mcArray,
857 nSelectedAnal, fAnalCuts,
858 array3Prong, mcHeader);
859 }
860 fCounter->StoreCandidates(aodEvent,nSelectedAnal,kFALSE);
861
862 PostData(1, fOutput);
863 PostData(2, fCounter);
864 PostData(4, fVariablesTreeSgn);
865 PostData(5, fVariablesTreeBkg);
866 PostData(6, fOutputKF);
867
868}
869
870//-------------------------------------------------------------------------------
871void AliAnalysisTaskSELc2V0bachelorTMVA::MakeAnalysisForLc2prK0S(TClonesArray *arrayLctopKos,
872 TClonesArray *mcArray,
873 Int_t &nSelectedAnal,
874 AliRDHFCutsLctoV0 *cutsAnal, TClonesArray *array3Prong,
875 AliAODMCHeader* aodheader){
876 //Lc prong needed to MatchToMC method
877
878 Int_t pdgCand = 4122;
879 Int_t pdgDgLctoV0bachelor[2]={2212, 310};
880 Int_t pdgDgV0toDaughters[2]={211, 211};
881
882 Int_t pdgDgLctopKpi[3]={2212, 321, 211};
883
884 // loop to search for candidates Lc->p+K+pi
885 Int_t n3Prong = array3Prong->GetEntriesFast();
886 Int_t nCascades= arrayLctopKos->GetEntriesFast();
887
888 //AliInfo(Form("\n\n\n\n3 prong candidates = %d, ncascades = %d \n\n\n\n\n", n3Prong, nCascades));
889 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
890 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
891 //Filling a control histogram with no cuts
892 if (fUseMCInfo) {
893
894 // find associated MC particle for Lc -> p+K+pi
895 Int_t mcLabel = d->MatchToMC(4122, mcArray, 3, pdgDgLctopKpi);
896 //Int_t mcLabelTemp = d->MatchToMC(4122, mcArray);
897 //Printf("mcLabel = %d, mcLabelTemp = %d", mcLabel, mcLabelTemp);
898 if (mcLabel >= 0) {
899
900 AliAODMCParticle *partLcpKpi = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
901 if(partLcpKpi){
902 Int_t pdgCode = partLcpKpi->GetPdgCode();
903 AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
904 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
905 fHistoLcpKpiBeforeCuts->Fill(1);
906
907 }
908 }
909 else {
910 //AliInfo(Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~", mcLabel));
911 fHistoLcpKpiBeforeCuts->Fill(0);
912 }
913 }
914 }
915
916 // loop over cascades to search for candidates Lc->p+K0S
917
918 for (Int_t iLctopK0s = 0; iLctopK0s < nCascades; iLctopK0s++) {
919
920 // Lc candidates and K0s from Lc
921 AliAODRecoCascadeHF* lcK0spr = dynamic_cast<AliAODRecoCascadeHF*>(arrayLctopKos->At(iLctopK0s));
922 if (!lcK0spr) {
923 AliDebug(2,Form("Cascade %d doens't exist, skipping",iLctopK0s));
924 continue;
925 }
926
927 //Filling a control histogram with no cuts
928 if (fUseMCInfo) {
929
930 Int_t pdgCode=-2;
931
932 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
933 fmcLabelLc = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
934 if (fmcLabelLc>=0) {
935 AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
936
937 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(fmcLabelLc));
938 if(partLc){
939 pdgCode = partLc->GetPdgCode();
940 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", fmcLabelLc, pdgCode));
941 pdgCode = TMath::Abs(pdgCode);
942 fHistoLcBeforeCuts->Fill(1);
943
944 }
945 }
946 else {
947 fHistoLcBeforeCuts->Fill(0);
948 pdgCode=-1;
949 }
950 }
951
952 //if (!lcK0spr->GetSecondaryVtx()) {
953 // AliInfo("No secondary vertex");
954 //continue;
955 //}
956
957 if (lcK0spr->GetNDaughters()!=2) {
958 AliDebug(2,Form("Cascade %d has not 2 daughters (nDaughters=%d)",iLctopK0s,lcK0spr->GetNDaughters()));
959 continue;
960 }
961
962 AliAODv0 * v0part = dynamic_cast<AliAODv0*>(lcK0spr->Getv0());
963 AliAODTrack * bachPart = dynamic_cast<AliAODTrack*>(lcK0spr->GetBachelor());
964 if (!v0part || !bachPart) {
965 AliDebug(2,Form("Cascade %d has no V0 or no bachelor object",iLctopK0s));
966 continue;
967 }
968
969
970 if (!v0part->GetSecondaryVtx()) {
971 AliDebug(2,Form("No secondary vertex for V0 by cascade %d",iLctopK0s));
972 continue;
973 }
974
975 if (v0part->GetNDaughters()!=2) {
976 AliDebug(2,Form("current V0 has not 2 daughters (onTheFly=%d, nDaughters=%d)",v0part->GetOnFlyStatus(),v0part->GetNDaughters()));
977 continue;
978 }
979
980 AliAODTrack * v0Pos = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0PositiveTrack());
981 AliAODTrack * v0Neg = dynamic_cast<AliAODTrack*>(lcK0spr->Getv0NegativeTrack());
982 if (!v0Neg || !v0Neg) {
983 AliDebug(2,Form("V0 by cascade %d has no V0positive of V0negative object",iLctopK0s));
984 continue;
985 }
986
987
988 if (v0Pos->Charge() == v0Neg->Charge()) {
989 AliDebug(2,Form("V0 by cascade %d has daughters with the same sign: IMPOSSIBLE!",iLctopK0s));
990 continue;
991 }
992
993 Int_t isLc = 0;
994
995 if (fUseMCInfo) {
996
997 Int_t pdgCode=-2;
998
999 // find associated MC particle for Lc -> p+K0 and K0S->pi+pi
1000 Int_t mcLabel = lcK0spr->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1001 if (mcLabel>=0) {
1002 AliDebug(2,Form(" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~cascade number %d (total cascade number = %d)", iLctopK0s, nCascades));
1003
1004 AliAODMCParticle *partLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
1005 if(partLc){
1006 pdgCode = partLc->GetPdgCode();
1007 if (pdgCode<0) AliDebug(2,Form(" ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ MClabel=%d ~~~~~~~~~~ pdgCode=%d", mcLabel, pdgCode));
1008 pdgCode = TMath::Abs(pdgCode);
1009 isLc = 1;
1010 fHistoLc->Fill(1);
1011
1012 }
1013 }
1014 else {
1015 fHistoLc->Fill(0);
1016 pdgCode=-1;
1017 }
1018 }
1019 AliDebug(2, Form("\n\n\n Analysing candidate %d\n", iLctopK0s));
1020 AliDebug(2, Form(">>>>>>>>>> Candidate is background, fFillOnlySgn = %d --> SKIPPING", fFillOnlySgn));
1021 if (!isLc) {
1022 if (fFillOnlySgn) { // if it is background, and we want only signal, we do not fill the tree
1023 continue;
1024 }
1025 else { // checking if we want to fill the background
1026 if (fKeepingOnlyHIJINGBkg){
1027 // we have decided to fill the background only when the candidate has the daugthers that all come from HIJING underlying event!
122f753d 1028 Bool_t isCandidateInjected = fUtils->HasCascadeCandidateAnyDaughInjected(lcK0spr, aodheader, mcArray);
0abe5247 1029 if (!isCandidateInjected){
1030 AliDebug(2, "The candidate is from HIJING (i.e. not injected), keeping it to fill background");
1031 fHistoBackground->Fill(1);
1032 }
1033 else {
1034 AliDebug(2, "The candidate is NOT from HIJING, we skip it when filling background");
1035 fHistoBackground->Fill(0);
1036 continue;
1037 }
1038 }
1039 }
1040 }
1041
1042 FillLc2pK0Sspectrum(lcK0spr, isLc, nSelectedAnal, cutsAnal, mcArray);
1043
1044 }
1045
1046 return;
1047
1048}
1049//________________________________________________________________________
1050void AliAnalysisTaskSELc2V0bachelorTMVA::FillLc2pK0Sspectrum(AliAODRecoCascadeHF *part,
1051 Int_t isLc,
1052 Int_t &nSelectedAnal,
1053 AliRDHFCutsLctoV0 *cutsAnal,
1054 TClonesArray *mcArray){
1055 //
1056 // Fill histos for Lc -> K0S+proton
1057 //
1058
1059 /*
1060 if (!part->GetOwnPrimaryVtx()) {
1061 //Printf("No primary vertex for Lc found!!");
1062 part->SetOwnPrimaryVtx(fVtx1);
1063 }
1064 else {
1065 //Printf("Yu-huuuu!!! primary vertex for Lc found!!");
1066 }
1067 */
1068 Double_t invmassLc = part->InvMassLctoK0sP();
1069
1070 AliAODv0 * v0part = part->Getv0();
1071 Bool_t onFlyV0 = v0part->GetOnFlyStatus(); // on-the-flight V0s
1072 if (onFlyV0){ // on-the-fly V0
1073 if (isLc) { // Lc
1074 fHistoLcOnTheFly->Fill(2.);
1075 }
1076 else { // not Lc
1077 fHistoLcOnTheFly->Fill(0.);
1078 }
1079 }
1080 else { // offline V0
1081 if (isLc) { // Lc
1082 fHistoLcOnTheFly->Fill(3.);
1083 }
1084 else { // not Lc
1085 fHistoLcOnTheFly->Fill(1.);
1086 }
1087 }
1088
1089 Double_t dcaV0 = v0part->GetDCA();
1090 Double_t invmassK0s = v0part->MassK0Short();
1091
1092 if ( (cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) {
1093 if (isLc) {
1094 fHistoFiducialAcceptance->Fill(3.);
1095 }
1096 else {
1097 fHistoFiducialAcceptance->Fill(1.);
1098 }
1099 }
1100 else {
1101 if (isLc) {
1102 fHistoFiducialAcceptance->Fill(2.);
1103 }
1104 else {
1105 fHistoFiducialAcceptance->Fill(0.);
1106 }
1107 }
1108
1109 Int_t isInV0window = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 2)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on V0 invMass
1110
1111 if (isInV0window == 0) {
1112 AliDebug(2, "No: The candidate has NOT passed the mass cuts!");
1113 if (isLc) Printf("SIGNAL candidate rejected");
1114 return;
1115 }
1116 else AliDebug(2, "Yes: The candidate has passed the mass cuts!");
1117
1118 // Printf("ciccio!!!!!!!");
1119
1120 Bool_t isInCascadeWindow = (((cutsAnal->IsSelectedSingleCut(part, AliRDHFCuts::kCandidate, 0)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // cut on Lc->p+K0S invMass
1121
1122 if (!isInCascadeWindow) {
1123 AliDebug(2, "No: The candidate has NOT passed the cascade window cuts!");
1124 if (isLc) Printf("SIGNAL candidate rejected");
1125 return;
1126 }
1127 else AliDebug(2, "Yes: The candidate has passed the cascade window cuts!");
1128
1129 Bool_t isCandidateSelectedCuts = (((cutsAnal->IsSelected(part, AliRDHFCuts::kCandidate)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr)); // kinematic/topological cuts
1130 if (!isCandidateSelectedCuts){
1131 AliDebug(2, "No: Analysis cuts kCandidate level NOT passed");
1132 if (isLc) Printf("SIGNAL candidate rejected");
1133 return;
1134 }
1135 else {
1136 AliDebug(2, "Yes: Analysis cuts kCandidate level passed");
1137 }
1138
1139 AliAODTrack *bachelor = (AliAODTrack*)part->GetBachelor();
1140 if (!bachelor) {
1141 AliDebug(2, Form("Very weird, the bachelor is not there... returning for this candidate"));
1142 return;
1143 }
1144
1145 //Bool_t isBachelorID = (((cutsAnal->IsSelected(part,AliRDHFCuts::kPID))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)); // ID x bachelor
1146 Double_t probTPCTOF[AliPID::kSPECIES]={0.};
1147
1148 UInt_t detUsed = fPIDCombined->ComputeProbabilities(bachelor, fPIDResponse, probTPCTOF);
1149 Double_t probProton = 0.;
1150 Double_t probPion = 0.;
1151 Double_t probKaon = 0.;
1152 if (detUsed == (UInt_t)fPIDCombined->GetDetectorMask() ) {
1153 probProton = probTPCTOF[AliPID::kProton];
1154 probPion = probTPCTOF[AliPID::kPion];
1155 probKaon = probTPCTOF[AliPID::kKaon];
1156 }
1157
1158 if ( !( !onFlyV0 || (onFlyV0 && fUseOnTheFlyV0) ) ) {
1159 AliDebug(2, "On-the-fly discarded");
1160 return;
1161 }
1162
1163 if ( (((cutsAnal->IsSelected(part,AliRDHFCuts::kAll))&(AliRDHFCutsLctoV0::kLcToK0Spr))==(AliRDHFCutsLctoV0::kLcToK0Spr)) ) {
1164 nSelectedAnal++;
1165 }
1166
1167 if ( !(cutsAnal->IsInFiducialAcceptance(part->Pt(),part->Y(4122))) ) return;
1168
1169 if ( !( ( (cutsAnal->IsSelected(part, AliRDHFCuts::kTracks)) & (AliRDHFCutsLctoV0::kLcToK0Spr)) == (AliRDHFCutsLctoV0::kLcToK0Spr) ) ) { // esd track cuts
1170 if (isLc) Printf("SIGNAL candidate rejected");
1171 AliDebug(2, "No: Analysis cuts kTracks level NOT passed");
1172 return;
1173 }
1174 else {
1175 AliDebug(2, "Yes: Analysis cuts kTracks level passed");
1176 }
1177
1178 Int_t pdgCand = 4122;
1179 Int_t pdgDgLctoV0bachelor[2]={211, 3122}; // case of wrong decay! Lc --> L + pi
1180 Int_t pdgDgV0toDaughters[2]={2212, 211}; // case of wrong decay! Lc --> L + pi
1181 Int_t isLc2LBarpi=0, isLc2Lpi=0;
1182 Int_t currentLabel = part->GetLabel();
1183 Int_t mcLabel = 0;
1184 if (fUseMCInfo) {
1185 mcLabel = part->MatchToMC(pdgCand, pdgDgLctoV0bachelor[1], pdgDgLctoV0bachelor, pdgDgV0toDaughters, mcArray, kTRUE);
1186 if (mcLabel>=0) {
1187 if (bachelor->Charge()==-1) isLc2LBarpi=1;
1188 if (bachelor->Charge()==+1) isLc2Lpi=1;
1189 }
1190 }
1191
1192 Int_t pdgDg2prong[2] = {211, 211};
1193 Int_t labelK0S = 0;
1194 Int_t isK0S = 0;
1195 if (fUseMCInfo) {
1196 labelK0S = v0part->MatchToMC(310, mcArray, 2, pdgDg2prong);
1197 if (labelK0S>=0) isK0S = 1;
1198 }
1199
1200 pdgDg2prong[0] = 211;
1201 pdgDg2prong[1] = 2212;
1202 Int_t isLambda = 0;
1203 Int_t isLambdaBar = 0;
1204 Int_t lambdaLabel = 0;
1205 if (fUseMCInfo) {
1206 lambdaLabel = v0part->MatchToMC(3122, mcArray, 2, pdgDg2prong);
1207 if (lambdaLabel>=0) {
1208 AliAODMCParticle *lambdaTrack = (AliAODMCParticle*)mcArray->At(lambdaLabel);
1209 if (lambdaTrack->GetPdgCode()==3122) isLambda = 1;
1210 else if (lambdaTrack->GetPdgCode()==-3122) isLambdaBar = 1;
1211 }
1212 }
1213
1214 pdgDg2prong[0] = 11;
1215 pdgDg2prong[1] = 11;
1216 Int_t isGamma = 0;
1217 Int_t gammaLabel = 0;
1218 if (fUseMCInfo) {
1219 gammaLabel = v0part->MatchToMC(22, mcArray, 2, pdgDg2prong);
1220 if (gammaLabel>=0) {
1221 AliAODMCParticle *gammaTrack = (AliAODMCParticle*)mcArray->At(gammaLabel);
1222 if (gammaTrack->GetPdgCode()==22) isGamma = 1;
1223 }
1224 }
1225
1226 Int_t pdgTemp = -1;
1227 if (currentLabel != -1){
1228 AliAODMCParticle *tempPart = (AliAODMCParticle*)mcArray->At(currentLabel);
1229 pdgTemp = tempPart->GetPdgCode();
1230 }
1231 if (isLc) AliDebug(2, Form("Signal: Candidate is a Lc in K0s+p"));
1232 else if (isLc2LBarpi) AliDebug(2, Form("Background: Candidate is a Lc in Lbar + pi"));
1233 else if (isLc2Lpi) AliDebug(2, Form("Background: Candidate is a Lc in L + pi"));
1234 else AliDebug(2, Form("Pure bkg: Candidate has pdg = %d", pdgTemp));
1235 if (isK0S) AliDebug(2, Form("V0 is a K0S"));
1236 else if (isLambda) AliDebug(2, Form("V0 is a Lambda"));
1237 else if (isLambdaBar) AliDebug(2, Form("V0 is a LambdaBar"));
1238 else if (isGamma) AliDebug(2, Form("V0 is a Gamma"));
1239 //else AliDebug(2, Form("V0 is something else!!"));
1240
1241 Double_t invmassLc2Lpi = part->InvMassLctoLambdaPi();
1242 Double_t invmassLambda = v0part->MassLambda();
1243 Double_t invmassLambdaBar = v0part->MassAntiLambda();
1244
1245 //Double_t nSigmaITSpr=-999.;
1246 Double_t nSigmaTPCpr=-999.;
1247 Double_t nSigmaTOFpr=-999.;
1248
1249 //Double_t nSigmaITSpi=-999.;
1250 Double_t nSigmaTPCpi=-999.;
1251 Double_t nSigmaTOFpi=-999.;
1252
1253 //Double_t nSigmaITSka=-999.;
1254 Double_t nSigmaTPCka=-999.;
1255 Double_t nSigmaTOFka=-999.;
1256
1257 /*
1258 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,4,nSigmaITSpr);
1259 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,4,nSigmaTPCpr);
1260 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,4,nSigmaTOFpr);
1261 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,2,nSigmaITSpi);
1262 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,2,nSigmaTPCpi);
1263 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,2,nSigmaTOFpi);
1264 cutsAnal->GetPidHF()->GetnSigmaITS(bachelor,3,nSigmaITSka);
1265 cutsAnal->GetPidHF()->GetnSigmaTPC(bachelor,3,nSigmaTPCka);
1266 cutsAnal->GetPidHF()->GetnSigmaTOF(bachelor,3,nSigmaTOFka);
1267 */
1268
1269 nSigmaTPCpi = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kPion));
1270 nSigmaTPCka = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kKaon));
1271 nSigmaTPCpr = fPIDResponse->NumberOfSigmasTPC(bachelor,(AliPID::kProton));
1272 nSigmaTOFpi = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kPion));
1273 nSigmaTOFka = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kKaon));
1274 nSigmaTOFpr = fPIDResponse->NumberOfSigmasTOF(bachelor,(AliPID::kProton));
1275
1276
1277 // Fill candidate variable Tree (track selection, V0 invMass selection)
1278 if (!onFlyV0 && isInV0window && isInCascadeWindow && part->CosV0PointingAngle()>0.99 && TMath::Abs(nSigmaTPCpr) <= 3 && v0part->Getd0Prong(0) < 20 && v0part->Getd0Prong(1) < 20) {
1279
1280 fCandidateVariables[0] = invmassLc;
1281 fCandidateVariables[1] = invmassLc2Lpi;
1282 fCandidateVariables[2] = invmassK0s;
1283 fCandidateVariables[3] = invmassLambda;
1284 fCandidateVariables[4] = invmassLambdaBar;
1285 fCandidateVariables[5] = part->CosV0PointingAngle();
1286 fCandidateVariables[6] = dcaV0;
1287 fCandidateVariables[7] = part->Getd0Prong(0);
1288 fCandidateVariables[8] = part->Getd0Prong(1);
1289 fCandidateVariables[9] = nSigmaTPCpr;
1290 fCandidateVariables[10] = nSigmaTPCpi;
1291 fCandidateVariables[11] = nSigmaTPCka;
1292 fCandidateVariables[12] = nSigmaTOFpr;
1293 fCandidateVariables[13] = nSigmaTOFpi;
1294 fCandidateVariables[14] = nSigmaTOFka;
1295 fCandidateVariables[15] = bachelor->Pt();
1296 AliAODTrack *v0pos = (AliAODTrack*)part->Getv0PositiveTrack();
1297 fCandidateVariables[16] = v0pos->Pt();
1298 AliAODTrack *v0neg = (AliAODTrack*)part->Getv0NegativeTrack();
1299 fCandidateVariables[17] = v0neg->Pt();
1300 fCandidateVariables[18] = v0part->Getd0Prong(0);
1301 fCandidateVariables[19] = v0part->Getd0Prong(1);
1302 fCandidateVariables[20] = v0part->Pt();
1303 fCandidateVariables[21] = v0part->InvMass2Prongs(0,1,11,11);
1304 fCandidateVariables[22] = part->Pt();
1305 fCandidateVariables[23] = probProton;
1306 fCandidateVariables[24] = part->Eta();
1307 fCandidateVariables[25] = v0pos->Eta();
1308 fCandidateVariables[26] = v0neg->Eta();
1309 fCandidateVariables[27] = probPion;
1310 fCandidateVariables[28] = probKaon;
1311 fCandidateVariables[29] = bachelor->Eta();
1312
1313 fCandidateVariables[30] = part->P();
1314 fCandidateVariables[31] = bachelor->P();
1315 fCandidateVariables[32] = v0part->P();
1316 fCandidateVariables[33] = v0pos->P();
1317 fCandidateVariables[34] = v0neg->P();
1318
1319 fCandidateVariables[35] = part->Y(4122);
1320 fCandidateVariables[36] = bachelor->Y(2212);
1321 fCandidateVariables[37] = v0part->Y(310);
1322 fCandidateVariables[38] = v0pos->Y(211);
1323 fCandidateVariables[39] = v0neg->Y(211);
1324
1325 fCandidateVariables[40] = v0part->Eta();
1326
1327 fCandidateVariables[41] = part->DecayLength();
1328 fCandidateVariables[42] = part->DecayLengthV0();
1329 fCandidateVariables[43] = part->Ct(4122);
1330 fCandidateVariables[44] = v0part->Ct(310, v0part->GetSecondaryVtx());
1331
1332 EBachelor bachCode = CheckBachelor(part, bachelor, mcArray);
1333 EK0S k0SCode = CheckK0S(part, v0part, mcArray);
1334
1335 fCandidateVariables[45] = bachCode;
1336 fCandidateVariables[46] = k0SCode;
1337
1338 Double_t V0KF[3] = {-999999, -999999, -999999};
1339 Double_t errV0KF[3] = {-999999, -999999, -999999};
1340 Double_t LcKF[3] = {-999999, -999999, -999999};
1341 Double_t errLcKF[3] = {-999999, -999999, -999999};
1342 Double_t distances[3] = {-999999, -999999, -999999};
1343 Double_t armPolKF[2] = {-999999, -999999};
1344
1345 if (fCallKFVertexing){
1346 Int_t kfResult = CallKFVertexing(part, v0part, bachelor, mcArray, &V0KF[0], &errV0KF[0], &LcKF[0], &errLcKF[0], &distances[0], &armPolKF[0]);
1347 AliDebug(2, Form("Result from KF = %d", kfResult));
1348 }
1349
1350 /*
1351 for (Int_t i = 0; i< 3; i++){
1352 Printf("i = %d, V0KF = %f, errV0KF = %f, LcKF = %f, errLcKF = %f", V0KF[i], errV0KF[i], LcKF[i], errLcKF[i]);
1353 }
1354 */
1355
1356 fCandidateVariables[47] = V0KF[0];
1357 fCandidateVariables[48] = V0KF[1];
1358 fCandidateVariables[49] = V0KF[2];
1359
1360 fCandidateVariables[50] = errV0KF[0];
1361 fCandidateVariables[51] = errV0KF[1];
1362 fCandidateVariables[52] = errV0KF[2];
1363
1364 fCandidateVariables[53] = LcKF[0];
1365 fCandidateVariables[54] = LcKF[1];
1366 fCandidateVariables[55] = LcKF[2];
1367
1368 fCandidateVariables[56] = errLcKF[0];
1369 fCandidateVariables[57] = errLcKF[1];
1370 fCandidateVariables[58] = errLcKF[2];
1371
1372 fCandidateVariables[59] = distances[0];
1373 fCandidateVariables[60] = distances[1];
1374 fCandidateVariables[61] = distances[2];
1375 fCandidateVariables[62] = armPolKF[0];
1376 fCandidateVariables[63] = armPolKF[1];
1377 fCandidateVariables[64] = v0part->AlphaV0();
1378 fCandidateVariables[65] = v0part->PtArmV0();
1379
1380 if (fUseMCInfo) {
1381 if (isLc){
1382 AliDebug(2, "Filling Sgn");
1383 fVariablesTreeSgn->Fill();
1384 fHistoCodesSgn->Fill(bachCode, k0SCode);
1385 }
1386 else {
1387 if (fFillOnlySgn == kFALSE){
1388 AliDebug(2, "Filling Bkg");
1389 fVariablesTreeBkg->Fill();
1390 fHistoCodesBkg->Fill(bachCode, k0SCode);
1391 }
1392 }
1393 }
1394 else {
1395 fVariablesTreeSgn->Fill();
1396 }
1397 }
1398
1399 return;
1400
1401}
1402
1403//________________________________________________________________________
1404Int_t AliAnalysisTaskSELc2V0bachelorTMVA::CallKFVertexing(AliAODRecoCascadeHF *cascade, AliAODv0* v0part, AliAODTrack* bach, TClonesArray *mcArray,
1405 Double_t* V0KF, Double_t* errV0KF, Double_t* LcKF, Double_t* errLcKF,
1406 Double_t* distances, Double_t* armPolKF) {
1407
1408 //
1409 // method to perform KF vertexing
1410 // elements: [0] = mass, [1] = DecayLength, [2] = lifeTime
1411 //
1412
1413 Int_t codeKFV0 = -1, codeKFLc = -1;
1414
1415 AliKFVertex primVtxCopy;
1416 Int_t nt = 0, ntcheck = 0;
1417 Double_t pos[3] = {0., 0., 0.};
1418
1419 fVtx1->GetXYZ(pos);
1420 Int_t contr = fVtx1->GetNContributors();
1421 Double_t covmatrix[6] = {0.};
1422 fVtx1->GetCovarianceMatrix(covmatrix);
1423 Double_t chi2 = fVtx1->GetChi2();
1424 AliESDVertex primaryESDVtxCopy(pos, covmatrix, chi2, contr, "Vertex");
1425
1426 // topological constraint
1427 primVtxCopy = AliKFVertex(primaryESDVtxCopy);
1428 nt = primaryESDVtxCopy.GetNContributors();
1429 ntcheck = nt;
1430
1431 Int_t pdg[2] = {211, -211};
1432 Int_t pdgLc[2] = {2212, 310};
1433
1434 Int_t pdgDgV0toDaughters[2] = {211, 211};
1435
1436 Int_t mcLabelV0 = v0part->MatchToMC(310, mcArray, 2, pdgDgV0toDaughters);
1437
1438 // the KF vertex for the V0 has to be built with the prongs of the V0!
1439 Bool_t isMCokV0 = kTRUE, isBkgV0 = kFALSE;
1440 AliKFParticle V0, positiveV0KF, negativeV0KF;;
1441 Int_t labelsv0daugh[2] = {-1, -1};
1442 for(Int_t ipr= 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1443 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1444 if(!aodTrack) {
1445 AliDebug(2, "No V0 daughters available");
1446 return -1;
1447 }
1448 labelsv0daugh[ipr] = TMath::Abs(aodTrack->GetLabel());
1449 if (labelsv0daugh[ipr] == -1) isBkgV0 = kTRUE;
1450
1451 //Printf("v0 daughter %d has label %d", ipr, labelsv0daugh[ipr]);
1452 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1453 if (aodTrack->Charge() > 0) { // assigning positive and negative track to KF V0 for Armenteros-Podolanski plot
1454 positiveV0KF = daughterKF;
1455 }
1456 else {
1457 negativeV0KF = daughterKF;
1458 }
1459 V0.AddDaughter(daughterKF);
1460 }
1461
1462 // Checking the quality of the KF V0 vertex
1463 if( V0.GetNDF() < 1 ) {
1464 //Printf("Number of degrees of freedom < 1, continuing");
1465 return -1;
1466 }
1467 if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) > 3. ) {
1468 //Printf("Chi2 per DOF too big, continuing");
1469 return -1;
1470 }
1471
1472 if(ftopoConstraint && nt > 0){
1473 for(Int_t ipr = 0; ipr < 2; ipr++){ // 0 is positive, 1 is negative
1474 AliAODTrack *aodTrack = (AliAODTrack*)v0part->GetDaughter(ipr);
1475 //* subtruct daughters from primary vertex
1476 if(!aodTrack->GetUsedForPrimVtxFit()) {
1477 //Printf("Track %d was not used for primary vertex, continuing", i);
1478 continue;
1479 }
1480 AliKFParticle daughterKF(*aodTrack, pdg[ipr]); // we assume that the PDG is correct
1481 primVtxCopy -= daughterKF;
1482 ntcheck--;
1483 }
1484 }
1485
1486 /*
1487 // Check V0 Chi^2 deviation from primary vertex // not needed for V0 for Lc decay!!
1488 if( V0.GetDeviationFromVertex( primVtxCopy ) >3. ) {
1489 //Printf("Deviation from vertex too big, continuing");
1490 return -1;
1491 }
1492 */
1493
1494 if(ftopoConstraint){
1495 if(ntcheck>0) { // number of constraints left with which the vertex was built, after removing V0 daughters
1496 //* Add V0 to primary vertex to improve the primary vertex resolution
1497 primVtxCopy += V0;
1498 V0.SetProductionVertex(primVtxCopy); // FIXME: is this correct here?
1499 }
1500 }
1501
1502 //* Check chi^2
1503 if( TMath::Sqrt( TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) {
1504 AliDebug(2, "Final Chi2 per DOF too big, continuing");
1505 return -1;
1506 }
1507
1508 //* Get V0 invariant mass
1509 Double_t massV0 = 999999, sigmaMassV0 = 999999;
1510 Int_t retMV0 = V0.GetMass( massV0, sigmaMassV0 );
1511 if( retMV0 ) {
1512 if (massV0 < 0) {
1513 codeKFV0 = 1; // Mass not ok
1514 if (sigmaMassV0 > 1e19) codeKFV0 = 5; // Mass and SigmaMass not ok
1515 }
1516 else if (sigmaMassV0 > 1e19) codeKFV0 = 2; // SigmaMass not ok
1517 }
1518 fHistoMassKFV0->Fill(massV0, sigmaMassV0);
1519
1520 //* Get V0 decayLength
1521 Double_t decayLengthV0 = 999999, sigmaDecayLengthV0 = 999999;
1522 Int_t retDLV0 = V0.GetDecayLength( decayLengthV0, sigmaDecayLengthV0 );
1523 if( retDLV0 ) {
1524 if (sigmaDecayLengthV0 > 1e19) {
1525 codeKFV0 = 3; // DecayLength not ok
1526 if (massV0 < 0) {
1527 codeKFV0 = 6; // DecayLength and Mass not ok
1528 if (sigmaMassV0 > 1e19) codeKFV0 = 11; // DecayLength and Mass and SigmaMass not ok
1529 }
1530 else if (sigmaMassV0 > 1e19) codeKFV0 = 8; // DecayLength and SigmaMass not ok
1531 }
1532 }
1533 fHistoDecayLengthKFV0->Fill(decayLengthV0, sigmaDecayLengthV0);
1534
1535 //* Get V0 life time
1536 Double_t lifeTimeV0 = 999999, sigmaLifeTimeV0 = 999999;
1537 Int_t retTLV0 = V0.GetLifeTime( lifeTimeV0, sigmaLifeTimeV0 );
1538 if( retTLV0 ) {
1539 if (sigmaLifeTimeV0 > 1e19) {
1540 codeKFV0 = 4; // LifeTime not ok
1541 if (sigmaDecayLengthV0 > 1e19) {
1542 codeKFV0 = 9; // LifeTime and DecayLength not ok
1543 if (massV0 < 0) {
1544 codeKFV0 = 14; // LifeTime and DecayLength and Mass not ok
1545 if (sigmaMassV0 > 1e19) codeKFV0 = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
1546 }
1547 else if (sigmaMassV0 > 1e19) codeKFV0 = 13; // LifeTime and DecayLength and SigmaMass not ok
1548 }
1549 else if (massV0 < 0) { // LifeTime and Mass and SigmaMass not ok
1550 codeKFV0 = 7; // LifeTime and Mass not ok
1551 if (sigmaMassV0 > 1e19) codeKFV0 = 12; // LifeTime and Mass and SigmaMass not ok
1552 }
1553 else if (sigmaMassV0 > 1e19) codeKFV0 = 10; // LifeTime and SigmaMass not ok
1554 }
1555 }
1556 fHistoLifeTimeKFV0->Fill(lifeTimeV0, sigmaLifeTimeV0);
1557
1558 if (codeKFV0 == -1) codeKFV0 = 0;
1559 fHistoKFV0->Fill(codeKFV0);
1560
1561 AliDebug(2, Form("V0: mass = %f, decay length = %f, life time = %f", massV0, decayLengthV0, lifeTimeV0 ));
1562
1563 fHistoMassV0All->Fill(massV0);
1564 fHistoDecayLengthV0All->Fill(decayLengthV0);
1565 fHistoLifeTimeV0All->Fill(lifeTimeV0);
1566
1567 Double_t qtAlphaV0[2] = {0., 0.};
1568 Double_t vtxV0KF[3] = {V0.GetX(), V0.GetY(), V0.GetZ()};
1569 positiveV0KF.TransportToPoint(vtxV0KF);
1570 negativeV0KF.TransportToPoint(vtxV0KF);
1571 V0.GetArmenterosPodolanski(positiveV0KF, negativeV0KF, qtAlphaV0);
1572 AliDebug(2, Form("Armenteros-Podolanski variables: alpha = %f, qt = %f", qtAlphaV0[1], qtAlphaV0[0]));
1573 fHistoArmenterosPodolanskiV0KF->Fill(qtAlphaV0[1], qtAlphaV0[0]);
1574 fHistoArmenterosPodolanskiV0AOD->Fill(v0part->AlphaV0(), v0part->PtArmV0());
1575 armPolKF[0] = qtAlphaV0[1];
1576 armPolKF[1] = qtAlphaV0[0];
1577
1578 // Checking MC info for V0
1579
1580 AliAODMCParticle *motherV0 = 0x0;
1581 AliAODMCParticle *daughv01 = 0x0;
1582 AliAODMCParticle *daughv02 = 0x0;
1583
1584 if (fUseMCInfo) {
1585 daughv01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[0]));
1586 daughv02 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsv0daugh[1]));
1587 if (!daughv01 && labelsv0daugh[0] > 0){
1588 AliDebug(2, "Could not access MC info for first daughter of V0, continuing");
1589 isMCokV0 = kFALSE;
1590 }
1591 if (!daughv02 && labelsv0daugh[1] > 0){
1592 AliDebug(2, "Could not access MC info for second daughter of V0, continuing");
1593 isMCokV0 = kFALSE;
1594 }
1595 if (isMCokV0){
1596 if( daughv01->GetMother() == daughv02->GetMother() && daughv01->GetMother()>=0 ){
1597 AliDebug(3, Form("The mother has label %d", daughv01->GetMother()));
1598 motherV0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01->GetMother()));
1599 if( motherV0 && TMath::Abs(motherV0->GetPdgCode()) != 21 ){ // These are all V0 that are truly V0, not only K0S, but no gluons
1600 if( motherV0->GetNDaughters() == 2 ){
1601 fHistoMassV0True->Fill(massV0);
1602 fHistoDecayLengthV0True->Fill(decayLengthV0);
1603 fHistoLifeTimeV0True->Fill(lifeTimeV0);
1604 fHistoMassV0TrueFromAOD->Fill(v0part->MassK0Short());
1605 if (TMath::Abs(motherV0->GetPdgCode()) == 310){ // These are true V0 that are also K0S
1606 fHistoMassV0TrueK0S->Fill(massV0);
1607 fHistoDecayLengthV0TrueK0S->Fill(decayLengthV0);
1608 fHistoLifeTimeV0TrueK0S->Fill(lifeTimeV0);
1609 fHistoMassV0TrueK0SFromAOD->Fill(v0part->MassK0Short());
1610 }
1611 }
1612 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()));
1613 }
1614 else if (!motherV0){
1615 AliDebug(3, "could not access MC info for mother, continuing");
1616 }
1617 else {
1618 AliDebug(3, "MC mother is a gluon, continuing");
1619 }
1620 }
1621 else {
1622 AliDebug(3, "Background V0!");
1623 isBkgV0 = kTRUE;
1624 }
1625 }
1626 }
1627
1628 AliDebug(2, Form("isMCokV0 = %d, isBkgV0 = %d", (Int_t)isMCokV0, (Int_t)isBkgV0));
1629
1630 // now use what we just try with the bachelor, to build the Lc
1631
1632 // topological constraint
1633 nt = primVtxCopy.GetNContributors();
1634 ntcheck = nt;
1635
1636 Bool_t isMCokLc = kTRUE, isBkgLc = kFALSE;
1637 AliKFParticle Lc;
1638 Int_t labelsLcdaugh[2] = {-1, -1};
1639 labelsLcdaugh[0] = TMath::Abs(bach->GetLabel());
1640 labelsLcdaugh[1] = mcLabelV0;
1641
1642 if (bach->Charge() < 0) pdgLc[0] = -pdgLc[0];
1643 AliKFParticle daughterKFLc(*bach, pdgLc[0]);
1644 Lc.AddDaughter(daughterKFLc);
1645 Lc.AddDaughter(V0);
1646 if( Lc.GetNDF() < 1 ) {
1647 AliDebug(2, Form("Lc: Number of degrees of freedom < 1 (%d), continuing", Lc.GetNDF()));
1648 return -1;
1649 }
1650 if( TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) >3. ) {
1651 AliDebug(2, Form("Lc: Chi2 per DOF too big, continuing (%f)", TMath::Sqrt(TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1652 return -1;
1653 }
1654
1655 if(ftopoConstraint && nt > 0){
1656 //* subtruct daughters from primary vertex
1657 if(!bach->GetUsedForPrimVtxFit()) {
1658 AliDebug(3, "Lc: Bachelor was not used for primary vertex, not subtracting it from primary vertex");
1659 }
1660 else{
1661 primVtxCopy -= daughterKFLc;
1662 ntcheck--;
1663 }
1664 /* the V0 was added above, so it is ok to remove it without checking
1665 if(!V0->GetUsedForPrimVtxFit()) {
1666 Printf("Lc: V0 was not used for primary vertex, continuing");
1667 continue;
1668 }
1669 */
1670 primVtxCopy -= V0;
1671 ntcheck--;
1672 }
1673
1674 //* Check V0 Chi^2 deviation from primary vertex
1675 if( Lc.GetDeviationFromVertex( primVtxCopy ) >3. ) {
1676 AliDebug(2, Form("Lc: Deviation from vertex too big, continuing (%f)", Lc.GetDeviationFromVertex( primVtxCopy )));
1677 return -1;
1678 }
1679
1680 if(ftopoConstraint){
1681 if(ntcheck>0) {
1682 //* Add Lc to primary vertex to improve the primary vertex resolution
1683 primVtxCopy += Lc;
1684 Lc.SetProductionVertex(primVtxCopy);
1685 }
1686 }
1687
1688 //* Check chi^2
1689 if( TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF())) >3. ) {
1690 AliDebug(2, Form("Lc: Final Chi2 per DOF too big, continuing (%f)", TMath::Sqrt( TMath::Abs(Lc.GetChi2()/Lc.GetNDF()))));
1691 return -1;
1692 }
1693
1694 //* Get Lc invariant mass
1695 Double_t massLc = 999999, sigmaMassLc= 999999;
1696 Int_t retMLc = Lc.GetMass( massLc, sigmaMassLc );
1697 if( retMLc ) {
1698 AliDebug(3, Form("----> Could not get mass (%e), and sigma(%e) for Lc, continuing", massLc, sigmaMassLc));
1699 if (massLc < 0) {
1700 codeKFLc = 1; // Mass not ok
1701 if (sigmaMassLc > 1e19) codeKFLc = 5; // Mass and SigmaMass not ok
1702 }
1703 else if (sigmaMassLc > 1e19) codeKFLc = 2; // SigmaMass not ok
1704 }
1705 fHistoMassKFLc->Fill(massLc, sigmaMassLc);
1706
1707 //* Get Lc Decay length
1708 Double_t decayLengthLc = 999999, sigmaDecayLengthLc = 999999;
1709 Int_t retDLLc = Lc.GetDecayLength( decayLengthLc, sigmaDecayLengthLc );
1710 if( retDLLc ) {
1711 AliDebug(3, "----> Lc: Could not get decay length, and sigma");
1712 if (sigmaDecayLengthLc > 1e19) {
1713 codeKFLc = 3; // DecayLength not ok
1714 if (massLc < 0) {
1715 codeKFLc = 6; // DecayLength and Mass not ok
1716 if (sigmaMassLc > 1e19) codeKFLc = 11; // DecayLength and Mass and SigmaMass not ok
1717 }
1718 else if (sigmaMassLc > 1e19) codeKFLc = 8; // DecayLength and SigmaMass not ok
1719 }
1720 }
1721 AliDebug(3, Form("retDLLc = %d, with decayLength = %f and error = %e", retDLLc, decayLengthLc, sigmaDecayLengthLc));
1722
1723 fHistoDecayLengthKFLc->Fill(decayLengthLc, sigmaDecayLengthLc);
1724
1725 //* Get Lc life time
1726 Double_t lifeTimeLc = 999999, sigmaLifeTimeLc = 999999;
1727 Int_t retTLLc = Lc.GetLifeTime( lifeTimeLc, sigmaLifeTimeLc );
1728 if( retTLLc ) {
1729 AliDebug(3, "----> Lc: Could not get lifeTime, and sigma");
1730 if (sigmaLifeTimeLc > 1e19) {
1731 codeKFLc = 4; // LifeTime not ok
1732 if (sigmaDecayLengthLc > 1e19) {
1733 codeKFLc = 9; // LifeTime and DecayLength not ok
1734 if (massLc < 0) {
1735 codeKFLc = 14; // LifeTime and DecayLength and Mass not ok
1736 if (sigmaMassLc > 1e19) codeKFLc = 15; // LifeTime and DecayLength and Mass and SigmaMass not ok
1737 }
1738 else if (sigmaMassLc > 1e19) codeKFLc = 13; // LifeTime and DecayLength and SigmaMass not ok
1739 }
1740 else if (massLc < 0) { // LifeTime and Mass and SigmaMass not ok
1741 codeKFLc = 7; // LifeTime and Mass not ok
1742 if (sigmaMassLc > 1e19) codeKFLc = 12; // LifeTime and Mass and SigmaMass not ok
1743 }
1744 else if (sigmaMassLc > 1e19) codeKFLc = 10; // LifeTime and SigmaMass not ok
1745 }
1746 }
1747
1748 fHistoLifeTimeKFLc->Fill(lifeTimeLc, sigmaLifeTimeLc);
1749
1750 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));
1751
1752 if (codeKFLc == -1) codeKFLc = 0;
1753 fHistoKFLc->Fill(codeKFLc);
1754
1755 fHistoKF->Fill(codeKFV0, codeKFLc);
1756
1757 // here we fill the histgrams for all the reconstructed KF vertices for the cascade
1758 fHistoMassLcAll->Fill(massLc);
1759 fHistoDecayLengthLcAll->Fill(decayLengthLc);
1760 fHistoLifeTimeLcAll->Fill(lifeTimeLc);
1761
1762 fHistoMassV0fromLcAll->Fill(massV0);
1763 fHistoDecayLengthV0fromLcAll->Fill(decayLengthV0);
1764 fHistoLifeTimeV0fromLcAll->Fill(lifeTimeV0);
1765
1766 Double_t xV0 = V0.GetX();
1767 Double_t yV0 = V0.GetY();
1768 Double_t zV0 = V0.GetZ();
1769
1770 Double_t xLc = Lc.GetX();
1771 Double_t yLc = Lc.GetY();
1772 Double_t zLc = Lc.GetZ();
1773
1774 Double_t xPrimVtx = primVtxCopy.GetX();
1775 Double_t yPrimVtx = primVtxCopy.GetY();
1776 Double_t zPrimVtx = primVtxCopy.GetZ();
1777
1778 Double_t distanceLcToPrimVtx = TMath::Sqrt((xPrimVtx - xLc) * (xPrimVtx - xLc) +
1779 (yPrimVtx - yLc) * (yPrimVtx - yLc) +
1780 (zPrimVtx - zLc) * (zPrimVtx - zLc));
1781
1782 Double_t distanceV0ToPrimVtx = TMath::Sqrt((xPrimVtx - xV0) * (xPrimVtx - xV0) +
1783 (yPrimVtx - yV0) * (yPrimVtx - yV0) +
1784 (zPrimVtx - zV0) * (zPrimVtx - zV0));
1785
1786 Double_t distanceV0ToLc = TMath::Sqrt((xLc - xV0)*(xLc - xV0) +
1787 (yLc - yV0)*(yLc - yV0) +
1788 (zLc - zV0)*(zLc - zV0));
1789
1790 //Printf("distanceLcToPrimVtx = %e, distanceV0ToPrimVtx= %f, distanceV0ToLc = %f", distanceLcToPrimVtx, distanceV0ToPrimVtx, distanceV0ToLc);
1791
1792 fHistoDistanceLcToPrimVtx->Fill(distanceLcToPrimVtx);
1793 fHistoDistanceV0ToPrimVtx->Fill(distanceV0ToPrimVtx);
1794 fHistoDistanceV0ToLc->Fill(distanceV0ToLc);
1795
1796 distances[0] = distanceLcToPrimVtx;
1797 distances[1] = distanceV0ToPrimVtx;
1798 distances[2] = distanceV0ToLc;
1799
1800 if (fUseMCInfo) {
1801
1802 AliAODMCParticle *daughv01Lc = 0x0;
1803 AliAODMCParticle *K0S = 0x0;
1804 AliAODMCParticle *daughv02Lc = 0x0;
1805
1806 if (labelsLcdaugh[0] >= 0) {
1807 // Printf("Getting Bachelor from label %d", labelsLcdaugh[1]);
1808 daughv01Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[0]));
1809 if (!daughv01Lc){
1810 AliDebug(3, "Could not access MC info for first daughter of Lc");
1811 isMCokLc = kFALSE;
1812 }
1813 else {
1814 AliDebug(2, Form("The bachelor has label = %d", daughv01Lc->GetLabel()));
1815 if (TMath::Abs(daughv01Lc->GetPdgCode()) != 2212) isBkgLc = kTRUE;
1816 }
1817 }
1818 else { // The bachelor is a fake
1819 isBkgLc = kTRUE;
1820 }
1821
1822 if (labelsLcdaugh[1] >= 0) {
1823 //Printf("Getting K0S");
1824 K0S = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelsLcdaugh[1]));
1825 if (!K0S) {
1826 AliDebug(3, "Could not access MC info for second daughter of Lc");
1827 isMCokLc = kFALSE;
1828 }
1829 else{
1830 if (TMath::Abs(K0S->GetPdgCode()) != 310) isBkgLc = kTRUE;
1831 }
1832 }
1833 else{
1834 AliDebug(2, "The K0S is not true --> it does not have a label, continuing...");
1835 isBkgLc = kTRUE;
1836 }
1837
1838 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
1839 if (isMCokLc) { // We can then access its MC info, and it might then be that also the Lc is a true Lc
1840 Int_t iK0 = K0S->GetMother();
1841 if (iK0 < 0) {
1842 Printf("The K0S has no mother... IMPOSSIBLE"); // the K0S MUST have a mother!
1843 }
1844 else { // The K0S has a mother
1845 daughv02Lc = dynamic_cast<AliAODMCParticle*>(mcArray->At(iK0));
1846 if (!daughv02Lc){
1847 AliDebug(3, "Could not access MC info for second daughter of Lc");
1848 }
1849 else { // we can access safely the K0S mother in the MC
1850 if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 ){ // This is a true cascade! bachelor and V0 come from the same mother
1851 //Printf("Lc: The mother has label %d", daughv01Lc->GetMother());
1852 AliAODMCParticle *motherLc = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughv01Lc->GetMother()));
1853 Int_t pdgMum = 0, pdgBach = 0, pdgV0 = 0;
1854 if (motherLc) pdgMum = motherLc->GetPdgCode();
1855 if (daughv01Lc) pdgBach = daughv01Lc->GetPdgCode();
1856 if (daughv02Lc) pdgV0 = daughv02Lc->GetPdgCode();
1857 AliDebug(2, Form("pdgMum = %d, pdgBach = %d, pdgV0 = %d", pdgMum, pdgBach, pdgV0));
1858
1859 if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 ){ // These are true cascades, we don't know yet if they are Lc
1860 fHistoMassLcTrue->Fill(massLc);
1861 fHistoDecayLengthLcTrue->Fill(decayLengthLc);
1862 fHistoLifeTimeLcTrue->Fill(lifeTimeLc);
1863 fHistoMassLcTrueFromAOD->Fill(cascade->InvMassLctoK0sP());
1864
1865 fHistoMassV0fromLcTrue->Fill(massV0);
1866 fHistoDecayLengthV0fromLcTrue->Fill(decayLengthV0);
1867 fHistoLifeTimeV0fromLcTrue->Fill(lifeTimeV0);
1868
1869 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...
1870 AliDebug(2, Form("IT IS SIGNAL!!! with label = %d", motherLc->GetLabel()));
1871
1872 fHistoArmenterosPodolanskiV0KFSgn->Fill(qtAlphaV0[1], qtAlphaV0[0]);
1873 fHistoArmenterosPodolanskiV0AODSgn->Fill(v0part->AlphaV0(), v0part->PtArmV0());
1874
1875 fHistoDistanceLcToPrimVtxSgn->Fill(distanceLcToPrimVtx);
1876 fHistoDistanceV0ToPrimVtxSgn->Fill(distanceV0ToPrimVtx);
1877 fHistoDistanceV0ToLcSgn->Fill(distanceV0ToLc);
1878
1879 fHistoMassLcSgn->Fill(massLc);
1880 fHistoMassLcSgnFromAOD->Fill(cascade->InvMassLctoK0sP());
1881
1882 fHistoDecayLengthLcSgn->Fill(decayLengthLc);
1883 fHistoLifeTimeLcSgn->Fill(lifeTimeLc);
1884
1885 fHistoMassV0fromLcSgn->Fill(massV0);
1886 fHistoDecayLengthV0fromLcSgn->Fill(decayLengthV0);
1887 fHistoLifeTimeV0fromLcSgn->Fill(lifeTimeV0);
1888 }
1889 else {
1890 isBkgLc = kTRUE; // it is not a Lc, since the pdg != 4122
1891 }
1892
1893 // if we got here, we can compare with MC information; this part is done also for cases where the cascade is not a Lc!
1894 // First, we compare the vtx of the Lc
1895 Double_t xLcMC = motherLc->Xv();
1896 Double_t yLcMC = motherLc->Yv();
1897 Double_t zLcMC = motherLc->Zv();
1898 //Printf("Got MC vtx for Lc");
1899 Double_t xProtonMC = daughv01Lc->Xv();
1900 Double_t yProtonMC = daughv01Lc->Yv();
1901 Double_t zProtonMC = daughv01Lc->Zv();
1902 //Printf("Got MC vtx for Proton");
1903
1904 Double_t vtxLcResidualToPrimVtx = TMath::Sqrt((xLcMC - xProtonMC) * (xLcMC - xProtonMC) +
1905 (yLcMC - yProtonMC) * (yLcMC - yProtonMC) +
1906 (zLcMC - zProtonMC) * (zLcMC - zProtonMC)) - distanceLcToPrimVtx;
1907
1908 // Then, we compare the vtx of the K0s
1909 Double_t xV0MC = motherV0->Xv(); //Production vertex of K0S --> Where Lc decays
1910 Double_t yV0MC = motherV0->Yv();
1911 Double_t zV0MC = motherV0->Zv();
1912 //Printf("Got MC vtx for V0");
1913 Double_t xPionMC = daughv01->Xv(); //Production vertex of Pion --> Where K0S decays
1914 Double_t yPionMC = daughv01->Yv();
1915 Double_t zPionMC = daughv01->Zv();
1916 //Printf("Got MC vtx for Pion");
1917
1918 Double_t vtxV0ResidualToLc = TMath::Sqrt((xV0MC - xPionMC) * (xV0MC - xPionMC) +
1919 (yV0MC - yPionMC) * (yV0MC - yPionMC) +
1920 (zV0MC - zPionMC) * (zV0MC - zPionMC)) - distanceV0ToLc;
1921
1922 // Then, the K0S vertex wrt primary vertex
1923
1924 Double_t vtxV0ResidualToPrimVtx = TMath::Sqrt((xPionMC - xLcMC) * (xPionMC - xLcMC) +
1925 (yPionMC - yLcMC) * (yPionMC - yLcMC) +
1926 (zPionMC - zLcMC) * (zPionMC - zLcMC)) - distanceV0ToPrimVtx;
1927
1928 fHistoVtxLcResidualToPrimVtx->Fill(vtxLcResidualToPrimVtx);
1929 fHistoVtxV0ResidualToLc->Fill(vtxV0ResidualToLc);
1930 fHistoVtxV0ResidualToPrimVtx->Fill(vtxV0ResidualToPrimVtx);
1931
1932 } //closing: if( motherLc && TMath::Abs(motherLc->GetPdgCode()) != 21 )
1933 else if (!motherLc){
1934 AliDebug(2, "We could not access MC info for Lc mother, so we did nothing");
1935 }
1936 else {
1937 AliDebug(2, "MC Lc mother is a gluon, so we did nothing");
1938 }
1939 } //closing: if( daughv01Lc->GetMother() == daughv02Lc->GetMother() && daughv01Lc->GetMother()>=0 )
1940 else {
1941 isBkgLc = kTRUE; // it cannot be a Lc, since the daughters do not have the same mother
1942 }
1943 } // closing: else { // we can access safely the K0S mother in the MC
1944 } // closing: else { // The K0S has a mother
1945 } // closing isMCLcok
1946 } // closing !isBkgLc
1947 } // closing fUseMCInfo
1948
1949 //Printf("retMV0 = %d, retMLc = %d", retMV0, retMLc);
1950 if ( retMV0 == 0 && retMLc == 0){
1951 V0KF[0] = massV0;
1952 errV0KF[0] = sigmaMassV0;
1953 V0KF[1] = decayLengthV0;
1954 errV0KF[1] = sigmaDecayLengthV0;
1955 V0KF[2] = lifeTimeV0;
1956 errV0KF[2] = sigmaLifeTimeV0;
1957 LcKF[0] = massLc;
1958 errLcKF[0] = sigmaMassLc;
1959 LcKF[1] = decayLengthLc;
1960 errLcKF[1] = sigmaDecayLengthLc;
1961 LcKF[2] = lifeTimeLc;
1962 errLcKF[2] = sigmaLifeTimeLc;
1963 }
1964
1965 return 1;
1966
1967}
1968//________________________________________________________________________
1969AliAnalysisTaskSELc2V0bachelorTMVA::EBachelor AliAnalysisTaskSELc2V0bachelorTMVA::CheckBachelor( AliAODRecoCascadeHF *part,
1970 AliAODTrack* bachelor,
1971 TClonesArray *mcArray ){
1972
1973 //Printf("In CheckBachelor");
1974
1975 // function to check if the bachelor is a p, if it is a p but not from Lc
1976 // to be filled for background candidates
1977
1978 Int_t label = bachelor->GetLabel();
1979 if (label == -1) {
1980 return kBachFake;
1981 }
1982
1983 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(TMath::Abs(label)));
1984 if (!mcpart) return kBachInvalid;
1985 Int_t pdg = mcpart->PdgCode();
1986 if (TMath::Abs(pdg) != 2212) {
1987 AliDebug(2, Form("Bachelor is not a p, but a particle with pdg code = %d", pdg));
1988 return kBachNoProton;
1989 }
1990 else { // it is a proton
1991 //Int_t labelLc = part->GetLabel();
1992 Int_t labelLc = FindLcLabel(part, mcArray);
1993 //Printf(">>>>>>>>>>>>> label for Lc = %d", labelLc);
1994 Int_t bachelorMotherLabel = mcpart->GetMother();
1995 //Printf(">>>>>>>>>>>>> label for bachelorMother = %d", bachelorMotherLabel);
1996 if (bachelorMotherLabel == -1) {
1997 return kBachPrimary;
1998 }
1999 AliAODMCParticle *bachelorMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(bachelorMotherLabel));
2000 if (!bachelorMother) return kBachInvalid;
2001 Int_t pdgMother = bachelorMother->PdgCode();
2002 if (TMath::Abs(pdgMother) != 4122) {
2003 AliDebug(2, Form("The proton does not come from a Lc, but from a particle with pdgcode = %d", pdgMother));
2004 return kBachNoLambdaMother;
2005 }
2006 else { // it comes from Lc
2007 if (labelLc != bachelorMotherLabel){
2008 //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));
2009 AliDebug(2, Form("The proton comes from a Lc, but it is not the candidate we are analyzing"));
2010 return kBachDifferentLambdaMother;
2011 }
2012 else { // it comes from the correct Lc
2013 AliDebug(2, Form("The proton comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2014 return kBachCorrectLambdaMother;
2015 }
2016 }
2017 }
2018
2019 return kBachInvalid;
2020
2021}
2022
2023//________________________________________________________________________
2024AliAnalysisTaskSELc2V0bachelorTMVA::EK0S AliAnalysisTaskSELc2V0bachelorTMVA::CheckK0S( AliAODRecoCascadeHF *part,
2025 AliAODv0* v0part,
2026 //AliAODTrack* v0part,
2027 TClonesArray *mcArray ){
2028
2029 // function to check if the K0Spart is a p, if it is a p but not from Lc
2030 // to be filled for background candidates
2031
2032 //Printf(" CheckK0S");
2033
2034 //Int_t labelMatchToMC = v0part->MatchToMC(310, mcArray);
2035 //Int_t label = v0part->GetLabel();
2036 Int_t labelFind = FindV0Label(v0part, mcArray);
2037 //Printf("\n >>>>>>>>>>>>> label for V0 = %d, from MatchToMC = %d, from Find = %d", label, labelMatchToMC, labelFind);
2038 if (labelFind == -1) {
2039 return kK0SFake;
2040 }
2041
2042 AliAODMCParticle *mcpart = dynamic_cast<AliAODMCParticle*>(mcArray->At(labelFind));
2043 if (!mcpart) return kK0SInvalid;
2044 Int_t pdg = mcpart->PdgCode();
2045 if (TMath::Abs(pdg) != 310) {
2046 AliDebug(2, Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2047 //AliInfo(Form("K0Spart is not a K0S, but a particle with pdg code = %d", pdg));
2048 return kK0SNoK0S;
2049 }
2050 else { // it is a K0S
2051 //Int_t labelLc = part->GetLabel();
2052 Int_t labelLc = FindLcLabel(part, mcArray);
2053 Int_t K0SpartMotherLabel = mcpart->GetMother();
2054 if (K0SpartMotherLabel == -1) {
2055 return kK0SWithoutMother;
2056 }
2057 AliAODMCParticle *K0SpartMother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0SpartMotherLabel)); // should be a K0 in case of signal
2058 if (!K0SpartMother) return kK0SInvalid;
2059 Int_t pdgMotherK0S = K0SpartMother->PdgCode();
2060 if (TMath::Abs(pdgMotherK0S) != 311) {
2061 AliDebug(2, Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2062 //AliInfo(Form("The K0S does not come from a K0, but from a particle with pdgcode = %d", pdgMotherK0S));
2063 return kK0SNotFromK0;
2064 }
2065 else { // the K0S comes from a K0
2066 Int_t K0MotherLabel = K0SpartMother->GetMother(); // mother of K0 --> Lc in case of signal
2067 if (K0MotherLabel == -1) {
2068 return kK0Primary;
2069 }
2070 AliAODMCParticle *K0Mother = dynamic_cast<AliAODMCParticle*>(mcArray->At(K0MotherLabel));
2071 if (!K0Mother) return kK0SInvalid;
2072 Int_t pdgK0Mother = K0Mother->PdgCode();
2073 if (TMath::Abs(pdgK0Mother) != 4122) { // the K0 does not come from a Lc
2074 AliDebug(2, Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2075 //AliInfo(Form("The K0 does not come from a Lc, but from a particle with pdgcode = %d", pdgK0Mother));
2076 return kK0NoLambdaMother;
2077 }
2078 else { // the K0 comes from Lc
2079 if (labelLc != K0MotherLabel){ // The K0 comes from a different Lc
2080 AliDebug(2, Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2081 //AliInfo(Form("The K0S comes from a Lc, but it is not the candidate we are analyzing"));
2082 return kK0DifferentLambdaMother;
2083 }
2084 else { // it comes from the correct Lc
2085 AliDebug(2, Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2086 //AliInfo(Form("The K0S comes from a Lc, and it is EXACTLY the candidate we are analyzing"));
2087 return kK0CorrectLambdaMother;
2088 }
2089 }
2090 }
2091 }
2092
2093 return kK0SInvalid;
2094
2095}
2096
2097//----------------------------------------------------------------------------
2098Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindV0Label(AliAODRecoDecay* v0part, TClonesArray *mcArray) const
2099{
2100
2101 //Printf(" FindV0Label");
2102
2103 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2104
2105 Int_t labMother[2]={-1, -1};
2106 AliAODMCParticle *part=0;
2107 AliAODMCParticle *mother=0;
2108 Int_t dgLabels = -1;
2109
2110 Int_t ndg = v0part->GetNDaughters();
2111 if (ndg != 2) {
2112 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2113 }
2114
2115 for(Int_t i = 0; i < 2; i++) {
2116 AliAODTrack *trk = (AliAODTrack*)v0part->GetDaughter(i);
2117 dgLabels = trk->GetLabel();
2118 if (dgLabels == -1) {
2119 //printf("daughter with negative label %d\n", dgLabels);
2120 return -1;
2121 }
2122 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2123 if (!part) {
2124 //printf("no MC particle\n");
2125 return -1;
2126 }
2127 labMother[i] = part->GetMother();
2128 if (labMother[i] != -1){
2129 mother = (AliAODMCParticle*)mcArray->At(labMother[i]);
2130 if(!mother) {
2131 //printf("no MC mother particle\n");
2132 return -1;
2133 }
2134 }
2135 else {
2136 return -1;
2137 }
2138 }
2139
2140 if (labMother[0] == labMother[1]) return labMother[0];
2141
2142 else return -1;
2143
2144}
2145
2146//----------------------------------------------------------------------------
2147Int_t AliAnalysisTaskSELc2V0bachelorTMVA::FindLcLabel(AliAODRecoCascadeHF* cascade, TClonesArray *mcArray) const
2148{
2149
2150 // finding the label of teh V0; inspired from AliAODRecoDecay::MatchToMC
2151
2152 //Printf(" FindLcLabel");
2153
2154 AliAODMCParticle *part=0;
2155 AliAODMCParticle *mother=0;
2156 AliAODMCParticle *grandmother=0;
2157 Int_t dgLabels = -1;
2158
2159 Int_t ndg = cascade->GetNDaughters();
2160 if (ndg != 2) {
2161 AliFatal(Form("IMPOSSIBLE!! There are %d daughters, but they should be 2!", ndg));
2162 }
2163
2164 // daughter 0 --> bachelor, daughter 1 --> V0
2165 AliAODTrack *trk = (AliAODTrack*)cascade->GetDaughter(0); // bachelor
2166 dgLabels = trk->GetLabel();
2167 if (dgLabels == -1 ) {
2168 //printf("daughter with negative label %d\n", dgLabels);
2169 return -1;
2170 }
2171 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2172 if (!part) {
2173 //printf("no MC particle\n");
2174 return -1;
2175 }
2176 Int_t labMotherBach = part->GetMother();
2177 if (labMotherBach == -1){
2178 return -1;
2179 }
2180 mother = (AliAODMCParticle*)mcArray->At(labMotherBach);
2181 if(!mother) {
2182 //printf("no MC mother particle\n");
2183 return -1;
2184 }
2185
2186 AliAODv0 *v0 = (AliAODv0*)cascade->GetDaughter(1); // V0
2187 dgLabels = FindV0Label(v0, mcArray);
2188 if (dgLabels == -1 ) {
2189 //printf("daughter with negative label (v0 was a fake) %d\n", dgLabels);
2190 return -1;
2191 }
2192 part = (AliAODMCParticle*)mcArray->At(TMath::Abs(dgLabels));
2193 if (!part) {
2194 //printf("no MC particle\n");
2195 return -1;
2196 }
2197 Int_t labMotherv0 = part->GetMother();
2198 if (labMotherv0 == -1){
2199 return -1;
2200 }
2201 mother = (AliAODMCParticle*)mcArray->At(labMotherv0);
2202 if(!mother) {
2203 //printf("no MC mother particle\n");
2204 return -1;
2205 }
2206 Int_t labGrandMotherv0 = mother->GetMother();
2207 if (labGrandMotherv0 == -1){
2208 return -1;
2209 }
2210 grandmother = (AliAODMCParticle*)mcArray->At(labGrandMotherv0);
2211 if(!grandmother) {
2212 //printf("no MC mother particle\n");
2213 return -1;
2214 }
2215
2216 //Printf("labMotherBach = %d, labMotherv0 = %d, labGrandMotherv0 = %d", labMotherBach, labMotherv0, labGrandMotherv0);
2217 if (labGrandMotherv0 == labMotherBach) return labMotherBach;
2218
2219 else return -1;
2220
2221}
0abe5247 2222
2223
2224
2225
2226
2227