dd86c11ce3120f4b08d843c8ccdbbcccd48cc556
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCutsDstoKKpi.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 Ds->KKpi
21 //
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <TDatabasePDG.h>
26 #include <Riostream.h>
27
28 #include "AliRDHFCutsDstoKKpi.h"
29 #include "AliAODRecoDecayHF3Prong.h"
30 #include "AliAODTrack.h"
31 #include "AliESDtrack.h"
32
33 using std::cout;
34 using std::endl;
35
36 ClassImp(AliRDHFCutsDstoKKpi)
37
38 //--------------------------------------------------------------------------
39 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const char* name) : 
40 AliRDHFCuts(name),
41 fPidOption(0),
42 fMaxPtStrongPid(0.),
43 fMaxPStrongPidK(0.),
44 fMaxPStrongPidpi(0.)
45 {
46   //
47   // Default Constructor
48   //
49   Int_t nvars=20;
50   SetNVars(nvars);
51   TString varNames[20]={"inv. mass [GeV]",
52                         "pTK [GeV/c]",
53                         "pTPi [GeV/c]",
54                         "d0K [cm]",
55                         "d0Pi [cm]",
56                         "dist12 [cm]",
57                         "sigmavert [cm]",
58                         "decLen [cm]",
59                         "ptMax [GeV/c]",
60                         "cosThetaPoint",
61                         "Sum d0^2 (cm^2)",
62                         "dca [cm]",
63                         "inv. mass (Mphi-MKK) [GeV]",
64                         "inv. mass (MKo*-MKpi) [GeV]",
65                         "Abs(CosineKpiPhiRFrame)^3",
66                         "CosPiDsLabFrame",
67                         "decLenXY [cm]"
68                         "NormdecLen",
69                         "NormdecLenXY [cm]",
70                         "cosThetaPointXY"};
71                         
72   Bool_t isUpperCut[20]={kTRUE,
73                          kFALSE,
74                          kFALSE,
75                          kFALSE,
76                          kFALSE,
77                          kFALSE,
78                          kTRUE,
79                          kFALSE,
80                          kFALSE,
81                          kFALSE,
82                          kFALSE,
83                          kTRUE,
84                          kTRUE,
85                          kTRUE,
86                          kFALSE,
87                          kTRUE,
88                          kFALSE,
89                          kFALSE,
90                          kFALSE,
91                          kFALSE};
92   SetVarNames(20,varNames,isUpperCut);
93   Bool_t forOpt[20]={kFALSE,
94                     kFALSE,
95                     kFALSE,
96                     kFALSE,
97                     kFALSE,
98                     kFALSE,
99                     kTRUE,
100                     kTRUE,
101                     kTRUE,
102                     kTRUE,
103                     kTRUE,
104                     kFALSE,
105                     kTRUE,
106                     kTRUE,
107                     kFALSE,
108                     kFALSE,
109                     kTRUE,
110                     kTRUE,
111                     kTRUE,
112                     kTRUE};
113                     
114   SetVarsForOpt(11,forOpt);
115   Float_t limits[2]={0,999999999.};
116   SetPtBins(2,limits);
117   if(fPidHF)delete fPidHF;
118   fPidHF=new AliAODPidHF();
119   Double_t plim[2]={0.6,0.8};
120   Double_t nsigma[5]={2.,1.,2.,3.,0.};
121   
122   fPidHF->SetPLimit(plim);
123   fPidHF->SetAsym(kTRUE);
124   fPidHF->SetSigma(nsigma);
125   fPidHF->SetMatch(1);
126   fPidHF->SetTPC(1);
127   fPidHF->SetTOF(1);
128   fPidHF->SetITS(0);
129   fPidHF->SetTRD(0);
130   fPidHF->SetCompat(kTRUE);
131
132 }
133 //--------------------------------------------------------------------------
134 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const AliRDHFCutsDstoKKpi &source) :
135   AliRDHFCuts(source),
136   fPidOption(source.fPidOption),
137   fMaxPtStrongPid(source.fMaxPtStrongPid),
138   fMaxPStrongPidK(source.fMaxPStrongPidK),
139   fMaxPStrongPidpi(source.fMaxPStrongPidpi)
140 {
141   //
142   // Copy constructor
143   //
144
145 }
146 //--------------------------------------------------------------------------
147 AliRDHFCutsDstoKKpi &AliRDHFCutsDstoKKpi::operator=(const AliRDHFCutsDstoKKpi &source)
148 {
149   //
150   // assignment operator
151   //
152   if(&source == this) return *this;
153
154   AliRDHFCuts::operator=(source);
155
156   fPidOption=source.fPidOption;
157   fMaxPtStrongPid=source.fMaxPtStrongPid;
158   fMaxPStrongPidK=source.fMaxPStrongPidK;
159   fMaxPStrongPidpi=source.fMaxPStrongPidpi;
160
161   return *this;
162 }
163
164
165 //---------------------------------------------------------------------------
166 void AliRDHFCutsDstoKKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
167   // 
168   // Fills in vars the values of the variables 
169   //
170
171   if(nvars!=fnVarsForOpt) {
172     printf("AliRDHFCutsDstoKKpi::GetCutsVarsForOpt: wrong number of variables\n");
173     return;
174   }
175
176   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
177   
178   //recalculate vertex w/o daughters
179   Bool_t cleanvtx=kFALSE;
180   AliAODVertex *origownvtx=0x0;
181   if(fRemoveDaughtersFromPrimary) {
182     if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
183     cleanvtx=kTRUE;
184     if(!RecalcOwnPrimaryVtx(dd,aod)) {
185       CleanOwnPrimaryVtx(dd,aod,origownvtx);
186       cleanvtx=kFALSE;
187     }
188   }
189
190   Int_t iter=-1;
191   if(fVarsForOpt[0]){
192     iter++;
193     if(TMath::Abs(pdgdaughters[0])==321){
194       vars[iter]=dd->InvMassDsKKpi();
195     }else{
196       vars[iter]=dd->InvMassDspiKK();
197     }
198   }
199   if(fVarsForOpt[1]){
200     iter++;
201     Float_t minPtDau=99999.;
202     for(Int_t iprong=0;iprong<3;iprong++){
203       if(TMath::Abs(pdgdaughters[iprong])==321 && 
204          dd->PtProng(iprong)<minPtDau) minPtDau=dd->PtProng(iprong);
205     }
206     vars[iter]=minPtDau;
207   }
208   if(fVarsForOpt[2]){
209     iter++;
210     for(Int_t iprong=0;iprong<3;iprong++){
211       if(TMath::Abs(pdgdaughters[iprong])==211) {
212         vars[iter]=dd->PtProng(iprong);
213       }
214     }
215   }
216   if(fVarsForOpt[3]){
217     iter++;
218     Float_t minImpParDau=99999.;
219     for(Int_t iprong=0;iprong<3;iprong++){
220       if(TMath::Abs(pdgdaughters[iprong])==321 &&
221          dd->Getd0Prong(iprong)<minImpParDau) minImpParDau=dd->Getd0Prong(iprong);
222     }
223     vars[iter]=minImpParDau;
224   }
225   if(fVarsForOpt[4]){
226     iter++;
227     for(Int_t iprong=0;iprong<3;iprong++){
228       if(TMath::Abs(pdgdaughters[iprong])==211) {
229         vars[iter]=dd->Getd0Prong(iprong);
230       }
231     }
232   }
233   if(fVarsForOpt[5]){
234     iter++;
235     Float_t minDistPair=TMath::Min(dd->GetDist12toPrim(),dd->GetDist23toPrim());
236     vars[iter]=minDistPair;
237   }
238   if(fVarsForOpt[6]){
239     iter++;
240     vars[iter]=dd->GetSigmaVert(aod);
241   }
242   if(fVarsForOpt[7]){
243     iter++;
244     vars[iter] = dd->DecayLength();
245   }
246   if(fVarsForOpt[8]){
247     iter++;
248     Float_t ptmax=0;
249     for(Int_t i=0;i<3;i++){
250       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
251     }
252     vars[iter]=ptmax;
253   }
254   if(fVarsForOpt[9]){
255     iter++;
256     vars[iter]=dd->CosPointingAngle();
257   }
258   if(fVarsForOpt[10]){
259     iter++;
260     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
261   }
262   if(fVarsForOpt[11]){
263     iter++;
264     Float_t maxDCA=0.;
265     for(Int_t i=0;i<3;i++){ 
266       if(d->GetDCA(i)>maxDCA) maxDCA=d->GetDCA(i);
267     }
268     vars[iter]=maxDCA;
269   }
270   if(fVarsForOpt[12]){
271     iter++;
272     Double_t mPDGPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
273     if(TMath::Abs(pdgdaughters[0])==321){
274       
275       Double_t phimass01=d->InvMass2Prongs(0,1,321,321);
276        vars[iter]=TMath::Abs(phimass01-mPDGPhi);
277        // vars[iter]=dd->InvMass2Prongs(0,1,321,321);
278     }else{
279       Double_t phimass12=d->InvMass2Prongs(1,2,321,321);
280        vars[iter]=TMath::Abs(phimass12-mPDGPhi);
281        // vars[iter]=dd->InvMass2Prongs(1,2,321,321);      
282     }
283   }
284   if(fVarsForOpt[13]){
285     iter++;
286     Double_t mPDGK0star = TDatabasePDG::Instance()->GetParticle(313)->Mass();
287     if(TMath::Abs(pdgdaughters[0])==321){
288       
289       Double_t mass12kpi=d->InvMass2Prongs(1,2,321,211);
290       vars[iter]=TMath::Abs(mass12kpi-mPDGK0star);
291       //              vars[iter]=dd->InvMass2Prongs(1,2,321,211);
292     }else{
293       Double_t mass01pik=d->InvMass2Prongs(0,1,211,321);
294       vars[iter]=TMath::Abs(mass01pik-mPDGK0star);
295       //        vars[iter]=dd->InvMass2Prongs(0,1,211,321);      
296     }
297   }
298   if(fVarsForOpt[14]){
299     iter++;
300     if(TMath::Abs(pdgdaughters[0])==321){
301       vars[iter]=dd->CosPiKPhiRFrameKKpi();
302     }else{
303       vars[iter]=dd->CosPiKPhiRFramepiKK();
304     }
305   }
306   if(fVarsForOpt[15]){
307     iter++;
308     if(TMath::Abs(pdgdaughters[0])==321){
309       vars[iter]=dd->CosPiDsLabFrameKKpi();
310     }else{
311       vars[iter]=dd->CosPiDsLabFramepiKK();
312     }
313   }
314   
315   if(fVarsForOpt[16]){
316     iter++;
317     vars[iter]=dd->DecayLengthXY();
318   }
319   
320   if(fVarsForOpt[17]){
321     iter++;
322     vars[iter]=dd->NormalizedDecayLength();
323   }
324   
325   if(fVarsForOpt[18]){
326     iter++;
327     vars[iter]=dd->NormalizedDecayLengthXY();
328   }
329   
330   if(fVarsForOpt[19]){
331     iter++;
332     vars[iter]=dd->CosPointingAngleXY();
333   }
334
335   if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx); 
336   return;
337 }
338 //---------------------------------------------------------------------------
339 Bool_t AliRDHFCutsDstoKKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
340 {
341   //
342   // Checking if Ds is in fiducial acceptance region 
343   //
344
345   if(pt > 5.) {
346     // applying cut for pt > 5 GeV
347     AliDebug(2,Form("pt of Ds = %f (> 5), cutting at |y| < 0.8",pt)); 
348     if (TMath::Abs(y) > 0.8) return kFALSE;
349     
350   } else {
351     // appliying smooth cut for pt < 5 GeV
352     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
353     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;         
354     AliDebug(2,Form("pt of Ds = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
355     if (y < minFiducialY || y > maxFiducialY) return kFALSE;    
356   }
357
358   return kTRUE;
359 }
360
361 //---------------------------------------------------------------------------
362 Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
363   // PID selection
364   // return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both 
365   Int_t retCode=3;
366   Bool_t okKKpi=kTRUE;
367   Bool_t okpiKK=kTRUE;
368   if(!fUsePID || !rd) return retCode;
369   if(!fPidHF){
370     AliWarning("AliAODPidHF not created!");
371     return retCode;
372   }
373
374   Double_t origCompatTOF=fPidHF->GetPCompatTOF();
375   Double_t origThreshTPC=fPidHF->GetPtThresholdTPC();
376   if(fPidOption==kStrong){
377     fPidHF->SetPCompatTOF(999999.);
378     fPidHF->SetPtThresholdTPC(999999.);
379   }
380
381   Int_t nKaons=0;
382   Int_t nNotKaons=0;
383   Int_t sign= rd->GetCharge(); 
384   for(Int_t iDaught=0; iDaught<3; iDaught++){
385     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(iDaught);
386     
387     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
388     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
389     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
390     
391     if(isProton>0 &&  isKaon<0  && isPion<0){
392       fPidHF->SetPCompatTOF(origCompatTOF);
393       fPidHF->SetPtThresholdTPC(origThreshTPC);
394       return 0;
395     }
396     if(sign!=track->Charge()){// must be kaon
397       if(isKaon<0){
398         fPidHF->SetPCompatTOF(origCompatTOF);
399         fPidHF->SetPtThresholdTPC(origThreshTPC);
400         return 0;
401       }
402       if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid && isKaon<=0){
403         fPidHF->SetPCompatTOF(origCompatTOF);
404         fPidHF->SetPtThresholdTPC(origThreshTPC);
405         return 0;
406       }
407       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
408         if(isKaon<=0 && track->P()<fMaxPStrongPidK) return 0;
409       }
410     }
411     
412     if(isKaon>0 && isPion<0) nKaons++;
413     if(isKaon<0) nNotKaons++;
414     if(iDaught==0){
415       if(isKaon<0) okKKpi=kFALSE;
416       if(isPion<0) okpiKK=kFALSE;      
417       if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid){
418         if(isKaon<=0) okKKpi=kFALSE;
419         if(isPion<=0) okpiKK=kFALSE;
420       }
421       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
422         if(isKaon<=0 && track->P()<fMaxPStrongPidK) okKKpi=kFALSE;
423         if(isPion<=0 && track->P()<fMaxPStrongPidpi) okpiKK=kFALSE;
424       }
425     }
426     else if(iDaught==2){
427       if(isKaon<0) okpiKK=kFALSE;
428       if(isPion<0) okKKpi=kFALSE;
429        if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid){
430          if(isKaon<=0) okpiKK=kFALSE;
431          if(isPion<=0) okKKpi=kFALSE;
432       }
433       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
434         if(isKaon<=0 && track->P()<fMaxPStrongPidK) okpiKK=kFALSE;  
435         if(isPion<=0 && track->P()<fMaxPStrongPidpi) okKKpi=kFALSE; 
436       }
437     }
438   }
439
440   fPidHF->SetPCompatTOF(origCompatTOF);
441   fPidHF->SetPtThresholdTPC(origThreshTPC);
442   
443   if(nKaons>2)return 0;
444   if(nNotKaons>1) return 0;
445   
446   if(!okKKpi) retCode-=1;
447   if(!okpiKK) retCode-=2;
448
449   return retCode;
450 }
451
452 //---------------------------------------------------------------------------
453 Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
454   //
455   // Apply selection
456   //
457
458   if(!fCutsRD){
459     cout<<"Cut matrix not inizialized. Exit..."<<endl;
460     return 0;
461   }
462   //PrintAll();
463   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
464
465   if(!d){
466     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
467     return 0;
468   }
469   
470   if(fKeepSignalMC) if(IsSignalMC(d,aod,431)) return 3;
471  
472   Double_t ptD=d->Pt();
473   if(ptD<fMinPtCand) return 0;
474   if(ptD>fMaxPtCand) return 0;
475
476   if(d->HasBadDaughters()) return 0;
477   
478
479   // selection on daughter tracks 
480   if(selectionLevel==AliRDHFCuts::kAll || 
481      selectionLevel==AliRDHFCuts::kTracks) {
482     if(!AreDaughtersSelected(d)) return 0;
483   }
484
485
486
487  
488   // selection on candidate
489   if(selectionLevel==AliRDHFCuts::kAll || 
490      selectionLevel==AliRDHFCuts::kCandidate) {
491     //recalculate vertex w/o daughters
492     AliAODVertex *origownvtx=0x0;
493     if(fRemoveDaughtersFromPrimary) {
494       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
495       if(!RecalcOwnPrimaryVtx(d,aod)) {
496         CleanOwnPrimaryVtx(d,aod,origownvtx);
497         return 0;
498       }
499     }
500
501     Int_t okDsKKpi=1;
502     Int_t okDspiKK=1;
503     Int_t okMassPhiKKpi=0;
504     Int_t okMassPhipiKK=0;
505     Int_t okMassK0starKKpi=0;
506     Int_t okMassK0starpiKK=0;
507     Int_t okDsPhiKKpi=0;
508     Int_t okDsPhipiKK=0;
509     Int_t okDsK0starKKpi=0;
510     Int_t okDsK0starpiKK=0;
511
512     Double_t pt=d->Pt();
513     Int_t ptbin=PtBin(pt);
514     if (ptbin==-1) {
515       CleanOwnPrimaryVtx(d,aod,origownvtx);
516       return 0;
517     }
518  
519     Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
520     Double_t mDsKKpi=d->InvMassDsKKpi();
521     Double_t mDspiKK=d->InvMassDspiKK();
522     if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
523     if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
524     if(!okDsKKpi && !okDspiKK){
525       CleanOwnPrimaryVtx(d,aod,origownvtx);
526       return 0;
527     }
528
529
530
531     // cuts on resonant decays (via Phi or K0*)
532     Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
533     Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
534     if(okDsKKpi){
535       Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
536       Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
537       if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhiKKpi=1;
538       if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starKKpi = 1;
539       if(!okMassPhiKKpi && !okMassK0starKKpi) okDsKKpi=0;
540       if(okMassPhiKKpi) okDsPhiKKpi=1;
541       if(okMassK0starKKpi) okDsK0starKKpi=1;
542     }
543     if(okDspiKK){
544       Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
545       Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
546       if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starpiKK = 1;
547       if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhipiKK=1;
548       if(!okMassPhipiKK && !okMassK0starpiKK) okDspiKK=0;
549       if(okMassPhipiKK) okDsPhipiKK=1;
550       if(okMassK0starpiKK) okDsK0starpiKK=1;
551     }
552     if(!okDsKKpi && !okDspiKK){
553       CleanOwnPrimaryVtx(d,aod,origownvtx);
554       return 0;
555     }
556
557     // Cuts on track pairs
558     for(Int_t i=0;i<3;i++){
559       if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]){
560         CleanOwnPrimaryVtx(d,aod,origownvtx);
561         return 0;
562       }
563     }
564     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] || 
565        d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]){
566       CleanOwnPrimaryVtx(d,aod,origownvtx);
567       return 0;
568     }
569
570
571
572     //single track
573     if(TMath::Abs(d->Pt2Prong(1)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] || 
574        TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]){
575       CleanOwnPrimaryVtx(d,aod,origownvtx);
576       return 0;
577     }
578
579     if(okDsKKpi){
580       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] || 
581          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
582       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] || 
583          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
584     }
585     if(okDspiKK){
586       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] || 
587          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
588       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] || 
589          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
590     }
591     if(!okDsKKpi && !okDspiKK){
592       CleanOwnPrimaryVtx(d,aod,origownvtx);
593       return 0;
594     }
595
596     // Cuts on candidate triplet
597
598
599     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]){
600       CleanOwnPrimaryVtx(d,aod,origownvtx); 
601       return 0;
602     }
603      
604     if(d->Pt2Prong(0)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] && 
605        d->Pt2Prong(1)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] && 
606        d->Pt2Prong(2)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)]) {
607       CleanOwnPrimaryVtx(d,aod,origownvtx); 
608       return 0;
609     }
610
611     if(d->DecayLength2()<fCutsRD[GetGlobalIndex(7,ptbin)]*fCutsRD[GetGlobalIndex(7,ptbin)]){
612       CleanOwnPrimaryVtx(d,aod,origownvtx);
613       return 0;
614     }
615
616
617     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
618     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)]){
619       CleanOwnPrimaryVtx(d,aod,origownvtx);
620       return 0;
621     }
622
623    
624     //sec vert
625     Double_t sigmavert=d->GetSigmaVert(aod);
626     if(sigmavert>fCutsRD[GetGlobalIndex(6,ptbin)]){
627       CleanOwnPrimaryVtx(d,aod,origownvtx);
628       return 0;
629     }
630     
631     // decay length XY
632     if(d->DecayLengthXY()<fCutsRD[GetGlobalIndex(16,ptbin)]){
633       CleanOwnPrimaryVtx(d,aod,origownvtx); 
634       return 0;
635     }
636     
637     //norm decay length
638     if(d->NormalizedDecayLength()<fCutsRD[GetGlobalIndex(17,ptbin)]){
639       CleanOwnPrimaryVtx(d,aod,origownvtx); 
640       return 0;
641     }
642     
643     //norm decay length XY
644     if(d->NormalizedDecayLengthXY()<fCutsRD[GetGlobalIndex(18,ptbin)]){
645       CleanOwnPrimaryVtx(d,aod,origownvtx); 
646       return 0;
647     }
648     
649     //cos pointing XY
650     if(d->CosPointingAngleXY()<fCutsRD[GetGlobalIndex(19,ptbin)]){
651       CleanOwnPrimaryVtx(d,aod,origownvtx); 
652       return 0;
653     }
654     
655
656     if(okDsKKpi){
657       Double_t cosPiKPhiRFKKpi=d->CosPiKPhiRFrameKKpi();
658       Double_t kincutPiKPhiKKpi=TMath::Abs(cosPiKPhiRFKKpi*cosPiKPhiRFKKpi*cosPiKPhiRFKKpi);
659       if(kincutPiKPhiKKpi<fCutsRD[GetGlobalIndex(14,ptbin)]) okDsKKpi=0;
660     }
661     if(okDspiKK){
662       Double_t cosPiKPhiRFpiKK=d->CosPiKPhiRFramepiKK();
663       Double_t kincutPiKPhipiKK=TMath::Abs(cosPiKPhiRFpiKK*cosPiKPhiRFpiKK*cosPiKPhiRFpiKK);
664       if(kincutPiKPhipiKK<fCutsRD[GetGlobalIndex(14,ptbin)]) okDspiKK=0;
665     }
666     if(!okDsKKpi && !okDspiKK){
667       CleanOwnPrimaryVtx(d,aod,origownvtx);
668       return 0;
669     }
670     
671     
672     
673     if(okDsKKpi){
674       Double_t cosPiDsLabFrameKKpi=d->CosPiDsLabFrameKKpi();
675       if(cosPiDsLabFrameKKpi>fCutsRD[GetGlobalIndex(15,ptbin)]) okDsKKpi=0;
676     }
677     if(okDspiKK){
678       Double_t cosPiDsLabFramepiKK=d->CosPiDsLabFramepiKK();
679       if(cosPiDsLabFramepiKK>fCutsRD[GetGlobalIndex(15,ptbin)]) okDspiKK=0;
680     }
681     if(!okDsKKpi && !okDspiKK){
682       CleanOwnPrimaryVtx(d,aod,origownvtx);
683       return 0;
684     }
685     
686      // unset recalculated primary vertex when not needed any more
687     CleanOwnPrimaryVtx(d,aod,origownvtx);
688       
689     
690
691     if(!okDsKKpi){
692       okDsPhiKKpi=0;
693       okDsK0starKKpi=0;
694     }
695     if(!okDspiKK){
696       okDsPhipiKK=0;
697       okDsK0starpiKK=0;
698     }
699
700     // PID selection
701     Int_t returnvaluePID=3;  
702     if(selectionLevel==AliRDHFCuts::kAll || 
703        selectionLevel==AliRDHFCuts::kCandidate ||     
704        selectionLevel==AliRDHFCuts::kPID) {
705       returnvaluePID = IsSelectedPID(d);
706       fIsSelectedPID=returnvaluePID;
707     }
708     if(returnvaluePID==0)return 0;
709
710     Bool_t okPidDsKKpi=returnvaluePID&1;
711     Bool_t okPidDspiKK=returnvaluePID&2;
712     if(!okPidDsKKpi){
713       okDsPhiKKpi=0;
714       okDsK0starKKpi=0;
715     }
716     if(!okPidDspiKK){
717       okDsPhipiKK=0;
718       okDsK0starpiKK=0;
719     }
720
721     if((okPidDsKKpi && okDsKKpi)||(okPidDspiKK && okDspiKK)){
722       Int_t returnvalue=0;
723       if(okDsKKpi) returnvalue+=1;
724       if(okDspiKK) returnvalue+=2;
725       if(okDsPhiKKpi) returnvalue+=4;
726       if(okDsPhipiKK) returnvalue+=8;
727       if(okDsK0starKKpi) returnvalue+=16;
728       if(okDsK0starpiKK) returnvalue+=32;
729       return returnvalue;
730     }else{
731       return 0;
732     }
733   }
734   return 15;
735
736 }
737
738 //--------------------------------------------------------------------------
739
740 UInt_t AliRDHFCutsDstoKKpi::GetPIDTrackTPCTOFBitMap(AliAODTrack *track) const{
741
742   UInt_t bitmap=0;
743
744   Double_t sigmaTPCPionHyp=-999.;
745   Double_t sigmaTPCKaonHyp=-999.;
746   Double_t sigmaTPCProtonHyp=-999.;
747   Double_t sigmaTOFPionHyp=-999.;
748   Double_t sigmaTOFKaonHyp=-999.;
749   Double_t sigmaTOFProtonHyp=-999.;
750   
751   Int_t oksigmaTPCPionHyp=fPidHF->GetnSigmaTPC(track,2,sigmaTPCPionHyp);
752   Int_t oksigmaTPCKaonHyp=fPidHF->GetnSigmaTPC(track,3,sigmaTPCKaonHyp);
753   Int_t oksigmaTPCProtonHyp=fPidHF->GetnSigmaTPC(track,4,sigmaTPCProtonHyp);
754   Int_t oksigmaTOFPionHyp=fPidHF->GetnSigmaTOF(track,2,sigmaTOFPionHyp);
755   Int_t oksigmaTOFKaonHyp=fPidHF->GetnSigmaTOF(track,3,sigmaTOFKaonHyp);
756   Int_t oksigmaTOFProtonHyp=fPidHF->GetnSigmaTOF(track,4,sigmaTOFProtonHyp);
757   
758   sigmaTPCPionHyp=TMath::Abs(sigmaTPCPionHyp);
759   sigmaTPCKaonHyp=TMath::Abs(sigmaTPCKaonHyp);
760   sigmaTPCProtonHyp=TMath::Abs(sigmaTPCProtonHyp);
761   sigmaTOFPionHyp=TMath::Abs(sigmaTOFPionHyp);
762   sigmaTOFKaonHyp=TMath::Abs(sigmaTOFKaonHyp);
763   sigmaTOFProtonHyp=TMath::Abs(sigmaTOFProtonHyp);
764
765  if (oksigmaTPCPionHyp && sigmaTPCPionHyp>0.){
766     if (sigmaTPCPionHyp<1.) bitmap+=1<<kTPCPionLess1;
767     else{
768       if (sigmaTPCPionHyp<2.) bitmap+=1<<kTPCPionMore1Less2;
769       else { 
770         if (sigmaTPCPionHyp<3.) bitmap+=1<<kTPCPionMore2Less3; 
771         else bitmap+=1<<kTPCPionMore3;
772       }
773     }
774   }
775   
776   if (oksigmaTPCKaonHyp && sigmaTPCKaonHyp>0.){
777     if (sigmaTPCKaonHyp<1.) bitmap+=1<<kTPCKaonLess1;
778     else{
779       if (sigmaTPCKaonHyp<2.) bitmap+=1<<kTPCKaonMore1Less2;
780       else { 
781         if (sigmaTPCKaonHyp<3.) bitmap+=1<<kTPCKaonMore2Less3; 
782         else bitmap+=1<<kTPCKaonMore3;
783       }
784     }
785   }
786   
787   if (oksigmaTPCProtonHyp && sigmaTPCProtonHyp>0.){
788     if (sigmaTPCProtonHyp<1.) bitmap+=1<<kTPCProtonLess1;
789     else{
790       if (sigmaTPCProtonHyp<2.) bitmap+=1<<kTPCProtonMore1Less2;
791       else { 
792         if (sigmaTPCProtonHyp<3.) bitmap+=1<<kTPCProtonMore2Less3; 
793         else bitmap+=1<<kTPCProtonMore3;
794       }
795     }
796   }
797   
798   if (oksigmaTOFPionHyp && sigmaTOFPionHyp>0.){
799     if (sigmaTOFPionHyp<1.) bitmap+=1<<kTOFPionLess1;
800     else{
801       if (sigmaTOFPionHyp<2.) bitmap+=1<<kTOFPionMore1Less2;
802       else { 
803         if (sigmaTOFPionHyp<3.) bitmap+=1<<kTOFPionMore2Less3; 
804         else bitmap+=1<<kTOFPionMore3;
805       }
806     }
807   }
808   
809   if (oksigmaTOFKaonHyp && sigmaTOFKaonHyp>0.){
810     if (sigmaTOFKaonHyp<1.) bitmap+=1<<kTOFKaonLess1;
811     else{
812       if (sigmaTOFKaonHyp<2.) bitmap+=1<<kTOFKaonMore1Less2;
813       else { 
814         if (sigmaTOFKaonHyp<3.) bitmap+=1<<kTOFKaonMore2Less3; 
815         else bitmap+=1<<kTOFKaonMore3;
816       }
817     }
818   }
819   
820   if (oksigmaTOFProtonHyp && sigmaTOFProtonHyp>0.){
821     if (sigmaTOFProtonHyp<1.) bitmap+=1<<kTOFProtonLess1;
822     else{
823       if (sigmaTOFProtonHyp<2.) bitmap+=1<<kTOFProtonMore1Less2;
824       else { 
825         if (sigmaTOFProtonHyp<3.) bitmap+=1<<kTOFProtonMore2Less3; 
826         else bitmap+=1<<kTOFProtonMore3;
827       }
828     }
829   }
830   
831   
832   
833   return bitmap;
834
835 }
836
837
838 //---------------------------------------------------------------------------
839
840 void AliRDHFCutsDstoKKpi::SetStandardCutsPP2010() {
841   //
842   //STANDARD CUTS USED FOR 2010 pp analysis 
843   //                                            
844   
845   SetName("DstoKKpiCutsStandard");
846   SetTitle("Standard Cuts for D+s analysis");
847   
848   // PILE UP REJECTION
849   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
850
851   // EVENT CUTS
852   SetMinVtxContr(1);
853
854   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
855   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
856   //default
857   esdTrackCuts->SetRequireTPCRefit(kTRUE);
858   esdTrackCuts->SetRequireITSRefit(kTRUE);
859   //esdTrackCuts->SetMinNClustersITS(4); // default is 5
860   esdTrackCuts->SetMinNClustersTPC(70);
861   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
862                                          AliESDtrackCuts::kAny); 
863   // default is kBoth, otherwise kAny
864   esdTrackCuts->SetMinDCAToVertexXY(0.);
865   esdTrackCuts->SetPtRange(0.3,1.e10);
866   
867   AddTrackCuts(esdTrackCuts);
868   
869  
870    
871   const Int_t nptbins=4;
872   Float_t ptbins[nptbins+1];
873   ptbins[0]=2.;
874   ptbins[1]=4.;
875   ptbins[2]=6.;
876   ptbins[3]=8.;
877   ptbins[4]=12.;
878   
879   const Int_t nvars=20;
880       
881   Float_t** anacutsval;
882   anacutsval=new Float_t*[nvars];
883   
884   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}  
885   for(Int_t ipt=0;ipt<nptbins;ipt++){
886   
887     anacutsval[0][ipt]=0.35;
888     anacutsval[1][ipt]=0.3;
889     anacutsval[2][ipt]=0.3;
890     anacutsval[3][ipt]=0.;
891     anacutsval[4][ipt]=0.;
892     anacutsval[5][ipt]=0.005;
893     anacutsval[8][ipt]=0.;
894     anacutsval[10][ipt]=0.;
895     anacutsval[11][ipt]=1000.0;
896     anacutsval[13][ipt]=0.1;
897     anacutsval[16][ipt]=0.;
898     anacutsval[17][ipt]=0.;
899     anacutsval[18][ipt]=0.;
900     anacutsval[19][ipt]=-1.;
901     
902     
903   }
904  
905   //sigmavertex
906
907   anacutsval[6][0]=0.020;
908   anacutsval[6][1]=0.030;
909   anacutsval[6][2]=0.030;
910   anacutsval[6][3]=0.060;
911    
912   //Decay length
913     
914   anacutsval[7][0]=0.035;
915   anacutsval[7][1]=0.035;
916   anacutsval[7][2]=0.040;
917   anacutsval[7][3]=0.040;
918  
919   //cosThetaPoint
920     
921   anacutsval[9][0]=0.94;
922   anacutsval[9][1]=0.94;
923   anacutsval[9][2]=0.94;
924   anacutsval[9][3]=0.94;
925    
926   //phi DeltaMass
927     
928   anacutsval[12][0]=0.0080;
929   anacutsval[12][1]=0.0050;
930   anacutsval[12][2]=0.0045;
931   anacutsval[12][3]=0.0090;
932       
933   //Kin1
934     
935   anacutsval[14][0]=0.10;
936   anacutsval[14][1]=0.05;
937   anacutsval[14][2]=0.0;
938   anacutsval[14][3]=0.05;
939       
940   //Kin2
941   
942   anacutsval[15][0]=0.95;
943   anacutsval[15][1]=0.95;
944   anacutsval[15][2]=1.;
945   anacutsval[15][3]=0.95;     
946   
947   fPidHF->SetOldPid(kTRUE);
948   SetUsePID(kTRUE); 
949   SetPidOption(1);
950   SetMaxPtStrongPid(9999.);
951   SetGlobalIndex(nvars,nptbins);
952   SetPtBins(nptbins+1,ptbins);
953   SetCuts(nvars,nptbins,anacutsval);
954   SetRemoveDaughtersFromPrim(kTRUE);
955   
956   PrintAll();
957
958   for(Int_t iic=0;iic<nvars;iic++){delete [] anacutsval[iic];}
959   delete [] anacutsval;
960   anacutsval=NULL;
961
962   return;
963 }
964