]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG3/vertexingHF/AliRDHFCutsDstoKKpi.cxx
Rejection of centrality outliers in AOD049 (Giacomo)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsDstoKKpi.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 Ds->KKpi
21 //
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <TDatabasePDG.h>
26 #include <Riostream.h>
27
28 #include "AliRDHFCutsDstoKKpi.h"
29 #include "AliAODRecoDecayHF3Prong.h"
30 #include "AliAODTrack.h"
31 #include "AliESDtrack.h"
32
33 ClassImp(AliRDHFCutsDstoKKpi)
34
35 //--------------------------------------------------------------------------
36 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const char* name) : 
37 AliRDHFCuts(name)
38 {
39   //
40   // Default Constructor
41   //
42   Int_t nvars=14;
43   SetNVars(nvars);
44   TString varNames[14]={"inv. mass [GeV]",   
45                         "pTK [GeV/c]",
46                         "pTPi [GeV/c]",
47                         "d0K [cm]",
48                         "d0Pi [cm]",
49                         "dist12 [cm]",
50                         "sigmavert [cm]",
51                         "decLen [cm]",
52                         "ptMax [GeV/c]",
53                         "cosThetaPoint",
54                         "Sum d0^2 (cm^2)",
55                         "dca [cm]",
56                         "inv. mass (Mphi-MKK) [GeV]",
57                         "inv. mass (MKo*-MKpi) [GeV]"};
58   Bool_t isUpperCut[14]={kTRUE,
59                          kFALSE,
60                          kFALSE,
61                          kFALSE,
62                          kFALSE,
63                          kFALSE,
64                          kTRUE,
65                          kFALSE,
66                          kFALSE,
67                          kFALSE,
68                          kFALSE,
69                          kTRUE,
70                          kTRUE,
71                          kTRUE};
72   SetVarNames(14,varNames,isUpperCut);
73   Bool_t forOpt[14]={kFALSE,
74                     kFALSE,
75                     kFALSE,
76                     kFALSE,
77                     kFALSE,
78                     kFALSE,
79                     kTRUE,
80                     kTRUE,
81                     kTRUE,
82                     kTRUE,
83                     kTRUE,
84                     kFALSE,
85                     kTRUE,
86                     kTRUE};
87   SetVarsForOpt(7,forOpt);
88   Float_t limits[2]={0,999999999.};
89   SetPtBins(2,limits);
90   if(fPidHF)delete fPidHF;
91   fPidHF=new AliAODPidHF();
92   Double_t plim[2]={0.6,0.8};
93   Double_t nsigma[5]={2.,1.,2.,3.,0.};
94   
95   fPidHF->SetPLimit(plim);
96   fPidHF->SetAsym(kTRUE);
97   fPidHF->SetSigma(nsigma);
98   fPidHF->SetMatch(1);
99   fPidHF->SetTPC(1);
100   fPidHF->SetTOF(1);
101   fPidHF->SetITS(0);
102   fPidHF->SetTRD(0);
103   fPidHF->SetCompat(kTRUE);
104
105 }
106 //--------------------------------------------------------------------------
107 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const AliRDHFCutsDstoKKpi &source) :
108   AliRDHFCuts(source)
109 {
110   //
111   // Copy constructor
112   //
113
114 }
115 //--------------------------------------------------------------------------
116 AliRDHFCutsDstoKKpi &AliRDHFCutsDstoKKpi::operator=(const AliRDHFCutsDstoKKpi &source)
117 {
118   //
119   // assignment operator
120   //
121   if(&source == this) return *this;
122
123   AliRDHFCuts::operator=(source);
124
125   return *this;
126 }
127
128
129 //---------------------------------------------------------------------------
130 void AliRDHFCutsDstoKKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
131   // 
132   // Fills in vars the values of the variables 
133   //
134
135   if(nvars!=fnVarsForOpt) {
136     printf("AliRDHFCutsDstoKKpi::GetCutsVarsForOpt: wrong number of variables\n");
137     return;
138   }
139
140   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
141  
142   Int_t iter=-1;
143   if(fVarsForOpt[0]){
144     iter++;
145     if(TMath::Abs(pdgdaughters[0])==321){
146       vars[iter]=dd->InvMassDsKKpi();
147     }else{
148       vars[iter]=dd->InvMassDspiKK();
149     }
150   }
151   if(fVarsForOpt[1]){
152     iter++;
153     Float_t minPtDau=99999.;
154     for(Int_t iprong=0;iprong<3;iprong++){
155       if(TMath::Abs(pdgdaughters[iprong])==321 && 
156          dd->PtProng(iprong)<minPtDau) minPtDau=dd->PtProng(iprong);
157     }
158     vars[iter]=minPtDau;
159   }
160   if(fVarsForOpt[2]){
161     iter++;
162     for(Int_t iprong=0;iprong<3;iprong++){
163       if(TMath::Abs(pdgdaughters[iprong])==211) {
164         vars[iter]=dd->PtProng(iprong);
165       }
166     }
167   }
168   if(fVarsForOpt[3]){
169     iter++;
170     Float_t minImpParDau=99999.;
171     for(Int_t iprong=0;iprong<3;iprong++){
172       if(TMath::Abs(pdgdaughters[iprong])==321 &&
173          dd->Getd0Prong(iprong)<minImpParDau) minImpParDau=dd->Getd0Prong(iprong);
174     }
175     vars[iter]=minImpParDau;
176   }
177   if(fVarsForOpt[4]){
178     iter++;
179     for(Int_t iprong=0;iprong<3;iprong++){
180       if(TMath::Abs(pdgdaughters[iprong])==211) {
181         vars[iter]=dd->Getd0Prong(iprong);
182       }
183     }
184   }
185   if(fVarsForOpt[5]){
186     iter++;
187     Float_t minDistPair=TMath::Min(dd->GetDist12toPrim(),dd->GetDist23toPrim());
188     vars[iter]=minDistPair;
189   }
190   if(fVarsForOpt[6]){
191     iter++;
192     vars[iter]=dd->GetSigmaVert();
193   }
194   if(fVarsForOpt[7]){
195     iter++;
196     vars[iter] = dd->DecayLength();
197   }
198   if(fVarsForOpt[8]){
199     iter++;
200     Float_t ptmax=0;
201     for(Int_t i=0;i<3;i++){
202       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
203     }
204     vars[iter]=ptmax;
205   }
206   if(fVarsForOpt[9]){
207     iter++;
208     vars[iter]=dd->CosPointingAngle();
209   }
210   if(fVarsForOpt[10]){
211     iter++;
212     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
213   }
214   if(fVarsForOpt[11]){
215     iter++;
216     Float_t maxDCA=0.;
217     for(Int_t i=0;i<3;i++){ 
218       if(d->GetDCA(i)>maxDCA) maxDCA=d->GetDCA(i);
219     }
220     vars[iter]=maxDCA;
221   }
222   if(fVarsForOpt[12]){
223     iter++;
224     Double_t mPDGPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
225     if(TMath::Abs(pdgdaughters[0])==321){
226       
227       Double_t phimass01=d->InvMass2Prongs(0,1,321,321);
228        vars[iter]=TMath::Abs(phimass01-mPDGPhi);
229        // vars[iter]=dd->InvMass2Prongs(0,1,321,321);
230     }else{
231       Double_t phimass12=d->InvMass2Prongs(1,2,321,321);
232        vars[iter]=TMath::Abs(phimass12-mPDGPhi);
233        // vars[iter]=dd->InvMass2Prongs(1,2,321,321);      
234     }
235   }
236   if(fVarsForOpt[13]){
237     iter++;
238     Double_t mPDGK0star = TDatabasePDG::Instance()->GetParticle(313)->Mass();
239     if(TMath::Abs(pdgdaughters[0])==321){
240       
241       Double_t mass12kpi=d->InvMass2Prongs(1,2,321,211);
242       vars[iter]=TMath::Abs(mass12kpi-mPDGK0star);
243       //              vars[iter]=dd->InvMass2Prongs(1,2,321,211);
244     }else{
245       Double_t mass01pik=d->InvMass2Prongs(0,1,211,321);
246       vars[iter]=TMath::Abs(mass01pik-mPDGK0star);
247       //        vars[iter]=dd->InvMass2Prongs(0,1,211,321);      
248     }
249   }
250
251   
252   return;
253 }
254 //---------------------------------------------------------------------------
255 Bool_t AliRDHFCutsDstoKKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
256 {
257   //
258   // Checking if Ds is in fiducial acceptance region 
259   //
260
261   if(pt > 5.) {
262     // applying cut for pt > 5 GeV
263     AliDebug(2,Form("pt of Ds = %f (> 5), cutting at |y| < 0.8",pt)); 
264     if (TMath::Abs(y) > 0.8) return kFALSE;
265     
266   } else {
267     // appliying smooth cut for pt < 5 GeV
268     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
269     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;         
270     AliDebug(2,Form("pt of Ds = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
271     if (y < minFiducialY || y > maxFiducialY) return kFALSE;    
272   }
273
274   return kTRUE;
275 }
276
277 //---------------------------------------------------------------------------
278 Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
279   // PID selection
280   // return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both 
281   Int_t retCode=3;
282   Bool_t okKKpi=kTRUE;
283   Bool_t okpiKK=kTRUE;
284   if(!fUsePID || !rd) return retCode;
285   if(!fPidHF){
286     AliWarning("AliAODPidHF not created!");
287     return retCode;
288   }
289   Int_t nKaons=0;
290   Int_t nNotKaons=0;
291   Int_t sign= rd->GetCharge(); 
292   for(Int_t iDaught=0; iDaught<3; iDaught++){
293     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(iDaught);
294     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
295     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
296     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
297     
298     if(isProton>0 &&  isKaon<0  && isPion<0) return 0;
299     if(sign!=track->Charge()){// must be kaon
300       if(isKaon<0) return 0;
301     }
302     if(isKaon>0 && isPion<0) nKaons++;
303     if(isKaon<0) nNotKaons++;
304     if(iDaught==0){
305       if(isKaon<0) okKKpi=kFALSE;
306       if(isPion<0) okpiKK=kFALSE;
307     }
308     else if(iDaught==2){
309       if(isKaon<0) okpiKK=kFALSE;
310       if(isPion<0) okKKpi=kFALSE;
311     }
312   }
313   
314   if(nKaons>2)return 0;
315   if(nNotKaons>1) return 0;
316   
317   if(!okKKpi) retCode-=1;
318   if(!okpiKK) retCode-=2;
319
320   return retCode;
321 }
322
323 //---------------------------------------------------------------------------
324 Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
325   //
326   // Apply selection
327   //
328
329   if(!fCutsRD){
330     cout<<"Cut matrix not inizialized. Exit..."<<endl;
331     return 0;
332   }
333   //PrintAll();
334   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
335
336   if(!d){
337     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
338     return 0;
339   }
340  
341   Double_t ptD=d->Pt();
342   if(ptD<fMinPtCand) return 0;
343   if(ptD>fMaxPtCand) return 0;
344
345   // selection on daughter tracks 
346   if(selectionLevel==AliRDHFCuts::kAll || 
347      selectionLevel==AliRDHFCuts::kTracks) {
348     if(!AreDaughtersSelected(d)) return 0;
349   }
350
351
352
353  
354   // selection on candidate
355   if(selectionLevel==AliRDHFCuts::kAll || 
356      selectionLevel==AliRDHFCuts::kCandidate) {
357     //recalculate vertex w/o daughters
358     AliAODVertex *origownvtx=0x0;
359     if(fRemoveDaughtersFromPrimary) {
360       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
361       if(!RecalcOwnPrimaryVtx(d,aod)) {
362         CleanOwnPrimaryVtx(d,aod,origownvtx);
363         return 0;
364       }
365     }
366
367     Int_t okDsKKpi=1;
368     Int_t okDspiKK=1;
369     Int_t okMassPhiKKpi=0;
370     Int_t okMassPhipiKK=0;
371     Int_t okMassK0starKKpi=0;
372     Int_t okMassK0starpiKK=0;
373     Int_t okDsPhiKKpi=0;
374     Int_t okDsPhipiKK=0;
375     Int_t okDsK0starKKpi=0;
376     Int_t okDsK0starpiKK=0;
377
378     Double_t pt=d->Pt();
379     Int_t ptbin=PtBin(pt);
380     if (ptbin==-1) {
381       CleanOwnPrimaryVtx(d,aod,origownvtx);
382       return 0;
383     }
384  
385     Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
386     Double_t mDsKKpi=d->InvMassDsKKpi();
387     Double_t mDspiKK=d->InvMassDspiKK();
388     if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
389     if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
390     if(!okDsKKpi && !okDspiKK){
391       CleanOwnPrimaryVtx(d,aod,origownvtx);
392       return 0;
393     }
394
395
396
397     // cuts on resonant decays (via Phi or K0*)
398     Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
399     Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
400     if(okDsKKpi){
401       Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
402       Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
403       if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhiKKpi=1;
404       if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starKKpi = 1;
405       if(!okMassPhiKKpi && !okMassK0starKKpi) okDsKKpi=0;
406       if(okMassPhiKKpi) okDsPhiKKpi=1;
407       if(okMassK0starKKpi) okDsK0starKKpi=1;
408     }
409     if(okDspiKK){
410       Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
411       Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
412       if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starpiKK = 1;
413       if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhipiKK=1;
414       if(!okMassPhipiKK && !okMassK0starpiKK) okDspiKK=0;
415       if(okMassPhipiKK) okDsPhipiKK=1;
416       if(okMassK0starpiKK) okDsK0starpiKK=1;
417     }
418     if(!okDsKKpi && !okDspiKK){
419       CleanOwnPrimaryVtx(d,aod,origownvtx);
420       return 0;
421     }
422
423     // Cuts on track pairs
424     for(Int_t i=0;i<3;i++){
425       if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]){
426         CleanOwnPrimaryVtx(d,aod,origownvtx);
427         return 0;
428       }
429     }
430     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] || 
431        d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]){
432       CleanOwnPrimaryVtx(d,aod,origownvtx);
433       return 0;
434     }
435
436
437
438     //single track
439     if(TMath::Abs(d->Pt2Prong(1)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] || 
440        TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]){
441       CleanOwnPrimaryVtx(d,aod,origownvtx);
442       return 0;
443     }
444
445     if(okDsKKpi){
446       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] || 
447          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
448       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] || 
449          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
450     }
451     if(okDspiKK){
452       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] || 
453          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
454       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] || 
455          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
456     }
457     if(!okDsKKpi && !okDspiKK){
458       CleanOwnPrimaryVtx(d,aod,origownvtx);
459       return 0;
460     }
461
462     // Cuts on candidate triplet
463
464     if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)]){
465       CleanOwnPrimaryVtx(d,aod,origownvtx);
466       return 0;
467     }
468
469     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]){
470       CleanOwnPrimaryVtx(d,aod,origownvtx); 
471       return 0;
472     }
473      
474     if(d->Pt2Prong(0)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] && 
475        d->Pt2Prong(1)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] && 
476        d->Pt2Prong(2)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)]) {
477       CleanOwnPrimaryVtx(d,aod,origownvtx); 
478       return 0;
479     }
480
481     if(d->DecayLength2()<fCutsRD[GetGlobalIndex(7,ptbin)]*fCutsRD[GetGlobalIndex(7,ptbin)]){
482       CleanOwnPrimaryVtx(d,aod,origownvtx);
483       return 0;
484     }
485
486
487     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
488     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)]){
489       CleanOwnPrimaryVtx(d,aod,origownvtx);
490       return 0;
491     }
492
493      // unset recalculated primary vertex when not needed any more
494     CleanOwnPrimaryVtx(d,aod,origownvtx);
495    
496
497     if(!okDsKKpi){
498       okDsPhiKKpi=0;
499       okDsK0starKKpi=0;
500     }
501     if(!okDspiKK){
502       okDsPhipiKK=0;
503       okDsK0starpiKK=0;
504     }
505
506     // PID selection
507     Int_t returnvaluePID=3;  
508     if(selectionLevel==AliRDHFCuts::kAll || 
509        selectionLevel==AliRDHFCuts::kCandidate ||     
510        selectionLevel==AliRDHFCuts::kPID) {
511       returnvaluePID = IsSelectedPID(d);
512       fIsSelectedPID=returnvaluePID;
513     }
514     if(returnvaluePID==0)return 0;
515
516     Bool_t okPidDsKKpi=returnvaluePID&1;
517     Bool_t okPidDspiKK=returnvaluePID&2;
518     if(!okPidDsKKpi){
519       okDsPhiKKpi=0;
520       okDsK0starKKpi=0;
521     }
522     if(!okPidDspiKK){
523       okDsPhipiKK=0;
524       okDsK0starpiKK=0;
525     }
526
527     if((okPidDsKKpi && okDsKKpi)||(okPidDspiKK && okDspiKK)){
528       Int_t returnvalue=0;
529       if(okDsKKpi) returnvalue+=1;
530       if(okDspiKK) returnvalue+=2;
531       if(okDsPhiKKpi) returnvalue+=4;
532       if(okDsPhipiKK) returnvalue+=8;
533       if(okDsK0starKKpi) returnvalue+=16;
534       if(okDsK0starpiKK) returnvalue+=32;
535       return returnvalue;
536     }else{
537       return 0;
538     }
539   }
540   return 15;
541
542 }
543 //---------------------------------------------------------------------------