Global fit of alignment with E field distortion map (440 parameters)
[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 Bool_t AliRDHFCutsDplustoKpipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
229 {
230   //
231   // Checking if Dplus is in fiducial acceptance region 
232   //
233
234   if(pt > 5.) {
235     // applying cut for pt > 5 GeV
236     AliDebug(2,Form("pt of D+ = %f (> 5), cutting at |y| < 0.8",pt)); 
237     if (TMath::Abs(y) > 0.8) return kFALSE;
238     
239   } else {
240     // appliying smooth cut for pt < 5 GeV
241     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
242     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;         
243     AliDebug(2,Form("pt of D+ = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
244     if (y < minFiducialY || y > maxFiducialY) return kFALSE;    
245   }
246
247   return kTRUE;
248 }
249
250 //---------------------------------------------------------------------------
251 Int_t AliRDHFCutsDplustoKpipi::IsSelectedPID(AliAODRecoDecayHF *rd)
252 {
253   //
254   // PID selection
255   // 
256   if(!fUsePID || !rd) return 1;
257   //if(fUsePID)printf("i am inside the pid \n");
258   Int_t nkaons=0;
259   Int_t nNotKaons=0;
260   Int_t sign= rd->GetCharge(); 
261   for(Int_t daught=0;daught<3;daught++){
262     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(daught);
263     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
264     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
265     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
266     
267     if(isProton>0 &&  isKaon<0  && isPion<0) return 0;
268     if(isKaon>0 && isPion<0) nkaons++;
269     if(isKaon<0) nNotKaons++;  
270     if(sign==track->Charge()){//pions
271       if(isPion<0)return 0;
272     }
273       else{//kaons
274         if(isKaon<0)return 0;
275       }
276     
277       
278   }
279   
280   if(nkaons>1)return 0;
281   if(nNotKaons==3)return 0;
282   
283   return 1;   
284 }
285
286
287
288 //---------------------------------------------------------------------------
289 Int_t AliRDHFCutsDplustoKpipi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
290   //
291   // Apply selection
292   //
293
294   if(!fCutsRD){
295     cout<<"Cut matrix not inizialized. Exit..."<<endl;
296     return 0;
297   }
298   //PrintAll();
299   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj; 
300
301   
302   if(!d){
303     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
304     return 0;
305   }
306
307
308  
309   // selection on daughter tracks 
310   if(selectionLevel==AliRDHFCuts::kAll || 
311      selectionLevel==AliRDHFCuts::kTracks) {
312     if(!AreDaughtersSelected(d)) return 0;
313   }
314   
315   // PID selection
316   Int_t returnvaluePID=1;  
317   Int_t returnvalueCuts=1;
318                                           
319
320   //if(selectionLevel==AliRDHFCuts::kAll || 
321   if(selectionLevel==AliRDHFCuts::kCandidate ||     
322      selectionLevel==AliRDHFCuts::kPID) {
323     returnvaluePID = IsSelectedPID(d);
324   }
325   if(returnvaluePID==0)return 0;
326   
327   
328   // selection on candidate
329   if(selectionLevel==AliRDHFCuts::kAll || 
330      selectionLevel==AliRDHFCuts::kCandidate) {
331     
332     //recalculate vertex w/o daughters
333     AliAODVertex *origownvtx=0x0;
334     AliAODVertex *recvtx=0x0;
335     if(fRemoveDaughtersFromPrimary) {
336       if(!aod) {
337         AliError("Can not remove daughters from vertex without AOD event");
338         return 0;
339       }
340       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
341       recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
342       if(!recvtx){
343         AliDebug(2,"Removal of daughter tracks failed");
344         //recvtx=d->GetPrimaryVtx();
345         if(origownvtx){
346           delete origownvtx;
347           origownvtx=NULL;
348         }
349         return 0;
350       }
351       //set recalculed primary vertex
352       d->SetOwnPrimaryVtx(recvtx);
353       delete recvtx; recvtx=NULL;
354     }
355
356     Double_t pt=d->Pt();
357     
358     Int_t ptbin=PtBin(pt);
359     if (ptbin==-1) {
360       if(origownvtx){
361         d->SetOwnPrimaryVtx(origownvtx);
362         delete origownvtx;
363         origownvtx=NULL;
364       }
365       else d->UnsetOwnPrimaryVtx();
366       return 0;
367     }
368     
369     Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
370     Double_t mDplus=d->InvMassDplus();
371     if(TMath::Abs(mDplus-mDplusPDG)>fCutsRD[GetGlobalIndex(0,ptbin)])returnvalueCuts=0;
372     //    if(d->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || d->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
373     if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] || TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)])returnvalueCuts=0;//Kaon
374     if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)])returnvalueCuts=0;//Pion1
375     if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)])returnvalueCuts=0;//Pion2
376
377     
378
379     //2track cuts
380     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]|| d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)])returnvalueCuts=0;
381     if(d->Getd0Prong(0)*d->Getd0Prong(1)<0. && d->Getd0Prong(2)*d->Getd0Prong(1)<0.)returnvalueCuts=0;
382     
383     //sec vert
384     if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)])returnvalueCuts=0;
385     
386     if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)])returnvalueCuts=0;
387     
388     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)])returnvalueCuts=0;
389     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)])returnvalueCuts=0;
390     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
391     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)])returnvalueCuts=0;
392     
393     //DCA
394     for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]) returnvalueCuts=0;
395     // unset recalculated primary vertex when not needed any more
396     if(origownvtx) {
397       d->SetOwnPrimaryVtx(origownvtx);
398       delete origownvtx;
399       origownvtx=NULL;
400     } else if(fRemoveDaughtersFromPrimary) {
401       d->UnsetOwnPrimaryVtx();
402       AliDebug(3,"delete new vertex\n");
403     }
404     
405     return returnvalueCuts;
406   }
407
408   return 1;
409 }
410 //---------------------------------------------------------------------------