Fix in the last caall to CleanOwnPrimaryVertex
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsDplustoKpipi.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 D+->Kpipi
21 //
22 // Author: R. Bala, bala@to.infn.it
23 //         G. Ortona, ortona@to.infn.it
24 /////////////////////////////////////////////////////////////
25
26 #include <TDatabasePDG.h>
27 #include <Riostream.h>
28
29 #include "AliRDHFCutsDplustoKpipi.h"
30 #include "AliAODRecoDecayHF3Prong.h"
31 #include "AliAODTrack.h"
32 #include "AliESDtrack.h"
33
34
35 ClassImp(AliRDHFCutsDplustoKpipi)
36
37 //--------------------------------------------------------------------------
38 AliRDHFCutsDplustoKpipi::AliRDHFCutsDplustoKpipi(const char* name) : 
39 AliRDHFCuts(name),
40   fUseStrongPid(0),
41   fMaxPtStrongPid(0.),
42   fUseImpParProdCorrCut(kFALSE)
43 {
44   //
45   // Default Constructor
46   //
47   Int_t nvars=14;
48   SetNVars(nvars);
49   TString varNames[14]={"inv. mass [GeV]",
50                         "pTK [GeV/c]",
51                         "pTPi [GeV/c]",
52                         "d0K [cm]   lower limit!",
53                         "d0Pi [cm]  lower limit!",
54                         "dist12 (cm)",
55                         "sigmavert (cm)",
56                         "dist prim-sec (cm)",
57                         "pM=Max{pT1,pT2,pT3} (GeV/c)",
58                         "cosThetaPoint",
59                         "Sum d0^2 (cm^2)",
60                         "dca cut (cm)",
61                         "dec len XY (cm)",
62                         "cosThetaPointXY"};
63   Bool_t isUpperCut[14]={kTRUE,
64                          kFALSE,
65                          kFALSE,
66                          kFALSE,
67                          kFALSE,
68                          kFALSE,
69                          kTRUE,
70                          kFALSE,
71                          kFALSE,
72                          kFALSE,
73                          kFALSE,
74                          kTRUE,
75                          kFALSE,
76                          kFALSE};
77   SetVarNames(nvars,varNames,isUpperCut);
78   Bool_t forOpt[14]={kFALSE,
79                      kFALSE,
80                      kFALSE,
81                      kFALSE,
82                      kFALSE,
83                      kFALSE,
84                      kTRUE,
85                      kTRUE,
86                      kTRUE,
87                      kTRUE,
88                      kTRUE,
89                      kFALSE,
90                      kTRUE,
91                      kTRUE};
92   SetVarsForOpt(7,forOpt);
93   Float_t limits[2]={0,999999999.};
94   SetPtBins(2,limits);
95   if(fPidHF)delete fPidHF;
96   fPidHF=new AliAODPidHF();
97   Double_t plim[2]={0.6,0.8};
98   Double_t nsigma[5]={2.,1.,2.,3.,0.};
99   
100   fPidHF->SetPLimit(plim);
101   fPidHF->SetAsym(kTRUE);
102   fPidHF->SetSigma(nsigma);
103   fPidHF->SetMatch(1);
104   fPidHF->SetTPC(1);
105   fPidHF->SetTOF(1);
106   fPidHF->SetITS(0);
107   fPidHF->SetTRD(0);
108   fPidHF->SetCompat(kTRUE);
109
110   
111 }
112
113
114
115
116
117
118
119
120 //--------------------------------------------------------------------------
121 AliRDHFCutsDplustoKpipi::AliRDHFCutsDplustoKpipi(const AliRDHFCutsDplustoKpipi &source) :
122   AliRDHFCuts(source),
123   fUseStrongPid(source.fUseStrongPid),
124   fMaxPtStrongPid(source.fMaxPtStrongPid),
125   fUseImpParProdCorrCut(source.fUseImpParProdCorrCut)
126 {
127   //
128   // Copy constructor
129   //
130
131 }
132 //--------------------------------------------------------------------------
133 AliRDHFCutsDplustoKpipi &AliRDHFCutsDplustoKpipi::operator=(const AliRDHFCutsDplustoKpipi &source)
134 {
135   //
136   // assignment operator
137   //
138   if(&source == this) return *this;
139
140   AliRDHFCuts::operator=(source);
141
142   return *this;
143 }
144 //
145
146
147 //---------------------------------------------------------------------------
148 void AliRDHFCutsDplustoKpipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
149   // 
150   // Fills in vars the values of the variables 
151   //
152
153
154   if(nvars!=fnVarsForOpt) {
155     printf("AliRDHFCutsDplustoKpipi::GetCutsVarsForOpt: wrong number of variables\n");
156     return;
157   }
158
159   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
160  
161   //recalculate vertex w/o daughters
162   Bool_t cleanvtx=kFALSE;
163   AliAODVertex *origownvtx=0x0;
164   if(fRemoveDaughtersFromPrimary) {
165     if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
166     cleanvtx=kTRUE;
167     if(!RecalcOwnPrimaryVtx(dd,aod)) {
168       CleanOwnPrimaryVtx(dd,aod,origownvtx);
169       cleanvtx=kFALSE;
170     }
171   }
172
173   Int_t iter=-1;
174   if(fVarsForOpt[0]){
175     iter++;
176     vars[iter]=dd->InvMassDplus();
177   }
178   if(fVarsForOpt[1]){
179     iter++;
180     for(Int_t iprong=0;iprong<3;iprong++){
181       if(TMath::Abs(pdgdaughters[iprong])==321) {
182         vars[iter]=dd->PtProng(iprong);
183       }
184     }
185   }
186   if(fVarsForOpt[2]){
187     iter++;
188     Float_t minPtDau=1000000.0;
189     for(Int_t iprong=0;iprong<3;iprong++){
190       if(TMath::Abs(pdgdaughters[iprong])==211) {
191         if(dd->PtProng(iprong)<minPtDau){
192           minPtDau=dd->PtProng(iprong);
193         }
194       }
195     }
196     vars[iter]=minPtDau;
197   }
198   if(fVarsForOpt[3]){
199     iter++;
200     for(Int_t iprong=0;iprong<3;iprong++){
201       if(TMath::Abs(pdgdaughters[iprong])==321) {
202         vars[iter]=dd->Getd0Prong(iprong);
203       }
204     }
205   }
206   if(fVarsForOpt[4]){
207     iter++;
208     Float_t minImpParDau=1000000.0;
209     for(Int_t iprong=0;iprong<3;iprong++){
210       if(TMath::Abs(pdgdaughters[iprong])==211) {
211         if(dd->Getd0Prong(iprong)<minImpParDau){
212           minImpParDau=dd->Getd0Prong(iprong);
213         }
214       }
215     }
216    vars[iter]=minImpParDau;
217   }
218   if(fVarsForOpt[5]){
219     iter++;
220     Float_t dist12 = dd->GetDist12toPrim();
221     Float_t dist23 = dd->GetDist23toPrim();
222     if(dist12<dist23)vars[iter]=dist12;
223     else vars[iter]=dist23;
224   }
225   if(fVarsForOpt[6]){
226     iter++;
227     vars[iter]=dd->GetSigmaVert();
228   }
229   if(fVarsForOpt[7]){
230     iter++;
231     vars[iter] = dd->DecayLength();
232   }
233   if(fVarsForOpt[8]){
234     iter++;
235     Float_t ptmax=0;
236     for(Int_t i=0;i<3;i++){
237       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
238     }
239     vars[iter]=ptmax;
240   }
241   if(fVarsForOpt[9]){
242     iter++;
243     vars[iter]=dd->CosPointingAngle();
244   }
245   if(fVarsForOpt[10]){
246     iter++;
247     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
248   }
249   if(fVarsForOpt[11]){
250     iter++;
251     Float_t maxDCA=0;
252     for(Int_t iprong=0;iprong<3;iprong++){
253       if(dd->GetDCA(iprong)<maxDCA){
254         maxDCA=dd->GetDCA(iprong);
255       }
256     }
257     vars[iter]=maxDCA;
258   }
259   if(fVarsForOpt[12]){
260     iter++;
261     vars[iter]=dd->NormalizedDecayLengthXY()*dd->P()/dd->Pt();
262   }
263   if(fVarsForOpt[13]){
264     iter++;
265     vars[iter]=dd->CosPointingAngleXY();
266   }
267
268   if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
269
270   return;
271 }
272 //---------------------------------------------------------------------------
273 Bool_t AliRDHFCutsDplustoKpipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
274 {
275   //
276   // Checking if Dplus is in fiducial acceptance region 
277   //
278
279   if(pt > 5.) {
280     // applying cut for pt > 5 GeV
281     AliDebug(2,Form("pt of D+ = %f (> 5), cutting at |y| < 0.8",pt)); 
282     if (TMath::Abs(y) > 0.8) return kFALSE;
283     
284   } else {
285     // appliying smooth cut for pt < 5 GeV
286     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
287     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;         
288     AliDebug(2,Form("pt of D+ = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
289     if (y < minFiducialY || y > maxFiducialY) return kFALSE;    
290   }
291
292   return kTRUE;
293 }
294
295 //---------------------------------------------------------------------------
296 Int_t AliRDHFCutsDplustoKpipi::IsSelectedPID(AliAODRecoDecayHF *rd)
297 {
298   //
299   // PID selection, returns 3 if accepted, 0 if not accepted
300   // 
301   if(!fUsePID || !rd) return 3;
302   //if(fUsePID)printf("i am inside the pid \n");
303   Int_t nkaons=0;
304   Int_t nNotKaons=0;
305   Int_t sign= rd->GetCharge(); 
306   for(Int_t daught=0;daught<3;daught++){
307     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(daught);
308     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
309     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
310     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
311     
312     if(isProton>0 &&  isKaon<0  && isPion<0) return 0;
313     if(isKaon>0 && isPion<0) nkaons++;
314     if(isKaon<0) nNotKaons++;  
315     if(sign==track->Charge()){//pions
316       if(isPion<0)return 0;
317       if(rd->Pt()<fMaxPtStrongPid && isPion<=0 && fUseStrongPid>1)return 0;
318     }
319     else{//kaons
320       if(isKaon<0)return 0;
321         if(rd->Pt()<fMaxPtStrongPid && isKaon<=0 && fUseStrongPid>0)return 0;
322     }
323   }
324   
325   if(nkaons>1)return 0;
326   if(nNotKaons==3)return 0;
327   
328   return 3;   
329 }
330
331
332 //---------------------------------------------------------------------------
333 Int_t AliRDHFCutsDplustoKpipi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
334   //
335   // Apply selection, returns 3 if accepted, 0 if not accepted
336   //
337
338
339   fIsSelectedCuts=0;
340   fIsSelectedPID=0;
341
342   if(!fCutsRD){
343     cout<<"Cut matrix not inizialized. Exit..."<<endl;
344     return 0;
345   }
346   //PrintAll();
347   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj; 
348
349   
350   if(!d){
351     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
352     return 0;
353   }
354
355   if(fKeepSignalMC) if(IsSignalMC(d,aod,411)) return 3;
356
357   // PID selection
358   Int_t returnvaluePID=3;
359   Int_t returnvalueCuts=3;
360
361   Double_t pt=d->Pt();
362   if(pt<fMinPtCand) return 0;
363   if(pt>fMaxPtCand) return 0;
364
365   
366   // selection on candidate
367   if(selectionLevel==AliRDHFCuts::kAll || 
368      selectionLevel==AliRDHFCuts::kCandidate) {
369     
370     //recalculate vertex w/o daughters
371     AliAODVertex *origownvtx=0x0;
372     if(fRemoveDaughtersFromPrimary) {
373       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
374       if(!RecalcOwnPrimaryVtx(d,aod)) {
375         CleanOwnPrimaryVtx(d,aod,origownvtx);
376         return 0;
377       }
378     }
379
380     Int_t ptbin=PtBin(pt);
381     if (ptbin==-1) {
382       CleanOwnPrimaryVtx(d,aod,origownvtx);
383       return 0;
384     }
385     
386     Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
387     Double_t mDplus=d->InvMassDplus();
388     if(TMath::Abs(mDplus-mDplusPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
389
390     //2track cuts
391     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]|| d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
392
393     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
394     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
395
396     if(fUseImpParProdCorrCut){
397       if(d->Getd0Prong(0)*d->Getd0Prong(1)<0. && d->Getd0Prong(2)*d->Getd0Prong(1)<0.) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
398     }
399
400     //sec vert
401     if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
402
403     //DCA
404     for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
405
406     if(d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] || TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}//Kaon
407
408     if(d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}//Pion1
409
410     if(d->Pt2Prong(2) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}//Pion2
411     
412     if(d->Pt2Prong(0)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] && d->Pt2Prong(1)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] && d->Pt2Prong(2)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
413
414     if(d->DecayLength2()<fCutsRD[GetGlobalIndex(7,ptbin)]*fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
415
416     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
417
418     if(d->NormalizedDecayLengthXY()*d->P()/pt<fCutsRD[GetGlobalIndex(12,ptbin)]){CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
419
420     if(d->CosPointingAngleXY()<fCutsRD[GetGlobalIndex(13,ptbin)]){CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
421
422     // unset recalculated primary vertex when not needed any more
423     CleanOwnPrimaryVtx(d,aod,origownvtx);
424     
425     fIsSelectedCuts=returnvalueCuts;
426
427     //if(!returnvalueCuts) return 0; // returnvalueCuts cannot be 0 here
428   }
429   
430   if(selectionLevel==AliRDHFCuts::kAll || 
431      selectionLevel==AliRDHFCuts::kCandidate ||     
432      selectionLevel==AliRDHFCuts::kPID) {
433     returnvaluePID = IsSelectedPID(d);
434     fIsSelectedPID=returnvaluePID;
435   }
436   if(returnvaluePID==0)return 0;
437
438   // selection on daughter tracks 
439   if(selectionLevel==AliRDHFCuts::kAll || 
440      selectionLevel==AliRDHFCuts::kTracks) {
441     if(!AreDaughtersSelected(d)) return 0;
442   }
443   
444  
445
446
447   return 3;
448 }
449
450
451
452
453 //---------------------------------------------------------------------------
454
455
456 void AliRDHFCutsDplustoKpipi::SetStandardCutsPP2010() {
457   //
458   //STANDARD CUTS USED FOR 2010 pp analysis 
459   //                                            
460   
461   SetName("DplustoKpipiCutsStandard");
462   SetTitle("Standard Cuts for D+ analysis");
463   
464   // PILE UP REJECTION
465   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
466
467   // EVENT CUTS
468   SetMinVtxContr(1);
469
470   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
471   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
472   //default
473   esdTrackCuts->SetRequireTPCRefit(kTRUE);
474   esdTrackCuts->SetRequireITSRefit(kTRUE);
475   //esdTrackCuts->SetMinNClustersITS(4); // default is 5
476   esdTrackCuts->SetMinNClustersTPC(70);
477   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
478                                          AliESDtrackCuts::kAny); 
479   // default is kBoth, otherwise kAny
480   esdTrackCuts->SetMinDCAToVertexXY(0.);
481   esdTrackCuts->SetPtRange(0.3,1.e10);
482   
483   AddTrackCuts(esdTrackCuts);
484   
485  
486   const Int_t nptbins =13;
487   const Int_t nvars=14;
488   Float_t ptbins[nptbins+1];
489   ptbins[0]=0.;
490   ptbins[1]=1;  
491   ptbins[2]=2.;
492   ptbins[3]=3.;
493   ptbins[4]=4.;
494   ptbins[5]=5.;
495   ptbins[6]=6.;
496   ptbins[7]=8.;
497   ptbins[8]=10.;
498   ptbins[9]=12.;
499   ptbins[10]=14.;
500   ptbins[11]=16.;
501   ptbins[12]=24.;
502   ptbins[13]=99999.;
503       
504     
505   Float_t** anacutsval;
506   anacutsval=new Float_t*[nvars];
507   
508   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
509
510   //Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
511   for(Int_t ipt=0;ipt<nptbins;ipt++){
512     anacutsval[0][ipt]=0.2;
513     anacutsval[3][ipt]=0.;
514     anacutsval[4][ipt]=0.;
515     anacutsval[5][ipt]=0.01;
516     anacutsval[11][ipt]=10000000000.;
517     anacutsval[12][ipt]=0.;
518     anacutsval[13][ipt]=0.;
519   }
520
521   anacutsval[1][0]=0.3;
522   anacutsval[1][1]=0.3;
523   anacutsval[1][2]=0.4; 
524   anacutsval[2][0]=0.3;
525   anacutsval[2][1]=0.3;
526   anacutsval[2][2]=0.4;  
527   for(Int_t ipt=3;ipt<nptbins;ipt++){
528     anacutsval[1][ipt]=0.4;
529     anacutsval[2][ipt]=0.4;
530   }
531   
532   anacutsval[6][0]=0.022100;
533   anacutsval[6][1]=0.022100;
534   anacutsval[6][2]=0.034;
535   anacutsval[6][3]=0.020667;
536   anacutsval[6][4]=0.020667;
537   anacutsval[6][5]=0.023333;
538     
539   
540   anacutsval[7][0]=0.08;
541   anacutsval[7][1]=0.08;
542   anacutsval[7][2]=0.09;  
543   anacutsval[7][3]=0.095;
544   anacutsval[7][4]=0.095;
545    
546   anacutsval[8][0]=0.5;
547   anacutsval[8][1]=0.5;
548   anacutsval[8][2]=1.0;
549   anacutsval[8][3]=0.5;
550   anacutsval[8][4]=0.5;
551      
552     
553   anacutsval[9][0]=0.95;
554   anacutsval[9][1]=0.95;
555   anacutsval[9][2]=0.95; 
556   anacutsval[9][3]=0.95; 
557   anacutsval[9][4]= 0.95;
558   anacutsval[9][5]=0.92;
559   anacutsval[9][6]=0.92;
560   anacutsval[9][7]=0.92;
561   anacutsval[9][8]=0.92;
562   anacutsval[9][9]=0.90;
563   anacutsval[9][10]=0.90;
564   anacutsval[9][11]=0.90;
565   anacutsval[9][12]=0.90;
566   
567   anacutsval[10][0]=0.0055;
568   anacutsval[10][1]=0.0055;
569   anacutsval[10][2]= 0.0028;
570   anacutsval[10][3]=0.000883;
571   anacutsval[10][4]=0.000883;
572
573   
574   for(Int_t ipt=5;ipt<nptbins;ipt++){
575     anacutsval[6][ipt]=0.02333;
576     anacutsval[7][ipt]=0.115;
577     anacutsval[8][ipt]=0.5;
578     anacutsval[10][ipt]=0.000883;
579     }   
580
581
582   SetGlobalIndex(nvars,nptbins);
583   SetPtBins(nptbins+1,ptbins);
584   SetCuts(nvars,nptbins,anacutsval);
585   SetUsePID(kTRUE);
586   SetRemoveDaughtersFromPrim(kTRUE);
587   
588   PrintAll();
589
590   for(Int_t iic=0;iic<nvars;iic++){delete [] anacutsval[iic];}
591   delete [] anacutsval;
592   anacutsval=NULL;
593
594   return;
595 }
596
597
598 void AliRDHFCutsDplustoKpipi::SetStandardCutsPbPb2010() {
599   //
600   //STANDARD CUTS USED FOR 2010 Pb Pb analysis.... not optimized yet 
601   //                                              
602   
603   SetName("DplustoKpipiCutsStandard");
604   SetTitle("Standard Cuts for D+ analysis in PbPb2010 run");
605   
606   // PILE UP REJECTION
607   //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
608
609   // EVENT CUTS
610   SetMinVtxContr(1);
611
612   
613   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
614   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
615   //default
616   esdTrackCuts->SetRequireTPCRefit(kTRUE);
617   esdTrackCuts->SetRequireITSRefit(kTRUE);
618   //esdTrackCuts->SetMinNClustersITS(4); // default is 5
619   esdTrackCuts->SetMinNClustersTPC(70);
620   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
621                                          AliESDtrackCuts::kAny); 
622   // default is kBoth, otherwise kAny
623   esdTrackCuts->SetMinDCAToVertexXY(0.);
624   esdTrackCuts->SetPtRange(0.8,1.e10);
625     
626   AddTrackCuts(esdTrackCuts);
627     
628   const Int_t nptbins=10;
629   Float_t* ptbins;
630   ptbins=new Float_t[nptbins+1];
631
632   ptbins[0]=0.;
633   ptbins[1]=1.;
634   ptbins[2]=2.;
635   ptbins[3]=3.;
636   ptbins[4]=4.;
637   ptbins[5]=5.;
638   ptbins[6]=6.;
639   ptbins[7]=8.;
640   ptbins[8]=12.;
641   ptbins[9]=16.;
642   ptbins[10]=24.;
643   const Int_t nvars=14;
644     
645   Float_t** anacutsval;
646   anacutsval=new Float_t*[nvars];
647   
648   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
649   //Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
650
651   for(Int_t ipt=0;ipt<nptbins;ipt++){
652     anacutsval[0][ipt]=0.2;
653     anacutsval[1][ipt]=0.8;
654     anacutsval[2][ipt]=0.8;
655     anacutsval[3][ipt]=0.;
656     anacutsval[4][ipt]=0.;
657     anacutsval[5][ipt]=0.01;
658     anacutsval[11][ipt]=10000000000.;
659     anacutsval[12][ipt]=0.;
660     anacutsval[13][ipt]=0.;
661   }
662   anacutsval[1][5]=0.9;
663
664   anacutsval[6][0]=0.022100;
665   anacutsval[6][1]=0.022100;
666   anacutsval[6][2]=0.034;
667   anacutsval[6][3]=0.020667;
668   anacutsval[6][4]=0.020667;
669   anacutsval[6][5]=0.023333;
670   
671   anacutsval[7][0]=0.08;
672   anacutsval[7][1]=0.08;
673   anacutsval[7][2]=0.17;  
674   anacutsval[7][3]=0.14;
675   anacutsval[7][4]=0.14;
676   anacutsval[7][5]=0.19;
677   
678   anacutsval[8][0]=0.8;
679   anacutsval[8][1]=0.8;
680   anacutsval[8][2]=1.1;
681   anacutsval[8][3]=0.5;
682   anacutsval[8][4]=0.5;
683   anacutsval[8][5]=0.5;
684   
685   anacutsval[9][0]=0.995;
686   anacutsval[9][1]=0.995;
687   anacutsval[9][2]=0.997; 
688   anacutsval[9][3]=0.998; 
689   anacutsval[9][4]=0.998;
690   anacutsval[9][5]=0.995;
691   
692   anacutsval[10][0]=0.0055;
693   anacutsval[10][1]=0.0055;
694   anacutsval[10][2]= 0.0028;
695   anacutsval[10][3]=0.000883;
696   anacutsval[10][4]=0.000883;
697   anacutsval[10][5]=0.000883;
698
699   anacutsval[12][5]=12.;
700   anacutsval[13][5]=0.998571;
701   anacutsval[12][6]=10.;
702   anacutsval[13][6]=0.997143;
703
704   for(Int_t ipt=6;ipt<nptbins;ipt++){
705     anacutsval[6][ipt]=0.02333;
706     anacutsval[7][ipt]=0.19;
707     anacutsval[8][ipt]=2.0;
708     anacutsval[9][ipt]=0.997;
709     anacutsval[10][ipt]=0.000883;
710   }   
711   anacutsval[7][6]=0.14;
712   anacutsval[9][6]=0.995;
713
714   SetPtBins(nptbins+1,ptbins);
715   SetCuts(nvars,nptbins,anacutsval);
716   SetUsePID(kTRUE);
717   SetMinCentrality(1E-10);
718   SetMaxCentrality(20.);
719   SetUseCentrality(AliRDHFCuts::kCentV0M);
720   SetRemoveDaughtersFromPrim(kFALSE);
721     
722   PrintAll();
723
724   for(Int_t iic=0;iic<nvars;iic++){delete [] anacutsval[iic];}
725   delete [] anacutsval;
726   anacutsval=NULL;
727
728   delete [] ptbins;
729   ptbins=NULL;
730
731   return;
732 }