Dplus task used AliAODPidHF via AliRDHFCutsDplustoKpipi (Renu, Francesco)
[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 /////////////////////////////////////////////////////////////
17 //
18 // Class for cuts on AOD reconstructed D+->Kpipi
19 //
20 // Author: R. Bala, bala@to.infn.it
21 //         G. Ortona, ortona@to.infn.it
22 /////////////////////////////////////////////////////////////
23
24 #include <TDatabasePDG.h>
25 #include <Riostream.h>
26
27 #include "AliRDHFCutsDplustoKpipi.h"
28 #include "AliAODRecoDecayHF3Prong.h"
29 #include "AliAODTrack.h"
30 #include "AliESDtrack.h"
31
32
33 ClassImp(AliRDHFCutsDplustoKpipi)
34
35 //--------------------------------------------------------------------------
36 AliRDHFCutsDplustoKpipi::AliRDHFCutsDplustoKpipi(const char* name) : 
37   AliRDHFCuts(name)
38 {
39   //
40   // Default Constructor
41   //
42   Int_t nvars=12;
43   SetNVars(nvars);
44   TString varNames[12]={"inv. mass [GeV]",
45                         "pTK [GeV/c]",
46                         "pTPi [GeV/c]",
47                         "d0K [cm]   lower limit!",
48                         "d0Pi [cm]  lower limit!",
49                         "dist12 (cm)",
50                         "sigmavert (cm)",
51                         "dist prim-sec (cm)",
52                         "pM=Max{pT1,pT2,pT3} (GeV/c)",
53                         "cosThetaPoint",
54                         "Sum d0^2 (cm^2)",
55                         "dca cut (cm)"};
56   Bool_t isUpperCut[12]={kTRUE,
57                          kFALSE,
58                          kFALSE,
59                          kFALSE,
60                          kFALSE,
61                          kFALSE,
62                          kTRUE,
63                          kFALSE,
64                          kFALSE,
65                          kFALSE,
66                          kFALSE,
67                          kTRUE};
68   SetVarNames(nvars,varNames,isUpperCut);
69   Bool_t forOpt[12]={kFALSE,
70                      kFALSE,
71                      kFALSE,
72                      kFALSE,
73                      kFALSE,
74                      kFALSE,
75                      kTRUE,
76                      kTRUE,
77                      kTRUE,
78                      kTRUE,
79                      kTRUE,
80                      kFALSE};
81   SetVarsForOpt(5,forOpt);
82   Float_t limits[2]={0,999999999.};
83   SetPtBins(2,limits);
84   if(fPidHF)delete fPidHF;
85   fPidHF=new AliAODPidHF();
86   Double_t plim[2]={0.6,0.8};
87   Double_t nsigma[5]={2.,1.,2.,3.,0.};
88   
89   fPidHF->SetPLimit(plim);
90   fPidHF->SetAsym(kTRUE);
91   fPidHF->SetSigma(nsigma);
92   fPidHF->SetMatch(1);
93   fPidHF->SetTPC(1);
94   fPidHF->SetTOF(1);
95   fPidHF->SetITS(0);
96   fPidHF->SetTRD(0);
97   fPidHF->SetCompat(kTRUE);
98
99   
100 }
101 //--------------------------------------------------------------------------
102 AliRDHFCutsDplustoKpipi::AliRDHFCutsDplustoKpipi(const AliRDHFCutsDplustoKpipi &source) :
103   AliRDHFCuts(source)
104 {
105   //
106   // Copy constructor
107   //
108
109 }
110 //--------------------------------------------------------------------------
111 AliRDHFCutsDplustoKpipi &AliRDHFCutsDplustoKpipi::operator=(const AliRDHFCutsDplustoKpipi &source)
112 {
113   //
114   // assignment operator
115   //
116   if(&source == this) return *this;
117
118   AliRDHFCuts::operator=(source);
119
120   return *this;
121 }
122 //
123
124
125 //---------------------------------------------------------------------------
126 void AliRDHFCutsDplustoKpipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
127   // 
128   // Fills in vars the values of the variables 
129   //
130
131
132   if(nvars!=fnVarsForOpt) {
133     printf("AliRDHFCutsDplustoKpipi::GetCutsVarsForOpt: wrong number of variables\n");
134     return;
135   }
136
137   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
138
139   Int_t iter=-1;
140   if(fVarsForOpt[0]){
141     iter++;
142     vars[iter]=dd->InvMassDplus();
143   }
144   if(fVarsForOpt[1]){
145     iter++;
146     for(Int_t iprong=0;iprong<3;iprong++){
147       if(TMath::Abs(pdgdaughters[iprong])==321) {
148         vars[iter]=dd->PtProng(iprong);
149       }
150     }
151   }
152   if(fVarsForOpt[2]){
153     iter++;
154     Float_t minPtDau=1000000.0;
155     for(Int_t iprong=0;iprong<3;iprong++){
156       if(TMath::Abs(pdgdaughters[iprong])==211) {
157         if(dd->PtProng(iprong)<minPtDau){
158           minPtDau=dd->PtProng(iprong);
159         }
160       }
161     }
162     vars[iter]=minPtDau;
163   }
164   if(fVarsForOpt[3]){
165     iter++;
166     for(Int_t iprong=0;iprong<3;iprong++){
167       if(TMath::Abs(pdgdaughters[iprong])==321) {
168         vars[iter]=dd->Getd0Prong(iprong);
169       }
170     }
171   }
172   if(fVarsForOpt[4]){
173     iter++;
174     Float_t minImpParDau=1000000.0;
175     for(Int_t iprong=0;iprong<3;iprong++){
176       if(TMath::Abs(pdgdaughters[iprong])==211) {
177         if(dd->Getd0Prong(iprong)<minImpParDau){
178           minImpParDau=dd->Getd0Prong(iprong);
179         }
180       }
181     }
182    vars[iter]=minImpParDau;
183   }
184   if(fVarsForOpt[5]){
185     iter++;
186     Float_t dist12 = dd->GetDist12toPrim();
187     Float_t dist23 = dd->GetDist23toPrim();
188     if(dist12<dist23)vars[iter]=dist12;
189     else vars[iter]=dist23;
190   }
191   if(fVarsForOpt[6]){
192     iter++;
193     vars[iter]=dd->GetSigmaVert();
194   }
195   if(fVarsForOpt[7]){
196     iter++;
197     vars[iter] = dd->DecayLength();
198   }
199   if(fVarsForOpt[8]){
200     iter++;
201     Float_t ptmax=0;
202     for(Int_t i=0;i<3;i++){
203       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
204     }
205     vars[iter]=ptmax;
206   }
207   if(fVarsForOpt[9]){
208     iter++;
209     vars[iter]=dd->CosPointingAngle();
210   }
211   if(fVarsForOpt[10]){
212     iter++;
213     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
214   }
215   if(fVarsForOpt[11]){
216     iter++;
217     Float_t maxDCA=0;
218     for(Int_t iprong=0;iprong<3;iprong++){
219       if(dd->GetDCA(iprong)<maxDCA){
220         maxDCA=dd->GetDCA(iprong);
221       }
222     }
223     vars[iter]=maxDCA;
224   }
225   return;
226 }
227
228 //---------------------------------------------------------------------------
229 Int_t AliRDHFCutsDplustoKpipi::IsSelectedPID(AliAODRecoDecayHF *rd)const {
230   //PID 
231   if(!fUsePID || !rd) return 1;
232   //if(fUsePID)printf("i am inside the pid \n");
233   Int_t nkaons=0;
234   Int_t nNotKaons=0;
235   Int_t sign= rd->GetCharge(); 
236   for(Int_t daught=0;daught<3;daught++){
237     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(daught);
238     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
239     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
240     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
241     
242     if(isProton>0 &&  isKaon<0  && isPion<0) return 0;
243     if(isKaon>0 && isPion<0) nkaons++;
244     if(isKaon<0) nNotKaons++;  
245     if(sign==track->Charge()){//pions
246       if(isPion<0)return 0;
247     }
248       else{//kaons
249         if(isKaon<0)return 0;
250       }
251     
252       
253   }
254   
255   if(nkaons>1)return 0;
256   if(nNotKaons==3)return 0;
257   
258   return 1;   
259 }
260
261
262
263 //---------------------------------------------------------------------------
264 Int_t AliRDHFCutsDplustoKpipi::IsSelected(TObject* obj,Int_t selectionLevel) {
265   //
266   // Apply selection
267   //
268
269   if(!fCutsRD){
270     cout<<"Cut matrix not inizialized. Exit..."<<endl;
271     return 0;
272   }
273   //PrintAll();
274   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj; 
275
276   
277   if(!d){
278     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
279     return 0;
280   }
281
282
283  
284   // selection on daughter tracks 
285   if(selectionLevel==AliRDHFCuts::kAll || 
286      selectionLevel==AliRDHFCuts::kTracks) {
287     if(!AreDaughtersSelected(d)) return 0;
288   }
289   
290   // PID selection
291   Int_t returnvaluePID=1;  
292                                           
293
294   //if(selectionLevel==AliRDHFCuts::kAll || 
295   if(selectionLevel==AliRDHFCuts::kCandidate ||     
296  selectionLevel==AliRDHFCuts::kPID) {
297     returnvaluePID = IsSelectedPID(d);
298   }
299   if(returnvaluePID==0)return 0;
300   
301   
302   // selection on candidate
303   if(selectionLevel==AliRDHFCuts::kAll || 
304      selectionLevel==AliRDHFCuts::kCandidate) {
305     
306     Double_t pt=d->Pt();
307     
308     Int_t ptbin=PtBin(pt);
309     
310     Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
311     Double_t mDplus=d->InvMassDplus();
312     if(TMath::Abs(mDplus-mDplusPDG)>fCutsRD[GetGlobalIndex(0,ptbin)])return 0;
313     //    if(d->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || d->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
314     if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] || TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)])return 0;//Kaon
315     if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)])return 0;//Pion1
316     if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)])return 0;//Pion2
317
318     
319
320   //2track cuts
321   if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]|| d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)])return 0;
322   if(d->Getd0Prong(0)*d->Getd0Prong(1)<0. && d->Getd0Prong(2)*d->Getd0Prong(1)<0.)return 0;
323
324   //sec vert
325   if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)])return 0;
326
327   if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)])return 0;
328
329   if(TMath::Abs(d->PtProng(0))<fCutsRD[GetGlobalIndex(8,ptbin)] && TMath::Abs(d->PtProng(1))<fCutsRD[GetGlobalIndex(8,ptbin)] && TMath::Abs(d->PtProng(2))<fCutsRD[GetGlobalIndex(8,ptbin)])return 0;
330   if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)])return 0;
331   Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
332   if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)])return 0;
333
334   //DCA
335   for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]) return 0;
336   
337   return 1;
338   }
339   return 1;
340 }
341 //---------------------------------------------------------------------------