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