Updates: finer mass bin, default cuts per system (Rossella, Marcel, Norman)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCutsLctopKpi.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->pKpi
21 //
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <TDatabasePDG.h>
26 #include <Riostream.h>
27 #include <AliAnalysisManager.h>
28 #include <AliInputEventHandler.h>
29 #include <AliPIDResponse.h>
30
31 #include "AliRDHFCutsLctopKpi.h"
32 #include "AliAODRecoDecayHF3Prong.h"
33 #include "AliRDHFCuts.h"
34 #include "AliAODTrack.h"
35 #include "AliESDtrack.h"
36 #include "AliKFParticle.h"
37 #include "AliESDVertex.h"
38
39 using std::cout;
40 using std::endl;
41
42 ClassImp(AliRDHFCutsLctopKpi)
43
44 //--------------------------------------------------------------------------
45 AliRDHFCutsLctopKpi::AliRDHFCutsLctopKpi(const char* name) : 
46 AliRDHFCuts(name),
47 fPidObjprot(0),
48 fPidObjpion(0),
49 fUseImpParProdCorrCut(kFALSE),
50 fPIDStrategy(kNSigma),
51 fCutsStrategy(kStandard),
52 fUseSpecialCut(kFALSE)
53 {
54   //
55   // Default Constructor
56   //
57   Int_t nvars=13;
58   SetNVars(nvars);
59   TString varNames[13]={"inv. mass [GeV]",
60                         "pTK [GeV/c]",
61                         "pTP [GeV/c]",
62                         "d0K [cm]   lower limit!",
63                         "d0Pi [cm]  lower limit!",
64                         "dist12 (cm)",
65                         "sigmavert (cm)",
66                         "dist prim-sec (cm)",
67                         "pM=Max{pT1,pT2,pT3} (GeV/c)",
68                         "cosThetaPoint",
69                         "Sum d0^2 (cm^2)",
70                         "dca cut (cm)",
71                         "cut on pTpion [GeV/c]"};
72   Bool_t isUpperCut[13]={kTRUE,
73                          kFALSE,
74                          kFALSE,
75                          kFALSE,
76                          kFALSE,
77                          kFALSE,
78                          kTRUE,
79                          kFALSE,
80                          kFALSE,
81                          kFALSE,
82                          kFALSE,
83                          kTRUE,
84                          kFALSE
85                          };
86   SetVarNames(nvars,varNames,isUpperCut);
87   Bool_t forOpt[13]={kFALSE,
88                      kTRUE,
89                      kTRUE,
90                      kFALSE,
91                      kFALSE,
92                      kFALSE,
93                      kFALSE,
94                      kTRUE,
95                      kFALSE,
96                      kFALSE,
97                      kFALSE,
98                      kFALSE,
99                      kTRUE};
100   SetVarsForOpt(4,forOpt);
101   Float_t limits[2]={0,999999999.};
102   SetPtBins(2,limits);
103   for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies)
104       fPIDThreshold[ispecies]=0.;
105 }
106 //--------------------------------------------------------------------------
107 AliRDHFCutsLctopKpi::AliRDHFCutsLctopKpi(const AliRDHFCutsLctopKpi &source) :
108   AliRDHFCuts(source),
109   fPidObjprot(0x0),
110   fPidObjpion(0x0),
111   fUseImpParProdCorrCut(source.fUseImpParProdCorrCut),
112   fPIDStrategy(source.fPIDStrategy),
113   fCutsStrategy(source.fCutsStrategy),
114   fUseSpecialCut(source.fUseSpecialCut)
115 {
116   //
117   // Copy constructor
118   //
119   if (source.fPidObjprot) fPidObjprot = new AliAODPidHF(*(source.fPidObjprot));
120   else fPidObjprot = new AliAODPidHF();
121   if (source.fPidObjpion) fPidObjpion = new AliAODPidHF(*(source.fPidObjpion));
122   else fPidObjpion = new AliAODPidHF();
123   memcpy(fPIDThreshold,source.fPIDThreshold,AliPID::kSPECIES*sizeof(Double_t));
124 }
125 //--------------------------------------------------------------------------
126 AliRDHFCutsLctopKpi &AliRDHFCutsLctopKpi::operator=(const AliRDHFCutsLctopKpi &source)
127 {
128   //
129   // assignment operator
130   //
131   if(this != &source) {
132     
133     AliRDHFCuts::operator=(source);
134     delete fPidObjprot;
135     fPidObjprot = new AliAODPidHF(*(source.fPidObjprot));
136     delete fPidObjpion;
137     fPidObjpion = new AliAODPidHF(*(source.fPidObjpion));
138     fPIDStrategy=source.fPIDStrategy;
139     fCutsStrategy=source.fCutsStrategy;
140     memcpy(fPIDThreshold,source.fPIDThreshold,AliPID::kSPECIES*sizeof(Double_t));
141   }
142     
143   return *this;
144 }
145 //---------------------------------------------------------------------------
146 AliRDHFCutsLctopKpi::~AliRDHFCutsLctopKpi() {
147  //
148  //  // Default Destructor
149  //   
150  if(fPidObjpion){
151   delete fPidObjpion;
152   fPidObjpion=0;
153  }
154  if(fPidObjprot){
155   delete fPidObjprot;
156   fPidObjprot=0;
157  }
158
159 }
160
161 //---------------------------------------------------------------------------
162 void AliRDHFCutsLctopKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters, AliAODEvent *aod) {
163   // 
164   // Fills in vars the values of the variables 
165   //
166
167   if(nvars!=fnVarsForOpt) {
168     printf("AliRDHFCutsLctopKpi::GetCutsVarsForOpt: wrong number of variables\n");
169     return;
170   }
171
172   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
173
174   Int_t iter=-1;
175   if(fVarsForOpt[0]){
176     iter++;
177     vars[iter]=dd->InvMassLcpKpi();
178   }
179   if(fVarsForOpt[1]){
180     iter++;
181     for(Int_t iprong=0;iprong<3;iprong++){
182       if(TMath::Abs(pdgdaughters[iprong])==321) {
183         vars[iter]=dd->PtProng(iprong);
184       }
185     }
186   }
187   if(fVarsForOpt[2]){
188     iter++;
189     for(Int_t iprong=0;iprong<3;iprong++){
190       if(TMath::Abs(pdgdaughters[iprong])==2212) {
191         vars[iter]=dd->PtProng(iprong);
192       }
193     }
194   }
195   if(fVarsForOpt[3]){
196     iter++;
197     for(Int_t iprong=0;iprong<3;iprong++){
198       if(TMath::Abs(pdgdaughters[iprong])==2212) {
199         vars[iter]=dd->Getd0Prong(iprong);
200       }
201     }
202   }
203   if(fVarsForOpt[4]){
204     iter++;
205     for(Int_t iprong=0;iprong<3;iprong++){
206       if(TMath::Abs(pdgdaughters[iprong])==211) {
207         vars[iter]=dd->Getd0Prong(iprong);
208       }
209     }
210   }
211   if(fVarsForOpt[5]){
212     iter++;
213     vars[iter]=dd->GetDist12toPrim();
214   }
215   if(fVarsForOpt[6]){
216     iter++;
217     vars[iter]=dd->GetSigmaVert(aod);
218   }
219   if(fVarsForOpt[7]){
220     iter++;
221     vars[iter] = dd->DecayLength();
222   }
223   if(fVarsForOpt[8]){
224     iter++;
225     Float_t ptmax=0;
226     for(Int_t i=0;i<3;i++){
227       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
228     }
229     vars[iter]=ptmax;
230   }
231   if(fVarsForOpt[9]){
232     iter++;
233     vars[iter]=dd->CosPointingAngle();
234   }
235   if(fVarsForOpt[10]){
236     iter++;
237     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
238   }
239   if(fVarsForOpt[11]){
240     iter++;
241     vars[iter]=dd->GetDCA();
242   }
243   if(fVarsForOpt[12]){
244     iter++;
245     for(Int_t iprong=0;iprong<3;iprong++){
246       if(TMath::Abs(pdgdaughters[iprong])==211) {
247         vars[iter]=dd->PtProng(iprong);
248       }
249     }
250   }
251
252   return;
253 }
254 //---------------------------------------------------------------------------
255 Int_t AliRDHFCutsLctopKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent *aod) {
256   //
257   // Apply selection
258   //
259
260   if(!fCutsRD){
261     AliError("Cut matrice not inizialized. Exit...\n");
262     return 0;
263   }
264   //PrintAll();
265   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
266
267   if(!d){
268     AliError("AliAODRecoDecayHF3Prong null \n");
269     return 0;
270   }
271
272
273   if(fKeepSignalMC) if(IsSignalMC(d,aod,4122)) return 3;
274
275   Int_t returnvalue=3;
276   Int_t returnvaluePID=3;
277
278   if(d->Pt()<fMinPtCand) return 0;
279   if(d->Pt()>fMaxPtCand) return 0;
280
281   if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
282
283
284   // selection on daughter tracks 
285   if(selectionLevel==AliRDHFCuts::kAll || 
286      selectionLevel==AliRDHFCuts::kTracks) {
287     if(!AreDaughtersSelected(d)) return 0;
288   }
289
290
291   // PID selection
292   if(selectionLevel==AliRDHFCuts::kAll ||
293      selectionLevel==AliRDHFCuts::kCandidate|| 
294      selectionLevel==AliRDHFCuts::kPID) {
295     switch (fPIDStrategy) {
296     case kNSigma:
297       returnvaluePID = IsSelectedPID(d);
298       break;
299      case kNSigmaPbPb:
300       returnvaluePID = IsSelectedNSigmaPbPb(d);
301       break;
302     case kCombined:
303       returnvaluePID = IsSelectedCombinedPID(d);
304       break;
305     case kCombinedSoft:
306       returnvaluePID = IsSelectedCombinedPIDSoft(d);
307       break;
308     case kNSigmaStrong:
309       returnvaluePID = IsSelectedPIDStrong(d);
310       break;
311     case kCombinedpPb:
312       returnvaluePID = IsSelectedCombinedPIDpPb(d);
313       break;
314     case kCombinedpPb2:
315       returnvaluePID = IsSelectedCombinedPIDpPb2(d);
316       break;
317     }
318     fIsSelectedPID=returnvaluePID;
319   }
320   //  if(fUsePID || selectionLevel==AliRDHFCuts::kPID) returnvaluePID = IsSelectedCombinedPID(d);   // to test!!
321   if(returnvaluePID==0) return 0;
322
323
324
325
326   // selection on candidate
327   if(selectionLevel==AliRDHFCuts::kAll || 
328      selectionLevel==AliRDHFCuts::kCandidate) {
329
330     Double_t pt=d->Pt();
331     
332     Int_t ptbin=PtBin(pt);
333     
334     Double_t mLcpKpi=0.,mLcpiKp=0.;
335     Int_t okLcpKpi=1,okLcpiKp=1;
336
337     Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
338
339     mLcpKpi=d->InvMassLcpKpi();
340     mLcpiKp=d->InvMassLcpiKp();
341
342     if(TMath::Abs(mLcpKpi-mLcPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okLcpKpi = 0;
343     if(TMath::Abs(mLcpiKp-mLcPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okLcpiKp = 0;
344     if(!okLcpKpi && !okLcpiKp) return 0;
345
346   switch (fCutsStrategy) {
347
348     case kStandard:
349     if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] || TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;//Kaon
350     if(d->Pt()>=3. && d->PProng(1)<0.55) return 0;
351     if(fUseSpecialCut) {
352       if(TMath::Abs(d->PtProng(0)) < TMath::Abs(d->PtProng(2)) )okLcpKpi=0;
353       if(TMath::Abs(d->PtProng(2)) < TMath::Abs(d->PtProng(0)) )okLcpiKp=0;
354     }
355     if((TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(12,ptbin)])) okLcpKpi=0;
356     if((TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(12,ptbin)]))okLcpiKp=0;
357     if(!okLcpKpi && !okLcpiKp) return 0;
358     //2track cuts
359     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]|| d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]) return 0;
360     if(d->GetDist12toPrim()>0.5) return 0;
361     if(d->GetDist23toPrim()>0.5) return 0;
362     if(fUseImpParProdCorrCut){
363       if(d->Getd0Prong(0)*d->Getd0Prong(1)<0. && d->Getd0Prong(2)*d->Getd0Prong(1)<0.) return 0;
364     }
365     //sec vert
366     if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
367     if(d->DecayLength()>0.5) return 0;
368
369   //  Double_t sumd0s=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
370   //  if(sumd0s<fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
371     if((d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
372     
373     if(TMath::Abs(d->PtProng(0))<fCutsRD[GetGlobalIndex(8,ptbin)] && TMath::Abs(d->PtProng(1))<fCutsRD[GetGlobalIndex(8,ptbin)] && TMath::Abs(d->PtProng(2))<fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
374     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]) return 0;
375     if(d->GetSigmaVert(aod)>fCutsRD[GetGlobalIndex(6,ptbin)]) return 0;
376     
377     //DCA
378     for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]) return 0;
379
380     break;
381
382    case kKF:
383     Int_t pdgs[3]={0,321,0};
384     Bool_t constraint=kFALSE;
385     if(fCutsRD[GetGlobalIndex(1,ptbin)]>0.) constraint=kTRUE;
386     Double_t field=aod->GetMagneticField();
387     if (returnvaluePID==1 || returnvaluePID==3){
388
389       pdgs[0]=2122;pdgs[2]=211;
390       AliKFParticle *lc1=ReconstructKF(d,pdgs,field,constraint);
391       if(!lc1){
392         okLcpKpi=0;
393       }else{
394         if(lc1->GetChi2()/lc1->GetNDF()>fCutsRD[GetGlobalIndex(2,ptbin)]) okLcpKpi=0;
395       }
396     } else if(returnvaluePID>=2){
397
398       pdgs[0]=211;pdgs[2]=2212;
399       AliKFParticle *lc2=ReconstructKF(d,pdgs,field,constraint);
400       if(!lc2){ 
401         okLcpiKp=0;
402       }else{
403         if(lc2->GetChi2()/lc2->GetNDF()>fCutsRD[GetGlobalIndex(2,ptbin)])okLcpiKp=0; 
404       }
405     }
406     break;
407
408    }
409
410     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
411     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
412     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
413    
414   }
415
416
417   Int_t returnvalueTot=CombinePIDCuts(returnvalue,returnvaluePID);
418   return returnvalueTot;
419 }
420 //---------------------------------------------------------------------------
421 Int_t AliRDHFCutsLctopKpi::IsSelectedPID(AliAODRecoDecayHF* obj) {
422
423
424     if(!fUsePID || !obj) return 3;
425     Int_t okLcpKpi=0,okLcpiKp=0;
426     Int_t returnvalue=0;
427     Bool_t isPeriodd=fPidHF->GetOnePad();
428     Bool_t isMC=fPidHF->GetMC();
429     Bool_t ispion0=kTRUE,ispion2=kTRUE;
430     Bool_t isproton0=kFALSE,isproton2=kFALSE;
431     Bool_t iskaon1=kFALSE;
432     if(isPeriodd) {
433      fPidObjprot->SetOnePad(kTRUE);
434      fPidObjpion->SetOnePad(kTRUE);
435     }
436     if(isMC) {
437      fPidObjprot->SetMC(kTRUE);
438      fPidObjpion->SetMC(kTRUE);
439     }
440
441    if(fPidObjprot->GetPidResponse()==0x0){
442       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
443       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
444       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
445       fPidObjprot->SetPidResponse(pidResp);
446     }
447     if(fPidObjpion->GetPidResponse()==0x0){
448       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
449       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
450       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
451       fPidObjpion->SetPidResponse(pidResp);
452     }
453     if(fPidHF->GetPidResponse()==0x0){
454       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
455       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
456       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
457       fPidHF->SetPidResponse(pidResp);
458     }
459
460     for(Int_t i=0;i<3;i++){
461      AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
462      if(!track) return 0;
463      // identify kaon
464      if(track->P()<0.55){
465       fPidHF->SetTOF(kFALSE);
466       fPidHF->SetTOFdecide(kFALSE);
467      }
468      if(i==1) {
469       Int_t isKaon=fPidHF->MakeRawPid(track,3); 
470       if(isKaon>=1) iskaon1=kTRUE;
471       if(track->P()<0.55){
472        fPidHF->SetTOF(kTRUE);
473        fPidHF->SetTOFdecide(kTRUE);
474       }
475       
476       if(!iskaon1) return 0;
477      
478      }else{
479      //pion or proton
480     if(track->P()<1.){
481       fPidObjprot->SetTOF(kFALSE);
482       fPidObjprot->SetTOFdecide(kFALSE);
483      }
484       
485      Int_t isProton=fPidObjprot->MakeRawPid(track,4);
486      
487
488      Int_t isPion=fPidObjpion->MakeRawPid(track,2);
489      
490      if(track->P()<1.){
491       fPidObjprot->SetTOF(kTRUE);
492       fPidObjprot->SetTOFdecide(kTRUE);
493      }
494      
495
496      if(i==0) {
497       if(isPion<0) ispion0=kFALSE;
498       if(isProton>=1) isproton0=kTRUE;
499
500      }
501       if(!ispion0 && !isproton0) return 0;
502      if(i==2) {
503       if(isPion<0) ispion2=kFALSE;
504       if(isProton>=1) isproton2=kTRUE;
505      }
506
507     }
508    }
509
510     if(ispion2 && isproton0 && iskaon1) okLcpKpi=1;
511     if(ispion0 && isproton2 && iskaon1) okLcpiKp=1;
512     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
513     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
514     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
515
516  return returnvalue;
517 }
518
519 //--------------------------------------------------------------------------
520
521 Int_t AliRDHFCutsLctopKpi::IsSelectedNSigmaPbPb(AliAODRecoDecayHF* obj) {
522
523
524   if(!fUsePID || !obj) return 3;
525   Int_t okLcpKpi=0,okLcpiKp=0;
526   Int_t returnvalue=0;
527   Bool_t isPeriodd=fPidHF->GetOnePad();
528   Bool_t isMC=fPidHF->GetMC();
529   Bool_t ispion0=kTRUE,ispion2=kTRUE;
530   Bool_t isproton0=kFALSE,isproton2=kFALSE;
531   Bool_t iskaon1=kFALSE;
532   if(isPeriodd) {
533     fPidObjprot->SetOnePad(kTRUE);
534     fPidObjpion->SetOnePad(kTRUE);
535   }
536   if(isMC) {
537     fPidObjprot->SetMC(kTRUE);
538     fPidObjpion->SetMC(kTRUE);
539   }
540
541   if(fPidObjprot->GetPidResponse()==0x0){
542     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
543     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
544     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
545     fPidObjprot->SetPidResponse(pidResp);
546   }
547   if(fPidObjpion->GetPidResponse()==0x0){
548     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
549     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
550     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
551     fPidObjpion->SetPidResponse(pidResp);
552   }
553   if(fPidHF->GetPidResponse()==0x0){
554     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
555     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
556     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
557     fPidHF->SetPidResponse(pidResp);
558   }
559
560   for(Int_t i=0;i<3;i++){
561     AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
562     if(!track) return 0;
563
564     if(i==1) {
565       //kaon
566       Int_t isKaon=fPidHF->MakeRawPid(track,3); 
567       if(isKaon>=1) iskaon1=kTRUE;
568       if(!iskaon1) return 0;
569     }
570
571     else {
572       //pion or proton
573       Int_t isProton=fPidObjprot->MakeRawPid(track,4);
574       Int_t isPion=fPidObjpion->MakeRawPid(track,2);
575
576       if(i==0) {
577         if(isPion<0) ispion0=kFALSE;
578         if(isProton>=1) isproton0=kTRUE;
579       }
580       if(!ispion0 && !isproton0) return 0;
581
582       if(i==2) {
583         if(isPion<0) ispion2=kFALSE;
584         if(isProton>=1) isproton2=kTRUE;
585       }
586     }
587   }
588
589   if(ispion2 && isproton0 && iskaon1) okLcpKpi=1;
590   if(ispion0 && isproton2 && iskaon1) okLcpiKp=1;
591   if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
592   if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
593   if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
594
595   return returnvalue;
596 }
597
598
599
600
601 //---------------------------------------------------------------------------
602 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPID(AliAODRecoDecayHF* obj) {
603     
604     if(!fUsePID || !obj) {return 3;}
605     Int_t okLcpKpi=0,okLcpiKp=0;
606     Int_t returnvalue=0;
607     Bool_t isPeriodd=fPidHF->GetOnePad();
608     Bool_t isMC=fPidHF->GetMC();
609
610     if(isPeriodd) {
611             fPidObjprot->SetOnePad(kTRUE);
612             fPidObjpion->SetOnePad(kTRUE);
613     }
614     if(isMC) {
615             fPidObjprot->SetMC(kTRUE);
616             fPidObjpion->SetMC(kTRUE);
617     }
618
619     AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
620     AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
621     AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
622     if (!track0 || !track1 || !track2) return 0;
623     Double_t prob0[AliPID::kSPECIES];
624     Double_t prob1[AliPID::kSPECIES];
625     Double_t prob2[AliPID::kSPECIES];
626     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
627     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
628     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
629
630    if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
631     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
632    if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
633
634     if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
635     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
636    if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
637
638     if(fPIDThreshold[AliPID::kPion]>0. && fPIDThreshold[AliPID::kKaon]>0. && fPIDThreshold[AliPID::kProton]>0.){
639     okLcpiKp=  (prob0[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ])
640              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
641              &&(prob2[AliPID::kProton]>fPIDThreshold[AliPID::kProton]);
642     okLcpKpi=  (prob0[AliPID::kProton]>fPIDThreshold[AliPID::kProton])
643              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
644              &&(prob2[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ]);
645    }else{ 
646                     //pion or proton
647                     
648                     
649     if(TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]){
650     if(TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) okLcpKpi = 1;  
651     if(TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) okLcpiKp = 1; 
652             }
653            }
654     
655     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
656     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
657     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
658     
659     return returnvalue;
660 }
661 //-----------------------
662 Int_t AliRDHFCutsLctopKpi::CombinePIDCuts(Int_t returnvalue, Int_t returnvaluePID) const {
663
664  Int_t returnvalueTot=0;
665  Int_t okLcpKpi=0,okLcpiKp=0;
666  if(returnvaluePID==1){
667    if(returnvalue==1 || returnvalue==3) okLcpKpi=1;
668  }
669  if(returnvaluePID==2){
670    if(returnvalue>=2) okLcpiKp=1;
671  }
672  if(returnvaluePID==3 && returnvalue>0){
673   if(returnvalue==1 || returnvalue==3) okLcpKpi=1;
674   if(returnvalue>=2) okLcpiKp=1;
675  } 
676
677  if(okLcpKpi) returnvalueTot=1; //cuts passed as Lc->pKpi
678  if(okLcpiKp) returnvalueTot=2; //cuts passed as Lc->piKp
679  if(okLcpKpi && okLcpiKp) returnvalueTot=3; //cuts passed as both pKpi and piKp
680  return returnvalueTot;
681 }
682 //----------------------------------
683 void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
684
685  SetName("LctopKpiProdCuts");
686  SetTitle("Production cuts for Lc analysis");
687
688  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
689  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
690  esdTrackCuts->SetRequireTPCRefit(kTRUE);
691  esdTrackCuts->SetMinNClustersTPC(70);
692  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
693                                           AliESDtrackCuts::kAny);
694  esdTrackCuts->SetRequireITSRefit(kTRUE);
695  esdTrackCuts->SetMinNClustersITS(4);
696  esdTrackCuts->SetMinDCAToVertexXY(0.);
697  esdTrackCuts->SetEtaRange(-0.8,0.8);
698  esdTrackCuts->SetPtRange(0.3,1.e10);
699  AddTrackCuts(esdTrackCuts);
700
701  const Int_t nptbins=8;
702  const Int_t nvars=13;
703  Float_t* ptbins;
704  ptbins=new Float_t[nptbins+1];
705  
706  ptbins[0]=1.;
707  ptbins[1]=2.;
708  ptbins[2]=3.;
709  ptbins[3]=4.;
710  ptbins[4]=5.;
711  ptbins[5]=6.;
712  ptbins[6]=8.;
713  ptbins[7]=10.;
714  ptbins[8]=20.;
715
716  SetGlobalIndex(nvars,nptbins);
717  SetPtBins(nptbins+1,ptbins);
718
719  Float_t** prodcutsval;
720  prodcutsval=new Float_t*[nvars];
721  for(Int_t iv=0;iv<nvars;iv++){
722   prodcutsval[iv]=new Float_t[nptbins];
723  }
724
725  for(Int_t ipt=0;ipt<nptbins;ipt++){
726   prodcutsval[0][ipt]=0.13;
727   prodcutsval[1][ipt]=0.4;
728   prodcutsval[2][ipt]=0.4;
729   prodcutsval[3][ipt]=0.;
730   prodcutsval[4][ipt]=0.;
731   prodcutsval[5][ipt]=0.;
732   prodcutsval[6][ipt]=0.06;
733   prodcutsval[7][ipt]=0.005;
734   prodcutsval[8][ipt]=0.;
735   prodcutsval[9][ipt]=0.;
736   prodcutsval[10][ipt]=0.;
737   prodcutsval[11][ipt]=0.05;
738   prodcutsval[12][ipt]=0.4;
739  }
740  SetCuts(nvars,nptbins,prodcutsval);
741
742  AliAODPidHF* pidObjK=new AliAODPidHF();
743  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
744  pidObjK->SetSigma(sigmasK);
745  pidObjK->SetAsym(kTRUE);
746  pidObjK->SetMatch(1);
747  pidObjK->SetTPC(kTRUE);
748  pidObjK->SetTOF(kTRUE);
749  pidObjK->SetITS(kTRUE);
750  Double_t plimK[2]={0.5,0.8};
751  pidObjK->SetPLimit(plimK,2);
752  pidObjK->SetTOFdecide(kTRUE);
753
754  SetPidHF(pidObjK);
755
756  AliAODPidHF* pidObjpi=new AliAODPidHF();
757  pidObjpi->SetTPC(kTRUE);
758  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
759  pidObjpi->SetSigma(sigmaspi);
760  pidObjpi->SetTOFdecide(kTRUE);
761  SetPidpion(pidObjpi);
762
763  AliAODPidHF* pidObjp=new AliAODPidHF();
764  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
765  pidObjp->SetSigma(sigmasp);
766  pidObjp->SetAsym(kTRUE);
767  pidObjp->SetMatch(1);
768  pidObjp->SetTPC(kTRUE);
769  pidObjp->SetTOF(kTRUE);
770  pidObjp->SetITS(kTRUE);
771  Double_t plimp[2]={1.,2.};
772  pidObjp->SetPLimit(plimp,2);
773  pidObjp->SetTOFdecide(kTRUE);
774
775  SetPidprot(pidObjp);
776
777  SetUsePID(kTRUE);
778  SetOptPileup(kTRUE);
779  
780  PrintAll();
781
782  for(Int_t iiv=0;iiv<nvars;iiv++){
783   delete [] prodcutsval[iiv];
784  }
785  delete [] prodcutsval;
786  prodcutsval=NULL;
787  delete [] ptbins;
788  ptbins=NULL;
789
790  delete pidObjK;
791  pidObjK=NULL;
792  delete pidObjpi;
793  pidObjpi=NULL;
794  delete pidObjp;
795  pidObjp=NULL;
796
797  return;
798 }
799 //------------------
800 void AliRDHFCutsLctopKpi::SetStandardCutsPbPb2010() {
801
802  SetName("LctopKpiProdCuts");
803  SetTitle("Production cuts for Lc analysis");
804
805  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
806
807  esdTrackCuts->SetRequireTPCRefit(kTRUE);
808  esdTrackCuts->SetMinNClustersTPC(70);
809  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
810                                           AliESDtrackCuts::kAny);
811  esdTrackCuts->SetRequireITSRefit(kTRUE);
812  esdTrackCuts->SetMinNClustersITS(4);
813  esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0100*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
814  esdTrackCuts->SetEtaRange(-0.8,0.8);
815  esdTrackCuts->SetMaxDCAToVertexXY(1.);
816  esdTrackCuts->SetMaxDCAToVertexZ(1.);
817  esdTrackCuts->SetPtRange(0.49,1.e10);
818  AddTrackCuts(esdTrackCuts);
819
820  const Int_t nptbins=8;
821  const Int_t nvars=13;
822  Float_t* ptbins;
823  ptbins=new Float_t[nptbins+1];
824  
825  ptbins[0]=1.;
826  ptbins[1]=2.;
827  ptbins[2]=3.;
828  ptbins[3]=4.;
829  ptbins[4]=5.;
830  ptbins[5]=6.;
831  ptbins[6]=8.;
832  ptbins[7]=10.;
833  ptbins[8]=99999.;
834
835
836  SetGlobalIndex(nvars,nptbins);
837  SetPtBins(nptbins+1,ptbins);
838
839  Float_t** prodcutsval;
840  prodcutsval=new Float_t*[nvars];
841  for(Int_t iv=0;iv<nvars;iv++){
842   prodcutsval[iv]=new Float_t[nptbins];
843  }
844
845  for(Int_t ipt=0;ipt<nptbins;ipt++){
846   prodcutsval[0][ipt]=0.13;
847   prodcutsval[1][ipt]=0.5;
848   prodcutsval[2][ipt]=0.6;
849   prodcutsval[3][ipt]=0.;
850   prodcutsval[4][ipt]=0.;
851   prodcutsval[5][ipt]=0.01;
852   prodcutsval[6][ipt]=0.04;
853   prodcutsval[7][ipt]=0.006;
854   prodcutsval[8][ipt]=0.8;
855   prodcutsval[9][ipt]=0.3;
856   prodcutsval[10][ipt]=0.;
857   prodcutsval[11][ipt]=0.05;
858   prodcutsval[12][ipt]=0.4;
859  }
860  SetCuts(nvars,nptbins,prodcutsval);
861
862  AliAODPidHF* pidObj=new AliAODPidHF();
863  pidObj->SetTPC(kTRUE);
864  pidObj->SetTOF(kTRUE);
865  SetPidHF(pidObj);
866  SetPidpion(pidObj);
867  SetPidprot(pidObj);
868
869
870  // bayesian pid
871  GetPidHF()->SetUseCombined(kTRUE);
872  GetPidHF()->SetUseDefaultPriors(kTRUE);
873  GetPidHF()->SetCombDetectors(AliAODPidHF::kTPCTOF);
874  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies){
875    SetPIDThreshold(static_cast<AliPID::EParticleType>(ispecies),0);
876  }
877  SetPIDStrategy(AliRDHFCutsLctopKpi::kCombinedpPb);
878  SetUsePID(kTRUE);
879
880
881  PrintAll();
882
883  for(Int_t iiv=0;iiv<nvars;iiv++){
884   delete [] prodcutsval[iiv];
885  }
886  delete [] prodcutsval;
887  prodcutsval=NULL;
888  delete [] ptbins;
889  ptbins=NULL;
890
891  delete pidObj;
892  pidObj=NULL;
893
894  return;
895 }
896 //------------------
897 AliKFParticle* AliRDHFCutsLctopKpi::ReconstructKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field,Bool_t constraint) const{
898   // Method to construct the KF particle from the candidate
899
900  const Int_t nprongs=d->GetNProngs();
901  if(nprongs<=0) return 0x0;
902
903  Int_t iprongs[nprongs];
904  for(Int_t i=0;i<nprongs;i++) iprongs[i]=i;
905
906  Double_t mass[2]={0.,0.};
907
908  AliKFParticle *decay=d->ApplyVertexingKF(iprongs,nprongs,pdgs,constraint,field,mass);
909  if(!decay) return 0x0;
910  AliESDVertex *vertexESD = new AliESDVertex(decay->Parameters(),
911                                             decay->CovarianceMatrix(),
912                                             decay->GetChi2(),
913                                             nprongs);
914  Double_t pos[3],cov[6],chi2perNDF;
915  vertexESD->GetXYZ(pos);
916  vertexESD->GetCovMatrix(cov);
917  chi2perNDF = vertexESD->GetChi2toNDF();
918  delete vertexESD; vertexESD=NULL;
919  AliAODVertex *vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
920  d->SetSecondaryVtx(vertexAOD);
921  return decay;
922 }
923
924 //------------------
925 void AliRDHFCutsLctopKpi::SetStandardCutsPbPb2011() {
926
927   // Default 2010 PbPb cut object
928   SetStandardCutsPbPb2010();
929
930   //
931   // Enable all 2011 PbPb run triggers
932   //  
933   SetTriggerClass("");
934   ResetMaskAndEnableMBTrigger();
935   EnableCentralTrigger();
936   EnableSemiCentralTrigger();
937 }
938 //-----------------
939
940 Bool_t AliRDHFCutsLctopKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
941 {
942   //
943   //  // Checking if Lc is in fiducial acceptance region 
944   //    //
945   //
946   if(fMaxRapidityCand>-998.){
947     if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
948     else return kTRUE;
949   }
950
951   if(pt > 5.) {
952     // applying cut for pt > 5 GeV
953    AliDebug(2,Form("pt of Lc = %f (> 5), cutting at |y| < 0.8",pt));
954    if (TMath::Abs(y) > 0.8) return kFALSE;
955   
956   } else {
957    // appliying smooth cut for pt < 5 GeV
958    Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
959    Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
960   AliDebug(2,Form("pt of Lc = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
961   if (y < minFiducialY || y > maxFiducialY) return kFALSE;
962  }
963   //
964   return kTRUE;
965 }
966 //--------------------------------------------------------
967 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPIDSoft(AliAODRecoDecayHF* obj) {
968  if(!fUsePID || !obj) {return 3;}
969  Int_t okLcpKpi=0,okLcpiKp=0;
970  Int_t returnvalue=0;
971  
972  AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
973  AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
974  AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
975  if (!track0 || !track1 || !track2) return 0;
976  Double_t prob0[AliPID::kSPECIES];
977  Double_t prob1[AliPID::kSPECIES];
978  Double_t prob2[AliPID::kSPECIES];
979
980  Bool_t isTOF0=fPidHF->CheckTOFPIDStatus((AliAODTrack*)obj->GetDaughter(0));
981  Bool_t isTOF1=fPidHF->CheckTOFPIDStatus((AliAODTrack*)obj->GetDaughter(1));
982  Bool_t isTOF2=fPidHF->CheckTOFPIDStatus((AliAODTrack*)obj->GetDaughter(2));
983
984 Bool_t isK1=kFALSE;
985  if(isTOF1){ //kaon
986   if(track1->P()<1.8) {
987     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
988     if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
989     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
990    
991   }else{
992     AliAODTrack *trackaod1=(AliAODTrack*)(obj->GetDaughter(1));
993     if(trackaod1->P()<0.55){
994       fPidHF->SetTOF(kFALSE);
995       fPidHF->SetTOFdecide(kFALSE);
996      }
997     Int_t isKaon=fPidHF->MakeRawPid(trackaod1,3);
998     if(isKaon>=1) isK1=kTRUE;
999     if(trackaod1->P()<0.55){
1000       fPidHF->SetTOF(kTRUE);
1001       fPidHF->SetTOFdecide(kTRUE);
1002      }
1003   }
1004  }else{
1005   if(track1->P()<0.8){
1006     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1007     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob0);
1008   }else{
1009     AliAODTrack *trackaod1=(AliAODTrack*)(obj->GetDaughter(1));
1010      if(trackaod1->P()<0.55){
1011       fPidHF->SetTOF(kFALSE);
1012       fPidHF->SetTOFdecide(kFALSE);
1013      }
1014     Int_t isKaon=fPidHF->MakeRawPid(trackaod1,3);
1015     if(isKaon>=1) isK1=kTRUE;
1016      if(trackaod1->P()<0.55){
1017       fPidHF->SetTOF(kTRUE);
1018       fPidHF->SetTOFdecide(kTRUE);
1019      }
1020   }
1021  }
1022
1023  Bool_t ispi0=kFALSE;
1024  Bool_t isp0=kFALSE;
1025  Bool_t ispi2=kFALSE;
1026  Bool_t isp2=kFALSE;
1027
1028  if(isTOF0){ //proton
1029   if(track0->P()<2.2) {
1030     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
1031     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1032     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1033   }else{
1034    AliAODTrack *trackaod0=(AliAODTrack*)(obj->GetDaughter(0));
1035    if(trackaod0->P()<1.){
1036       fPidObjprot->SetTOF(kFALSE);
1037       fPidObjprot->SetTOFdecide(kFALSE);
1038    }
1039    Int_t isProton=fPidObjprot->MakeRawPid(trackaod0,4);
1040    if(isProton>=1) isp0=kTRUE;
1041    if(trackaod0->P()<1.){
1042       fPidObjprot->SetTOF(kTRUE);
1043       fPidObjprot->SetTOFdecide(kTRUE);
1044    }
1045   }
1046  }else{
1047    if(track0->P()<1.2){
1048     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1049     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1050    }else{
1051     AliAODTrack *trackaod0=(AliAODTrack*)(obj->GetDaughter(0));
1052     if(trackaod0->P()<1.){
1053       fPidObjprot->SetTOF(kFALSE);
1054       fPidObjprot->SetTOFdecide(kFALSE);
1055     }
1056     Int_t isProton=fPidObjprot->MakeRawPid(trackaod0,4);
1057     if(isProton>=1) isp0=kTRUE;
1058     if(trackaod0->P()<1.){
1059       fPidObjprot->SetTOF(kTRUE);
1060       fPidObjprot->SetTOFdecide(kTRUE);
1061     }
1062    }
1063  }
1064
1065  if(isTOF2){ //proton
1066   if(track2->P()<2.2) {
1067     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
1068     if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1069     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1070   }else{
1071     AliAODTrack *trackaod2=(AliAODTrack*)(obj->GetDaughter(2));
1072     if(trackaod2->P()<1.){
1073       fPidObjprot->SetTOF(kFALSE);
1074       fPidObjprot->SetTOFdecide(kFALSE);
1075     }
1076    Int_t isProton=fPidObjprot->MakeRawPid(trackaod2,4);
1077    if(isProton>=1) isp2=kTRUE;
1078    if(trackaod2->P()<1.){
1079       fPidObjprot->SetTOF(kTRUE);
1080       fPidObjprot->SetTOFdecide(kTRUE);
1081     }
1082   }
1083  }else{
1084    if(track2->P()<1.2){
1085     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1086     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1087    }else{
1088     AliAODTrack *trackaod2=(AliAODTrack*)(obj->GetDaughter(2));
1089      if(trackaod2->P()<1.){
1090       fPidObjprot->SetTOF(kFALSE);
1091       fPidObjprot->SetTOFdecide(kFALSE);
1092     }
1093     Int_t isProton=fPidObjprot->MakeRawPid(trackaod2,4);
1094     if(isProton>=1) isp2=kTRUE;
1095     if(trackaod2->P()<1.){
1096       fPidObjprot->SetTOF(kTRUE);
1097       fPidObjprot->SetTOFdecide(kTRUE);
1098     }
1099    }
1100  }
1101   AliAODTrack *trackaod2=(AliAODTrack*)(obj->GetDaughter(2));
1102   if(fPidObjpion->MakeRawPid(trackaod2,2)>=1)ispi2=kTRUE;
1103   AliAODTrack *trackaod0=(AliAODTrack*)(obj->GetDaughter(2));
1104   if(fPidObjpion->MakeRawPid(trackaod0,2)>=1)ispi0=kTRUE;
1105
1106   if(TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]){
1107
1108     if(TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) okLcpKpi = 1;
1109     if(TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) okLcpiKp = 1;
1110     }
1111
1112    if(!isK1 && TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]) isK1=kTRUE;
1113     if(!ispi0 && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) ispi0=kTRUE;
1114     if(!ispi2 && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) ispi2=kTRUE;
1115     if(!isp0 && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton]) isp0=kTRUE;
1116     if(!isp2 && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton]) isp2=kTRUE;
1117     if(isK1 && ispi0 && isp2) okLcpiKp = 1;
1118     if(isK1 && isp0 && ispi2) okLcpKpi = 1;
1119
1120    if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
1121     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
1122     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
1123
1124     return returnvalue;
1125  
1126
1127 }
1128 //----------------------------------------------------------
1129 Int_t AliRDHFCutsLctopKpi::IsSelectedPIDStrong(AliAODRecoDecayHF* obj) {
1130
1131
1132     if(!fUsePID || !obj) return 3;
1133     Int_t okLcpKpi=0,okLcpiKp=0;
1134     Int_t returnvalue=0;
1135     Bool_t isPeriodd=fPidHF->GetOnePad();
1136     Bool_t isMC=fPidHF->GetMC();
1137     Bool_t ispion0=kTRUE,ispion2=kTRUE;
1138     Bool_t isproton0=kFALSE,isproton2=kFALSE;
1139     Bool_t iskaon1=kFALSE;
1140     if(isPeriodd) {
1141      fPidObjprot->SetOnePad(kTRUE);
1142      fPidObjpion->SetOnePad(kTRUE);
1143     }
1144     if(isMC) {
1145      fPidObjprot->SetMC(kTRUE);
1146      fPidObjpion->SetMC(kTRUE);
1147     }
1148
1149     for(Int_t i=0;i<3;i++){
1150      AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
1151      if(!track) return 0;
1152      // identify kaon
1153      if(i==1) {
1154       Int_t isKaon=fPidHF->MakeRawPid(track,3);
1155       if(isKaon>=1) {
1156        iskaon1=kTRUE;
1157       if(fPidHF->MakeRawPid(track,2)>=1) iskaon1=kFALSE;
1158       }
1159       if(!iskaon1) return 0;
1160      
1161      }else{
1162      //pion or proton
1163      
1164      Int_t isProton=fPidObjprot->MakeRawPid(track,4);
1165      if(isProton>=1){
1166       if(fPidHF->MakeRawPid(track,2)>=1) isProton=-1;
1167       if(fPidHF->MakeRawPid(track,3)>=1) isProton=-1;
1168      }
1169
1170      Int_t isPion=fPidObjpion->MakeRawPid(track,2);
1171      if(fPidHF->MakeRawPid(track,3)>=1) isPion=-1;
1172      if(fPidObjprot->MakeRawPid(track,4)>=1) isPion=-1;
1173
1174
1175      if(i==0) {
1176       if(isPion<0) ispion0=kFALSE;
1177       if(isProton>=1) isproton0=kTRUE;
1178
1179      }
1180       if(!ispion0 && !isproton0) return 0;
1181      if(i==2) {
1182       if(isPion<0) ispion2=kFALSE;
1183       if(isProton>=1) isproton2=kTRUE;
1184      }
1185
1186     }
1187    }
1188
1189     if(ispion2 && isproton0 && iskaon1) okLcpKpi=1;
1190     if(ispion0 && isproton2 && iskaon1) okLcpiKp=1;
1191     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
1192     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
1193     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
1194
1195  return returnvalue;
1196 }
1197 //--------------------
1198 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPIDpPb(AliAODRecoDecayHF* obj) {
1199     
1200     if(!fUsePID || !obj) {return 3;}
1201     Int_t okLcpKpi=0,okLcpiKp=0;
1202     Int_t returnvalue=0;
1203    
1204     Bool_t isMC=fPidHF->GetMC();
1205
1206     
1207     if(isMC) {
1208             fPidObjprot->SetMC(kTRUE);
1209             fPidObjpion->SetMC(kTRUE);
1210     }
1211
1212     AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
1213     AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
1214     AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
1215     if (!track0 || !track1 || !track2) return 0;
1216     Double_t prob0[AliPID::kSPECIES];
1217     Double_t prob1[AliPID::kSPECIES];
1218     Double_t prob2[AliPID::kSPECIES];
1219     
1220     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1221     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
1222     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1223   
1224
1225     if(fPIDThreshold[AliPID::kPion]>0. && fPIDThreshold[AliPID::kKaon]>0. && fPIDThreshold[AliPID::kProton]>0.){
1226     okLcpiKp=  (prob0[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ])
1227              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
1228              &&(prob2[AliPID::kProton]>fPIDThreshold[AliPID::kProton]);
1229     okLcpKpi=  (prob0[AliPID::kProton]>fPIDThreshold[AliPID::kProton])
1230              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
1231              &&(prob2[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ]);
1232    }else{ 
1233                     //pion or proton
1234                     
1235                     
1236     if(TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]){
1237     if(TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) okLcpKpi = 1;  
1238     if(TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) okLcpiKp = 1; 
1239             }
1240            }
1241     
1242     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
1243     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
1244     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
1245     
1246     return returnvalue;
1247 }
1248 //-----------------------
1249 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPIDpPb2(AliAODRecoDecayHF* obj) {
1250
1251   Int_t returnvalue =0;
1252   Double_t thresholdK =0.7;
1253   Double_t thresholdPi =0.3; //!
1254   Double_t thresholdPr =0.7;
1255
1256   AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
1257   AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
1258   AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
1259
1260   if (!track0 || !track1 || !track2) return 0;
1261   Double_t prob0[AliPID::kSPECIES];
1262   Double_t prob1[AliPID::kSPECIES];
1263   Double_t prob2[AliPID::kSPECIES];
1264
1265   fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1266   fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
1267   fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1268
1269   if(prob1[AliPID::kKaon]>thresholdK){
1270     if(TMath::MaxElement(AliPID::kSPECIES,prob0)>TMath::MaxElement(AliPID::kSPECIES,prob2)){
1271       if(((prob0[AliPID::kPion  ]>prob0[AliPID::kProton  ])&& prob0[AliPID::kPion]>thresholdPi) && prob2[AliPID::kProton]>thresholdPr) returnvalue=2;//piKp     
1272       else if(((prob0[AliPID::kProton  ]>prob0[AliPID::kPion  ])&& prob0[AliPID::kProton]>thresholdPr) && prob2[AliPID::kPion]>thresholdPi)returnvalue =1;//pKpi
1273     }
1274
1275     else if(TMath::MaxElement(AliPID::kSPECIES,prob0)<TMath::MaxElement(AliPID::kSPECIES,prob2)){
1276       if(((prob2[AliPID::kPion  ]>prob2[AliPID::kProton  ])&& prob2[AliPID::kPion]>thresholdPi) && prob0[AliPID::kProton]>thresholdPr) returnvalue=1; //pKpi    
1277       else if(((prob2[AliPID::kProton  ]>prob2[AliPID::kPion  ])&& prob2[AliPID::kProton]>thresholdPr) && prob0[AliPID::kPion]>thresholdPi)returnvalue =2;      //piKp
1278     }
1279
1280   }
1281   return returnvalue;
1282
1283 }
1284 //------------------------------------------------------
1285 void AliRDHFCutsLctopKpi::SetStandardCutsPPb2013() {
1286
1287  SetName("LctopKpiProdCuts");
1288  SetTitle("Production cuts for Lc analysis");
1289
1290  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1291  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1292  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1293  esdTrackCuts->SetMinNClustersTPC(70);
1294  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1295                                           AliESDtrackCuts::kAny);
1296  esdTrackCuts->SetRequireITSRefit(kTRUE);
1297  esdTrackCuts->SetMinNClustersITS(4);
1298  esdTrackCuts->SetMinDCAToVertexXY(0.);
1299  esdTrackCuts->SetEtaRange(-0.8,0.8);
1300  esdTrackCuts->SetPtRange(0.3,1.e10);
1301  AddTrackCuts(esdTrackCuts);
1302  
1303  const Int_t nvars=13;
1304  const Int_t nptbins=9;
1305  Float_t* ptbins;
1306  ptbins=new Float_t[nptbins+1];
1307  ptbins[0]=0.;
1308  ptbins[1]=1.;
1309  ptbins[2]=2.;
1310  ptbins[3]=3.;
1311  ptbins[4]=4.;
1312  ptbins[5]=5.;
1313  ptbins[6]=6.;
1314  ptbins[7]=8.;
1315  ptbins[8]=10.;  
1316  ptbins[9]=99999.;
1317
1318  SetGlobalIndex(nvars,nptbins);
1319  SetPtBins(nptbins+1,ptbins);
1320
1321  Float_t** prodcutsval;
1322  prodcutsval=new Float_t*[nvars];
1323  for(Int_t iv=0;iv<nvars;iv++){
1324   prodcutsval[iv]=new Float_t[nptbins];
1325  }
1326
1327  for(Int_t ipt=0;ipt<nptbins;ipt++){
1328     prodcutsval[0][ipt]=0.13;
1329     prodcutsval[1][ipt]=0.4;
1330     prodcutsval[2][ipt]=0.4;
1331     prodcutsval[3][ipt]=0.;
1332     prodcutsval[4][ipt]=0.;
1333     prodcutsval[5][ipt]=0.;
1334     prodcutsval[6][ipt]=0.06;
1335     prodcutsval[7][ipt]=0.;
1336     prodcutsval[8][ipt]=0.;
1337     prodcutsval[9][ipt]=0.;
1338     prodcutsval[10][ipt]=0.;
1339     prodcutsval[11][ipt]=0.05;
1340     prodcutsval[12][ipt]=0.4;
1341  }
1342  SetCuts(nvars,nptbins,prodcutsval);
1343
1344  AliAODPidHF* pidObjK=new AliAODPidHF();
1345  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
1346  pidObjK->SetSigma(sigmasK);
1347  pidObjK->SetAsym(kTRUE);
1348  pidObjK->SetMatch(1);
1349  pidObjK->SetTPC(kTRUE);
1350  pidObjK->SetTOF(kTRUE);
1351  pidObjK->SetITS(kTRUE);
1352  Double_t plimK[2]={0.5,0.8};
1353  pidObjK->SetPLimit(plimK,2);
1354  pidObjK->SetTOFdecide(kTRUE);
1355  SetPidHF(pidObjK);
1356
1357  AliAODPidHF* pidObjpi=new AliAODPidHF();
1358  pidObjpi->SetTPC(kTRUE);
1359  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
1360  pidObjpi->SetSigma(sigmaspi);
1361  pidObjpi->SetTOFdecide(kTRUE);
1362  SetPidpion(pidObjpi);
1363
1364  AliAODPidHF* pidObjp=new AliAODPidHF();
1365  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
1366  pidObjp->SetSigma(sigmasp);
1367  pidObjp->SetAsym(kTRUE);
1368  pidObjp->SetMatch(1);
1369  pidObjp->SetTPC(kTRUE);
1370  pidObjp->SetTOF(kTRUE);
1371  pidObjp->SetITS(kTRUE);
1372  Double_t plimp[2]={1.,2.};
1373  pidObjp->SetPLimit(plimp,2);
1374  pidObjp->SetTOFdecide(kTRUE);
1375
1376  SetUsePID(kTRUE);
1377
1378  PrintAll();
1379
1380  for(Int_t iiv=0;iiv<nvars;iiv++){
1381   delete [] prodcutsval[iiv];
1382  }
1383  delete [] prodcutsval;
1384  prodcutsval=NULL;
1385  delete [] ptbins;
1386  ptbins=NULL;
1387
1388  delete pidObjK;
1389  pidObjK=NULL;
1390  delete pidObjpi;
1391  pidObjpi=NULL;
1392  delete pidObjp;
1393  pidObjp=NULL;
1394
1395  return;
1396 }
1397