1fe5d9167f36b7ec7afc3c251a71fc58209c536f
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCutsDStartoKpipi.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 DStar->Kpipi
21 //
22 // Author: A.Grelli, alessandro.grelli@uu.nl
23 //
24 // PID method implemented by   Y.Wang, yifei@physi.uni-heidelberg.de
25 //           
26 /////////////////////////////////////////////////////////////
27
28 #include <TDatabasePDG.h>
29 #include <Riostream.h>
30
31 #include "AliAODRecoDecayHF2Prong.h"
32 #include "AliAODRecoCascadeHF.h"
33 #include "AliRDHFCutsD0toKpi.h"
34 #include "AliRDHFCutsDStartoKpipi.h"
35 #include "AliAODTrack.h"
36 #include "AliESDtrack.h"
37 #include "AliAODPid.h"
38 #include "AliTPCPIDResponse.h"
39 #include "AliAODVertex.h"
40 #include "AliESDVertex.h"
41
42 ClassImp(AliRDHFCutsDStartoKpipi)
43
44 //--------------------------------------------------------------------------
45 AliRDHFCutsDStartoKpipi::AliRDHFCutsDStartoKpipi(const char* name) : 
46   AliRDHFCuts(name),
47   fTrackCutsSoftPi(0)
48 {
49   //
50   // Default Constructor
51   //
52   
53   Int_t nvars=16;
54   SetNVars(nvars);
55   TString varNames[16]={
56     "inv. mass [GeV]",   
57     "dca [cm]",
58     "cosThetaStar", 
59     "pTK [GeV/c]",
60     "pTPi [GeV/c]",
61     "d0K [cm]",
62     "d0Pi [cm]",
63     "d0d0 [cm^2]",
64     "cosThetaPoint",
65     "inv. mass half width of D* [GeV]",
66     "half width of (M_Kpipi-M_D0) [GeV]",
67     "PtMin of pi_s [GeV/c]",
68     "PtMax of pi_s [GeV/c]",
69     "theta, angle between the pi_s and decay plane of the D0 [rad]",
70     "|cosThetaPointXY|", 
71     "NormDecayLenghtXY"};
72   Bool_t isUpperCut[16]={
73     kTRUE,
74     kTRUE,
75     kTRUE,
76     kFALSE,
77     kFALSE,
78     kTRUE,
79     kTRUE,
80     kTRUE,
81     kFALSE,
82     kTRUE,
83     kTRUE,
84     kTRUE,
85     kTRUE,
86     kFALSE,
87     kFALSE, 
88     kFALSE};
89   SetVarNames(nvars,varNames,isUpperCut);
90   Bool_t forOpt[16]={
91     kFALSE,
92     kTRUE,
93     kTRUE,
94     kFALSE,
95     kFALSE,
96     kFALSE,
97     kFALSE,
98     kTRUE,
99     kTRUE,
100     kFALSE,
101     kTRUE,
102     kFALSE,
103     kFALSE,
104     kFALSE,
105     kFALSE,
106     kFALSE};
107   SetVarsForOpt(5,forOpt);
108   Float_t limits[2]={0,999999999.};
109   SetPtBins(2,limits);
110 }
111 //--------------------------------------------------------------------------
112 AliRDHFCutsDStartoKpipi::AliRDHFCutsDStartoKpipi(const AliRDHFCutsDStartoKpipi &source) :
113   AliRDHFCuts(source),
114   fTrackCutsSoftPi(0)
115 {
116   //
117   // Copy constructor
118   //
119   
120   if(source.GetTrackCutsSoftPi()) AddTrackCutsSoftPi(source.GetTrackCutsSoftPi());
121   
122 }
123 //--------------------------------------------------------------------------
124 AliRDHFCutsDStartoKpipi &AliRDHFCutsDStartoKpipi::operator=(const AliRDHFCutsDStartoKpipi &source)
125 {
126   //
127   // assignment operator
128   //
129   if(&source == this) return *this;
130
131   AliRDHFCuts::operator=(source);
132   if(source.GetTrackCutsSoftPi()) {
133     delete fTrackCutsSoftPi;
134     fTrackCutsSoftPi = new AliESDtrackCuts(*(source.GetTrackCutsSoftPi()));
135   }
136
137   return *this;
138 }
139
140
141 //---------------------------------------------------------------------------
142 void AliRDHFCutsDStartoKpipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
143   // 
144   // Fills in vars the values of the variables 
145   //
146   if(nvars!=fnVarsForOpt) {
147     printf("AliRDHFCutsDStartoKpipi::GetCutsVarsForOpt: wrong number of variables\n");
148     return;
149   }
150   
151  
152   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)d;
153
154   AliAODTrack *softPi = (AliAODTrack*)dstarD0pi->GetBachelor();
155
156   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
157   
158    Int_t iter=-1;
159   if(fVarsForOpt[0]){
160     iter++;
161     if(TMath::Abs(pdgdaughters[0])==211) {
162       vars[iter]=dd->InvMassD0();
163     } else {
164       vars[iter]=dd->InvMassD0bar();
165     }
166   }
167   if(fVarsForOpt[1]){
168     iter++;
169     vars[iter]=dd->GetDCA();
170   }
171   if(fVarsForOpt[2]){
172     iter++;
173     if(TMath::Abs(pdgdaughters[0])==211) {
174       vars[iter] = dd->CosThetaStarD0();
175     } else {
176       vars[iter] = dd->CosThetaStarD0bar();
177     }
178   }
179   if(fVarsForOpt[3]){
180     iter++;
181    if(TMath::Abs(pdgdaughters[0])==321) {
182      vars[iter]=dd->PtProng(0);
183    }
184    else{
185      vars[iter]=dd->PtProng(1);
186    }
187   }
188   if(fVarsForOpt[4]){
189     iter++;
190    if(TMath::Abs(pdgdaughters[0])==211) {
191      vars[iter]=dd->PtProng(0);
192    }
193    else{
194      vars[iter]=dd->PtProng(1);
195    }
196   }
197   if(fVarsForOpt[5]){
198     iter++;
199     if(TMath::Abs(pdgdaughters[0])==321) {
200      vars[iter]=dd->Getd0Prong(0);
201    }
202    else{
203      vars[iter]=dd->Getd0Prong(1);
204    }
205   }
206   if(fVarsForOpt[6]){
207     iter++;
208      if(TMath::Abs(pdgdaughters[0])==211) {
209      vars[iter]=dd->Getd0Prong(0);
210    }
211    else{
212      vars[iter]=dd->Getd0Prong(1);
213    }
214   }
215   if(fVarsForOpt[7]){
216     iter++;
217     vars[iter]= dd->Prodd0d0();
218   }
219   if(fVarsForOpt[8]){
220     iter++;
221     vars[iter]=dd->CosPointingAngle();
222   }
223   if(fVarsForOpt[9]){
224     iter++;
225     vars[iter]=dstarD0pi->InvMassDstarKpipi();
226   }
227   if(fVarsForOpt[10]){
228     iter++;
229     vars[iter]=dstarD0pi->DeltaInvMass();
230   }
231   if(fVarsForOpt[11]){
232     iter++;
233     vars[iter] = softPi->Pt();
234   }
235   if(fVarsForOpt[12]){
236     iter++;
237     vars[iter] = softPi->Pt();
238   }
239   if(fVarsForOpt[13]){
240     iter++;
241     vars[iter] =dstarD0pi->AngleD0dkpPisoft();
242   }
243   if(fVarsForOpt[14]){
244     iter++;
245     vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
246   }
247   if(fVarsForOpt[15]){
248     iter++;
249     vars[iter]=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
250   }
251  
252   return;
253 }
254 //---------------------------------------------------------------------------
255 Int_t AliRDHFCutsDStartoKpipi::IsSelected(TObject* obj,Int_t selectionLevel) {
256   //
257   // Apply selection for D*.
258   //
259   if(!fCutsRD){
260     cout<<"Cut matrice not inizialized. Exit..."<<endl;
261     return 0;
262   }
263   
264   AliAODRecoCascadeHF* d = (AliAODRecoCascadeHF*)obj;
265   if(!d){
266     cout<<"AliAODRecoCascadeHF null"<<endl;
267     return 0;
268   }
269   
270   Double_t ptD=d->Pt();
271   if(ptD<fMinPtCand) return 0;
272   if(ptD>fMaxPtCand) return 0;
273   
274
275   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)d->Get2Prong();  
276   if(!dd){
277     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
278     return 0;
279   }
280
281   if(dd->HasBadDaughters()) return 0;
282
283   AliAODTrack *b = (AliAODTrack*)d->GetBachelor();
284   if(fTrackCutsSoftPi && fTrackCutsSoftPi->GetRequireTPCRefit()){
285     if(!(b->TestFilterMask(BIT(4)))) return 0;
286   }
287   
288   Int_t returnvalue=1;
289   Int_t returnvaluePID=3;
290
291
292   // selection on candidate
293   if(selectionLevel==AliRDHFCuts::kAll || 
294      selectionLevel==AliRDHFCuts::kCandidate) {
295     
296     Double_t pt=d->Pt();
297     Int_t ptbin=PtBin(pt);
298  
299     // DStarMass and D0mass
300     Double_t mDSPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
301     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
302     // delta mass PDG
303     Double_t deltaPDG = mDSPDG-mD0PDG;
304    
305     // Half width DStar mass
306     if(TMath::Abs(mDSPDG - (d->InvMassDstarKpipi()))>fCutsRD[GetGlobalIndex(9,ptbin)]) return 0;
307     // Half width Delta mass
308     
309     if(TMath::Abs(deltaPDG-(d->DeltaInvMass())) > fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
310     
311     // cut on soft pion pt
312     if(b->Pt() < fCutsRD[GetGlobalIndex(11,ptbin)] || b->Pt() > fCutsRD[GetGlobalIndex(12,ptbin)]) return 0;
313     // cut on the angle between D0 decay plane and soft pion
314     if(d->AngleD0dkpPisoft() > fCutsRD[GetGlobalIndex(13,ptbin)]) return 0;
315   
316     // select D0 that passes D* cuts
317     returnvalue = IsD0FromDStarSelected(pt,dd,selectionLevel);
318     if((b->Charge()==+1 && returnvalue==2) || (b->Charge()==-1 && returnvalue==1)) return 0; 
319     
320   }
321
322   // selection on PID 
323   if(selectionLevel==AliRDHFCuts::kAll || 
324      selectionLevel==AliRDHFCuts::kCandidate ||
325      selectionLevel==AliRDHFCuts::kPID) {
326     returnvaluePID = IsSelectedPID(d);
327   }
328   if(returnvaluePID!=3) returnvalue =0;
329
330
331   // selection on daughter tracks 
332   if(selectionLevel==AliRDHFCuts::kAll || 
333      selectionLevel==AliRDHFCuts::kTracks) {
334     if(!AreDaughtersSelected(dd)) return 0;
335     if(fTrackCutsSoftPi) {
336       AliAODVertex *vAOD = d->GetPrimaryVtx();
337       Double_t pos[3],cov[6];
338       vAOD->GetXYZ(pos);
339       vAOD->GetCovarianceMatrix(cov);
340       const AliESDVertex vESD(pos,cov,100.,100);
341       if(!IsDaughterSelected(b,&vESD,fTrackCutsSoftPi)) return 0;
342     }
343   }
344   
345   return returnvalue;
346
347 }
348 //_________________________________________________________________________________________________
349 Int_t AliRDHFCutsDStartoKpipi::IsD0FromDStarSelected(Double_t pt, TObject* obj,Int_t selectionLevel) const {
350   //
351   // Apply selection for D0 from D*. The selection in on D0 prongs
352   //
353   
354   if(!fCutsRD){
355     cout<<"Cut matrice not inizialized. Exit..."<<endl;
356     return 0;
357   }
358   
359   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)obj;
360   
361   if(!dd){
362     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
363     return 0;
364   }
365
366   // selection on daughter tracks is done in IsSelected()
367   
368   Int_t returnvalue=1;
369   
370   // selection on candidate
371   if(selectionLevel==AliRDHFCuts::kAll || 
372      selectionLevel==AliRDHFCuts::kCandidate) {
373     
374     // D0 mass
375     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
376     // delta mass PDG
377  
378     Int_t ptbin=PtBin(pt);
379     
380     Double_t mD0,mD0bar,ctsD0,ctsD0bar;
381   
382     Int_t okD0     =0;
383     Int_t okD0bar  =0;
384     okD0=1; okD0bar=1;
385
386     if(dd->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
387     if(dd->PtProng(0) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(1) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
388     
389     if(!okD0 && !okD0bar) return 0;
390     
391     if(TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] || 
392        TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
393     if(TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
394        TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
395     if(!okD0 && !okD0bar) return 0;
396     
397     if(dd->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;
398     
399     dd->InvMassD0(mD0,mD0bar);
400     if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
401     if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
402     if(!okD0 && !okD0bar) return 0;
403     
404     dd->CosThetaStarD0(ctsD0,ctsD0bar);
405     if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
406     if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
407     if(!okD0 && !okD0bar) return 0;
408     
409     if(dd->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
410     
411     if(dd->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
412
413     if(TMath::Abs(dd->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(14,ptbin)])  return 0;
414         
415     Double_t normalDecayLengXY=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
416     if (normalDecayLengXY < fCutsRD[GetGlobalIndex(15, ptbin)]) return 0;
417
418     if (okD0) returnvalue=1; //cuts passed as D0
419     if (okD0bar) returnvalue=2; //cuts passed as D0bar
420     if (okD0 && okD0bar) returnvalue=3; //both
421   }
422  
423   return returnvalue;
424 }
425 //----------------------------------------------------------------------------------
426 Bool_t AliRDHFCutsDStartoKpipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
427 {
428   //
429   // D* fiducial acceptance region 
430   //
431
432   if(pt > 5.) {
433     // applying cut for pt > 5 GeV
434     AliDebug(4,Form("pt of D* = %f (> 5), cutting at |y| < 0.8\n",pt)); 
435     if (TMath::Abs(y) > 0.8){
436       return kFALSE;
437     }
438   } else {    
439     // appliying smooth cut for pt < 5 GeV
440     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
441     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;         
442     AliDebug(2,Form("pt of D* = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
443     if (y < minFiducialY || y > maxFiducialY){
444       return kFALSE;
445     }
446   }
447     
448   return kTRUE;
449 }
450
451 //_______________________________________________________________________________-
452 Int_t AliRDHFCutsDStartoKpipi::IsSelectedPID(AliAODRecoDecayHF* obj)
453 {
454   //
455   // PID method, n signa approach default
456   //
457   
458   if(!fUsePID) return 3;
459   
460   AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)obj;
461   if(!dstar){
462     cout<<"AliAODRecoCascadeHF null"<<endl;
463     return 0;
464   }  
465   AliAODRecoDecayHF2Prong* d0 = (AliAODRecoDecayHF2Prong*)dstar->Get2Prong();  
466   if(!d0){
467     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
468     return 0;
469   }
470
471   //  here the PID
472   AliAODTrack *pos = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(0);
473   AliAODTrack *neg = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(1);
474
475   if (dstar->Charge()>0){
476     if(!SelectPID(pos,2)) return 0;//pion+
477     if(!SelectPID(neg,3)) return 0;//kaon-
478   }else{
479     if(!SelectPID(pos,3)) return 0;//kaon+
480     if(!SelectPID(neg,2)) return 0;//pion-
481   }
482
483   return 3;
484 }
485
486 //_______________________________________________________________________________-
487 Int_t AliRDHFCutsDStartoKpipi::SelectPID(AliAODTrack *track, Int_t type)
488 {
489   //
490   //  here the PID
491     
492   Bool_t isParticle=kTRUE;
493
494   if(fPidHF->GetMatch()==1){//n-sigma
495     Bool_t TPCon=TMath::Abs(2)>1e-4?kTRUE:kFALSE;
496     Bool_t TOFon=TMath::Abs(3)>1e-4?kTRUE:kFALSE;
497     
498     Bool_t isTPC=kTRUE;
499     Bool_t isTOF=kTRUE;
500
501     if (TPCon){//TPC
502       if(fPidHF->CheckStatus(track,"TPC")){
503         if(type==2) isTPC=fPidHF->IsPionRaw(track,"TPC");
504         if(type==3) isTPC=fPidHF->IsKaonRaw(track,"TPC");
505       }
506     }
507     if (TOFon){//TOF
508       if(fPidHF->CheckStatus(track,"TOF")){
509         if(type==2) isTOF=fPidHF->IsPionRaw(track,"TOF");
510         if(type==3) isTOF=fPidHF->IsKaonRaw(track,"TOF");
511       }
512     }
513     isParticle = isTPC&&isTOF;
514   }
515   
516   if(fPidHF->GetMatch()==2){//bayesian
517     //Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
518     Double_t prob[5]={1.,1.,1.,1.,1.};
519     
520     //fPidHF->SetPriors(priors);
521     fPidHF->BayesianProbability(track,prob);
522     
523     Double_t max=0.;
524     Int_t k=-1;
525     for (Int_t i=0; i<5; i++) {
526       if (prob[i]>max) {k=i; max=prob[i];}
527     }
528     isParticle = Bool_t(k==type);
529   }
530   
531   return isParticle;
532   
533 }
534 //__________________________________________________________________________________-
535 void  AliRDHFCutsDStartoKpipi::SetStandardCutsPP2010() {
536   //
537   //STANDARD CUTS USED FOR 2010 pp analysis 
538   //                                           
539   // Need to be updated for the final cut version
540   //
541
542   SetName("DStartoD0piCutsStandard");
543   SetTitle("Standard Cuts for D* analysis");
544   
545   // PILE UP REJECTION
546   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
547
548   // EVENT CUTS
549   SetMinVtxContr(1);
550   
551   // CUTS ON SINGLE TRACKS
552   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
553   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
554   esdTrackCuts->SetRequireTPCRefit(kTRUE);
555   esdTrackCuts->SetRequireITSRefit(kTRUE);
556   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
557   esdTrackCuts->SetMinDCAToVertexXY(0.);
558   esdTrackCuts->SetEtaRange(-0.8,0.8);
559   esdTrackCuts->SetPtRange(0.3,1.e10);
560   
561   // CUTS on SOFT PION
562   AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts();
563   esdSoftPicuts->SetRequireSigmaToVertex(kFALSE);
564   esdSoftPicuts->SetRequireTPCRefit(kFALSE);
565   esdSoftPicuts->SetRequireITSRefit(kFALSE);
566   esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
567                                           AliESDtrackCuts::kAny); 
568   esdSoftPicuts->SetPtRange(0.0,1.e10);
569
570   AddTrackCuts(esdTrackCuts);
571   AddTrackCutsSoftPi(esdSoftPicuts);
572
573   const Int_t nptbins =13;
574   const Double_t ptmax = 9999.;
575   const Int_t nvars=16;
576   Float_t ptbins[nptbins+1];
577   ptbins[0]=0.;
578   ptbins[1]=0.5;        
579   ptbins[2]=1.;
580   ptbins[3]=2.;
581   ptbins[4]=3.;
582   ptbins[5]=4.;
583   ptbins[6]=5.;
584   ptbins[7]=6.;
585   ptbins[8]=7.;
586   ptbins[9]=8.;
587   ptbins[10]=12.;
588   ptbins[11]=16.;
589   ptbins[12]=24.;
590   ptbins[13]=ptmax;
591
592   SetGlobalIndex(nvars,nptbins);
593   SetPtBins(nptbins+1,ptbins);
594   
595   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-2000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* pt<0.5*/
596                                                   {0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-16000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* 0.5<pt<1*/
597                                                   {0.7,400.*1E-4,0.8,0.7,0.7,400.*1E-4,400.*1E-4,-36000.*1E-8,0.82,0.3,0.1,0.05,100,0.5,-1.,0.},/* 1<pt<2 */
598                                                   {0.7,200.*1E-4,0.8,0.7,0.7,800.*1E-4,800.*1E-4,-16000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 2<pt<3 */
599                                                   {0.7,500.*1E-4,0.8,1.0,1.0,420.*1E-4,560.*1E-4,-6500.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 3<pt<4 */
600                                                   {0.7,800.*1E-4,0.9,1.2,1.2,700.*1E-4,700.*1E-4,1000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 4<pt<5 */
601                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,800.*1E-4,800.*1E-4,50000.*1E-8,0.8,0.3,0.1,0.05,100,0.5,-1.,0.},/* 5<pt<6 */
602                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 6<pt<7 */
603                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 7<pt<8 */
604                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,600000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 8<pt<12 */
605                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 12<pt<16 */
606                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 16<pt<20 */
607                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.}};/* pt>24 */
608   
609   
610   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
611   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
612   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
613   
614   for (Int_t ibin=0;ibin<nptbins;ibin++){
615     for (Int_t ivar = 0; ivar<nvars; ivar++){
616       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];      
617     }
618   }
619   
620   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
621
622   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
623   delete [] cutsMatrixTransposeStand;
624   cutsMatrixTransposeStand=NULL;
625
626   // PID SETTINGS FOR D* analysis
627   AliAODPidHF* pidObj=new AliAODPidHF();
628   //pidObj->SetName("pid4DSatr");
629   Int_t mode=1;
630   Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
631   pidObj->SetPriors(priors);
632   pidObj->SetMatch(mode);
633   pidObj->SetSigma(0,2); // TPC
634   pidObj->SetSigma(3,3); // TOF
635   pidObj->SetTPC(kTRUE);
636   pidObj->SetTOF(kTRUE);
637   
638   SetPidHF(pidObj);
639   SetUsePID(kTRUE);
640
641   PrintAll();
642
643   delete pidObj;
644   pidObj=NULL;
645
646   return;
647 }
648 //_____________________________________________________________________________-
649 void  AliRDHFCutsDStartoKpipi::SetStandardCutsPbPb2010(){  
650   //
651   // TEMPORARY, WORK IN PROGRESS ... BUT WORKING! 
652   //
653   //  Lead Lead
654   //
655
656   SetName("DStartoD0piCutsStandard");
657   SetTitle("Standard Cuts for D* analysis in PbPb 2010");
658
659   // EVENT CUTS
660   SetMinVtxContr(1);
661   
662   // CUTS ON SINGLE TRACKS
663   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
664   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
665   esdTrackCuts->SetRequireTPCRefit(kTRUE);
666   esdTrackCuts->SetRequireITSRefit(kTRUE);
667   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
668   esdTrackCuts->SetMinDCAToVertexXY(0.);
669   esdTrackCuts->SetEtaRange(-0.8,0.8);
670   esdTrackCuts->SetPtRange(0.3,1.e10);
671   
672   // CUTS on SOFT PION
673   AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts();
674   esdSoftPicuts->SetRequireSigmaToVertex(kFALSE);
675   esdSoftPicuts->SetRequireTPCRefit(kTRUE);
676   esdSoftPicuts->SetRequireITSRefit(kTRUE);
677   esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
678                                           AliESDtrackCuts::kAny); //test d0 asimmetry
679   esdSoftPicuts->SetPtRange(0.25,5);
680
681   AddTrackCuts(esdTrackCuts);
682   AddTrackCutsSoftPi(esdSoftPicuts);
683
684   const Int_t nptbins =13;
685   const Double_t ptmax = 9999.;
686   const Int_t nvars=16;
687   Float_t ptbins[nptbins+1];
688   ptbins[0]=0.;
689   ptbins[1]=0.5;        
690   ptbins[2]=1.;
691   ptbins[3]=2.;
692   ptbins[4]=3.;
693   ptbins[5]=4.;
694   ptbins[6]=5.;
695   ptbins[7]=6.;
696   ptbins[8]=7.;
697   ptbins[9]=8.;
698   ptbins[10]=12.;
699   ptbins[11]=16.;
700   ptbins[12]=24.;
701   ptbins[13]=ptmax;
702
703   SetGlobalIndex(nvars,nptbins);
704   SetPtBins(nptbins+1,ptbins);
705   
706   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-2000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* pt<0.5*/
707                                                   {0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-16000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* 0.5<pt<1*/
708                                                   {0.7,400.*1E-4,0.8,0.7,0.7,800.*1E-4,800.*1E-4,-36000.*1E-8,0.82,0.3,0.1,0.05,100,0.5,-1.,0.},/* 1<pt<2 */
709                                                   {0.7,200.*1E-4,0.8,0.7,0.7,800.*1E-4,800.*1E-4,-16000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 2<pt<3 */
710                                                   {0.7,500.*1E-4,0.8,1.0,1.0,420.*1E-4,560.*1E-4,-6500.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 3<pt<4 */
711                                                   {0.7,800.*1E-4,0.9,1.2,1.2,700.*1E-4,700.*1E-4,1000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 4<pt<5 */
712                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,800.*1E-4,800.*1E-4,50000.*1E-8,0.8,0.3,0.1,0.05,100,0.5,-1.,0.},/* 5<pt<6 */
713                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 6<pt<7 */
714                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 7<pt<8 */
715                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,600000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 8<pt<12 */
716                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 12<pt<16 */
717                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 16<pt<24 */
718                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.}};/* pt>24 */
719   
720   
721   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
722   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
723   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
724   
725   for (Int_t ibin=0;ibin<nptbins;ibin++){
726     for (Int_t ivar = 0; ivar<nvars; ivar++){
727       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];      
728     }
729   }
730   
731   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
732
733   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
734   delete [] cutsMatrixTransposeStand;
735   cutsMatrixTransposeStand=NULL;
736   
737   // PID SETTINGS
738   AliAODPidHF* pidObj=new AliAODPidHF();
739   // pidObj->SetName("pid4DSatr");
740   Int_t mode=1;
741   Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
742   pidObj->SetPriors(priors);
743   pidObj->SetMatch(mode);
744   pidObj->SetSigma(0,2); // TPC
745   pidObj->SetSigma(3,3); // TOF
746   pidObj->SetTPC(kTRUE);
747   pidObj->SetTOF(kTRUE);
748   
749   SetPidHF(pidObj);
750   SetUsePID(kTRUE);
751
752   PrintAll();
753
754   delete pidObj;
755   pidObj=NULL;
756
757   return;
758
759 }
760
761 //_____________________________________________________________________________
762 void  AliRDHFCutsDStartoKpipi::SetStandardCutsPbPb2011(){
763
764   // Default 2010 PbPb cut object
765   SetStandardCutsPbPb2010();
766
767   // Enable all 2011 PbPb run triggers
768   //  
769   SetTriggerClass("");
770   ResetMaskAndEnableMBTrigger();
771   EnableCentralTrigger();
772   EnableSemiCentralTrigger();
773
774 }