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