Added IsInFiducualAcceptance() (Alessandro)
[u/mrichter/AliRoot.git] / PWG3 / 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 /////////////////////////////////////////////////////////////
17 //
18 // Class for cuts on AOD reconstructed DStar->Kpipi
19 //
20 // Author: A.Grelli, alessandro.grelli@uu.nl
21 /////////////////////////////////////////////////////////////
22
23 #include <TDatabasePDG.h>
24 #include <Riostream.h>
25
26 #include "AliAODRecoDecayHF2Prong.h"
27 #include "AliAODRecoCascadeHF.h"
28 #include "AliRDHFCutsD0toKpi.h"
29 #include "AliRDHFCutsDStartoKpipi.h"
30 #include "AliAODTrack.h"
31 #include "AliESDtrack.h"
32 #include "AliAODVertex.h"
33 #include "AliESDVertex.h"
34
35
36
37 ClassImp(AliRDHFCutsDStartoKpipi)
38
39 //--------------------------------------------------------------------------
40 AliRDHFCutsDStartoKpipi::AliRDHFCutsDStartoKpipi(const char* name) : 
41 AliRDHFCuts(name),
42 fTrackCutsSoftPi(0)
43 {
44   //
45   // Default Constructor
46   //
47  
48   Int_t nvars=14;
49   SetNVars(nvars);
50   TString varNames[14]={
51     "inv. mass [GeV]",   
52     "dca [cm]",
53     "cosThetaStar", 
54     "pTK [GeV/c]",
55     "pTPi [GeV/c]",
56     "d0K [cm]",
57     "d0Pi [cm]",
58     "d0d0 [cm^2]",
59     "cosThetaPoint",
60     "inv. mass half width of D* [GeV]",
61     "half width of (M_Kpipi-M_D0) [GeV]",
62     "PtMin of pi_s [GeV/c]",
63     "PtMax of pi_s [GeV/c]",
64     "theta, angle between the pi_s and decay plane of the D0 [rad]"};
65   Bool_t isUpperCut[14]={
66     kTRUE,
67     kTRUE,
68     kTRUE,
69     kFALSE,
70     kFALSE,
71     kTRUE,
72     kTRUE,
73     kTRUE,
74     kFALSE,
75     kTRUE,
76     kTRUE,
77     kTRUE,
78     kTRUE,
79     kFALSE};
80   SetVarNames(nvars,varNames,isUpperCut);
81   Bool_t forOpt[14]={
82     kFALSE,
83     kTRUE,
84     kTRUE,
85     kFALSE,
86     kFALSE,
87     kFALSE,
88     kFALSE,
89     kTRUE,
90     kTRUE,
91     kFALSE,
92     kTRUE,
93     kFALSE,
94     kFALSE,
95     kFALSE};
96   SetVarsForOpt(5,forOpt);
97   Float_t limits[2]={0,999999999.};
98   SetPtBins(2,limits);
99 }
100 //--------------------------------------------------------------------------
101 AliRDHFCutsDStartoKpipi::AliRDHFCutsDStartoKpipi(const AliRDHFCutsDStartoKpipi &source) :
102   AliRDHFCuts(source),
103   fTrackCutsSoftPi(0)
104 {
105   //
106   // Copy constructor
107   //
108   
109   if(source.GetTrackCutsSoftPi()) AddTrackCutsSoftPi(source.GetTrackCutsSoftPi());
110   
111 }
112 //--------------------------------------------------------------------------
113 AliRDHFCutsDStartoKpipi &AliRDHFCutsDStartoKpipi::operator=(const AliRDHFCutsDStartoKpipi &source)
114 {
115   //
116   // assignment operator
117   //
118   if(&source == this) return *this;
119
120   AliRDHFCuts::operator=(source);
121
122   if(source.GetTrackCutsSoftPi()) AddTrackCutsSoftPi(source.GetTrackCutsSoftPi());
123   return *this;
124 }
125
126
127 //---------------------------------------------------------------------------
128 void AliRDHFCutsDStartoKpipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
129   // 
130   // Fills in vars the values of the variables 
131   //
132   if(nvars!=fnVarsForOpt) {
133     printf("AliRDHFCutsDStartoKpipi::GetCutsVarsForOpt: wrong number of variables\n");
134     return;
135   }
136   
137  
138   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)d;
139
140   AliAODTrack *softPi = (AliAODTrack*)dstarD0pi->GetBachelor();
141
142   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
143   
144    Int_t iter=-1;
145   if(fVarsForOpt[0]){
146     iter++;
147     if(TMath::Abs(pdgdaughters[0])==211) {
148       vars[iter]=dd->InvMassD0();
149     } else {
150       vars[iter]=dd->InvMassD0bar();
151     }
152   }
153   if(fVarsForOpt[1]){
154     iter++;
155     vars[iter]=dd->GetDCA();
156   }
157   if(fVarsForOpt[2]){
158     iter++;
159     if(TMath::Abs(pdgdaughters[0])==211) {
160       vars[iter] = dd->CosThetaStarD0();
161     } else {
162       vars[iter] = dd->CosThetaStarD0bar();
163     }
164   }
165   if(fVarsForOpt[3]){
166     iter++;
167    if(TMath::Abs(pdgdaughters[0])==321) {
168      vars[iter]=dd->PtProng(0);
169    }
170    else{
171      vars[iter]=dd->PtProng(1);
172    }
173   }
174   if(fVarsForOpt[4]){
175     iter++;
176    if(TMath::Abs(pdgdaughters[0])==211) {
177      vars[iter]=dd->PtProng(0);
178    }
179    else{
180      vars[iter]=dd->PtProng(1);
181    }
182   }
183   if(fVarsForOpt[5]){
184     iter++;
185     if(TMath::Abs(pdgdaughters[0])==321) {
186      vars[iter]=dd->Getd0Prong(0);
187    }
188    else{
189      vars[iter]=dd->Getd0Prong(1);
190    }
191   }
192   if(fVarsForOpt[6]){
193     iter++;
194      if(TMath::Abs(pdgdaughters[0])==211) {
195      vars[iter]=dd->Getd0Prong(0);
196    }
197    else{
198      vars[iter]=dd->Getd0Prong(1);
199    }
200   }
201   if(fVarsForOpt[7]){
202     iter++;
203     vars[iter]= dd->Prodd0d0();
204   }
205   if(fVarsForOpt[8]){
206     iter++;
207     vars[iter]=dd->CosPointingAngle();
208   }
209   if(fVarsForOpt[9]){
210     iter++;
211     vars[iter]=dstarD0pi->InvMassDstarKpipi();
212   }
213   if(fVarsForOpt[10]){
214     iter++;
215     vars[iter]=dstarD0pi->DeltaInvMass();
216   }
217   if(fVarsForOpt[11]){
218     iter++;
219     vars[iter] = softPi->Pt();
220   }
221   if(fVarsForOpt[12]){
222     iter++;
223     vars[iter] = softPi->Pt();
224   }
225   if(fVarsForOpt[13]){
226     iter++;
227     vars[iter] =dstarD0pi->AngleD0dkpPisoft();
228   }
229  
230   return;
231 }
232 //---------------------------------------------------------------------------
233 Int_t AliRDHFCutsDStartoKpipi::IsSelected(TObject* obj,Int_t selectionLevel) {
234   //
235   // Apply selection for D*.
236   //
237   if(!fCutsRD){
238     cout<<"Cut matrice not inizialized. Exit..."<<endl;
239     return 0;
240   }
241   
242   AliAODRecoCascadeHF* d = (AliAODRecoCascadeHF*)obj;
243   if(!d){
244     cout<<"AliAODRecoCascadeHF null"<<endl;
245     return 0;
246   }
247   
248   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)d->Get2Prong();  
249   if(!dd){
250     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
251     return 0;
252   }
253
254   AliAODTrack *b = (AliAODTrack*)d->GetBachelor();
255
256   // selection on daughter tracks 
257   if(selectionLevel==AliRDHFCuts::kAll || 
258      selectionLevel==AliRDHFCuts::kTracks) {
259     if(!AreDaughtersSelected(dd)) return 0;
260     if(fTrackCutsSoftPi) {
261       AliAODVertex *vAOD = d->GetPrimaryVtx();
262       Double_t pos[3],cov[6];
263       vAOD->GetXYZ(pos);
264       vAOD->GetCovarianceMatrix(cov);
265       const AliESDVertex vESD(pos,cov,100.,100);
266       if(!IsDaughterSelected(b,&vESD,fTrackCutsSoftPi)) return 0;
267     }
268   }
269   
270   Int_t returnvalue=1;
271   
272   // selection on candidate
273   if(selectionLevel==AliRDHFCuts::kAll || 
274      selectionLevel==AliRDHFCuts::kCandidate) {
275     
276     Double_t pt=d->Pt();
277     Int_t ptbin=PtBin(pt);
278
279     // select D0 that passes D* cuts
280     returnvalue = IsD0FromDStarSelected(pt,dd,selectionLevel);
281     //if(retunvalue==0) return 0;
282
283     if((b->Charge()==+1 && returnvalue==2) || (b->Charge()==-1 && returnvalue==1)) return 0; 
284
285     // DStarMass and D0mass
286     Double_t mDSPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
287     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
288     // delta mass PDG
289     Double_t deltaPDG = mDSPDG-mD0PDG;
290    
291     // Half width DStar mass
292     if(TMath::Abs(mDSPDG - (d->InvMassDstarKpipi()))>fCutsRD[GetGlobalIndex(9,ptbin)]) return 0;
293     // Half width Delta mass
294     
295     if(TMath::Abs(deltaPDG-(d->DeltaInvMass())) > fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
296     
297     // cut on soft pion pt
298     if(b->Pt() < fCutsRD[GetGlobalIndex(11,ptbin)] || b->Pt() > fCutsRD[GetGlobalIndex(12,ptbin)]) return 0;
299     // cut on the angle between D0 decay plane and soft pion
300     if(d->AngleD0dkpPisoft() > fCutsRD[GetGlobalIndex(13,ptbin)]) return 0;
301   
302   }
303
304   return returnvalue;
305 }
306 //_________________________________________________________________________________________________
307 Int_t AliRDHFCutsDStartoKpipi::IsD0FromDStarSelected(Double_t pt, TObject* obj,Int_t selectionLevel) const {
308   //
309   // Apply selection for D0 from D*. The selection in on D0 prongs
310   //
311   
312   if(!fCutsRD){
313     cout<<"Cut matrice not inizialized. Exit..."<<endl;
314     return 0;
315   }
316   
317   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)obj;
318   
319   if(!dd){
320     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
321     return 0;
322   }
323
324   // selection on daughter tracks is done in IsSelected()
325   
326   Int_t returnvalue=1;
327   
328   // selection on candidate
329   if(selectionLevel==AliRDHFCuts::kAll || 
330      selectionLevel==AliRDHFCuts::kCandidate) {
331     
332     // D0 mass
333     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
334     // delta mass PDG
335  
336     Int_t ptbin=PtBin(pt);
337     
338     Double_t mD0,mD0bar,ctsD0,ctsD0bar;
339   
340     Int_t okD0     =0;
341     Int_t okD0bar  =0;
342     okD0=1; okD0bar=1;
343
344     if(dd->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
345     if(dd->PtProng(0) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(1) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
346     
347     if(!okD0 && !okD0bar) return 0;
348     
349     if(TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] || 
350        TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
351     if(TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
352        TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
353     if(!okD0 && !okD0bar) return 0;
354     
355     if(dd->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;
356     
357     dd->InvMassD0(mD0,mD0bar);
358     if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
359     if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
360     if(!okD0 && !okD0bar) return 0;
361     
362     dd->CosThetaStarD0(ctsD0,ctsD0bar);
363     if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
364     if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
365     if(!okD0 && !okD0bar) return 0;
366     
367     if(dd->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
368     
369     if(dd->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
370
371     if (okD0) returnvalue=1; //cuts passed as D0
372     if (okD0bar) returnvalue=2; //cuts passed as D0bar
373     if (okD0 && okD0bar) returnvalue=3; //both
374   }
375
376   //if((Charge()==+1 && !okD0) || (Charge()==-1 && !okD0bar)) return kFALSE; 
377   return returnvalue;
378 }
379 //----------------------------------------------------------------------------------
380 Bool_t AliRDHFCutsDStartoKpipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
381 {
382   //
383   // D* fiducial acceptance region 
384   //
385
386   if(pt > 5.) {
387     // applying cut for pt > 5 GeV
388     AliDebug(4,Form("pt of D* = %f (> 5), cutting at |y| < 0.8\n",pt)); 
389     if (TMath::Abs(y) > 0.8){
390       return kFALSE;
391     }
392   } else {
393     // appliying smooth cut for pt < 5 GeV
394     Double_t maxFiducialY = -0.151/15*pt*pt+1.9/15*pt+0.4; 
395     Double_t minFiducialY = 0.151/15*pt*pt-1.9/15*pt-0.4;               
396     AliDebug(4,Form("pt of D* = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
397     if (y < minFiducialY || y > maxFiducialY){
398       return kFALSE;
399     }
400   }
401
402   return kTRUE;
403 }