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