Updates in LC->Kos+proton code (Annalisa)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCutsLctoV0.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2010, 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  * appear 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$ */
17
18 /////////////////////////////////////////////////////////////
19 //
20 // Class for cuts on AOD reconstructed Lc->V0+X
21 //
22 // Modified by A.De Caro - decaro@sa.infn.it
23 //
24 /////////////////////////////////////////////////////////////
25
26 #include <Riostream.h>
27
28 #include <TDatabasePDG.h>
29 #include <TMath.h>
30
31 #include "AliAnalysisManager.h"
32 #include "AliInputEventHandler.h"
33 #include "AliPIDResponse.h"
34 #include "AliRDHFCutsLctoV0.h"
35 #include "AliAODRecoCascadeHF.h"
36 #include "AliAODTrack.h"
37 #include "AliESDtrack.h"
38 #include "AliESDVertex.h"
39 #include "AliAODVertex.h"
40 #include "AliAODv0.h"
41 #include "AliESDv0.h"
42
43 using std::cout;
44 using std::endl;
45
46 ClassImp(AliRDHFCutsLctoV0)
47
48 //--------------------------------------------------------------------------
49 AliRDHFCutsLctoV0::AliRDHFCutsLctoV0(const char* name, Short_t /*v0channel*/) :
50 AliRDHFCuts(name),
51   fPidSelectionFlag(0),
52   fV0daughtersCuts(0),
53   fV0Type(0),
54   fHighPtCut(2.5),
55   fLowPtCut(1.0)
56 {
57   //
58   // Default Constructor
59   //
60
61   const Int_t nvars=17;
62   SetNVars(nvars);
63   TString varNames[nvars]={"Lc inv. mass if K0S [GeV/c2]",                   //  0
64                            "Lc inv. mass if Lambda [GeV/c2]",                //  1
65                            "K0S inv. mass [GeV/c2]",                         //  2
66                            "Lambda/LambdaBar inv. mass[GeV/c2]",             //  3
67                            "pT min bachelor track [GeV/c]",                  //  4
68                            "pT min V0-positive track [GeV/c]",               //  5
69                            "pT min V0-negative track [GeV/c]",               //  6
70                            "dca cascade (prong-to-prong) cut [cm]",          //  7
71                            "dca V0 (prong-to-prong) cut [number of sigmas]", //  8
72                            "V0 cosPA min wrt PV",                            //  9
73                            "d0 max bachelor wrt PV [cm]",                    // 10
74                            "d0 max V0 wrt PV [cm]",                          // 11
75                            "mass K0S veto [GeV/c2]",                         // 12
76                            "mass Lambda/LambdaBar veto [GeV/c2]",            // 13
77                            "mass Gamma veto [GeV/c2]",                       // 14
78                            "pT min V0 track [GeV/c]",                        // 15
79                            "V0 type"                                         // 16
80   };
81
82   Bool_t isUpperCut[nvars]={kTRUE,  //  0
83                             kTRUE,  //  1
84                             kTRUE,  //  2
85                             kTRUE,  //  3
86                             kFALSE, //  4
87                             kFALSE, //  5
88                             kFALSE, //  6
89                             kTRUE,  //  7
90                             kTRUE,  //  8
91                             kFALSE, //  9
92                             kTRUE,  // 10
93                             kTRUE,  // 11
94                             kFALSE, // 12
95                             kFALSE, // 13
96                             kFALSE, // 14
97                             kFALSE, // 15
98                             kFALSE  // 16
99   };
100   SetVarNames(nvars,varNames,isUpperCut);
101   Bool_t forOpt[nvars]={kFALSE, //  0
102                         kFALSE, //  1
103                         kTRUE,  //  2
104                         kTRUE,  //  3
105                         kTRUE,  //  4
106                         kTRUE,  //  5
107                         kTRUE,  //  6
108                         kTRUE,  //  7
109                         kTRUE,  //  8
110                         kTRUE,  //  9
111                         kTRUE,  // 10
112                         kTRUE,  // 11
113                         kTRUE,  // 12
114                         kTRUE,  // 13
115                         kTRUE,  // 14
116                         kTRUE,  // 15
117                         kFALSE  // 16
118   };
119   SetVarsForOpt(nvars,forOpt);
120
121   Float_t limits[2]={0,999999999.};
122   SetPtBins(2,limits);
123
124   /*
125     switch (v0channel) {
126     case 0:
127     fV0channel = 0x0001;
128     break;
129     case 1:
130     fV0channel = 0x0002;
131     break;
132     case 2:
133     fV0channel = 0x0004;
134     break;
135     }
136   */
137
138 }
139 //--------------------------------------------------------------------------
140 AliRDHFCutsLctoV0::AliRDHFCutsLctoV0(const AliRDHFCutsLctoV0 &source) :
141   AliRDHFCuts(source),
142   fPidSelectionFlag(source.fPidSelectionFlag),
143   fV0daughtersCuts(0),
144   fV0Type(source.fV0Type),
145   fHighPtCut(source.fHighPtCut),
146   fLowPtCut(source.fLowPtCut)
147   /*fV0channel(source.fV0channel)*/
148 {
149   //
150   // Copy constructor
151   //
152
153   if (source.fV0daughtersCuts) AddTrackCutsV0daughters(source.fV0daughtersCuts);
154   else fV0daughtersCuts = new AliESDtrackCuts();
155
156 }
157 //--------------------------------------------------------------------------
158 AliRDHFCutsLctoV0 &AliRDHFCutsLctoV0::operator=(const AliRDHFCutsLctoV0 &source)
159 {
160   //
161   // assignment operator
162   //
163
164   if (this != &source) {
165
166     AliRDHFCuts::operator=(source);
167     fPidSelectionFlag = source.fPidSelectionFlag;
168
169     delete fV0daughtersCuts;
170     fV0daughtersCuts = new AliESDtrackCuts(*(source.fV0daughtersCuts));
171
172     fV0Type  = source.fV0Type;
173
174     fHighPtCut = source.fHighPtCut;
175     fLowPtCut = source.fLowPtCut;
176
177   }
178
179   return *this;
180 }
181
182
183 //---------------------------------------------------------------------------
184 AliRDHFCutsLctoV0::~AliRDHFCutsLctoV0() {
185   //
186   //  Default Destructor
187   //
188
189   if (fV0daughtersCuts) {
190     delete fV0daughtersCuts;
191     fV0daughtersCuts=0;
192   }
193
194 }
195
196 //---------------------------------------------------------------------------
197 void AliRDHFCutsLctoV0::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
198   //
199   // Fills in vars the values of the variables
200   //
201
202   if (pdgdaughters[0]==-9999) return; // dummy
203
204   if (nvars!=fnVarsForOpt) {
205     printf("AliRDHFCutsLctoV0::GetCutVarsForOpt: wrong number of variables\n");
206     return;
207   }
208
209   Double_t mLcPDG  = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
210   Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
211   Double_t mLPDG   = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
212
213   AliAODRecoCascadeHF *dd = (AliAODRecoCascadeHF*)d;
214
215   // Get the v0 and all daughter tracks
216   AliAODTrack *bachelorTrack = (AliAODTrack*)dd->GetBachelor();
217   AliAODv0 *v0 = (AliAODv0*)dd->Getv0();
218   AliAODTrack *v0positiveTrack = (AliAODTrack*)dd->Getv0PositiveTrack();
219   AliAODTrack *v0negativeTrack = (AliAODTrack*)dd->Getv0NegativeTrack();
220
221   Int_t iter=-1;
222   // cut on cascade mass, if K0S + p
223   if (fVarsForOpt[0]) {
224     iter++;
225     vars[iter]=TMath::Abs(dd->InvMassLctoK0sP()-mLcPDG);
226   }
227   // cut on cascade mass, if Lambda/LambdaBar + pi
228   if (fVarsForOpt[1]) {
229     iter++;
230     vars[iter]=TMath::Abs(dd->InvMassLctoLambdaPi()-mLcPDG);
231   }
232
233   // cut on V0 mass if K0S
234   if (fVarsForOpt[2]) {
235     iter++;
236     vars[iter]=TMath::Abs(v0->MassK0Short()-mk0sPDG);
237   }
238
239   // cut on V0 mass if Lambda/LambdaBar
240   if (fVarsForOpt[3]) {
241
242     if (bachelorTrack->Charge()==1) {
243       iter++;
244       vars[iter]=TMath::Abs(v0->MassLambda()-mLPDG);
245     } else if (bachelorTrack->Charge()==-1) {
246       iter++;
247       vars[iter]=TMath::Abs(v0->MassAntiLambda()-mLPDG);
248     }
249
250   }
251
252   // cut bachelor min pt
253   if (fVarsForOpt[4]) {
254     iter++;
255     vars[iter]=bachelorTrack->Pt();
256   }
257
258   // cut on V0-positive min pt
259   if (fVarsForOpt[5]) {
260     iter++;
261     vars[iter]=v0positiveTrack->Pt();
262   }
263
264   // cut on V0-negative min pt
265   if (fVarsForOpt[6]) {
266     iter++;
267     vars[iter]=v0negativeTrack->Pt();
268   }
269
270   // cut on cascade dca (prong-to-prong)
271   if (fVarsForOpt[7]) {
272     iter++;
273     vars[iter]=dd->GetDCA(); // prong-to-prong cascade DCA
274   }
275
276   // cut on V0 dca (prong-to-prong)
277   if (fVarsForOpt[8]) {
278     iter++;
279     vars[iter]=v0->GetDCA(); // prong-to-prong V0 DCA
280   }
281
282   // cut on V0 cosPA wrt PV
283   if (fVarsForOpt[9]) {
284     iter++;
285     vars[iter]=dd->CosV0PointingAngle(); // cosine of V0 pointing angle wrt primary vertex
286   }
287
288   // cut on bachelor transverse impact parameter wrt PV
289   if (fVarsForOpt[10]) {
290     iter++;
291     vars[iter]=dd->Getd0Prong(0); // bachelor transverse impact parameter wrt primary vertex
292   }
293
294   // cut on V0 transverse impact parameter wrt PV
295   if (fVarsForOpt[11]) {
296     iter++;
297     vars[iter]=dd->Getd0Prong(1); // V0 transverse impact parameter wrt primary vertex
298   }
299
300   // cut on K0S invariant mass veto
301   if (fVarsForOpt[12]) {
302     iter++;
303     vars[iter]=TMath::Abs(v0->MassK0Short()-mk0sPDG); // K0S invariant mass veto
304   }
305
306   // cut on Lambda/LambdaBar invariant mass veto
307   if (fVarsForOpt[13]) {
308
309     if (bachelorTrack->Charge()==1) {
310       iter++;
311       vars[iter]=TMath::Abs(v0->MassLambda()-mLPDG);
312     } else if (bachelorTrack->Charge()==-1) {
313       iter++;
314       vars[iter]=TMath::Abs(v0->MassAntiLambda()-mLPDG);
315     }
316
317   }
318
319   // cut on gamma invariant mass veto
320   if (fVarsForOpt[14]) {
321     iter++;
322     vars[iter]= v0->InvMass2Prongs(0,1,11,11); // gamma invariant mass veto
323   }
324
325
326   // cut on V0 pT min
327   if (fVarsForOpt[15]) {
328     iter++;
329     vars[iter]= v0->Pt(); // V0 pT min
330   }
331
332
333   return;
334 }
335 //---------------------------------------------------------------------------
336 Int_t AliRDHFCutsLctoV0::IsSelected(TObject* obj,Int_t selectionLevel) {
337   //
338   // Apply selection
339   //
340
341   if (!fCutsRD) {
342     AliFatal("Cut matrice not inizialized. Exit...");
343     return 0;
344   }
345
346   AliAODRecoCascadeHF* d = (AliAODRecoCascadeHF*)obj;
347   if (!d) {
348     AliDebug(2,"AliAODRecoCascadeHF null");
349     return 0;
350   }
351
352   if (!d->GetSecondaryVtx()) {
353     AliDebug(2,"No secondary vertex for cascade");
354     return 0;
355   }
356
357   if (d->GetNDaughters()!=2) {
358     AliDebug(2,Form("No 2 daughters for current cascade (nDaughters=%d)",d->GetNDaughters()));
359     return 0;
360   }
361
362   AliAODv0 * v0 = dynamic_cast<AliAODv0*>(d->Getv0());
363   AliAODTrack * bachelorTrack = dynamic_cast<AliAODTrack*>(d->GetBachelor());
364   if (!v0 || !bachelorTrack) {
365     AliDebug(2,"No V0 or no bachelor for current cascade");
366     return 0;
367   }
368
369   if (bachelorTrack->GetID()<0) {
370     AliDebug(2,Form("Bachelor has negative ID %d",bachelorTrack->GetID()));
371     return 0;
372   }
373
374   if (!v0->GetSecondaryVtx()) {
375     AliDebug(2,"No secondary vertex for V0 by cascade");
376     return 0;
377   }
378
379   if (v0->GetNDaughters()!=2) {
380     AliDebug(2,Form("No 2 daughters for V0 of current cascade (onTheFly=%d, nDaughters=%d)",v0->GetOnFlyStatus(),v0->GetNDaughters()));
381     return 0;
382   }
383
384
385   // Get the V0 daughter tracks
386   AliAODTrack *v0positiveTrack = dynamic_cast<AliAODTrack*>(d->Getv0PositiveTrack());
387   AliAODTrack *v0negativeTrack = dynamic_cast<AliAODTrack*>(d->Getv0NegativeTrack());
388   if (!v0positiveTrack || !v0negativeTrack ) {
389     AliDebug(2,"No V0 daughters' objects");
390     return 0;
391   }
392
393   if (v0positiveTrack->GetID()<0 || v0negativeTrack->GetID()<0) {
394     AliDebug(2,Form("At least one of V0 daughters has negative ID %d %d",v0positiveTrack->GetID(),v0negativeTrack->GetID()));
395     return 0;
396   }
397
398   //if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
399   if ( fUseTrackSelectionWithFilterBits && !(bachelorTrack->TestFilterMask(BIT(4))) ) {
400     AliDebug(2,"Check on the bachelor FilterBit: no BIT(4). Candidate rejected.");
401     return 0;
402   }
403
404   Int_t returnvalueTrack = 7;
405
406   // selection on daughter tracks
407   if (selectionLevel==AliRDHFCuts::kAll ||
408       selectionLevel==AliRDHFCuts::kTracks) {
409
410     if (!AreLctoV0DaughtersSelected(d)) return 0;
411
412   }
413
414   Bool_t okLck0sp=kTRUE, okLcLpi=kTRUE, okLcLBarpi=kTRUE;
415
416   // selection on candidate
417   if (selectionLevel==AliRDHFCuts::kAll ||
418       selectionLevel==AliRDHFCuts::kCandidate) {
419
420     Double_t pt = d->Pt();
421     Int_t ptbin = PtBin(pt);
422
423     Double_t mLcPDG  = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
424     Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
425     Double_t mLPDG   = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
426
427     // K0S + p
428     Double_t mk0s    = v0->MassK0Short();
429     Double_t mLck0sp = d->InvMassLctoK0sP();
430
431     // Lambda + pi
432     Double_t mlambda  = v0->MassLambda();
433     Double_t malambda = v0->MassAntiLambda();
434     Double_t mLcLpi   = d->InvMassLctoLambdaPi();
435
436     Bool_t okK0spipi=kTRUE, okLppi=kTRUE, okLBarpip=kTRUE;
437     Bool_t isNotK0S = kTRUE, isNotLambda = kTRUE, isNotLambdaBar = kTRUE, isNotGamma = kTRUE;
438
439     // cut on Lc mass with K0S+p hypothesis
440     if (TMath::Abs(mLck0sp-mLcPDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) {
441       okLck0sp = kFALSE;
442       AliDebug(4,Form(" cascade mass is %2.2e and does not correspond to Lambda_c into K0S+p cut",mLck0sp));
443     }
444
445     // cuts on the V0 mass: K0S case
446     if (TMath::Abs(mk0s-mk0sPDG) > fCutsRD[GetGlobalIndex(2,ptbin)]) {
447       okK0spipi = kFALSE;
448       AliDebug(4,Form(" V0 mass is %2.2e and does not correspond to K0S cut",mk0s));
449     }
450
451     // cut on Lc mass with Lambda+pi hypothesis
452     if (TMath::Abs(mLcLpi-mLcPDG) > fCutsRD[GetGlobalIndex(1,ptbin)]) {
453       okLcLpi = kFALSE;
454       okLcLBarpi = kFALSE;
455       AliDebug(4,Form(" cascade mass is %2.2e and does not correspond to Lambda_c into Lambda+pi cut",mLcLpi));
456     }
457
458     // cuts on the V0 mass: Lambda/LambdaBar case
459     //if ( !(bachelorTrack->Charge()==+1 && TMath::Abs(mlambda-mLPDG) <= fCutsRD[GetGlobalIndex(3,ptbin)] ) ) {
460     if ( TMath::Abs(mlambda-mLPDG) > fCutsRD[GetGlobalIndex(3,ptbin)] ) {
461       okLppi = kFALSE;
462       AliDebug(4,Form(" V0 mass is %2.2e and does not correspond to LambdaBar cut",mlambda));
463     }
464
465     //if ( !(bachelorTrack->Charge()==-1 && TMath::Abs(malambda-mLPDG) <= fCutsRD[GetGlobalIndex(3,ptbin)] ) ) {
466     if ( TMath::Abs(malambda-mLPDG) > fCutsRD[GetGlobalIndex(3,ptbin)] ) {
467       okLBarpip = kFALSE;
468       AliDebug(4,Form(" V0 mass is %2.2e and does not correspond to LambdaBar cut",malambda));
469     }
470
471     // cut on K0S invariant mass veto
472     if (TMath::Abs(v0->MassK0Short()-mk0sPDG) < fCutsRD[GetGlobalIndex(12,ptbin)]) { // K0S invariant mass veto
473       AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
474       //return 0;
475       isNotK0S=kFALSE;
476     }
477
478     // cut on Lambda/LambdaBar invariant mass veto
479     if (TMath::Abs(v0->MassLambda()-mLPDG) < fCutsRD[GetGlobalIndex(13,ptbin)]) { // Lambda invariant mass veto
480       AliDebug(4,Form(" veto on Lambda invariant mass doesn't pass the cut"));
481       isNotLambda=kFALSE;
482       //return 0;
483     }
484     if (TMath::Abs(v0->MassAntiLambda()-mLPDG) < fCutsRD[GetGlobalIndex(13,ptbin)] ) { // LambdaBar invariant mass veto
485       AliDebug(4,Form(" veto on LambdaBar invariant mass doesn't pass the cut"));
486       isNotLambdaBar=kFALSE;
487       //return 0;
488     }
489
490     // cut on gamma invariant mass veto
491     if (v0->InvMass2Prongs(0,1,11,11) < fCutsRD[GetGlobalIndex(14,ptbin)]) { // K0S invariant mass veto
492       AliDebug(4,Form(" veto on gamma invariant mass doesn't pass the cut"));
493       isNotGamma=kFALSE;
494       //return 0;
495     }
496
497     okLck0sp   = okLck0sp   && okK0spipi && isNotLambda && isNotLambdaBar && isNotGamma;
498     okLcLpi    = okLcLpi    && okLppi    && isNotK0S    && isNotLambdaBar && isNotGamma;
499     okLcLBarpi = okLcLBarpi && okLBarpip && isNotK0S    && isNotLambda    && isNotGamma;
500
501     if (!okLck0sp && !okLcLpi && !okLcLBarpi) return 0;
502
503     // cuts on the minimum pt of the tracks
504     if (TMath::Abs(bachelorTrack->Pt()) < fCutsRD[GetGlobalIndex(4,ptbin)]) {
505       AliDebug(4,Form(" bachelor track Pt=%2.2e > %2.2e",bachelorTrack->Pt(),fCutsRD[GetGlobalIndex(4,ptbin)]));
506       return 0;
507     }
508     if (TMath::Abs(v0positiveTrack->Pt()) < fCutsRD[GetGlobalIndex(5,ptbin)]) {
509       AliDebug(4,Form(" V0-positive track Pt=%2.2e > %2.2e",v0positiveTrack->Pt(),fCutsRD[GetGlobalIndex(5,ptbin)]));
510       return 0;
511     }
512     if (TMath::Abs(v0negativeTrack->Pt()) < fCutsRD[GetGlobalIndex(6,ptbin)]) {
513       AliDebug(4,Form(" V0-negative track Pt=%2.2e > %2.2e",v0negativeTrack->Pt(),fCutsRD[GetGlobalIndex(6,ptbin)]));
514       return 0;
515     }
516
517     // cut on cascade dca (prong-to-prong)
518     if ( TMath::Abs(d->GetDCA()) > fCutsRD[GetGlobalIndex(7,ptbin)] ) { // prong-to-prong cascade DCA
519       AliDebug(4,Form(" cascade tracks DCA don't pass the cut"));
520       return 0;
521     }
522
523     // cut on V0 dca (prong-to-prong)
524     if ( TMath::Abs(v0->GetDCA()) > fCutsRD[GetGlobalIndex(8,ptbin)] ) { // prong-to-prong V0 DCA
525       AliDebug(4,Form(" V0 DCA don't pass the cut"));
526       return 0;
527     }
528
529     // cut on V0 cosine of pointing angle wrt PV
530     if (d->CosV0PointingAngle() < fCutsRD[GetGlobalIndex(9,ptbin)]) { // cosine of V0 pointing angle wrt primary vertex
531       AliDebug(4,Form(" V0 cosine of pointing angle doesn't pass the cut"));
532       return 0;
533     }
534
535     // cut on bachelor transverse impact parameter wrt PV
536     if (TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(10,ptbin)]) { // bachelor transverse impact parameter wrt PV
537       AliDebug(4,Form(" bachelor transverse impact parameter doesn't pass the cut"));
538       return 0;
539     }
540
541     // cut on V0 transverse impact parameter wrt PV
542     if (TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(11,ptbin)]) { // V0 transverse impact parameter wrt PV
543       AliDebug(4,Form(" V0 transverse impact parameter doesn't pass the cut"));
544       return 0;
545     }
546
547     // cut on V0 pT min
548     if (v0->Pt() < fCutsRD[GetGlobalIndex(15,ptbin)]) { // V0 pT min
549       AliDebug(4,Form(" V0 track Pt=%2.2e > %2.2e",v0->Pt(),fCutsRD[GetGlobalIndex(15,ptbin)]));
550       return 0;
551     }
552
553   }
554
555   Int_t returnvalue = okLck0sp+2*okLcLBarpi+4*okLcLpi;
556   /*
557     retvalue case
558     1  Lc->K0S + p
559     2  Lc->LambdaBar + pi
560     3  Lc->K0S + p AND Lc->LambdaBar + pi
561     4  Lc->Lambda + pi
562     5  Lc->K0S + p AND Lc->Lambda + pi
563     6  Lc->LambdaBar + pi AND Lc->Lambda + pi
564     7  Lc->K0S + p AND Lc->LambdaBar + pi AND Lc->Lambda + pi
565   */
566
567   Int_t returnvaluePID = 7;
568
569   // selection on candidate
570   if (selectionLevel==AliRDHFCuts::kAll ||
571       selectionLevel==AliRDHFCuts::kCandidate ||
572       selectionLevel==AliRDHFCuts::kPID )
573     returnvaluePID = IsSelectedPID(d);
574
575   //if (fUsePID && returnvaluePID==0) return 0;
576
577   Int_t returnvalueTot = 0;
578   if ( fUsePID )
579     returnvalueTot = CombineCuts(returnvalueTrack,returnvalue,returnvaluePID);
580   else
581     returnvalueTot = CombineCuts(returnvalueTrack,returnvalue,7);
582
583   return returnvalueTot;
584
585 }
586 //---------------------------------------------------------------------------
587 Int_t AliRDHFCutsLctoV0::IsSelectedPID(AliAODRecoDecayHF* obj) {
588
589   // fPidHF -> PID object for bachelor
590
591   if (!fUsePID || !obj) {
592     AliDebug(2,"PID selection inactive. Candidate accepted.");
593     return 7; // all hypothesis are valid
594   }
595
596   AliAODRecoCascadeHF *objD = (AliAODRecoCascadeHF*)obj;
597
598   AliAODTrack *bachelor = (AliAODTrack*)objD->GetBachelor();
599   AliAODTrack *v0Pos = (AliAODTrack*)objD->Getv0PositiveTrack();
600   AliAODTrack *v0Neg = (AliAODTrack*)objD->Getv0NegativeTrack();
601
602   if (!bachelor || !v0Pos || !v0Neg) return 0;
603
604   Bool_t okLcK0Sp = kTRUE; // K0S case
605   Bool_t okLcLambdaBarPi = kTRUE; // LambdaBar case
606   Bool_t okLcLambdaPi = kTRUE; // Lambda case
607
608   CheckPID(bachelor,v0Neg,v0Pos,okLcK0Sp,okLcLambdaBarPi,okLcLambdaPi);
609
610   Int_t returnvalue = okLcK0Sp+2*okLcLambdaBarPi+4*okLcLambdaPi;
611
612   return returnvalue;
613 }
614 //-----------------------
615 void AliRDHFCutsLctoV0::CheckPID(AliAODTrack *bachelor,
616                                  AliAODTrack * /*v0Neg*/, AliAODTrack * /*v0Pos*/,
617                                  Bool_t &isBachelorID1, Bool_t &isBachelorID2, Bool_t &isBachelorID4) {
618   // identification strategy
619
620   Int_t idxIDbyTOF = -1;
621   Int_t idxIDbyTPC = -1;
622
623   Int_t tpcID = -1;
624   Int_t tofID = -1;
625   Double_t nTPCsigmas = -999;
626   Double_t nTOFsigmas = -999;
627
628   Bool_t trackIDByTOF = -1;
629   Bool_t trackIDByTPC = -1;
630
631   Bool_t dummy = kFALSE;
632
633   switch (fPidSelectionFlag) {
634
635   case 0:
636
637     // identify bachelor
638     idxIDbyTOF = fPidHF->ApplyPidTOFRaw(bachelor,4);
639     idxIDbyTPC = fPidHF->ApplyPidTPCRaw(bachelor,4);
640     isBachelorID1 = (idxIDbyTOF==4) && (idxIDbyTPC==4); // K0S case
641
642     idxIDbyTOF = fPidHF->ApplyPidTOFRaw(bachelor,2);
643     idxIDbyTPC = fPidHF->ApplyPidTPCRaw(bachelor,2);
644     isBachelorID2 = (idxIDbyTOF==2) && (idxIDbyTPC==2); // LambdaBar case
645
646     isBachelorID4 = isBachelorID2; // Lambda case
647
648     break;
649
650   case 1:
651
652     // identify bachelor
653     idxIDbyTOF = fPidHF->ApplyPidTOFRaw(bachelor,4);
654     idxIDbyTPC = fPidHF->ApplyPidTPCRaw(bachelor,4);
655     dummy = ( !(fPidHF->CheckTOFPIDStatus(bachelor)) && (idxIDbyTPC==4) &&
656               fPidHF->IsExcluded(bachelor,2,2.,"TPC") && fPidHF->IsExcluded(bachelor,3,2.,"TPC") ); // K0S case
657     isBachelorID1 = ( (idxIDbyTOF==4) || dummy );
658
659     idxIDbyTOF = fPidHF->ApplyPidTOFRaw(bachelor,2);
660     idxIDbyTPC = fPidHF->ApplyPidTPCRaw(bachelor,2);
661     dummy = ( !(fPidHF->CheckTOFPIDStatus(bachelor)) && (idxIDbyTPC==2) &&
662               fPidHF->IsExcluded(bachelor,3,2.,"TPC") && fPidHF->IsExcluded(bachelor,4,2.,"TPC") ); // LambdaBar case
663     isBachelorID2 = ( (idxIDbyTOF==2) || dummy );
664
665     isBachelorID4 = isBachelorID2; // Lambda case
666
667     break;
668
669   case 2:
670
671     // identify bachelor
672     nTOFsigmas = -999;
673     tofID = fPidHF->GetnSigmaTOF(bachelor,4,nTOFsigmas);
674     nTPCsigmas = -999;
675     tpcID = fPidHF->GetnSigmaTPC(bachelor,4,nTPCsigmas);
676     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
677                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
678     trackIDByTPC = ( (tpcID==1) && ( (bachelor->P()<fLowPtCut && TMath::Abs(nTPCsigmas)<2) ||
679                                      (bachelor->P()>=fLowPtCut && TMath::Abs(nTPCsigmas)<3) ) );
680     isBachelorID1 = (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && trackIDByTPC && trackIDByTOF); // K0S case
681
682     nTOFsigmas = -999;
683     tofID = fPidHF->GetnSigmaTOF(bachelor,2,nTOFsigmas);
684     nTPCsigmas = -999;
685     tpcID = fPidHF->GetnSigmaTPC(bachelor,2,nTPCsigmas);
686     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
687                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
688     trackIDByTPC = ( (tpcID==1) && ( (bachelor->P()<fLowPtCut && TMath::Abs(nTPCsigmas)<2) ||
689                                      (bachelor->P()>=fLowPtCut && TMath::Abs(nTPCsigmas)<3) ) );
690     isBachelorID2 = (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && trackIDByTPC && trackIDByTOF); // LambdaBar case
691
692     isBachelorID4 = isBachelorID2; // Lambda case
693
694     break;
695
696   case 3:
697
698     // identify bachelor
699     nTOFsigmas = -999;
700     tofID = fPidHF->GetnSigmaTOF(bachelor,4,nTOFsigmas);
701     nTPCsigmas = -999;
702     tpcID = fPidHF->GetnSigmaTPC(bachelor,4,nTPCsigmas);
703     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
704                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
705     trackIDByTPC = ( (tpcID==1) && ( (bachelor->P()<fLowPtCut && TMath::Abs(nTPCsigmas)<2) ||
706                                      (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTPCsigmas)<3) ||
707                                      (bachelor->P()>=fHighPtCut && nTPCsigmas>-3 && nTPCsigmas<2) ) );
708     isBachelorID1 = (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && trackIDByTPC && trackIDByTOF); // K0S case
709
710     nTOFsigmas = -999;
711     tofID = fPidHF->GetnSigmaTOF(bachelor,2,nTOFsigmas);
712     nTPCsigmas = -999;
713     tpcID = fPidHF->GetnSigmaTPC(bachelor,2,nTPCsigmas);
714     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
715                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
716     trackIDByTPC = ( (tpcID==1) && ( (bachelor->P()<fLowPtCut && TMath::Abs(nTPCsigmas)<2) ||
717                                      (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTPCsigmas)<3) ||
718                                      (bachelor->P()>=fHighPtCut && nTPCsigmas>-3 && nTPCsigmas<2) ) );
719     isBachelorID2 = (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && trackIDByTPC && trackIDByTOF); // LambdaBar case
720
721     isBachelorID4 = isBachelorID2; // Lambda case
722
723     break;
724
725   case 4:
726
727     // identify bachelor
728     nTOFsigmas = -999;
729     tofID = fPidHF->GetnSigmaTOF(bachelor,4,nTOFsigmas);
730     nTPCsigmas = -999;
731     tpcID = fPidHF->GetnSigmaTPC(bachelor,4,nTPCsigmas);
732     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
733                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
734     trackIDByTPC = ( (tpcID==1) && (TMath::Abs(nTPCsigmas)<2) );
735     isBachelorID1 = ( (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && trackIDByTOF) ); // K0S case
736
737     nTOFsigmas = -999;
738     tofID = fPidHF->GetnSigmaTOF(bachelor,2,nTOFsigmas);
739     nTPCsigmas = -999;
740     tpcID = fPidHF->GetnSigmaTPC(bachelor,2,nTPCsigmas);
741     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
742                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
743     trackIDByTPC = ( (tpcID==1) && (TMath::Abs(nTPCsigmas)<2) );
744     isBachelorID2 = ( (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && trackIDByTOF) ); // LambdaBar case
745
746     isBachelorID4 = isBachelorID2; // Lambda case
747
748     break;
749
750   case 5:
751
752     // identify bachelor
753     nTOFsigmas = -999;
754     tofID = fPidHF->GetnSigmaTOF(bachelor,4,nTOFsigmas);
755     nTPCsigmas = -999;
756     tpcID = fPidHF->GetnSigmaTPC(bachelor,4,nTPCsigmas);
757     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
758                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
759     trackIDByTPC = ( (tpcID==1) && ( (bachelor->P()<fLowPtCut && TMath::Abs(nTPCsigmas)<2) || (bachelor->P()>=fHighPtCut && !trackIDByTOF && TMath::Abs(nTPCsigmas)<3) ) );
760     isBachelorID1 = (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && trackIDByTOF) || (bachelor->P()>=fHighPtCut && (trackIDByTOF || trackIDByTPC) ); // K0S case
761
762     nTOFsigmas = -999;
763     tofID = fPidHF->GetnSigmaTOF(bachelor,2,nTOFsigmas);
764     nTPCsigmas = -999;
765     tpcID = fPidHF->GetnSigmaTPC(bachelor,2,nTPCsigmas);
766     trackIDByTOF = ( (tofID==1) && ( (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && TMath::Abs(nTOFsigmas)<3) ||
767                                      (bachelor->P()>=fHighPtCut && nTOFsigmas>-2 && nTOFsigmas<3) ) );
768     trackIDByTPC = ( (tpcID==1) && ( (bachelor->P()<fLowPtCut && TMath::Abs(nTPCsigmas)<2) || (bachelor->P()>=fHighPtCut && !trackIDByTOF && TMath::Abs(nTPCsigmas)<3) ) );
769     isBachelorID2 = (bachelor->P()<fLowPtCut && trackIDByTPC) || (bachelor->P()>=fLowPtCut && bachelor->P()<fHighPtCut && trackIDByTOF) || (bachelor->P()>=fHighPtCut && (trackIDByTOF || trackIDByTPC) ); // LambdaBar case
770
771     isBachelorID4 = isBachelorID2; // Lambda case
772
773     break;
774   }
775
776 }
777 //----------------
778 Int_t AliRDHFCutsLctoV0::CombineCuts(Int_t returnvalueTrack, Int_t returnvalue, Int_t returnvaluePID) const {
779   //
780   // combine track selection, topological cuts and PID
781   //
782
783   Int_t returnvalueTot=returnvalueTrack&returnvalue;
784   returnvalueTot=returnvalueTot&returnvaluePID;
785
786   return returnvalueTot;
787 }
788
789 //----------------------------------
790 Int_t AliRDHFCutsLctoV0::IsSelectedSingleCut(TObject* obj, Int_t selectionLevel, Int_t cutIndex) {
791   //
792   // Apply selection on single cut
793   //
794
795   if (!fCutsRD) {
796     AliDebug(2,"Cut matrice not inizialized. Exit...");
797     return 0;
798   }
799
800   AliAODRecoCascadeHF* d = (AliAODRecoCascadeHF*)obj;
801   if (!d) {
802     AliDebug(2,"AliAODRecoCascadeHF null");
803     return 0;
804   }
805
806   if (!d->GetSecondaryVtx()) {
807     AliDebug(2,"No secondary vertex for cascade");
808     return 0;
809   }
810
811   if (d->GetNDaughters()!=2) {
812     AliDebug(2,Form("No 2 daughters for current cascade (nDaughters=%d)",d->GetNDaughters()));
813     return 0;
814   }
815
816   AliAODv0 * v0 = dynamic_cast<AliAODv0*>(d->Getv0());
817   AliAODTrack * bachelorTrack = dynamic_cast<AliAODTrack*>(d->GetBachelor());
818   if (!v0 || !bachelorTrack) {
819     AliDebug(2,"No V0 or no bachelor for current cascade");
820     return 0;
821   }
822
823   if (bachelorTrack->GetID()<0) {
824     AliDebug(2,Form("Bachelor has negative ID %d",bachelorTrack->GetID()));
825     return 0;
826   }
827
828   if (!v0->GetSecondaryVtx()) {
829     AliDebug(2,"No secondary vertex for V0 by cascade");
830     return 0;
831   }
832
833   if (v0->GetNDaughters()!=2) {
834     AliDebug(2,Form("No 2 daughters for V0 of current cascade (onTheFly=%d, nDaughters=%d)",v0->GetOnFlyStatus(),v0->GetNDaughters()));
835     return 0;
836   }
837
838
839   // Get the V0 daughter tracks
840   AliAODTrack *v0positiveTrack = dynamic_cast<AliAODTrack*>(d->Getv0PositiveTrack());
841   AliAODTrack *v0negativeTrack = dynamic_cast<AliAODTrack*>(d->Getv0NegativeTrack());
842   if (!v0positiveTrack || !v0negativeTrack ) {
843     AliDebug(2,"No V0 daughters' objects");
844     return 0;
845   }
846
847   if (v0positiveTrack->GetID()<0 || v0negativeTrack->GetID()<0) {
848     AliDebug(2,Form("At least one of V0 daughters has negative ID %d %d",v0positiveTrack->GetID(),v0negativeTrack->GetID()));
849     return 0;
850   }
851
852   //if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
853   if ( fUseTrackSelectionWithFilterBits && !(bachelorTrack->TestFilterMask(BIT(4))) ) {
854     AliDebug(2,"Check on the bachelor FilterBit: no BIT(4). Candidate rejected.");
855     return 0;
856   }
857
858
859   // selection on daughter tracks
860   if (selectionLevel==AliRDHFCuts::kAll ||
861       selectionLevel==AliRDHFCuts::kTracks) {
862
863     if (!AreLctoV0DaughtersSelected(d)) return 0;
864
865   }
866
867   Bool_t okLck0sp=kFALSE, okLcLpi=kFALSE, okLcLBarpi=kFALSE;
868
869   // selection on candidate
870   if (selectionLevel==AliRDHFCuts::kAll ||
871       selectionLevel==AliRDHFCuts::kCandidate) {
872
873     Double_t pt = d->Pt();
874     Int_t ptbin = PtBin(pt);
875
876     Double_t mLcPDG  = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
877     Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
878     Double_t mLPDG   = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
879
880     // K0S + p
881     Double_t mk0s    = v0->MassK0Short();
882     Double_t mLck0sp = d->InvMassLctoK0sP();
883
884     // Lambda + pi
885     Double_t mlambda  = v0->MassLambda();
886     Double_t malambda = v0->MassAntiLambda();
887     Double_t mLcLpi   = d->InvMassLctoLambdaPi();
888
889     switch (cutIndex) {
890     case 0:
891       // cut on Lc mass with K0S+p hypothesis
892       okLck0sp   = TMath::Abs(mLck0sp-mLcPDG)<=fCutsRD[GetGlobalIndex(0,ptbin)];
893       okLcLpi    = kFALSE;
894       okLcLBarpi = kFALSE;
895       break;
896     case 1:
897       // cut on Lc mass with Lambda+pi hypothesis
898       okLck0sp   = kFALSE;
899       okLcLpi    = TMath::Abs(mLcLpi-mLcPDG)<=fCutsRD[GetGlobalIndex(1,ptbin)];
900       okLcLBarpi = okLcLpi;
901       break;
902     case 2:
903       // cuts on the V0 mass: K0S case
904       okLck0sp   = TMath::Abs(mk0s-mk0sPDG)<=fCutsRD[GetGlobalIndex(2,ptbin)];
905       okLcLpi    = kFALSE;
906       okLcLBarpi = kFALSE;
907       break;
908     case 3:
909       // cuts on the V0 mass: Lambda/LambdaBar case
910       okLck0sp   = kFALSE;
911       okLcLpi    = TMath::Abs(mlambda-mLPDG)<=fCutsRD[GetGlobalIndex(3,ptbin)];
912       //okLcLpi    = okLcLpi && (bachelorTrack->Charge()==+1);
913       okLcLBarpi = TMath::Abs(malambda-mLPDG)<=fCutsRD[GetGlobalIndex(3,ptbin)];
914       //okLcLBarpi = okLcLBarpi && (bachelorTrack->Charge()==-1);
915       break;
916     case 4:
917       // cuts on the minimum pt of bachelor
918       okLck0sp   = TMath::Abs(bachelorTrack->Pt())>=fCutsRD[GetGlobalIndex(4,ptbin)];
919       okLcLpi    = okLck0sp;
920       okLcLBarpi = okLck0sp;
921       break;
922     case 5:
923       // cuts on the minimum pt of positive V0daughter
924       okLck0sp   = TMath::Abs(v0positiveTrack->Pt())>=fCutsRD[GetGlobalIndex(5,ptbin)];
925       okLcLpi    = okLck0sp;
926       okLcLBarpi = okLck0sp;
927       break;
928     case 6:
929       // cuts on the minimum pt of negative V0daughter
930       okLck0sp   = TMath::Abs(v0negativeTrack->Pt())>=fCutsRD[GetGlobalIndex(6,ptbin)];
931       okLcLpi    = okLck0sp;
932       okLcLBarpi = okLck0sp;
933       break;
934     case 7:
935       // cut on cascade dca
936       okLck0sp   = TMath::Abs(d->GetDCA())<=fCutsRD[GetGlobalIndex(7,ptbin)];
937       okLcLpi    = okLck0sp;
938       okLcLBarpi = okLck0sp;
939       break;
940     case 8:
941       // cut on V0 dca
942       okLck0sp   = TMath::Abs(v0->GetDCA())<=fCutsRD[GetGlobalIndex(8,ptbin)];
943       okLcLpi    = okLck0sp;
944       okLcLBarpi = okLck0sp;
945       break;
946     case 9:
947       // cut on V0 cosine of pointing angle wrt PV
948       okLck0sp   = d->CosV0PointingAngle()>=fCutsRD[GetGlobalIndex(9,ptbin)];
949       okLcLpi    = okLck0sp;
950       okLcLBarpi = okLck0sp;
951       break;
952     case 10:
953       // cut on bachelor transverse impact parameter wrt PV
954       okLck0sp   = TMath::Abs(d->Getd0Prong(0))<=fCutsRD[GetGlobalIndex(10,ptbin)];
955       okLcLpi    = okLck0sp;
956       okLcLBarpi = okLck0sp;
957       break;
958     case 11:
959       // cut on V0 transverse impact parameter wrt PV
960       okLck0sp   = TMath::Abs(d->Getd0Prong(1))<=fCutsRD[GetGlobalIndex(11,ptbin)];
961       okLcLpi    = okLck0sp;
962       okLcLBarpi = okLck0sp;
963       break;
964     case 12:
965       // cut on K0S invariant mass veto
966       okLcLpi    = TMath::Abs(mk0s-mk0sPDG)>=fCutsRD[GetGlobalIndex(12,ptbin)];
967       okLcLBarpi = TMath::Abs(mk0s-mk0sPDG)>=fCutsRD[GetGlobalIndex(12,ptbin)];
968       break;
969     case 13:
970       // cut on Lambda/LambdaBar invariant mass veto
971       okLck0sp   = (TMath::Abs(mlambda-mLPDG)>=fCutsRD[GetGlobalIndex(13,ptbin)] &&
972                     TMath::Abs(malambda-mLPDG)>=fCutsRD[GetGlobalIndex(13,ptbin)]);
973       break;
974     case 14:
975       // cut on gamma invariant mass veto
976       okLck0sp   = v0->InvMass2Prongs(0,1,11,11)>=fCutsRD[GetGlobalIndex(14,ptbin)];
977       okLcLpi    = okLck0sp;
978       okLcLBarpi = okLck0sp;
979       break;
980     case 15:
981       // cut on V0 pT min
982       okLck0sp   = v0->Pt()>=fCutsRD[GetGlobalIndex(15,ptbin)];
983       okLcLpi    = okLck0sp;
984       okLcLBarpi = okLck0sp;
985       break;
986     }
987   }
988
989   Int_t returnvalue = okLck0sp+2*okLcLBarpi+4*okLcLpi;
990   /*
991     retvalue case
992     1  Lc->K0S + p
993     2  Lc->LambdaBar + pi
994     3  Lc->K0S + p AND Lc->LambdaBar + pi
995     4  Lc->Lambda + pi
996     5  Lc->K0S + p AND Lc->Lambda + pi
997     6  Lc->LambdaBar + pi AND Lc->Lambda + pi
998     7  Lc->K0S + p AND Lc->LambdaBar + pi AND Lc->Lambda + pi
999   */
1000
1001
1002   /*
1003     Int_t returnvaluePID = 7;
1004
1005     // selection on PID
1006     if (selectionLevel==AliRDHFCuts::kAll ||
1007     selectionLevel==AliRDHFCuts::kCandidate ||
1008     selectionLevel==AliRDHFCuts::kPID )
1009     returnvaluePID = IsSelectedPID(d);
1010   */
1011
1012   Int_t returnvalueTot = 0;
1013   //if ( fUsePID )
1014   //returnvalueTot = CombineCuts(returnvalue,returnvaluePID);
1015   //else
1016   returnvalueTot = returnvalue;
1017
1018   return returnvalueTot;
1019
1020 }
1021 //----------------------------------
1022 void AliRDHFCutsLctoV0::SetStandardCutsPP2010() {
1023
1024   SetName("LctoV0ProductionCuts");
1025   SetTitle("Production cuts for Lc->V0+bachelor analysis");
1026
1027   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1028   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1029   //default
1030   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1031   esdTrackCuts->SetRequireITSRefit(kTRUE);
1032   esdTrackCuts->SetMinNClustersITS(0);//(4); // default is 5
1033   esdTrackCuts->SetMinNClustersTPC(70);
1034   //esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1035   //                                       AliESDtrackCuts::kAny);
1036   // default is kBoth, otherwise kAny
1037   esdTrackCuts->SetMinDCAToVertexXY(0.);
1038   esdTrackCuts->SetPtRange(0.3,1.e10);
1039   //esdTrackCuts->SetEtaRange(-0.8,+0.8);
1040   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
1041   AddTrackCuts(esdTrackCuts);
1042
1043
1044   AliESDtrackCuts* esdTrackCutsV0daughters=new AliESDtrackCuts();
1045   esdTrackCutsV0daughters->SetRequireSigmaToVertex(kFALSE);
1046   //default
1047   esdTrackCutsV0daughters->SetRequireTPCRefit(kTRUE);
1048   esdTrackCutsV0daughters->SetRequireITSRefit(kFALSE);//(kTRUE);
1049   esdTrackCutsV0daughters->SetMinNClustersITS(0);//(4); // default is 5
1050   esdTrackCutsV0daughters->SetMinNClustersTPC(70);
1051   //esdTrackCutsV0daughters->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1052   //                                                  AliESDtrackCuts::kAny);
1053   // default is kBoth, otherwise kAny
1054   esdTrackCutsV0daughters->SetMinDCAToVertexXY(0.);
1055   esdTrackCutsV0daughters->SetPtRange(0.,1.e10);
1056   esdTrackCutsV0daughters->SetAcceptKinkDaughters(kFALSE);
1057   AddTrackCutsV0daughters(esdTrackCutsV0daughters);
1058
1059   const Int_t nptbins=1;
1060   Float_t* ptbins;
1061   ptbins=new Float_t[nptbins+1];
1062   ptbins[0]=0.;
1063   ptbins[1]=99999999.;
1064
1065   SetPtBins(nptbins+1,ptbins);
1066   SetPtBins(nptbins+1,ptbins);
1067
1068   const Int_t nvars=17;
1069
1070   Float_t** prodcutsval;
1071   prodcutsval=new Float_t*[nvars];
1072   for(Int_t ic=0;ic<nvars;ic++){prodcutsval[ic]=new Float_t[nptbins];}
1073   for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
1074     prodcutsval[0][ipt2]=1.;    // inv. mass if K0S [GeV/c2]
1075     prodcutsval[1][ipt2]=1.;    // inv. mass if Lambda [GeV/c2]
1076     prodcutsval[2][ipt2]=0.05;  // inv. mass V0 if K0S [GeV/c2]
1077     prodcutsval[3][ipt2]=0.05;  // inv. mass V0 if Lambda [GeV/c2]
1078     prodcutsval[4][ipt2]=0.3;   // pT min bachelor track [GeV/c] // AOD by construction
1079     prodcutsval[5][ipt2]=0.;    // pT min V0-positive track [GeV/c]
1080     prodcutsval[6][ipt2]=0.;    // pT min V0-negative track [GeV/c]
1081     prodcutsval[7][ipt2]=1000.; // dca cascade cut [cm]
1082     prodcutsval[8][ipt2]=1.5;   // dca V0 cut [nSigma] // it's 1.5 x offline V0s
1083     prodcutsval[9][ipt2]=-1.;   // cosPA V0 cut // it's 0.90 x offline V0s at reconstruction level, 0.99 at filtering level
1084     prodcutsval[10][ipt2]=3.;   // d0 max bachelor wrt PV [cm]
1085     prodcutsval[11][ipt2]=1000.;// d0 max V0 wrt PV [cm]
1086     prodcutsval[12][ipt2]=0.;   // mass K0S veto [GeV/c2]
1087     prodcutsval[13][ipt2]=0.;   // mass Lambda/LambdaBar veto [GeV/c2]
1088     prodcutsval[14][ipt2]=0.;   // mass Gamma veto [GeV/c2]
1089     prodcutsval[15][ipt2]=0.;   // pT min V0 track [GeV/c]
1090     prodcutsval[16][ipt2]=0.;   // V0 type cut
1091   }
1092   SetCuts(nvars,nptbins,prodcutsval);
1093
1094   SetGlobalIndex(nvars,nptbins);
1095   SetPtBins(nptbins+1,ptbins);
1096
1097
1098   //pid settings
1099   //1. bachelor: default one
1100   AliAODPidHF* pidObjBachelor = new AliAODPidHF();
1101   Double_t sigmasBac[5]={3.,1.,1.,3.,3.}; // 0, 1(A), 2(A) -> TPC; 3 -> TOF; 4 -> ITS
1102   pidObjBachelor->SetSigma(sigmasBac);
1103   pidObjBachelor->SetAsym(kFALSE);
1104   pidObjBachelor->SetMatch(1);
1105   pidObjBachelor->SetTPC(kTRUE);
1106   pidObjBachelor->SetTOF(kTRUE);
1107   pidObjBachelor->SetTOFdecide(kFALSE);
1108   SetPidHF(pidObjBachelor);
1109
1110   SetUsePID(kFALSE);//(kTRUE);
1111
1112   //PrintAll();
1113
1114   for(Int_t iiv=0;iiv<nvars;iiv++){
1115     delete [] prodcutsval[iiv];
1116   }
1117   delete [] prodcutsval;
1118   prodcutsval=NULL;
1119   delete [] ptbins;
1120   ptbins=NULL;
1121
1122
1123   delete pidObjBachelor;
1124   pidObjBachelor=NULL;
1125
1126   return;
1127 }
1128 //------------------
1129 void AliRDHFCutsLctoV0::SetStandardCutsPbPb2010() {
1130
1131   SetName("LctoV0ProductionCuts");
1132   SetTitle("Production cuts for Lc->V0+bachelor analysis");
1133
1134   SetStandardCutsPP2010();
1135
1136   return;
1137 }
1138 //------------------
1139 void AliRDHFCutsLctoV0::SetStandardCutsPbPb2011() {
1140
1141   // Default 2010 PbPb cut object
1142   SetStandardCutsPbPb2010();
1143
1144   //
1145   // Enable all 2011 PbPb run triggers
1146   //
1147   SetTriggerClass("");
1148   ResetMaskAndEnableMBTrigger();
1149   EnableCentralTrigger();
1150   EnableSemiCentralTrigger();
1151
1152 }
1153 //-----------------------
1154 Int_t AliRDHFCutsLctoV0::GetV0Type(){
1155
1156   const Int_t nvars = this->GetNVars() ;
1157   //Float_t *vArray =GetCuts();
1158   //fV0Type = vArray[nvars-1];
1159   fV0Type = (this->GetCuts())[nvars-1];
1160   //this->GetCuts(vArray);
1161   TString *sVarNames=GetVarNames();
1162
1163   if(sVarNames[nvars-1].Contains("V0 type")) return (Int_t)fV0Type;
1164   else {AliInfo("AliRDHFCutsLctoV0 Last variable is not the V0 type!!!"); return -999;}
1165 }
1166
1167 //---------------------------------------------------------------------------
1168 void AliRDHFCutsLctoV0::PrintAll() const {
1169   //
1170   // print all cuts values
1171   // 
1172
1173   printf("Minimum vtx type %d\n",fMinVtxType);
1174   printf("Minimum vtx contr %d\n",fMinVtxContr);
1175   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
1176   printf("Min SPD mult %d\n",fMinSPDMultiplicity);
1177   printf("Use PID %d (PID selection flag = %d) OldPid=%d\n",(Int_t)fUsePID,(Int_t)fPidSelectionFlag,fPidHF ? fPidHF->GetOldPid() : -1);
1178   printf("High value for pT %f\n",fHighPtCut);
1179   printf("Low and high values for pT cuts: %f %f\n",fLowPtCut,fHighPtCut);
1180   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
1181   printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
1182   printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
1183   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
1184   printf("UseTrackSelectionWithFilterBits: %s\n",fUseTrackSelectionWithFilterBits ? "Yes" : "No");
1185   printf("Reject kink: %s\n",fKinkReject ? "Yes" : "No");
1186   if(fOptPileup==1) printf(" -- Reject pileup event");
1187   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
1188   if(fUseCentrality>0) {
1189     TString estimator="";
1190     if(fUseCentrality==1) estimator = "V0";
1191     if(fUseCentrality==2) estimator = "Tracks";
1192     if(fUseCentrality==3) estimator = "Tracklets";
1193     if(fUseCentrality==4) estimator = "SPD clusters outer"; 
1194     printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
1195   }
1196   if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
1197
1198   if(fVarNames){
1199     cout<<"Array of variables"<<endl;
1200     for(Int_t iv=0;iv<fnVars;iv++){
1201       cout<<fVarNames[iv]<<"\t";
1202     }
1203     cout<<endl;
1204   }
1205   if(fVarsForOpt){
1206     cout<<"Array of optimization"<<endl;
1207     for(Int_t iv=0;iv<fnVars;iv++){
1208       cout<<fVarsForOpt[iv]<<"\t";
1209     }
1210     cout<<endl;
1211   }
1212   if(fIsUpperCut){
1213     cout<<"Array of upper/lower cut"<<endl;
1214     for(Int_t iv=0;iv<fnVars;iv++){
1215       cout<<fIsUpperCut[iv]<<"\t";
1216     }
1217     cout<<endl;
1218   }
1219   if(fPtBinLimits){
1220     cout<<"Array of ptbin limits"<<endl;
1221     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
1222       cout<<fPtBinLimits[ib]<<"\t";
1223     }
1224     cout<<endl;
1225   }
1226   if(fCutsRD){
1227     cout<<"Matrix of cuts"<<endl;
1228     for(Int_t iv=0;iv<fnVars;iv++){
1229       for(Int_t ib=0;ib<fnPtBins;ib++){
1230         cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
1231       } 
1232       cout<<endl;
1233     }
1234     cout<<endl;
1235   }
1236
1237   if (fTrackCuts) {
1238     Float_t eta1=0, eta2=0; fTrackCuts->GetEtaRange(eta1,eta2);
1239     cout << " etaRange for Bachelor: [" << eta1 << "," << eta2 << "]\n";
1240   }
1241   if (fV0daughtersCuts) {
1242     Float_t eta3=0, eta4=0; fV0daughtersCuts->GetEtaRange(eta3,eta4);
1243     cout << " etaRange for V0daughters: [" << eta3 << "," << eta4 << "]\n";
1244   }
1245   return;
1246
1247 }
1248
1249 //-------------------------
1250 Bool_t AliRDHFCutsLctoV0::IsInFiducialAcceptance(Double_t pt, Double_t y) const
1251 {
1252   //
1253   //
1254   // Checking if Lc is in fiducial acceptance region
1255   //
1256   //
1257
1258   if(fMaxRapidityCand>-998.){
1259     if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
1260     else return kTRUE;
1261   }
1262
1263   if(pt > 5.) {
1264     // applying cut for pt > 5 GeV
1265     AliDebug(2,Form("pt of Lambda_c = %f (> 5), cutting at |y| < 0.8",pt));
1266     if (TMath::Abs(y) > 0.8) return kFALSE;
1267
1268   } else {
1269     // appliying smooth cut for pt < 5 GeV
1270     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
1271     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
1272     AliDebug(2,Form("pt of Lambda_c = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
1273     if (y < minFiducialY || y > maxFiducialY) return kFALSE;
1274   }
1275   //
1276   return kTRUE;
1277 }
1278 //---------------------------------------------------------------------------
1279 Bool_t AliRDHFCutsLctoV0::AreLctoV0DaughtersSelected(AliAODRecoDecayHF *dd) const{
1280   //
1281   // Daughter tracks selection
1282   //
1283
1284   AliAODRecoCascadeHF* d = (AliAODRecoCascadeHF*)dd;
1285   if (!d) {
1286     AliDebug(2,"AliAODRecoCascadeHF null");
1287     return kFALSE;
1288   }
1289
1290   if (!fTrackCuts) {
1291     AliFatal("Cut object is not defined for bachelor. Candidate accepted.");
1292     return kFALSE;
1293   }
1294
1295   AliAODTrack * bachelorTrack = dynamic_cast<AliAODTrack*>(d->GetBachelor());
1296   if (!bachelorTrack) return kFALSE;
1297
1298   if (fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst) {
1299     if(!bachelorTrack->HasPointOnITSLayer(0)) return kFALSE;
1300   }
1301
1302   if (fKinkReject != (!(fTrackCuts->GetAcceptKinkDaughters())) ) {
1303     AliError(Form("Not compatible setting: fKinkReject=%1d - fTrackCuts->GetAcceptKinkDaughters()=%1d",fKinkReject, fTrackCuts->GetAcceptKinkDaughters()));
1304     return kFALSE;
1305   }
1306
1307   AliAODVertex *vAOD = d->GetPrimaryVtx();
1308   Double_t pos[3]; vAOD->GetXYZ(pos);
1309   Double_t cov[6]; vAOD->GetCovarianceMatrix(cov);
1310   const AliESDVertex vESD(pos,cov,100.,100);
1311
1312   if (!IsDaughterSelected(bachelorTrack,&vESD,fTrackCuts)) return kFALSE;
1313
1314   if (!fV0daughtersCuts) {
1315     AliFatal("Cut object is not defined for V0daughters. Candidate accepted.");
1316     return kFALSE;
1317   }
1318
1319   AliAODv0 * v0 = dynamic_cast<AliAODv0*>(d->Getv0());
1320   if (!v0) return kFALSE;
1321   AliAODTrack *v0positiveTrack = dynamic_cast<AliAODTrack*>(d->Getv0PositiveTrack());
1322   if (!v0positiveTrack) return kFALSE;
1323   AliAODTrack *v0negativeTrack = dynamic_cast<AliAODTrack*>(d->Getv0NegativeTrack());
1324   if (!v0negativeTrack) return kFALSE;
1325
1326
1327   Float_t etaMin=0, etaMax=0; fV0daughtersCuts->GetEtaRange(etaMin,etaMax);
1328   if ( (v0positiveTrack->Eta()<=etaMin || v0positiveTrack->Eta()>=etaMax) ||
1329        (v0negativeTrack->Eta()<=etaMin || v0negativeTrack->Eta()>=etaMax) ) return kFALSE;
1330   Float_t ptMin=0, ptMax=0; fV0daughtersCuts->GetPtRange(ptMin,ptMax);
1331   if ( (v0positiveTrack->Pt()<=ptMin || v0positiveTrack->Pt()>=ptMax) ||
1332        (v0negativeTrack->Pt()<=ptMin || v0negativeTrack->Pt()>=ptMax) ) return kFALSE;
1333
1334   // Condition on nTPCclusters
1335   if (fV0daughtersCuts->GetMinNClusterTPC()>0) {
1336     if ( ( ( v0positiveTrack->GetTPCClusterInfo(2,1) ) < fV0daughtersCuts->GetMinNClusterTPC() ) || 
1337          ( ( v0negativeTrack->GetTPCClusterInfo(2,1) ) < fV0daughtersCuts->GetMinNClusterTPC() ) ) return kFALSE;
1338   }
1339
1340   // kTPCrefit status
1341   if (v0->GetOnFlyStatus()==kFALSE) { // only for offline V0s
1342     if (fV0daughtersCuts->GetRequireTPCRefit()) {
1343       if( !(v0positiveTrack->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1344       if( !(v0negativeTrack->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1345     }
1346   }
1347   // kink condition
1348   if (!fV0daughtersCuts->GetAcceptKinkDaughters()) {
1349     AliAODVertex *maybeKinkPos = (AliAODVertex*)v0positiveTrack->GetProdVertex();
1350     AliAODVertex *maybeKinkNeg = (AliAODVertex*)v0negativeTrack->GetProdVertex();
1351     if (maybeKinkPos->GetType()==AliAODVertex::kKink ||
1352         maybeKinkNeg->GetType()==AliAODVertex::kKink) return kFALSE;
1353   }
1354   // Findable clusters > 0 condition - from V0 analysis
1355   //if( v0positiveTrack->GetTPCNclsF()<=0 || v0negativeTrack->GetTPCNclsF()<=0 ) return kFALSE;
1356   /*
1357     Float_t lPosTrackCrossedRows = v0positiveTrack->GetTPCClusterInfo(2,1);
1358     Float_t lNegTrackCrossedRows = v0positiveTrack->GetTPCClusterInfo(2,1);
1359     fTreeVariableLeastNbrCrossedRows = (Int_t) lPosTrackCrossedRows;
1360     if( lNegTrackCrossedRows < fTreeVariableLeastNbrCrossedRows )
1361     fTreeVariableLeastNbrCrossedRows = (Int_t) lNegTrackCrossedRows;
1362     //Compute ratio Crossed Rows / Findable clusters
1363     //Note: above test avoids division by zero!
1364     Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF()));
1365     Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF()));
1366     fTreeVariableLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
1367     if( lNegTrackCrossedRowsOverFindable < fTreeVariableLeastRatioCrossedRowsOverFindable )
1368     fTreeVariableLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
1369     //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
1370     if ( fTreeVariableLeastRatioCrossedRowsOverFindable < 0.8 ) return kFALSE;
1371   */
1372
1373   return kTRUE;
1374
1375 }