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