1db74200e389744b0672fe7ec0062086bfa0c7b9
[u/mrichter/AliRoot.git] / ANALYSIS / 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 //
19 // Note: the two decay tracks are labelled: 0 (positive track)
20 //                                          1 (negative track)
21 //
22 //            Origin: A. Dainese    andrea.dainese@pd.infn.it            
23 //----------------------------------------------------------------------------
24 #include <Riostream.h>
25 #include <TH1.h>
26 #include <TH2.h>
27 #include <TCanvas.h>
28 #include <TPaveLabel.h>
29 #include <TVector3.h>
30
31 #include "AliD0toKpi.h"
32
33 ClassImp(AliD0toKpi)
34
35 //----------------------------------------------------------------------------
36 AliD0toKpi::AliD0toKpi() {
37   // Default constructor
38   
39   fSignal = kFALSE;
40
41   fEvent = 0;
42
43   fTrkNum[0] = 0;
44   fTrkNum[1] = 0;
45
46   fV1x = 0.;
47   fV1y = 0.;
48   fV1z = 0.;
49   fV2x = 0.;
50   fV2y = 0.;
51   fV2z = 0.;
52   fDCA = 0.;
53
54   fPx[0] = 0.;
55   fPy[0] = 0.;
56   fPz[0] = 0.;
57   fPx[1] = 0.;
58   fPy[1] = 0.;
59   fPz[1] = 0.;
60
61   fd0[0] = 0.;
62   fd0[1] = 0.;
63
64   fPdg[0] = 0;
65   fPdg[1] = 0;
66   fMum[0] = 0;
67   fMum[1] = 0;
68
69   fTagPi[0] = 0.;
70   fTagPi[1] = 0.;
71   fTagKa[0] = 0.;
72   fTagKa[1] = 0.;
73   fTagNid[0] = 0.;
74   fTagNid[1] = 0.;
75
76   fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=0;
77
78 }
79 //----------------------------------------------------------------------------
80 AliD0toKpi::AliD0toKpi(Int_t ev,Int_t trkNum[2],
81                        Double_t v1[3],Double_t v2[3], 
82                        Double_t dca,
83                        Double_t mom[6],Double_t d0[2]) {
84   // Constructor
85
86   fSignal = kFALSE;
87
88   fEvent = ev;
89   fTrkNum[0] = trkNum[0];
90   fTrkNum[1] = trkNum[1];
91
92   fV1x = v1[0];
93   fV1y = v1[1];
94   fV1z = v1[2];
95   fV2x = v2[0];
96   fV2y = v2[1];
97   fV2z = v2[2];
98   fDCA = dca;
99
100   fPx[0] = mom[0];
101   fPy[0] = mom[1];
102   fPz[0] = mom[2];
103   fPx[1] = mom[3];
104   fPy[1] = mom[4];
105   fPz[1] = mom[5];
106
107   fd0[0] = d0[0];
108   fd0[1] = d0[1];
109
110   fPdg[0] = 0;
111   fPdg[1] = 0;
112   fMum[0] = 0;
113   fMum[1] = 0;
114
115   fTagPi[0]  = 0.;
116   fTagPi[1]  = 0.;
117   fTagKa[0]  = 0.;
118   fTagKa[1]  = 0.;
119   fTagNid[0] = 0.;
120   fTagNid[1] = 0.;
121
122   fWgtAD0=fWgtAD0bar=fWgtBD0=fWgtBD0bar=fWgtCD0=fWgtCD0bar=fWgtDD0=fWgtDD0bar=0;
123 }
124 //----------------------------------------------------------------------------
125 AliD0toKpi::~AliD0toKpi() {}
126 //____________________________________________________________________________
127 AliD0toKpi::AliD0toKpi( const AliD0toKpi& d0toKpi):TObject(d0toKpi) {
128   // dummy copy constructor
129 }
130 //----------------------------------------------------------------------------
131 void AliD0toKpi::ApplyPID(TString pidScheme) {
132
133   const char *tofparampbpb = strstr(pidScheme.Data(),"TOFparam_PbPb");
134   const char *tofparampp   = strstr(pidScheme.Data(),"TOFparam_pp");
135
136   if((tofparampbpb || tofparampp) && fPdg[0]==0) {
137     printf("AliD0toKpi::ApplyPID :\n Warning: TOF parameterized PID can be used only for simulation!\n"); 
138     return;
139   }
140
141   if(tofparampbpb) {
142     // tagging of the positive track
143     if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13 
144        || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
145       fTagPi[0]  = LinearInterpolation(PChild(0),kPiBins_PbPb,kPiBinWidth_PbPb,kPiTagPi_PbPb);
146       fTagNid[0] = 1.-fTagPi[0];
147       fTagKa[0]   = 0.;
148       fTagPr[0]   = 0.;
149     } 
150     if(TMath::Abs(fPdg[0])==321) { // kaon
151       fTagKa[0]   = LinearInterpolation(PChild(0),kKBins_PbPb,kKBinWidth_PbPb,kKTagK_PbPb);
152       fTagNid[0] = LinearInterpolation(PChild(0),kKBins_PbPb,kKBinWidth_PbPb,kKTagNid_PbPb);
153       if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
154       fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
155       fTagPr[0] = 0.;
156     } 
157     if(TMath::Abs(fPdg[0])==2212) { // proton
158       fTagPr[0]  = LinearInterpolation(PChild(0),kPBins_PbPb,kPBinWidth_PbPb,kPTagP_PbPb);
159       fTagNid[0] = LinearInterpolation(PChild(0),kPBins_PbPb,kPBinWidth_PbPb,kPTagNid_PbPb);
160       if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
161       fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
162       fTagKa[0]   = 0.;
163     } 
164     // tagging of the negative track
165     if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13 
166        || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
167       fTagPi[1]  = LinearInterpolation(PChild(1),kPiBins_PbPb,kPiBinWidth_PbPb,kPiTagPi_PbPb);
168       fTagNid[1] = 1.-fTagPi[1];
169       fTagKa[1]   = 0.;
170       fTagPr[1]   = 0.;
171     } 
172     if(TMath::Abs(fPdg[1])==321) { // kaon
173       fTagKa[1]   = LinearInterpolation(PChild(1),kKBins_PbPb,kKBinWidth_PbPb,kKTagK_PbPb);
174       fTagNid[1] = LinearInterpolation(PChild(1),kKBins_PbPb,kKBinWidth_PbPb,kKTagNid_PbPb);
175       if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
176       fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
177       fTagPr[1] = 0.;
178     } 
179     if(TMath::Abs(fPdg[1])==2212) { // proton
180       fTagPr[1]  = LinearInterpolation(PChild(1),kPBins_PbPb,kPBinWidth_PbPb,kPTagP_PbPb);
181       fTagNid[1] = LinearInterpolation(PChild(1),kPBins_PbPb,kPBinWidth_PbPb,kPTagNid_PbPb);
182       if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
183       fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
184       fTagKa[1]   = 0.;
185     } 
186   }
187
188
189   if(tofparampp) {
190     // tagging of the positive track
191     if(TMath::Abs(fPdg[0])==211 || TMath::Abs(fPdg[0])==13 
192        || TMath::Abs(fPdg[0])==11) { // pion,muon,electron
193       fTagPi[0]  = LinearInterpolation(PChild(0),kPiBins_pp,kPiBinWidth_pp,kPiTagPi_pp);
194       fTagNid[0] = 1.-fTagPi[0];
195       fTagKa[0]   = 0.;
196       fTagPr[0]   = 0.;
197     } 
198     if(TMath::Abs(fPdg[0])==321) { // kaon
199       fTagKa[0]   = LinearInterpolation(PChild(0),kKBins_pp,kKBinWidth_pp,kKTagK_pp);
200       fTagNid[0] = LinearInterpolation(PChild(0),kKBins_pp,kKBinWidth_pp,kKTagNid_pp);
201       if((fTagNid[0]+fTagKa[0])>1.) fTagNid[0] = 1.-fTagKa[0];
202       fTagPi[0] = 1.-fTagNid[0]-fTagKa[0];
203       fTagPr[0] = 0.;
204     } 
205     if(TMath::Abs(fPdg[0])==2212) { // proton
206       fTagPr[0]  = LinearInterpolation(PChild(0),kPBins_pp,kPBinWidth_pp,kPTagP_pp);
207       fTagNid[0] = LinearInterpolation(PChild(0),kPBins_pp,kPBinWidth_pp,kPTagNid_pp);
208       if((fTagNid[0]+fTagPr[0])>1.) fTagNid[0] = 1.-fTagPr[0];
209       fTagPi[0] = 1.-fTagNid[0]-fTagPr[0];
210       fTagKa[0]   = 0.;
211     } 
212     // tagging of the negative track
213     if(TMath::Abs(fPdg[1])==211 || TMath::Abs(fPdg[1])==13 
214        || TMath::Abs(fPdg[1])==11) { // pion,muon,electron
215       fTagPi[1]  = LinearInterpolation(PChild(1),kPiBins_pp,kPiBinWidth_pp,kPiTagPi_pp);
216       fTagNid[1] = 1.-fTagPi[1];
217       fTagKa[1]   = 0.;
218       fTagPr[1]   = 0.;
219     } 
220     if(TMath::Abs(fPdg[1])==321) { // kaon
221       fTagKa[1]   = LinearInterpolation(PChild(1),kKBins_pp,kKBinWidth_pp,kKTagK_pp);
222       fTagNid[1] = LinearInterpolation(PChild(1),kKBins_pp,kKBinWidth_pp,kKTagNid_pp);
223       if((fTagNid[1]+fTagKa[1])>1.) fTagNid[1] = 1.-fTagKa[1];
224       fTagPi[1] = 1.-fTagNid[1]-fTagKa[1];
225       fTagPr[1] = 0.;
226     } 
227     if(TMath::Abs(fPdg[1])==2212) { // proton
228       fTagPr[1]  = LinearInterpolation(PChild(1),kPBins_pp,kPBinWidth_pp,kPTagP_pp);
229       fTagNid[1] = LinearInterpolation(PChild(1),kPBins_pp,kPBinWidth_pp,kPTagNid_pp);
230       if((fTagNid[1]+fTagPr[1])>1.) fTagNid[1] = 1.-fTagPr[1];
231       fTagPi[1] = 1.-fTagNid[1]-fTagPr[1];
232       fTagKa[1]   = 0.;
233     } 
234   }
235
236   return;
237 }
238 //----------------------------------------------------------------------------
239 Double_t AliD0toKpi::ChildrenRelAngle() const {
240   // relative angle between K and pi
241
242   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
243   TVector3 mom1(fPx[1],fPy[1],fPz[1]);
244
245   Double_t angle = mom0.Angle(mom1);
246
247   return angle; 
248 }
249 //----------------------------------------------------------------------------
250 void AliD0toKpi::ComputeWgts() {
251   // calculate the weights for PID
252
253
254   // assignement of the weights from PID
255   fWgtAD0    = fTagKa[1]*(fTagPi[0]+fTagNid[0]);
256   fWgtAD0bar = fTagKa[0]*(fTagPi[1]+fTagNid[1]);
257   fWgtBD0    = fTagPi[0]*fTagNid[1];
258   fWgtBD0bar = fTagPi[1]*fTagNid[0];
259   fWgtCD0    = fTagNid[0]*fTagNid[1];
260   fWgtCD0bar = fTagNid[0]*fTagNid[1];
261   fWgtDD0    = 1.-fWgtAD0-fWgtBD0-fWgtCD0;
262   fWgtDD0bar = 1.-fWgtAD0bar-fWgtBD0bar-fWgtCD0bar;
263
264   /*
265   cerr<<fWgtAD0<<"  "<<fWgtAD0bar<<endl;
266   cerr<<fWgtBD0<<"  "<<fWgtBD0bar<<endl;
267   cerr<<fWgtCD0<<"  "<<fWgtCD0bar<<endl;
268
269   if(fWgtAD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
270   if(fWgtAD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
271   if(fWgtBD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
272   if(fWgtBD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
273   if(fWgtCD0<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
274   if(fWgtCD0bar<0.) cerr<<"AliD0toKpi::ComputeWgts()  Negative weight!!!\n";
275   */
276
277   return;
278 }
279 //----------------------------------------------------------------------------
280 void AliD0toKpi::CorrectWgt4BR(Double_t factor) {
281   // correct weights of background from charm 
282
283   fWgtAD0    *= factor;
284   fWgtAD0bar *= factor;
285   fWgtBD0    *= factor;
286   fWgtBD0bar *= factor;
287   fWgtCD0    *= factor;
288   fWgtCD0bar *= factor;
289   fWgtDD0    *= factor;
290   fWgtDD0bar *= factor;
291
292   return;
293 }
294 //----------------------------------------------------------------------------
295 Double_t AliD0toKpi::CosPointing() const {
296   // cosine of pointing angle in space
297
298   TVector3 mom(Px(),Py(),Pz());
299   TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
300
301   Double_t pta = mom.Angle(flight);
302
303   return TMath::Cos(pta); 
304 }
305 //----------------------------------------------------------------------------
306 Double_t AliD0toKpi::CosPointingXY() const {
307   // cosine of pointing angle in transverse plane
308
309   TVector3 momXY(Px(),Py(),0.);
310   TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
311
312   Double_t ptaXY = momXY.Angle(flightXY);
313
314   return TMath::Cos(ptaXY); 
315 }
316 //----------------------------------------------------------------------------
317 void AliD0toKpi::CosThetaStar(Double_t &ctsD0,Double_t &ctsD0bar) const {
318   // cosine of decay angle in the D0 rest frame
319
320   Double_t pStar = TMath::Sqrt(TMath::Power(kMD0*kMD0-kMK*kMK-kMPi*kMPi,2.)-4.*kMK*kMK*kMPi*kMPi)/(2.*kMD0);
321
322   Double_t beta = P()/Energy();
323   Double_t gamma = Energy()/kMD0;
324
325   ctsD0 = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
326   // if(ctsD0 > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0<<"!\n"; }
327   // if(ctsD0 < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0<<"!\n"; }
328
329   ctsD0bar = (Ql(0)/gamma-beta*TMath::Sqrt(pStar*pStar+kMK*kMK))/pStar;
330   // if(ctsD0bar > 1.)  { cerr<<"AliD0toKpi::CosThetaStar: > 1 "<<ctsD0bar<<"!\n"; }
331   // if(ctsD0bar < -1.) { cerr<<"AliD0toKpi::CosThetaStar: < -1 "<<ctsD0bar<<"!\n";} 
332
333   return;
334 }
335 //----------------------------------------------------------------------------
336 Double_t AliD0toKpi::Eta() const {
337   // pseudorapidity of the D0
338
339   Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
340   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
341   return eta;
342 }
343 //----------------------------------------------------------------------------
344 Double_t AliD0toKpi::EtaChild(Int_t child) const {
345   // pseudorapidity of the decay tracks
346
347   Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
348   Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
349   return eta;
350 }
351 //----------------------------------------------------------------------------
352 void AliD0toKpi::GetWgts(Double_t &WgtD0,Double_t &WgtD0bar,TString sample) 
353   const {
354   // returns the weights for pid
355
356   const char *sampleA = strstr(sample.Data(),"A");
357   const char *sampleB = strstr(sample.Data(),"B");
358   const char *sampleC = strstr(sample.Data(),"C");
359   const char *sampleD = strstr(sample.Data(),"D");
360   const char *sampleABCD = strstr(sample.Data(),"ABCD");
361   const char *sampleABC = strstr(sample.Data(),"ABC");
362   const char *sampleBC = strstr(sample.Data(),"BC");
363
364   if(sampleA) { WgtD0 = fWgtAD0;  WgtD0bar = fWgtAD0bar; }
365   if(sampleB) { WgtD0 = fWgtBD0;  WgtD0bar = fWgtBD0bar; }
366   if(sampleC) { WgtD0 = fWgtCD0;  WgtD0bar = fWgtCD0bar; }
367   if(sampleD) { WgtD0 = fWgtDD0;  WgtD0bar = fWgtDD0bar; }
368   if(sampleABCD) { 
369     WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0+fWgtDD0; 
370     WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar+fWgtDD0bar; 
371   }
372   if(sampleABC) { 
373     WgtD0    = fWgtAD0+fWgtBD0+fWgtCD0; 
374     WgtD0bar = fWgtAD0bar+fWgtBD0bar+fWgtCD0bar; 
375   }
376   if(sampleBC) { 
377     WgtD0    = fWgtBD0+fWgtCD0; 
378     WgtD0bar = fWgtBD0bar+fWgtCD0bar; 
379   }
380
381
382   if(fSignal) {
383     if(fMum[0]==421)  WgtD0bar = 0.;
384     if(fMum[0]==-421) WgtD0 = 0.; 
385   }
386
387   return;
388 }
389 //----------------------------------------------------------------------------
390 void AliD0toKpi::InvMass(Double_t &mD0,Double_t &mD0bar) const {
391   // invariant mass as D0 and as D0bar
392
393   Double_t energy[2];
394
395   // D0 -> K- Pi+
396   energy[1] = TMath::Sqrt(kMK*kMK+PChild(1)*PChild(1));
397   energy[0] = TMath::Sqrt(kMPi*kMPi+PChild(0)*PChild(0));
398
399   mD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
400     
401
402   // D0bar -> K+ Pi-
403   energy[0] = TMath::Sqrt(kMK*kMK+PChild(0)*PChild(0));
404   energy[1] = TMath::Sqrt(kMPi*kMPi+PChild(1)*PChild(1));
405
406   mD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
407     
408   return;
409
410 }
411 //----------------------------------------------------------------------------
412 Double_t AliD0toKpi::Ql(Int_t child) const {
413   // longitudinal momentum of decay tracks w.r.t. to D0 momentum
414
415   Double_t qL;
416   TVector3 mom(fPx[child],fPy[child],fPz[child]);
417   TVector3 momD(Px(),Py(),Pz());
418
419   qL = mom.Dot(momD)/momD.Mag();
420
421   return qL ;
422 }
423 //----------------------------------------------------------------------------
424 Double_t AliD0toKpi::Qt() const {
425   // transverse momentum of decay tracks w.r.t. to D0 momentum  
426
427   TVector3 mom0(fPx[0],fPy[0],fPz[0]);
428   TVector3 momD(Px(),Py(),Pz());
429
430   return mom0.Perp(momD);
431 }
432 //----------------------------------------------------------------------------
433 Bool_t AliD0toKpi::Select(const Double_t* cuts,Int_t& okD0,Int_t& okD0bar) 
434   const {
435 //
436 // This function compares the D0 with a set of cuts:
437 //
438 // cuts[0] = inv. mass half width [GeV]   
439 // cuts[1] = dca [micron]
440 // cuts[2] = cosThetaStar 
441 // cuts[3] = pTK [GeV/c]
442 // cuts[4] = pTPi [GeV/c]
443 // cuts[5] = d0K [micron]   upper limit!
444 // cuts[6] = d0Pi [micron]  upper limit!
445 // cuts[7] = d0d0 [micron^2]
446 // cuts[8] = cosThetaPoint
447 //
448 // If the D0/D0bar doesn't pass the cuts it sets the weights to 0
449 // If neither D0 nor D0bar pass the cuts return kFALSE
450 //
451   Double_t mD0,mD0bar,ctsD0,ctsD0bar;
452   okD0=1; okD0bar=1;
453
454   if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okD0 = 0;
455   if(PtChild(0) < cuts[3] || PtChild(1) < cuts[4]) okD0bar = 0;
456   if(!okD0 && !okD0bar) return kFALSE;
457
458   if(TMath::Abs(Getd0Child(1)) > cuts[5] || 
459      TMath::Abs(Getd0Child(0)) > cuts[6]) okD0 = 0;
460   if(TMath::Abs(Getd0Child(0)) > cuts[6] ||
461      TMath::Abs(Getd0Child(1)) > cuts[5]) okD0bar = 0;
462   if(!okD0 && !okD0bar) return kFALSE;
463
464   if(GetDCA() > cuts[1]) { okD0 = okD0bar = 0; return kFALSE; }
465
466   InvMass(mD0,mD0bar);
467   if(TMath::Abs(mD0-kMD0)    > cuts[0]) okD0 = 0;
468   if(TMath::Abs(mD0bar-kMD0) > cuts[0]) okD0bar = 0;
469   if(!okD0 && !okD0bar) return kFALSE;
470
471   CosThetaStar(ctsD0,ctsD0bar);
472   if(TMath::Abs(ctsD0)    > cuts[2]) okD0 = 0;
473   if(TMath::Abs(ctsD0bar) > cuts[2]) okD0bar = 0;
474   if(!okD0 && !okD0bar) return kFALSE;
475
476   if(ProdImpParams() > cuts[7]) { okD0 = okD0bar = 0; return kFALSE; }
477
478   if(CosPointing()   < cuts[8]) { okD0 = okD0bar = 0; return kFALSE; }
479
480   return kTRUE;
481 }
482 //-----------------------------------------------------------------------------
483 void AliD0toKpi::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
484   // Set combined PID detector response probabilities
485
486   fPIDrespEl[0] = resp0[0];
487   fPIDrespEl[1] = resp1[0];
488   fPIDrespMu[0] = resp0[1];
489   fPIDrespMu[1] = resp1[1];
490   fPIDrespPi[0] = resp0[2];
491   fPIDrespPi[1] = resp1[2];
492   fPIDrespKa[0] = resp0[3];
493   fPIDrespKa[1] = resp1[3];
494   fPIDrespPr[0] = resp0[4];
495   fPIDrespPr[1] = resp1[4];
496
497   return;
498
499 //-----------------------------------------------------------------------------
500 void AliD0toKpi::DrawPIDinTOF(TString pidScheme) const {
501   // Draw parameterized PID probabilities in TOF
502
503   const char *tofparampbpb = strstr(pidScheme.Data(),"TOFparam_PbPb");
504   const char *tofparampp = strstr(pidScheme.Data(),"TOFparam_pp");
505
506   TH2F* framePi = new TH2F("framePi","Tag probabilities for PIONS",2,0,2.5,2,0,1);
507   framePi->SetXTitle("p [GeV/c]"); 
508   framePi->SetStats(0);
509   TH2F* frameK = new TH2F("frameK","Tag probabilities for KAONS",2,0,2.5,2,0,1);
510   frameK->SetXTitle("p [GeV/c]");
511   frameK->SetStats(0);
512   TH2F* frameP = new TH2F("frameP","Tag probabilities for PROTONS",2,0,4.5,2,0,1);
513   frameP->SetXTitle("p [GeV/c]");
514   frameP->SetStats(0);
515
516   TH1F* hPiPi = new TH1F("hPiPi","Tag probabilities for PIONS",kPiBins_PbPb,0,2.5);
517   TH1F* hPiNid = new TH1F("hPiNid","Tag probabilities for PIONS",kPiBins_PbPb,0,2.5);
518
519   TH1F* hKK = new TH1F("hKK","Tag probabilities for KAONS",kKBins_PbPb,0,2.5);
520   TH1F* hKNid = new TH1F("hKNid","Tag probabilities for KAONS",kKBins_PbPb,0,2.5);
521   TH1F* hKPi = new TH1F("hKPi","Tag probabilities for KAONS",kKBins_PbPb,0,2.5);
522
523   TH1F* hPP = new TH1F("hPP","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5);
524   TH1F* hPNid = new TH1F("hPNid","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5);
525   TH1F* hPPi = new TH1F("hPPi","Tag probabilities for PROTONS",kPBins_PbPb,0,4.5);
526
527
528   if(tofparampbpb) {
529
530     for(Int_t i=1; i<=kPiBins_PbPb; i++) {
531       hPiPi->SetBinContent(i,kPiTagPi_PbPb[i-1]);
532       hPiNid->SetBinContent(i,kPiTagPi_PbPb[i-1]+kPiTagNid_PbPb[i-1]);
533       
534       hKK->SetBinContent(i,kKTagK_PbPb[i-1]);
535       hKPi->SetBinContent(i,kKTagK_PbPb[i-1]+kKTagPi_PbPb[i-1]);
536       hKNid->SetBinContent(i,kKTagK_PbPb[i-1]+kKTagPi_PbPb[i-1]+kKTagNid_PbPb[i-1]);
537     }
538     for(Int_t i=1; i<=kPBins_PbPb; i++) {    
539       hPP->SetBinContent(i,kPTagP_PbPb[i-1]);
540       hPPi->SetBinContent(i,kPTagP_PbPb[i-1]+kPTagPi_PbPb[i-1]);
541       hPNid->SetBinContent(i,kPTagP_PbPb[i-1]+kPTagPi_PbPb[i-1]+kPTagNid_PbPb[i-1]);
542     }
543
544   } else if(tofparampp) {
545
546     for(Int_t i=1; i<=kPiBins_pp; i++) {
547       hPiPi->SetBinContent(i,kPiTagPi_pp[i-1]);
548       hPiNid->SetBinContent(i,kPiTagPi_pp[i-1]+kPiTagNid_pp[i-1]);
549       
550       hKK->SetBinContent(i,kKTagK_pp[i-1]);
551       hKPi->SetBinContent(i,kKTagK_pp[i-1]+kKTagPi_pp[i-1]);
552       hKNid->SetBinContent(i,kKTagK_pp[i-1]+kKTagPi_pp[i-1]+kKTagNid_pp[i-1]);
553     }
554     for(Int_t i=1; i<=kPBins_pp; i++) {    
555       hPP->SetBinContent(i,kPTagP_pp[i-1]);
556       hPPi->SetBinContent(i,kPTagP_pp[i-1]+kPTagPi_pp[i-1]);
557       hPNid->SetBinContent(i,kPTagP_pp[i-1]+kPTagPi_pp[i-1]+kPTagNid_pp[i-1]);
558     }
559
560   } 
561
562
563   TCanvas* c = new TCanvas("c","Parameterized PID in TOF",0,0,1000,400);
564   c->Divide(3,1);
565   c->cd(1);
566   framePi->Draw();
567   hPiNid->SetFillColor(18); hPiNid->Draw("same");
568   hPiPi->SetFillColor(4); hPiPi->Draw("same");
569   TPaveLabel* pav1 = new TPaveLabel(1,.2,1.4,.3,"#pi");
570   pav1->SetBorderSize(0);
571   pav1->Draw("same");
572   TPaveLabel* pav2 = new TPaveLabel(1,.8,1.8,.9,"non-id");
573   pav2->SetBorderSize(0);
574   pav2->Draw("same");
575
576   c->cd(2);
577   frameK->Draw();
578   hKNid->SetFillColor(18); hKNid->Draw("same");
579   hKPi->SetFillColor(4); hKPi->Draw("same");
580   hKK->SetFillColor(7); hKK->Draw("same");
581   TPaveLabel* pav3 = new TPaveLabel(1,.2,1.5,.3,"K");
582   pav3->SetBorderSize(0);
583   pav3->Draw("same");
584   TPaveLabel* pav4 = new TPaveLabel(1,.8,1.8,.9,"non-id");
585   pav4->SetBorderSize(0);
586   pav4->Draw("same");
587   TPaveLabel* pav5 = new TPaveLabel(.4,.5,.8,.6,"#pi");
588   pav5->SetBorderSize(0);
589   pav5->Draw("same");
590
591   c->cd(3);
592   frameP->Draw();
593   hPNid->SetFillColor(18); hPNid->Draw("same");
594   hPPi->SetFillColor(4); hPPi->Draw("same");
595   hPP->SetFillColor(3); hPP->Draw("same");
596   TPaveLabel* pav6 = new TPaveLabel(1,.2,1.5,.3,"p");
597   pav6->SetBorderSize(0);
598   pav6->Draw("same");
599   TPaveLabel* pav7 = new TPaveLabel(1,.8,2.6,.9,"non-id");
600   pav7->SetBorderSize(0);
601   pav7->Draw("same");
602   TPaveLabel* pav8 = new TPaveLabel(.2,.5,1,.6,"#pi");
603   pav8->SetBorderSize(0);
604   pav8->Draw("same");
605
606
607   return;
608 }
609 //----------------------------------------------------------------------------
610 Double_t AliD0toKpi::LinearInterpolation(Double_t p,Int_t nBins,Double_t Bin,
611                                          const Double_t *values) const {
612   // a linear interpolation method
613
614   Double_t value=0; 
615   Double_t slope;
616
617   if(p<0.5*Bin) {
618     value = values[0];
619   } else if(p>=(nBins-0.5)*Bin) {
620     slope = (2*values[nBins-1]-values[nBins-2]-values[nBins-3])/Bin/2;
621     value = values[nBins-2]+slope*(p-Bin*(nBins-1.5));
622   } else {
623     for(Int_t i=0; i<nBins; i++) {
624       if(p<(i+0.5)*Bin) {
625         slope = (values[i]-values[i-1])/Bin;
626         value = values[i-1]+slope*(p-Bin*(i-0.5));
627         break;
628       }
629     }
630   }
631
632   if(value<0.) value=0.;
633   if(value>1.) value=1.;
634
635   return value;
636 }
637 //----------------------------------------------------------------------------
638
639
640
641
642
643
644 /*
645 //____________________________________________________________________________
646 void AliD0toKpi::SetPtWgts4pp() {
647   // Correct pt distribution in order to reproduce MNR pt slope
648   // (for pp generated with PYTHIA min. bias)
649
650   if(TMath::Abs(fMum[0]) != 421 && TMath::Abs(fMum[1]) != 421 &&
651      TMath::Abs(fMum[0]) != 411 && TMath::Abs(fMum[1]) != 411) return;
652
653   Double_t ptWgt = 1.;
654   ptWgt = 2.05-0.47*Pt()+0.02*Pt()*Pt();
655   if(Pt() >= 5.) ptWgt = 0.56*TMath::Exp(-0.12*Pt());
656
657   fWgtAD0    *= ptWgt;
658   fWgtAD0bar *= ptWgt;
659   fWgtBD0    *= ptWgt;
660   fWgtBD0bar *= ptWgt;
661   fWgtCD0    *= ptWgt;
662   fWgtCD0bar *= ptWgt;
663   fWgtDD0    *= ptWgt;
664   fWgtDD0bar *= ptWgt;
665
666   return;
667 }
668 //____________________________________________________________________________
669 */
670
671
672
673
674
675
676
677
678