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