Removing PWG3TrackExtrapInterface from PWG3base library (Gines)
[u/mrichter/AliRoot.git] / PWG3 / AliD0toKpi.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 //               Implementation of the D0toKpi class
18 //                  for pp and PbPb interactions
19 // Note: the two decay tracks are labelled: 0 (positive track)
20 //                                          1 (negative track)
21 //            Origin: A. Dainese    andrea.dainese@lnl.infn.it            
22 //----------------------------------------------------------------------------
23
24 #include <TH1.h>
25 #include <TH2.h>
26 #include <TCanvas.h>
27 #include <TPaveLabel.h>
28 #include <TVector3.h>
29
30 #include "AliD0toKpi.h"
31
32 ClassImp(AliD0toKpi)
33
34 //----------------------------------------------------------------------------
35 AliD0toKpi::AliD0toKpi():
36 fSignal(kFALSE),
37 fEvent(0),
38 fV1x(0.),
39 fV1y(0.),
40 fV1z(0.),
41 fV2x(0.),
42 fV2y(0.),
43 fV2z(0.),
44 fDCA(0.),
45 fWgtAD0(0.),
46 fWgtAD0bar(0.),
47 fWgtBD0(0.),
48 fWgtBD0bar(0.),
49 fWgtCD0(0.),
50 fWgtCD0bar(0.),
51 fWgtDD0(0.),
52 fWgtDD0bar(0.)
53 {
54   // Default constructor
55
56   fTrkNum[0] = 0;
57   fTrkNum[1] = 0;
58
59   fPx[0] = 0.;
60   fPy[0] = 0.;
61   fPz[0] = 0.;
62   fPx[1] = 0.;
63   fPy[1] = 0.;
64   fPz[1] = 0.;
65
66   fd0[0] = 0.;
67   fd0[1] = 0.;
68
69   fPdg[0] = 0;
70   fPdg[1] = 0;
71   fMum[0] = 0;
72   fMum[1] = 0;
73
74   fTagPi[0] = 0.;
75   fTagPi[1] = 0.;
76   fTagKa[0] = 0.;
77   fTagKa[1] = 0.;
78   fTagNid[0] = 0.;
79   fTagNid[1] = 0.;
80
81 }
82 //----------------------------------------------------------------------------
83 AliD0toKpi::AliD0toKpi(Int_t ev,Int_t trkNum[2],
84                        Double_t v1[3],Double_t v2[3], 
85                        Double_t dca,
86                        Double_t mom[6],Double_t d0[2]): 
87 fSignal(kFALSE),
88 fEvent(ev),
89 fV1x(v1[0]),
90 fV1y(v1[1]),
91 fV1z(v1[2]),
92 fV2x(v2[0]),
93 fV2y(v2[1]),
94 fV2z(v2[2]),
95 fDCA(dca),
96 fWgtAD0(0.),
97 fWgtAD0bar(0.),
98 fWgtBD0(0.),
99 fWgtBD0bar(0.),
100 fWgtCD0(0.),
101 fWgtCD0bar(0.),
102 fWgtDD0(0.),
103 fWgtDD0bar(0.)
104 {
105   // Constructor
106
107   fTrkNum[0] = trkNum[0];
108   fTrkNum[1] = trkNum[1];
109
110   fPx[0] = mom[0];
111   fPy[0] = mom[1];
112   fPz[0] = mom[2];
113   fPx[1] = mom[3];
114   fPy[1] = mom[4];
115   fPz[1] = mom[5];
116
117   fd0[0] = d0[0];
118   fd0[1] = d0[1];
119
120   fPdg[0] = 0;
121   fPdg[1] = 0;
122   fMum[0] = 0;
123   fMum[1] = 0;
124
125   fTagPi[0]  = 0.;
126   fTagPi[1]  = 0.;
127   fTagKa[0]  = 0.;
128   fTagKa[1]  = 0.;
129   fTagNid[0] = 0.;
130   fTagNid[1] = 0.;
131 }
132 //----------------------------------------------------------------------------
133 AliD0toKpi::~AliD0toKpi() {}
134 //----------------------------------------------------------------------------
135 void AliD0toKpi::ApplyPID(TString pidScheme) {
136   // Applies particle identification
137
138   if((!pidScheme.CompareTo("TOFparamPbPb") || !pidScheme.CompareTo("TOFparamPP")) && fPdg[0]==0) {
139     printf("AliD0toKpi::ApplyPID :\n Warning: TOF parameterized PID can be used only for simulation!\n"); 
140     return;
141   }
142
143   if(!pidScheme.CompareTo("TOFparamPbPb")) {
144     // tagging of the positive track
145     if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13 
146        || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
147       fTagPi[0]  = LinearInterpolation(PChild(0),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
148       fTagNid[0] = 1.-fTagPi[0];
149       fTagKa[0]   = 0.;
150       fTagPr[0]   = 0.;
151     } 
152     if(TMath::Abs(fPdg[0])==321) { // kaon
153       fTagKa[0]   = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
154       fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
155       if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
156       fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
157       fTagPr[0] = 0.;
158     } 
159     if(TMath::Abs(fPdg[0])==2212) { // proton
160       fTagPr[0]  = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
161       fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
162       if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
163       fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
164       fTagKa[0]   = 0.;
165     } 
166     // tagging of the negative track
167     if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13 
168        || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
169       fTagPi[1]  = LinearInterpolation(PChild(1),kPiBinsPbPb,kPiBinWidthPbPb,kPiTagPiPbPb);
170       fTagNid[1] = 1.-fTagPi[1];
171       fTagKa[1]   = 0.;
172       fTagPr[1]   = 0.;
173     } 
174     if(TMath::Abs(fPdg[1])==321) { // kaon
175       fTagKa[1]   = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagKPbPb);
176       fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPbPb,kKBinWidthPbPb,kKTagNidPbPb);
177       if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
178       fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
179       fTagPr[1] = 0.;
180     } 
181     if(TMath::Abs(fPdg[1])==2212) { // proton
182       fTagPr[1]  = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagPPbPb);
183       fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPbPb,kPBinWidthPbPb,kPTagNidPbPb);
184       if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
185       fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
186       fTagKa[1]   = 0.;
187     } 
188   }
189
190
191   if(!pidScheme.CompareTo("TOFparamPP")) {
192     // tagging of the positive track
193     if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13 
194        || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
195       fTagPi[0]  = LinearInterpolation(PChild(0),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
196       fTagNid[0] = 1.-fTagPi[0];
197       fTagKa[0]   = 0.;
198       fTagPr[0]   = 0.;
199     } 
200     if(TMath::Abs(fPdg[0])==321) { // kaon
201       fTagKa[0]   = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagKPP);
202       fTagNid[0] = LinearInterpolation(PChild(0),kKBinsPP,kKBinWidthPP,kKTagNidPP);
203       if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
204       fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
205       fTagPr[0] = 0.;
206     } 
207     if(TMath::Abs(fPdg[0])==2212) { // proton
208       fTagPr[0]  = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagPPP);
209       fTagNid[0] = LinearInterpolation(PChild(0),kPBinsPP,kPBinWidthPP,kPTagNidPP);
210       if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
211       fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
212       fTagKa[0]   = 0.;
213     } 
214     // tagging of the negative track
215     if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13 
216        || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
217       fTagPi[1]  = LinearInterpolation(PChild(1),kPiBinsPP,kPiBinWidthPP,kPiTagPiPP);
218       fTagNid[1] = 1.-fTagPi[1];
219       fTagKa[1]   = 0.;
220       fTagPr[1]   = 0.;
221     } 
222     if(TMath::Abs(fPdg[1])==321) { // kaon
223       fTagKa[1]   = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagKPP);
224       fTagNid[1] = LinearInterpolation(PChild(1),kKBinsPP,kKBinWidthPP,kKTagNidPP);
225       if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
226       fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
227       fTagPr[1] = 0.;
228     } 
229     if(TMath::Abs(fPdg[1])==2212) { // proton
230       fTagPr[1]  = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagPPP);
231       fTagNid[1] = LinearInterpolation(PChild(1),kPBinsPP,kPBinWidthPP,kPTagNidPP);
232       if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
233       fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
234       fTagKa[1]   = 0.;
235     } 
236   }
237
238   return;
239 }
240 //----------------------------------------------------------------------------
241 Double_t AliD0toKpi::ChildrenRelAngle() const {
242   // relative angle between K and pi
243
244   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
245   TVector3 mom1(fPx[1],fPy[1],fPz[1]);
246
247   Double_t angle = mom0.Angle(mom1);
248
249   return angle; 
250 }
251 //----------------------------------------------------------------------------
252 void AliD0toKpi::ComputeWgts() {
253   // calculate the weights for PID
254
255
256   // assignement of the weights from PID
257   fWgtAD0    = fTagKa[1]*(fTagPi[0]+fTagNid[0]);
258   fWgtAD0bar = fTagKa[0]*(fTagPi[1]+fTagNid[1]);
259   fWgtBD0    = fTagPi[0]*fTagNid[1];
260   fWgtBD0bar = fTagPi[1]*fTagNid[0];
261   fWgtCD0    = fTagNid[0]*fTagNid[1];
262   fWgtCD0bar = fTagNid[0]*fTagNid[1];
263   fWgtDD0    = 1.-fWgtAD0-fWgtBD0-fWgtCD0;
264   fWgtDD0bar = 1.-fWgtAD0bar-fWgtBD0bar-fWgtCD0bar;
265
266   /*  
267   cerr<<fWgtAD0<<"  "<<fWgtAD0bar<<endl;
268   cerr<<fWgtBD0<<"  "<<fWgtBD0bar<<endl;
269   cerr<<fWgtCD0<<"  "<<fWgtCD0bar<<endl;
270   cerr<<fWgtDD0<<"  "<<fWgtDD0bar<<endl;
271   */
272   /*
273   if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
274   if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
275   if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
276   if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
277   if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
278   if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
279   */
280
281   return;
282 }
283 //----------------------------------------------------------------------------
284 void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
285   // correct weights of background from charm 
286
287   fWgtAD0    *= factor;
288   fWgtAD0bar *= factor;
289   fWgtBD0    *= factor;
290   fWgtBD0bar *= factor;
291   fWgtCD0    *= factor;
292   fWgtCD0bar *= factor;
293   fWgtDD0    *= factor;
294   fWgtDD0bar *= factor;
295
296   return;
297 }
298 //----------------------------------------------------------------------------
299 Double_t AliD0toKpi::CosPointing() const {
300   // cosine of pointing angle in space
301
302   TVector3 mom(Px(),Py(),Pz());
303   TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
304
305   Double_t pta = mom.Angle(flight);
306
307   return TMath::Cos(pta); 
308 }
309 //----------------------------------------------------------------------------
310 Double_t AliD0toKpi::CosPointingXY() const {
311   // cosine of pointing angle in transverse plane
312
313   TVector3 momXY(Px(),Py(),0.);
314   TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
315
316   Double_t ptaXY = momXY.Angle(flightXY);
317
318   return TMath::Cos(ptaXY); 
319 }
320 //----------------------------------------------------------------------------
321 void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
322   // cosine of decay angle in the D0 rest frame
323
324   Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
325
326   Double_t beta = P()/Energy();
327   Double_t gamma = Energy()/kMD0;
328
329   ctsD0 = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
330   // if(ctsD0 > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
331   // if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }
332
333   ctsD0bar = (Ql(0)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
334   // if(ctsD0bar > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
335   // if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";} 
336
337   return;
338 }
339 //----------------------------------------------------------------------------
340 Double_t AliD0toKpi::Eta() const {
341   // pseudorapidity of the D0
342
343   Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
344   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
345   return eta;
346 }
347 //----------------------------------------------------------------------------
348 Double_t AliD0toKpi::EtaChild(Int_t child) const {
349   // pseudorapidity of the decay tracks
350
351   Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
352   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
353   return eta;
354 }
355 //----------------------------------------------------------------------------
356 void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample) 
357   const {
358   // returns the weights for pid
359
360   if(!sample.CompareTo("A")) { WgtD0 = fWgtAD0;  WgtD0bar = fWgtAD0bar; }
361   if(!sample.CompareTo("B")) { WgtD0 = fWgtBD0;  WgtD0bar = fWgtBD0bar; }
362   if(!sample.CompareTo("C")) { WgtD0 = fWgtCD0;  WgtD0bar = fWgtCD0bar; }
363   if(!sample.CompareTo("D")) { WgtD0 = fWgtDD0;  WgtD0bar = fWgtDD0bar; }
364   if(!sample.CompareTo("ABCD")) { 
365     WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0+fWgtDD0; 
366     WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar+fWgtDD0bar; 
367   }
368   if(!sample.CompareTo("ABC")) { 
369     WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0; 
370     WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar; 
371   }
372   if(!sample.CompareTo("BC")) { 
373     WgtD0    = fWgtBD0+fWgtCD0; 
374     WgtD0bar = fWgtBD0bar+fWgtCD0bar; 
375   }
376
377   return;
378 }
379 //----------------------------------------------------------------------------
380 Double_t AliD0toKpi::ImpPar() const {
381   // D0 impact parameter in the bending plane
382   
383     Double_t k = -(fV2x-fV1x)*Px()-(fV2y-fV1y)*Py();
384     k /= Pt()*Pt();
385     Double_t dx = fV2x-fV1x+k*Px();
386     Double_t dy = fV2y-fV1y+k*Py();
387     Double_t absDD = TMath::Sqrt(dx*dx+dy*dy);
388     TVector3 mom(Px(),Py(),Pz());
389     TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
390     TVector3 cross = mom.Cross(flight);
391     return (cross.Z()>0. ? absDD : -absDD);
392 }
393 //----------------------------------------------------------------------------
394 void AliD0toKpi::InvMass(Double_t &mD0,Double_t &mD0bar) const {
395   // invariant mass as D0 and as D0bar
396
397   Double_t energy[2];
398
399   // D0 -> K- Pi+
400   energy[1] = TMath::Sqrt(kMK*kMK+PChild(1)*PChild(1));
401   energy[0] = TMath::Sqrt(kMPi*kMPi+PChild(0)*PChild(0));
402
403   mD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
404     
405
406   // D0bar -> K+ Pi-
407   energy[0] = TMath::Sqrt(kMK*kMK+PChild(0)*PChild(0));
408   energy[1] = TMath::Sqrt(kMPi*kMPi+PChild(1)*PChild(1));
409
410   mD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
411     
412   return;
413
414 }
415 //----------------------------------------------------------------------------
416 Double_t AliD0toKpi::Ql(Int_t child) const {
417   // longitudinal momentum of decay tracks w.r.t. to D0 momentum
418
419   Double_t qL;
420   TVector3 mom(fPx[child],fPy[child],fPz[child]);
421   TVector3 momD(Px(),Py(),Pz());
422
423   qL = mom.Dot(momD)/momD.Mag();
424
425   return qL ;
426 }
427 //----------------------------------------------------------------------------
428 Double_t AliD0toKpi::Qt() const {
429   // transverse momentum of decay tracks w.r.t. to D0 momentum  
430
431   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
432   TVector3 momD(Px(),Py(),Pz());
433
434   return mom0.Perp(momD);
435 }
436 //----------------------------------------------------------------------------
437 Bool_t AliD0toKpi::Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar) 
438   const {
439 //
440 // This function compares the D0 with a set of cuts:
441 //
442 // cuts[0] = inv. mass half width [GeV]   
443 // cuts[1] = dca [micron]
444 // cuts[2] = cosThetaStar 
445 // cuts[3] = pTK [GeV/c]
446 // cuts[4] = pTPi [GeV/c]
447 // cuts[5] = d0K [micron]   upper limit!
448 // cuts[6] = d0Pi [micron]  upper limit!
449 // cuts[7] = d0d0 [micron^2]
450 // cuts[8] = cosThetaPoint
451 //
452 // If the D0/D0bar doesn't pass the cuts it sets the weights to 0
453 // If neither D0 nor D0bar pass the cuts return kFALSE
454 //
455   Double_t mD0,mD0bar,ctsD0,ctsD0bar;
456   okD0=1; okD0bar=1;
457
458   if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okD0 = 0;
459   if(PtChild(0) < cuts[3] || PtChild(1) < cuts[4]) okD0bar = 0;
460   if(!okD0 && !okD0bar) return kFALSE;
461
462   if(TMath::Abs(Getd0Child(1)) > cuts[5] || 
463      TMath::Abs(Getd0Child(0)) > cuts[6]) okD0 = 0;
464   if(TMath::Abs(Getd0Child(0)) > cuts[6] ||
465      TMath::Abs(Getd0Child(1)) > cuts[5]) okD0bar = 0;
466   if(!okD0 && !okD0bar) return kFALSE;
467
468   if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }
469
470   InvMass(mD0,mD0bar);
471   if(TMath::Abs(mD0-kMD0)    > cuts[0]) okD0 = 0;
472   if(TMath::Abs(mD0bar-kMD0) > cuts[0]) okD0bar = 0;
473   if(!okD0 && !okD0bar) return kFALSE;
474
475   CosThetaStar(ctsD0,ctsD0bar);
476   if(TMath::Abs(ctsD0)    > cuts[2]) okD0 = 0;
477   if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
478   if(!okD0 && !okD0bar) return kFALSE;
479
480   if(ProdImpParams() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }
481
482   if(CosPointing()   < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }
483
484   return kTRUE;
485 }
486 //-----------------------------------------------------------------------------
487 void AliD0toKpi::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
488   // Set combined PID detector response probabilities
489
490   fPIDrespEl[0] = resp0[0];
491   fPIDrespEl[1] = resp1[0];
492   fPIDrespMu[0] = resp0[1];
493   fPIDrespMu[1] = resp1[1];
494   fPIDrespPi[0] = resp0[2];
495   fPIDrespPi[1] = resp1[2];
496   fPIDrespKa[0] = resp0[3];
497   fPIDrespKa[1] = resp1[3];
498   fPIDrespPr[0] = resp0[4];
499   fPIDrespPr[1] = resp1[4];
500
501   return;
502
503 //----------------------------------------------------------------------------
504 Double_t AliD0toKpi::LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,
505                                          const Double_t *values) const {
506   // a linear interpolation method
507
508   Double_t value=0; 
509   Double_t slope;
510
511   if(p<0.5*Bin) {
512     value = values[0];
513   } else if(p>=(nBins-0.5)*Bin) {
514     slope = (2*values[nBins-1]-values[nBins-2]-values[nBins-3])/Bin/2;
515     value = values[nBins-2]+slope*(p-Bin*(nBins-1.5));
516   } else {
517     for(Int_t i=0; i<nBins; i++) {
518       if(p<(i+0.5)*Bin) {
519         slope = (values[i]-values[i-1])/Bin;
520         value = values[i-1]+slope*(p-Bin*(i-0.5));
521         break;
522       }
523     }
524   }
525
526   if(value<0.) value=0.;
527   if(value>1.) value=1.;
528
529   return value;
530 }
531 //----------------------------------------------------------------------------
532
533
534
535
536
537
538
539
540