Updates in PbPb cuts (Andrea Rossi)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsD0toKpi.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 D0->Kpi
21 //
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <TDatabasePDG.h>
26 #include <Riostream.h>
27
28 #include "AliRDHFCutsD0toKpi.h"
29 #include "AliAODRecoDecayHF2Prong.h"
30 #include "AliAODTrack.h"
31 #include "AliESDtrack.h"
32 #include "AliAODPid.h"
33 #include "AliTPCPIDResponse.h"
34 #include "AliAODVertex.h"
35 #include "AliKFVertex.h"
36 #include "AliKFParticle.h"
37
38 ClassImp(AliRDHFCutsD0toKpi)
39
40 //--------------------------------------------------------------------------
41 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) : 
42 AliRDHFCuts(name),
43 fUseSpecialCuts(kTRUE),
44 fLowPt(kTRUE),
45 fDefaultPID(kFALSE),
46 fUseKF(kFALSE)
47 {
48   //
49   // Default Constructor
50   //
51   Int_t nvars=11;
52   SetNVars(nvars);
53   TString varNames[11]={"inv. mass [GeV]",   
54                         "dca [cm]",
55                         "cosThetaStar", 
56                         "pTK [GeV/c]",
57                         "pTPi [GeV/c]",
58                         "d0K [cm]",
59                         "d0Pi [cm]",
60                         "d0d0 [cm^2]",
61                         "cosThetaPoint",
62                         "|cosThetaPointXY|", 
63                         "NormDecayLenghtXY"};
64   Bool_t isUpperCut[11]={kTRUE,
65                          kTRUE,
66                          kTRUE,
67                          kFALSE,
68                          kFALSE,
69                          kTRUE,
70                          kTRUE,
71                          kTRUE,
72                          kFALSE,
73                          kFALSE, 
74                          kFALSE};
75   SetVarNames(nvars,varNames,isUpperCut);
76   Bool_t forOpt[11]={kFALSE,
77                      kTRUE,
78                      kTRUE,
79                      kFALSE,
80                      kFALSE,
81                      kFALSE,
82                      kFALSE,
83                      kTRUE,
84                      kTRUE,
85                      kFALSE,
86                      kFALSE};
87   SetVarsForOpt(4,forOpt);
88   Float_t limits[2]={0,999999999.};
89   SetPtBins(2,limits);
90
91 }
92 //--------------------------------------------------------------------------
93 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
94   AliRDHFCuts(source),   
95   fUseSpecialCuts(source.fUseSpecialCuts),
96   fLowPt(source.fLowPt),
97   fDefaultPID(source.fDefaultPID),
98   fUseKF(source.fUseKF)
99 {
100   //
101   // Copy constructor
102   //
103
104 }
105 //--------------------------------------------------------------------------
106 AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)
107 {
108   //
109   // assignment operator
110   //
111   if(&source == this) return *this;
112
113   AliRDHFCuts::operator=(source); 
114   fUseSpecialCuts=source.fUseSpecialCuts;
115   fLowPt=source.fLowPt;
116   fDefaultPID=source.fDefaultPID;
117
118   return *this;
119 }
120
121
122 //---------------------------------------------------------------------------
123 void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
124   // 
125   // Fills in vars the values of the variables 
126   //
127
128   if(nvars!=fnVarsForOpt) {
129     printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
130     return;
131   }
132
133   AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
134  
135   //recalculate vertex w/o daughters
136   Bool_t cleanvtx=kFALSE;
137   AliAODVertex *origownvtx=0x0;
138   if(fRemoveDaughtersFromPrimary) {
139     if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
140     cleanvtx=kTRUE;
141     if(!RecalcOwnPrimaryVtx(dd,aod)) {
142       CleanOwnPrimaryVtx(dd,aod,origownvtx);
143       cleanvtx=kFALSE;
144     }
145   }
146
147   Int_t iter=-1;
148   if(fVarsForOpt[0]){
149     iter++;
150     if(TMath::Abs(pdgdaughters[0])==211) {
151       vars[iter]=dd->InvMassD0();
152     } else {
153       vars[iter]=dd->InvMassD0bar();
154     }
155   }
156   if(fVarsForOpt[1]){
157     iter++;
158     vars[iter]=dd->GetDCA();
159   }
160   if(fVarsForOpt[2]){
161     iter++;
162     if(TMath::Abs(pdgdaughters[0])==211) {
163       vars[iter] = dd->CosThetaStarD0();
164     } else {
165       vars[iter] = dd->CosThetaStarD0bar();
166     }
167   }
168   if(fVarsForOpt[3]){
169     iter++;
170    if(TMath::Abs(pdgdaughters[0])==321) {
171      vars[iter]=dd->PtProng(0);
172    }
173    else{
174      vars[iter]=dd->PtProng(1);
175    }
176   }
177   if(fVarsForOpt[4]){
178     iter++;
179    if(TMath::Abs(pdgdaughters[0])==211) {
180      vars[iter]=dd->PtProng(0);
181    }
182    else{
183      vars[iter]=dd->PtProng(1);
184    }
185   }
186   if(fVarsForOpt[5]){
187     iter++;
188     if(TMath::Abs(pdgdaughters[0])==321) {
189      vars[iter]=dd->Getd0Prong(0);
190    }
191    else{
192      vars[iter]=dd->Getd0Prong(1);
193    }
194   }
195   if(fVarsForOpt[6]){
196     iter++;
197      if(TMath::Abs(pdgdaughters[0])==211) {
198      vars[iter]=dd->Getd0Prong(0);
199    }
200    else{
201      vars[iter]=dd->Getd0Prong(1);
202    }
203   }
204   if(fVarsForOpt[7]){
205     iter++;
206     vars[iter]= dd->Prodd0d0();
207   }
208   if(fVarsForOpt[8]){
209     iter++;
210     vars[iter]=dd->CosPointingAngle();
211   }
212   
213   if(fVarsForOpt[9]){
214                 iter++;
215                 vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
216         }
217   
218    if(fVarsForOpt[10]){
219                 iter++;
220            vars[iter]=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
221         }
222
223    if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
224
225   return;
226 }
227 //---------------------------------------------------------------------------
228 Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
229   //
230   // Apply selection
231   //
232  
233
234   fIsSelectedCuts=0;
235   fIsSelectedPID=0;
236
237   if(!fCutsRD){
238     cout<<"Cut matrice not inizialized. Exit..."<<endl;
239     return 0;
240   }
241   //PrintAll();
242   AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
243
244   if(!d){
245     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
246     return 0;
247   }
248
249   if(fKeepSignalMC) if(IsSignalMC(d,aod,421)) return 3;
250
251   Double_t ptD=d->Pt();
252   if(ptD<fMinPtCand) return 0;
253   if(ptD>fMaxPtCand) return 0;
254  
255   // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
256   Int_t returnvaluePID=3;
257   Int_t returnvalueCuts=3;
258
259   // selection on candidate
260   if(selectionLevel==AliRDHFCuts::kAll || 
261      selectionLevel==AliRDHFCuts::kCandidate) {
262
263     if(!fUseKF) {
264
265       //recalculate vertex w/o daughters
266       AliAODVertex *origownvtx=0x0;
267       if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
268         if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
269         if(!RecalcOwnPrimaryVtx(d,aod)) { 
270           CleanOwnPrimaryVtx(d,aod,origownvtx);
271           return 0;
272         }
273       }
274
275       if(fUseMCVertex) {
276         if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
277         if(!SetMCPrimaryVtx(d,aod)) {
278           CleanOwnPrimaryVtx(d,aod,origownvtx);
279           return 0;
280         }
281       }
282       
283       Double_t pt=d->Pt();
284       
285       Int_t okD0=0,okD0bar=0;
286       
287       Int_t ptbin=PtBin(pt);
288       if (ptbin==-1) {
289         CleanOwnPrimaryVtx(d,aod,origownvtx);
290         return 0;
291       }
292
293       Double_t mD0,mD0bar,ctsD0,ctsD0bar;
294       okD0=1; okD0bar=1;
295       
296       Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
297
298       d->InvMassD0(mD0,mD0bar);
299       if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
300       if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)])  okD0bar = 0;
301       if(!okD0 && !okD0bar)  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
302
303       if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
304     
305       
306       if(d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
307       if(d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
308       if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
309       
310       
311       if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] || 
312          TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
313       if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
314          TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
315       if(!okD0 && !okD0bar)  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
316       
317       if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
318       
319     
320       d->CosThetaStarD0(ctsD0,ctsD0bar);
321       if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0; 
322       if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
323       if(!okD0 && !okD0bar)   {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
324     
325       if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
326
327       
328       if(TMath::Abs(d->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(9,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
329         
330       Double_t normalDecayLengXY=(d->NormalizedDecayLengthXY()*(d->P()/d->Pt()));
331       if (normalDecayLengXY < fCutsRD[GetGlobalIndex(10, ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
332       
333       if (returnvalueCuts!=0) {
334         if (okD0) returnvalueCuts=1; //cuts passed as D0
335         if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
336         if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
337       }
338
339       // call special cuts
340       Int_t special=1;
341       if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
342       if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
343
344       // unset recalculated primary vertex when not needed any more
345       CleanOwnPrimaryVtx(d,aod,origownvtx);
346
347     } else {
348       // go to selection with Kalman vertexing, if requested
349       returnvalueCuts = IsSelectedKF(d,aod);
350     }
351
352     fIsSelectedCuts=returnvalueCuts;
353     if(!returnvalueCuts) return 0;
354   }
355
356
357   // selection on PID 
358   if(selectionLevel==AliRDHFCuts::kAll || 
359      selectionLevel==AliRDHFCuts::kCandidate ||
360      selectionLevel==AliRDHFCuts::kPID) {
361     returnvaluePID = IsSelectedPID(d);
362     fIsSelectedPID=returnvaluePID;
363     if(!returnvaluePID) return 0;
364   }
365
366   Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
367
368   if(!returnvalueComb) return 0;
369
370   // selection on daughter tracks 
371   if(selectionLevel==AliRDHFCuts::kAll || 
372      selectionLevel==AliRDHFCuts::kTracks) {
373     if(!AreDaughtersSelected(d)) return 0;
374   }
375  
376   //  cout<<"Pid = "<<returnvaluePID<<endl;
377   return returnvalueComb;
378 }
379
380 //------------------------------------------------------------------------------------------
381 Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
382                                        AliAODEvent* aod) const {
383   //
384   // Apply selection using KF-vertexing
385   //
386         
387   AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
388   AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1); 
389    
390   if(!track0 || !track1) {
391     cout<<"one or two D0 daughters missing!"<<endl;
392     return 0;
393   }
394
395   // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
396   Int_t returnvalueCuts=3;
397          
398   // candidate selection with AliKF
399   AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
400   
401   Int_t okD0=0,okD0bar=0;
402   okD0=1; okD0bar=1;
403   
404   // convert tracks into AliKFParticles
405   
406   AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
407   AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
408   AliKFParticle posPiKF(*track0,211); // pos pion kandidate
409   AliKFParticle posKKF(*track0,321); // pos kaon kandidate
410   
411   // build D0 candidates
412   
413   AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
414   AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
415   
416   // create kf primary vertices
417   
418   AliAODVertex *vtx1 = aod->GetPrimaryVertex();
419   AliKFVertex primVtx(*vtx1); 
420   AliKFVertex aprimVtx(*vtx1);
421   
422   if(primVtx.GetNContributors()<=0) okD0 = 0;
423   if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
424   if(!okD0 && !okD0bar) returnvalueCuts=0;
425         
426   // calculate mass
427         
428   Double_t d0mass = d0c.GetMass();
429   Double_t ad0mass = ad0c.GetMass();
430         
431   // calculate P of D0 and D0bar
432   Double_t d0P = d0c.GetP();
433   Double_t d0Px = d0c.GetPx();
434   Double_t d0Py = d0c.GetPy();
435   Double_t d0Pz = d0c.GetPz();
436   Double_t ad0P = ad0c.GetP(); 
437   Double_t ad0Px = ad0c.GetPx();
438   Double_t ad0Py = ad0c.GetPy();
439   Double_t ad0Pz = ad0c.GetPz();
440   
441   //calculate Pt of D0 and D0bar
442         
443   Double_t pt=d0c.GetPt(); 
444   Double_t apt=ad0c.GetPt();
445         
446   // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
447   
448   if(track0->GetUsedForPrimVtxFit()) {
449     primVtx -= posPiKF; 
450     aprimVtx -= posKKF;
451   }
452   
453   if(track1->GetUsedForPrimVtxFit()) { 
454     primVtx -= negKKF; 
455     aprimVtx -= negPiKF;
456   }
457   
458   primVtx += d0c;
459   aprimVtx += ad0c;
460   
461   if(primVtx.GetNContributors()<=0) okD0 = 0;
462   if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
463   if(!okD0 && !okD0bar) returnvalueCuts=0;
464   
465   //calculate cut variables
466   
467   // calculate impact params of daughters w.r.t recalculated vertices
468   
469   Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
470   Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
471   Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
472   Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
473         
474   // calculate Product of Impact Params
475         
476   Double_t prodParam = impactPi*impactKa;
477   Double_t aprodParam = aimpactPi*aimpactKa;
478         
479   // calculate cosine of pointing angles
480         
481   TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
482   TVector3 fline(d0c.GetX()-primVtx.GetX(),
483                  d0c.GetY()-primVtx.GetY(),
484                  d0c.GetZ()-primVtx.GetZ());
485   Double_t pta = mom.Angle(fline);
486   Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
487   
488   TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
489   TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
490                   ad0c.GetY()-aprimVtx.GetY(),
491                   ad0c.GetZ()-aprimVtx.GetZ());
492   Double_t apta = amom.Angle(afline);
493   Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
494   
495   // calculate P of Pions at Decay Position of D0 and D0bar candidates
496   negKKF.TransportToParticle(d0c);
497   posPiKF.TransportToParticle(d0c);
498   posKKF.TransportToParticle(ad0c);
499   negPiKF.TransportToParticle(ad0c);
500   
501   Double_t pxPi =  posPiKF.GetPx();
502   Double_t pyPi =  posPiKF.GetPy();
503   Double_t pzPi =  posPiKF.GetPz();
504   Double_t ptPi =  posPiKF.GetPt();
505   
506   Double_t apxPi =  negPiKF.GetPx();
507   Double_t apyPi =  negPiKF.GetPy();
508   Double_t apzPi =  negPiKF.GetPz();
509   Double_t aptPi =  negPiKF.GetPt();
510   
511   // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
512   
513   Double_t ptK =  negKKF.GetPt();
514   Double_t aptK =  posKKF.GetPt();
515         
516   //calculate cos(thetastar)
517   Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
518   Double_t massp[2];
519   massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
520   massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
521   Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
522                                -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
523   
524   // cos(thetastar) for D0 and Pion
525   
526   Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
527   Double_t beta = d0P/d0E;
528   Double_t gamma = d0E/massvtx;
529   TVector3 momPi(pxPi,pyPi,pzPi);
530   TVector3 momTot(d0Px,d0Py,d0Pz);
531   Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
532   Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;      
533         
534   // cos(thetastar) for D0bar and Pion  
535   
536   Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
537   Double_t abeta = ad0P/ad0E;
538   Double_t agamma = ad0E/massvtx;
539   TVector3 amomPi(apxPi,apyPi,apzPi);
540   TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
541   Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
542   Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;  
543   
544   // calculate reduced Chi2 for the full D0 fit
545   d0c.SetProductionVertex(primVtx);
546   ad0c.SetProductionVertex(aprimVtx);
547   negKKF.SetProductionVertex(d0c);
548   posPiKF.SetProductionVertex(d0c);
549   posKKF.SetProductionVertex(ad0c);
550   negPiKF.SetProductionVertex(ad0c);
551   d0c.TransportToProductionVertex();
552   ad0c.TransportToProductionVertex();
553         
554   // calculate the decay length
555   Double_t decayLengthD0 = d0c.GetDecayLength();
556   Double_t adecayLengthD0 = ad0c.GetDecayLength();
557   
558   Double_t chi2D0 = 50.;
559   if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
560     chi2D0 = d0c.GetChi2()/d0c.GetNDF();
561   }
562   
563   Double_t achi2D0 = 50.;
564   if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
565     achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
566   }
567         
568   // Get the Pt-bins
569   Int_t ptbin=PtBin(pt);
570   Int_t aptbin=PtBin(apt);
571
572   if(ptbin < 0) okD0 = 0;
573   if(aptbin < 0) okD0bar = 0;
574   if(!okD0 && !okD0bar) returnvalueCuts=0;
575   
576   if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
577   if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
578   if(!okD0 && !okD0bar) returnvalueCuts=0;
579   
580   
581   if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] || 
582      TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
583   
584   if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
585      TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
586   
587   if(!okD0 && !okD0bar)  returnvalueCuts=0;
588   
589   // for the moment via the standard method due to bug in AliKF
590   if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
591   if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
592   if(!okD0 && !okD0bar) returnvalueCuts=0;
593     
594     
595   if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
596   if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)])  okD0bar = 0;
597   if(!okD0 && !okD0bar)  returnvalueCuts=0;
598   
599   
600   if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0; 
601   if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
602   if(!okD0 && !okD0bar)   returnvalueCuts=0;
603   
604   if(prodParam  > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
605   if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
606   if(!okD0 && !okD0bar)  returnvalueCuts=0;
607     
608   if(cosP  < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0; 
609   if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
610   if(!okD0 && !okD0bar)  returnvalueCuts=0;
611         
612   if(chi2D0  > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0; 
613   if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
614   if(!okD0 && !okD0bar)  returnvalueCuts=0;
615         
616   if(decayLengthD0  < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0; 
617   if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
618   if(!okD0 && !okD0bar)  returnvalueCuts=0;
619     
620   if(returnvalueCuts!=0) {
621     if(okD0) returnvalueCuts=1; //cuts passed as D0
622     if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
623     if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
624   }
625
626   return returnvalueCuts;  
627 }
628
629 //---------------------------------------------------------------------------
630
631 Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
632 {
633   //
634   // Checking if D0 is in fiducial acceptance region 
635   //
636
637   if(pt > 5.) {
638     // applying cut for pt > 5 GeV
639     AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt)); 
640     if (TMath::Abs(y) > 0.8){
641       return kFALSE;
642     }
643   } else {
644     // appliying smooth cut for pt < 5 GeV
645     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
646     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;         
647     AliDebug(2,Form("pt of D0 = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
648     if (y < minFiducialY || y > maxFiducialY){
649       return kFALSE;
650     }
651   }
652
653   return kTRUE;
654 }
655 //---------------------------------------------------------------------------
656 Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d) 
657 {
658   // ############################################################
659   //
660   // Apply PID selection
661   //
662   //
663   // ############################################################
664
665   if(!fUsePID) return 3;
666   if(fDefaultPID) return IsSelectedPIDdefault(d);
667   fWhyRejection=0;
668   Int_t isD0D0barPID[2]={1,2};
669   Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value; 
670   //                                                                                                 same for prong 2
671   //                                               values convention -1 = discarded 
672   //                                                                  0 = not identified (but compatible) || No PID (->hasPID flag)
673   //                                                                  1 = identified
674   // PID search:   pion (TPC) or not K (TOF), Kaon hypothesis for both 
675   // Initial hypothesis: unknwon (but compatible) 
676   combinedPID[0][0]=0;  // prima figlia, Kaon
677   combinedPID[0][1]=0;  // prima figlia, pione
678   combinedPID[1][0]=0;  // seconda figlia, Kaon
679   combinedPID[1][1]=0;  // seconda figlia, pion
680   
681   Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
682   Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
683   for(Int_t daught=0;daught<2;daught++){
684     //Loop con prongs
685     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
686     if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0; 
687     
688     if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
689       checkPIDInfo[daught]=kFALSE; 
690       continue;
691     }
692
693     // identify kaon
694     combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
695
696     // identify pion
697
698     if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
699      combinedPID[daught][1]=0;
700     }else{
701       fPidHF->SetTOF(kFALSE);
702       combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
703       fPidHF->SetTOF(kTRUE);
704       fPidHF->SetCompat(kTRUE);
705      }
706
707
708     if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
709       isD0D0barPID[0]=0;
710       isD0D0barPID[1]=0;
711     }
712     else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
713       if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
714       else isD0D0barPID[0]=0;// if K+ D0 excluded
715     }
716     /*    else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
717           isD0D0barPID[0]=0;
718           isD0D0barPID[1]=0;
719           }
720     */
721     else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){ 
722       if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
723       else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
724     }
725     else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
726       if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
727       else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
728     }
729
730     if(fLowPt && d->Pt()<2.){
731      Double_t sigmaTPC[3]={3.,2.,0.};
732      fPidHF->SetSigmaForTPC(sigmaTPC);
733     // identify kaon
734     combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
735
736     Double_t ptProng=aodtrack->P();
737
738     if(ptProng<0.6){
739      fPidHF->SetCompat(kFALSE);
740      combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
741      fPidHF->SetCompat(kTRUE);
742     }
743
744     if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
745      combinedPID[daught][1]=0;
746     }else{
747       fPidHF->SetTOF(kFALSE);
748       Double_t sigmaTPCpi[3]={3.,3.,0.};
749       fPidHF->SetSigmaForTPC(sigmaTPCpi);
750       combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
751       fPidHF->SetTOF(kTRUE);
752        if(ptProng<0.8){
753         Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
754         if(isTPCpion){
755          combinedPID[daught][1]=1;
756         }else{
757          combinedPID[daught][1]=-1;
758         }
759       }
760     }
761
762    }
763    fPidHF->SetSigmaForTPC(sigma_tmp);
764   }// END OF LOOP ON DAUGHTERS
765
766    if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
767     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
768     return 0;
769    }
770
771
772   // FURTHER PID REQUEST (both daughter info is needed)
773   if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
774     fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
775     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
776     return 0;
777   }
778
779   if(fLowPt && d->Pt()<2.){
780     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
781     if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
782       fWhyRejection=32;// reject cases where the Kaon is not identified
783       return 0;
784     }
785   }
786     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
787
788   //  cout<<"Why? "<<fWhyRejection<<endl;  
789   return isD0D0barPID[0]+isD0D0barPID[1];
790 }
791 //---------------------------------------------------------------------------
792 Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d) 
793 {
794   // ############################################################
795   //
796   // Apply PID selection
797   //
798   //
799   // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
800   //
801   // d must be a AliAODRecoDecayHF2Prong object
802   // returns 0 if both D0 and D0bar are rejectecd
803   //         1 if D0 is accepted while D0bar is rejected
804   //         2 if D0bar is accepted while D0 is rejected
805   //         3 if both are accepted
806   // fWhyRejection variable (not returned for the moment, print it if needed)
807   //               keeps some information on why a candidate has been 
808   //               rejected according to the following (unfriendly?) scheme 
809   //             if more rejection cases are considered interesting, just add numbers
810   //
811   //      TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant) 
812   //              from 20 to 30: "detector" selection (PID acceptance)                                             
813   //                                                  26: TPC refit
814   //                                                  27: ITS refit
815   //                                                  28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
816   //
817   //              from 30 to 40: PID selection
818   //                                                  31: no Kaon compatible tracks found between daughters
819   //                                                  32: no Kaon identified tracks found (strong sel. at low momenta)
820   //                                                  33: both mass hypotheses are rejected 
821   //                  
822   // ############################################################
823
824   if(!fUsePID) return 3;
825   fWhyRejection=0;
826   Int_t isD0D0barPID[2]={1,2};
827   Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
828   Double_t tofSig,times[5];// used fot TOF pid
829   Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
830   Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
831   Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value; 
832   //                                                                                                 same for prong 2
833   //                                               values convention -1 = discarded 
834   //                                                                  0 = not identified (but compatible) || No PID (->hasPID flag)
835   //                                                                  1 = identified
836   // PID search:   pion (TPC) or not K (TOF), Kaon hypothesis for both 
837   // Initial hypothesis: unknwon (but compatible) 
838   isKaonPionTOF[0][0]=0;
839   isKaonPionTOF[0][1]=0;
840   isKaonPionTOF[1][0]=0;
841   isKaonPionTOF[1][1]=0;
842   
843   isKaonPionTPC[0][0]=0;
844   isKaonPionTPC[0][1]=0;
845   isKaonPionTPC[1][0]=0;
846   isKaonPionTPC[1][1]=0;
847   
848   combinedPID[0][0]=0;
849   combinedPID[0][1]=0;
850   combinedPID[1][0]=0;
851   combinedPID[1][1]=0;
852   
853   
854  
855   for(Int_t daught=0;daught<2;daught++){
856     //Loop con prongs
857     
858     // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
859
860     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught); 
861    
862     if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
863       fWhyRejection=26;
864       return 0;
865     } 
866     if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
867       fWhyRejection=27;
868       return 0;
869     } 
870     
871     AliAODPid *pid=aodtrack->GetDetPid();
872     if(!pid) {
873       //delete esdtrack;
874       hasPID[daught]--;
875       continue;
876     }
877   
878     // ########### Step 1- Check of TPC and TOF response ####################
879
880     Double_t ptrack=aodtrack->P();
881     //#################### TPC PID #######################
882      if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
883        // NO TPC PID INFO FOR THIS TRACK 
884        hasPID[daught]--;
885      }
886      else {
887        static AliTPCPIDResponse theTPCpid;
888        AliAODPid *pidObj = aodtrack->GetDetPid();
889        Double_t ptProng=pidObj->GetTPCmomentum();
890        nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
891        nsigmaTPCK =  theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
892        //if(ptrack<0.6){
893        if(ptProng<0.6){
894          if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
895          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
896          if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
897          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
898        }
899        //else if(ptrack<.8){
900        else if(ptProng<.8){
901          if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
902          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
903          if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
904          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
905        }     
906        else {
907          //     if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
908          if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
909          //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
910          if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
911        }
912      }
913     
914     
915     // ##### TOF PID: do not ask nothing for pion/protons ############
916      if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
917        // NO TOF PID INFO FOR THIS TRACK      
918        hasPID[daught]--;
919      }
920      else{
921        tofSig=pid->GetTOFsignal(); 
922        pid->GetIntegratedTimes(times);
923        if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
924        if(TMath::Abs(tofSig-times[3])>3.*160.){
925          isKaonPionTOF[daught][0]=-1;
926        }
927        else {    
928          if(ptrack<1.5){
929            isKaonPionTOF[daught][0]=1;
930          }
931        }
932      }
933      
934      //######### Step 2: COMBINE TOF and TPC PID ###############
935      // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
936      combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
937      combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
938      
939      
940      //######### Step 3:   USE PID INFO     
941      
942      if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
943        isD0D0barPID[0]=0;
944        isD0D0barPID[1]=0;
945      }
946      else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-), if k for both TPC and TOF -> is K
947        if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
948        else isD0D0barPID[0]=0;// if K+ D0 excluded
949      }
950      else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-) and k- only for TPC or TOF -> reject
951        isD0D0barPID[0]=0;
952        isD0D0barPID[1]=0;
953      }
954      else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
955        if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
956        else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
957      }
958      else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
959        if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
960       else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
961      }
962      
963      // ##########  ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification      ###############################
964      // ########## more tolerant criteria for single particle ID-> more selective criteria for D0   ##############################
965      // ###############                     NOT OPTIMIZED YET                                  ###################################
966      if(d->Pt()<2.){
967        isKaonPionTPC[daught][0]=0;
968        isKaonPionTPC[daught][1]=0;
969        AliAODPid *pidObj = aodtrack->GetDetPid();
970        Double_t ptProng=pidObj->GetTPCmomentum();
971        //if(ptrack<0.6){
972        if(ptProng<0.6){
973          if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
974          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
975          if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
976          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
977      }
978        //else if(ptrack<.8){
979        else if(ptProng<.8){
980          if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
981          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
982          if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
983          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
984        }     
985        else {
986          if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
987          if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
988        }
989      }
990      
991   }// END OF LOOP ON DAUGHTERS
992   
993   // FURTHER PID REQUEST (both daughter info is needed)
994   if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
995     fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
996     return 0;
997   }
998   else if(hasPID[0]==0&&hasPID[1]==0){
999     fWhyRejection=28;// reject cases in which no PID info is available  
1000     return 0;
1001   }
1002   if(d->Pt()<2.){
1003     // request of K identification at low D0 pt
1004     combinedPID[0][0]=0;
1005     combinedPID[0][1]=0;
1006     combinedPID[1][0]=0;
1007     combinedPID[1][1]=0;
1008     
1009     combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
1010     combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
1011     combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
1012     combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
1013     
1014     if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
1015       fWhyRejection=32;// reject cases where the Kaon is not identified
1016       return 0;
1017     }
1018   }
1019
1020   //  cout<<"Why? "<<fWhyRejection<<endl;  
1021   return isD0D0barPID[0]+isD0D0barPID[1];
1022 }
1023
1024
1025
1026 //---------------------------------------------------------------------------
1027 Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
1028                                                  Int_t selectionvalCand,
1029                                                  Int_t selectionvalPID) const
1030 {
1031   //
1032   // This method combines the tracks, PID and cuts selection results
1033   //
1034   if(selectionvalTrack==0) return 0;
1035
1036   Int_t returnvalue;
1037
1038   switch(selectionvalPID) {
1039   case 0:
1040     returnvalue=0;
1041     break;
1042   case 1:
1043     returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1044     break;
1045   case 2:
1046     returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1047     break;
1048   case 3:
1049     returnvalue=selectionvalCand;
1050     break;
1051   default:
1052     returnvalue=0;
1053     break;
1054   }
1055
1056   return returnvalue;
1057 }
1058 //----------------------------------------------------------------------------
1059 Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
1060 {
1061   //
1062   // Note: this method is temporary
1063   // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1064   //
1065
1066   //apply cuts
1067
1068   Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1069   // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1070   // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution), 
1071
1072   Int_t returnvalue=3; //cut passed
1073   for(Int_t i=0;i<2/*prongs*/;i++){
1074     if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1075   }
1076   if(d->DecayLength2()<decLengthCut*decLengthCut)  return 0; //decLengthCut not passed
1077   if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut)  return 0; //decLengthCut not passed
1078         
1079   return returnvalue;
1080 }
1081
1082 //----------------------------------------------
1083 void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
1084 {
1085   //switch on candidate selection via AliKFparticle
1086   if(!useKF) return;
1087   if(useKF){
1088     fUseKF=useKF;
1089     Int_t nvarsKF=11;
1090     SetNVars(nvarsKF);
1091     TString varNamesKF[11]={"inv. mass [GeV]",   
1092                             "dca [cm]",
1093                             "cosThetaStar", 
1094                             "pTK [GeV/c]",
1095                             "pTPi [GeV/c]",
1096                             "d0K [cm]",
1097                             "d0Pi [cm]",
1098                             "d0d0 [cm^2]",
1099                             "cosThetaPoint"
1100                             "DecayLength[cm]",
1101                             "RedChi2"};
1102     Bool_t isUpperCutKF[11]={kTRUE,
1103                              kTRUE,
1104                              kTRUE,
1105                              kFALSE,
1106                              kFALSE,
1107                              kTRUE,
1108                              kTRUE,
1109                              kTRUE,
1110                              kFALSE,
1111                              kFALSE,
1112                              kTRUE};
1113     SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1114     SetGlobalIndex();
1115     Bool_t forOpt[11]={kFALSE,
1116                     kTRUE,
1117                     kTRUE,
1118                     kFALSE,
1119                     kFALSE,
1120                     kFALSE,
1121                     kFALSE,
1122                     kTRUE,
1123                     kTRUE,
1124                     kFALSE,
1125                     kFALSE};
1126     SetVarsForOpt(4,forOpt);
1127   }
1128   return;
1129 }
1130
1131
1132 void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
1133   //
1134   //STANDARD CUTS USED FOR 2010 pp analysis 
1135   //dca cut will be enlarged soon to 400 micron
1136   //
1137   
1138   SetName("D0toKpiCutsStandard");
1139   SetTitle("Standard Cuts for D0 analysis");
1140   
1141   // PILE UP REJECTION
1142   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1143
1144   // EVENT CUTS
1145   SetMinVtxContr(1);
1146
1147   
1148   // TRACKS ON SINGLE TRACKS
1149   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1150   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1151   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1152   esdTrackCuts->SetRequireITSRefit(kTRUE);
1153   //  esdTrackCuts->SetMinNClustersITS(4);
1154   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1155   esdTrackCuts->SetMinDCAToVertexXY(0.);
1156   esdTrackCuts->SetEtaRange(-0.8,0.8);
1157   esdTrackCuts->SetPtRange(0.3,1.e10);
1158   
1159   AddTrackCuts(esdTrackCuts);
1160   
1161   const Int_t nptbins =14;
1162   const Double_t ptmax = 9999.;
1163   const Int_t nvars=11;
1164   Float_t ptbins[nptbins+1];
1165   ptbins[0]=0.;
1166   ptbins[1]=0.5;        
1167   ptbins[2]=1.;
1168   ptbins[3]=2.;
1169   ptbins[4]=3.;
1170   ptbins[5]=4.;
1171   ptbins[6]=5.;
1172   ptbins[7]=6.;
1173   ptbins[8]=7.;
1174   ptbins[9]=8.;
1175   ptbins[10]=12.;
1176   ptbins[11]=16.;
1177   ptbins[12]=20.;
1178   ptbins[13]=24.;
1179   ptbins[14]=ptmax;
1180
1181   SetGlobalIndex(nvars,nptbins);
1182   SetPtBins(nptbins+1,ptbins);
1183   
1184   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* pt<0.5*/
1185                                                   {0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* 0.5<pt<1*/
1186                                                   {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.80,0.,0.},/* 1<pt<2 */
1187                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,0.},/* 2<pt<3 */
1188                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 3<pt<4 */
1189                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 4<pt<5 */
1190                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 5<pt<6 */
1191                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */
1192                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-7000.*1E-8,0.85,0.,0.},/* 7<pt<8 */
1193                                                   {0.400,300.*1E-4,0.9,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.,0.},/* 8<pt<12 */
1194                                                   {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,10000.*1E-8,0.85,0.,0.},/* 12<pt<16 */
1195                                                   {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 16<pt<20 */
1196                                                   {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 20<pt<24 */
1197                                                   {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.}};/* pt>24 */
1198   
1199   
1200   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1201   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1202   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1203   
1204   for (Int_t ibin=0;ibin<nptbins;ibin++){
1205     for (Int_t ivar = 0; ivar<nvars; ivar++){
1206       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];      
1207     }
1208   }
1209   
1210   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1211   SetUseSpecialCuts(kTRUE);
1212   SetRemoveDaughtersFromPrim(kTRUE);
1213   
1214   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1215   delete [] cutsMatrixTransposeStand;
1216   cutsMatrixTransposeStand=NULL;
1217
1218   // PID SETTINGS
1219   AliAODPidHF* pidObj=new AliAODPidHF();
1220   //pidObj->SetName("pid4D0");
1221   Int_t mode=1;
1222   const Int_t nlims=2;
1223   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1224   Bool_t compat=kTRUE; //effective only for this mode
1225   Bool_t asym=kTRUE;
1226   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1227   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1228   pidObj->SetMatch(mode);
1229   pidObj->SetPLimit(plims,nlims);
1230   pidObj->SetSigma(sigmas);
1231   pidObj->SetCompat(compat);
1232   pidObj->SetTPC(kTRUE);
1233   pidObj->SetTOF(kTRUE);
1234   pidObj->SetPCompatTOF(1.5);
1235   pidObj->SetSigmaForTPCCompat(3.);
1236   pidObj->SetSigmaForTOFCompat(3.);
1237
1238   SetPidHF(pidObj);
1239   SetUsePID(kTRUE);
1240   SetUseDefaultPID(kFALSE);
1241
1242
1243   PrintAll();
1244
1245   delete pidObj;
1246   pidObj=NULL;
1247
1248   return;
1249
1250 }
1251
1252
1253 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1254   //
1255   //PRELIMINARY CUTS USED FOR 2010 PbPb analysis
1256   //... EVOLVING SOON 
1257   // 
1258   
1259   SetName("D0toKpiCutsStandard");
1260   SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1261   
1262   // PILE UP REJECTION
1263   //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1264   // CENTRALITY SELECTION
1265   SetMinCentrality(0.);
1266   SetMaxCentrality(80.);
1267   SetUseCentrality(AliRDHFCuts::kCentV0M);
1268
1269
1270   // EVENT CUTS
1271   SetMinVtxContr(1);
1272   // MAX Z-VERTEX CUT
1273   SetMaxVtxZ(10.);
1274   
1275   // TRACKS ON SINGLE TRACKS
1276   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1277   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1278   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1279   esdTrackCuts->SetRequireITSRefit(kTRUE);
1280   //  esdTrackCuts->SetMinNClustersITS(4);
1281   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1282   esdTrackCuts->SetMinDCAToVertexXY(0.);
1283   esdTrackCuts->SetEtaRange(-0.8,0.8);
1284   esdTrackCuts->SetPtRange(0.8,1.e10);
1285
1286   esdTrackCuts->SetMaxDCAToVertexXY(1.);  
1287   esdTrackCuts->SetMaxDCAToVertexZ(1.);
1288   esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");  
1289
1290
1291   AddTrackCuts(esdTrackCuts);
1292
1293   // CANDIDATE CUTS  
1294   const Int_t nptbins =13;
1295   const Double_t ptmax = 9999.;
1296   const Int_t nvars=11;
1297   Float_t ptbins[nptbins+1];
1298   ptbins[0]=0.;
1299   ptbins[1]=0.5;        
1300   ptbins[2]=1.;
1301   ptbins[3]=2.;
1302   ptbins[4]=3.;
1303   ptbins[5]=4.;
1304   ptbins[6]=5.;
1305   ptbins[7]=6.;
1306   ptbins[8]=8.;
1307   ptbins[9]=12.;
1308   ptbins[10]=16.;
1309   ptbins[11]=20.;
1310   ptbins[12]=24.;
1311   ptbins[13]=ptmax;
1312
1313   SetGlobalIndex(nvars,nptbins);
1314   SetPtBins(nptbins+1,ptbins);
1315   SetMinPtCandidate(2.);
1316
1317   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.85,0.,5.},/* pt<0.5*/
1318                                                   {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.9,0.,5.},/* 0.5<pt<1*/
1319                                                   {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.,5.},/* 1<pt<2 */
1320                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.95,0.998,5.},/* 2<pt<3 */
1321                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.998,5.},/* 3<pt<4 */
1322                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.998,5.},/* 4<pt<5 */
1323                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
1324                                                   {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
1325                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.998,5.},/* 8<pt<12 */
1326                                                   {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,5.},/* 12<pt<16 */
1327                                                   {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.998,5.},/* 16<pt<20 */
1328                                                   {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.998,5.},/* 20<pt<24 */
1329                                                   {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.998,5.}};/* pt>24 */
1330   
1331   
1332   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1333   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1334   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1335   
1336   for (Int_t ibin=0;ibin<nptbins;ibin++){
1337     for (Int_t ivar = 0; ivar<nvars; ivar++){
1338       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];      
1339     }
1340   }
1341   
1342   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1343   SetUseSpecialCuts(kTRUE);
1344   SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1345   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1346   delete [] cutsMatrixTransposeStand;
1347   cutsMatrixTransposeStand=NULL;
1348   
1349   // PID SETTINGS
1350   AliAODPidHF* pidObj=new AliAODPidHF();
1351   //pidObj->SetName("pid4D0");
1352   Int_t mode=1;
1353   const Int_t nlims=2;
1354   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1355   Bool_t compat=kTRUE; //effective only for this mode
1356   Bool_t asym=kTRUE;
1357   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1358   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1359   pidObj->SetMatch(mode);
1360   pidObj->SetPLimit(plims,nlims);
1361   pidObj->SetSigma(sigmas);
1362   pidObj->SetCompat(compat);
1363   pidObj->SetTPC(kTRUE);
1364   pidObj->SetTOF(kTRUE);
1365   pidObj->SetPCompatTOF(2.);
1366   pidObj->SetSigmaForTPCCompat(3.);
1367   pidObj->SetSigmaForTOFCompat(3.);  
1368
1369
1370   SetPidHF(pidObj);
1371   SetUsePID(kTRUE);
1372   SetUseDefaultPID(kFALSE);
1373
1374
1375   PrintAll();
1376
1377
1378   delete pidObj;
1379   pidObj=NULL;
1380
1381   return;
1382
1383 }