96c0ca4f9c5a16a98cc3448c6c428357948e99a5
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibV0.cxx
1
2 #include <TROOT.h>
3 #include <TChain.h>
4 #include <TFile.h>
5 #include <TH3F.h>
6 #include <TH2F.h>
7 //
8 #include "TParticle.h"
9 #include "TDatabasePDG.h"
10 #include "AliRunLoader.h"
11 #include "AliStack.h"
12
13
14
15 #include <TPDGCode.h>
16 #include <TStyle.h>
17 #include "TLinearFitter.h"
18 #include "TMatrixD.h"
19 #include "TTreeStream.h"
20 #include "TF1.h"
21
22
23
24 #include "AliMagF.h"
25 #include "AliTracker.h"
26 #include "AliESDEvent.h"
27 #include "AliESDtrack.h"
28 #include "AliESDfriend.h"
29 #include "AliESDfriendTrack.h" 
30 #include "AliMathBase.h" 
31 #include "AliTPCseed.h"
32 #include "AliTPCclusterMI.h"
33
34 #include "AliKFParticle.h"
35 #include "AliKFVertex.h"
36
37 #include "AliTrackPointArray.h"
38 #include "TCint.h"
39 #include "AliTPCcalibV0.h"
40 #include "AliV0.h"
41
42
43
44
45
46 ClassImp(AliTPCcalibV0)
47
48
49 AliTPCcalibV0::AliTPCcalibV0() : 
50    AliTPCcalibBase(),
51    fStack(0),
52    fESD(0),
53    fPdg(0),
54    fParticles(0),
55    fV0s(0),
56    fGammas(0),
57    fV0type(0),
58    fTPCdEdx(0),
59    fTPCdEdxPi(0),
60    fTPCdEdxEl(0),
61    fTPCdEdxP(0)
62 {
63   fPdg = new TDatabasePDG;       
64   // create output histograms
65   fTPCdEdx   = new TH2F("TPCdEdX",  "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
66   BinLogX(fTPCdEdx); 
67   fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
68   fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);  
69   fTPCdEdxP  = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);        
70   
71   
72 }   
73
74 AliTPCcalibV0::~AliTPCcalibV0(){
75   //
76   //
77   //
78 }
79
80
81
82
83
84 void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){
85   //
86   //
87   //
88   fESD = esd;
89   AliKFParticle::SetField(esd->GetMagneticField());
90   MakeV0s();
91   if (stack) {
92     fStack = stack;
93     MakeMC();
94   }else{
95     fStack =0;
96   }
97 }
98
99 void AliTPCcalibV0::MakeMC(){
100   //
101   // MC comparison 
102   //    1. Select interesting particles
103   //    2. Assign the recosntructed particles 
104   //
105   //1. Select interesting particles
106   const Float_t kMinP   = 0.2;
107   const Float_t kMinPt  = 0.1;
108   const Float_t kMaxR   = 0.5;
109   const Float_t kMaxTan = 1.2;
110   const Float_t kMaxRad = 150;
111   //
112   if (!fParticles) fParticles = new TObjArray;
113   TParticle *part=0;
114   //  
115   Int_t entries = fStack->GetNtrack();
116   for (Int_t ipart=0; ipart<entries; ipart++){
117     part = fStack->Particle(ipart);
118     if (!part) continue;
119     if (part->P()<kMinP) continue;
120     if (part->R()>kMaxR) continue;
121     if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue;
122     Bool_t isInteresting = kFALSE;
123     if (part->GetPdgCode()==22) isInteresting =kTRUE;
124     if (part->GetPdgCode()==310) isInteresting =kTRUE;
125     if (part->GetPdgCode()==111) isInteresting =kTRUE;
126     if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE;
127
128     //
129     if (!isInteresting) continue;    
130     fParticles->AddLast(new TParticle(*part));
131   }
132   if (fParticles->GetEntries()<1) {
133     return;
134   }
135   //
136   //
137   //
138   Int_t sentries=fParticles->GetEntries();;
139   for (Int_t ipart=0; ipart<sentries; ipart++){
140     part = (TParticle*)fParticles->At(ipart);
141     TParticle *p0 = 0;
142     TParticle *p1 = 0;
143
144     Int_t nold =0;
145     Int_t nnew =0;
146     Int_t id0  = part->GetDaughter(0);
147     Int_t id1  = part->GetDaughter(1);    
148     if (id0>=fStack->GetNtrack() ) id0*=-1;
149     if (id1>=fStack->GetNtrack() ) id1*=-1;
150     Bool_t findable = kTRUE;
151     if (id0<0 || id1<0) findable = kFALSE;
152     Int_t charge =0; 
153     if (findable){
154       p0 = fStack->Particle(id0);
155       if (p0->R()>kMaxRad) findable = kFALSE;
156       if (p0->Pt()<kMinPt) findable = kFALSE;
157       if (p0->Vz()>250) findable= kFALSE;
158       if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
159       if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE;
160       else
161         if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++;
162           
163       p1 = fStack->Particle(id1);
164       if (p1->R()>kMaxRad) findable = kFALSE;
165       if (p1->Pt()<kMinPt) findable = kFALSE;
166       if (TMath::Abs(p1->Vz())>250) findable= kFALSE;
167       if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
168       if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE;
169       else
170         if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++;
171                           
172     }
173     //   (*fDebugStream)<<"MC0"<<
174     //       "P.="<<part<<
175     //       "findable="<<findable<<
176     //       "id0="<<id0<<
177     //       "id1="<<id1<<
178     //       "\n";
179     if (!findable) continue;
180     Float_t minpt = TMath::Min(p0->Pt(), p1->Pt());
181     Int_t type=-1;
182     
183     //
184     // 
185     AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
186     for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
187       AliESDv0 * v0 = fESD->GetV0(ivertex);
188       // select coresponding track
189       AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
190       if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue;
191       AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
192       if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue;
193       TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel()));
194       TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel()));
195       //
196       //
197       if ( v0->GetOnFlyStatus()) nnew++;
198       if (!v0->GetOnFlyStatus()) nold++;
199       if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 &&  TMath::Abs(pp->GetPdgCode())==11) 
200         type =1;
201       if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 &&  TMath::Abs(pp->GetPdgCode())==211) 
202         type =0;
203       if (part->GetPdgCode()==3122){
204         if (TMath::Abs(pn->GetPdgCode())==210 ) type=2;
205         else type=3;
206       }
207       AliKFParticle *v0kf       = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
208       v0kf->SetProductionVertex( primVtx );
209       AliKFParticle *v0kfc       = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
210       v0kfc->SetProductionVertex( primVtx );
211       v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass());
212       Float_t chi2 = v0kf->GetChi2();
213       Float_t chi2C = v0kf->GetChi2();
214       //
215       //
216       TTreeSRedirector *cstream = GetDebugStreamer();
217       if (cstream){
218         (*cstream)<<"MCRC"<<
219           "P.="<<part<<
220           "type="<<type<<
221           "chi2="<<chi2<<
222           "chi2C="<<chi2C<<
223           "minpt="<<minpt<<
224           "id0="<<id0<<
225           "id1="<<id1<<
226           "Pn.="<<pn<<
227           "Pp.="<<pp<<
228           "tn.="<<trackN<<
229           "tp.="<<trackP<<
230           "nold.="<<nold<<
231           "nnew.="<<nnew<<
232           "v0.="<<v0<<
233           "v0kf.="<<v0kf<<
234           "v0kfc.="<<v0kfc<<
235           "\n";     
236         delete v0kf;
237         delete v0kfc; 
238         //
239       }
240       
241       if (cstream){
242         (*cstream)<<"MC"<<
243           "P.="<<part<<
244           "charge="<<charge<<
245           "type="<<type<<
246           "minpt="<<minpt<<
247           "id0="<<id0<<
248           "id1="<<id1<<
249           "P0.="<<p0<<
250           "P1.="<<p1<<
251           "nold="<<nold<<
252           "nnew="<<nnew<<
253           "\n";
254       }
255     }
256     fParticles->Delete(); 
257   }
258 }
259
260
261
262 void AliTPCcalibV0::MakeV0s(){
263   //
264   //
265   //
266   const Int_t kMinCluster=40;
267   const Float_t kMinR    =0;
268   if (! fV0s) fV0s = new TObjArray(10);
269   fV0s->Clear();
270   //
271   // Old V0 finder
272   //
273   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
274     AliESDv0 * v0 = fESD->GetV0(ivertex);
275     if (v0->GetOnFlyStatus()) continue;
276     fV0s->AddLast(v0);
277   }
278   ProcessV0(0);
279   fV0s->Clear(0);
280   //
281   // MI V0 finder
282   //
283   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
284     AliESDv0 * v0 = fESD->GetV0(ivertex);
285     if (!v0->GetOnFlyStatus()) continue;
286     fV0s->AddLast(v0);
287   }
288   ProcessV0(1);
289   fV0s->Clear();
290   return;
291   //
292   // combinatorial
293   //
294   Int_t ntracks = fESD->GetNumberOfTracks();
295   for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
296     AliESDtrack * track0 = fESD->GetTrack(itrack0);
297     if (track0->GetSign()>0) continue;
298     if ( track0->GetTPCNcls()<kMinCluster) continue;
299     if (track0->GetKinkIndex(0)>0) continue;    
300     //
301     for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
302       AliESDtrack * track1 = fESD->GetTrack(itrack1);
303       if (track1->GetSign()<0) continue;
304       if ( track1->GetTPCNcls()<kMinCluster) continue;
305       if (track1->GetKinkIndex(0)>0) continue;
306       //
307       //      AliExternalTrackParam param0(*track0);
308       // AliExternalTrackParam param1(*track1);
309       AliV0 vertex;
310       vertex.SetParamN(*track0);
311       vertex.SetParamP(*track1);
312       Float_t xyz[3];
313       xyz[0] = fESD->GetPrimaryVertex()->GetXv();
314       xyz[1] = fESD->GetPrimaryVertex()->GetYv();
315       xyz[2] = fESD->GetPrimaryVertex()->GetZv();
316       vertex.Update(xyz);
317       if (vertex.GetRr()<kMinR) continue;
318       if (vertex.GetDcaV0Daughters()>1.) continue;
319       if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue;
320       // if (vertex.GetPointAngle()<0.9) continue;
321       vertex.SetIndex(0,itrack0);
322       vertex.SetIndex(1,itrack1);      
323       fV0s->AddLast(new AliV0(vertex));
324     }
325   }
326   ProcessV0(2);
327   for (Int_t i=0;i<fV0s->GetEntries(); i++) delete fV0s->At(i);
328   fV0s->Clear();
329 }
330
331
332
333
334
335
336 // void AliTPCcalibV0::ProcessV0(Int_t ftype){
337 //   //
338 //   //
339 //   const Double_t ktimeK0     = 2.684;
340 //   const Double_t ktimeLambda = 7.89; 
341   
342   
343 //   if (! fGammas) fGammas = new TObjArray(10);
344 //   fGammas->Clear();
345 //   Int_t nV0s  = fV0s->GetEntries();
346 //   if (nV0s==0) return;
347 //   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
348 //   //
349 //   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
350 //     AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
351 //     AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
352 //     AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
353 //     // 
354 //     // 
355 //     //
356 //     AliKFParticle *v0K0       = Fit(primVtx,v0,211,211);
357 //     AliKFParticle *v0Gamma    = Fit(primVtx,v0,11,-11);
358 //     AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211);
359 //     AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212);
360 //     //Set production vertex
361 //     v0K0->SetProductionVertex( primVtx );
362 //     v0Gamma->SetProductionVertex( primVtx );
363 //     v0Lambda42->SetProductionVertex( primVtx );
364 //     v0Lambda24->SetProductionVertex( primVtx );
365 //     Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
366 //     v0K0->GetMass( massK0,massSigma);
367 //     v0Gamma->GetMass( massGamma,massSigma);
368 //     v0Lambda42->GetMass( massLambda42,massSigma);
369 //     v0Lambda24->GetMass( massLambda24,massSigma);
370 //     Float_t chi2K0       = v0K0->GetChi2()/v0K0->GetNDF();
371 //     Float_t chi2Gamma    = v0Gamma->GetChi2()/v0Gamma->GetNDF();
372 //     Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
373 //     Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
374 //     //
375 //     // Mass Contrained params
376 //     //
377 //     AliKFParticle *v0K0C       = Fit(primVtx,v0,211,211);
378 //     AliKFParticle *v0GammaC    = Fit(primVtx,v0,11,-11);
379 //     AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211);
380 //     AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212);
381 //     //   
382 //     v0K0C->SetProductionVertex( primVtx );
383 //     v0GammaC->SetProductionVertex( primVtx );
384 //     v0Lambda42C->SetProductionVertex( primVtx );
385 //     v0Lambda24C->SetProductionVertex( primVtx );
386
387 //     v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
388 //     v0GammaC->SetMassConstraint(0);
389 //     v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
390 //     v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
391 //     //    
392 //     Double_t timeK0, sigmaTimeK0;  
393 //     Double_t timeLambda42, sigmaTimeLambda42;  
394 //     Double_t timeLambda24, sigmaTimeLambda24;  
395 //     v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
396 //     //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
397 //     v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
398 //     v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
399     
400
401 //     //
402 //     Float_t chi2K0C       = v0K0C->GetChi2()/v0K0C->GetNDF();
403 //     if (chi2K0C<0) chi2K0C=100;
404 //     Float_t chi2GammaC    = v0GammaC->GetChi2()/v0GammaC->GetNDF();
405 //     if (chi2GammaC<0) chi2GammaC=100;
406 //     Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
407 //     if (chi2Lambda42C<0) chi2Lambda42C=100;
408 //     Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
409 //     if (chi2Lambda24C<0) chi2Lambda24C=100;
410 //     //
411 //     Float_t  minChi2C=99;
412 //     Int_t   type   =-1;
413 //     if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
414 //     if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
415 //     if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
416 //     if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
417 //     Float_t  minChi2=99;
418 //     Int_t   type0   =-1;
419 //     if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
420 //     if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
421 //     if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
422 //     if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
423 //     Float_t betaGammaP  = trackN->GetP()/fPdg->GetParticle(-2212)->Mass(); 
424 //     Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass();
425 //     Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass();
426 //     Float_t dedxTeorP = BetheBlochAleph(betaGammaP);
427 //     Float_t dedxTeorPi = BetheBlochAleph(betaGammaPi);;
428 //     Float_t dedxTeorEl = BetheBlochAleph(betaGammaEl);;
429 //     //
430 //     //
431 //     if (minChi2>50) continue;
432 //     (*fDebugStream)<<"V0"<<
433 //       "ftype="<<ftype<<
434 //       "v0.="<<v0<<
435 //       "trackN.="<<trackN<<
436 //       "trackP.="<<trackP<<
437 //       //
438 //       "dedxTeorP="<<dedxTeorP<<
439 //       "dedxTeorPi="<<dedxTeorPi<<
440 //       "dedxTeorEl="<<dedxTeorEl<<
441 //       //
442 //       "type="<<type<<
443 //       "chi2C="<<minChi2C<<
444 //       "v0K0.="<<v0K0<<
445 //       "v0Gamma.="<<v0Gamma<<
446 //       "v0Lambda42.="<<v0Lambda42<<
447 //       "v0Lambda24.="<<v0Lambda24<<
448 //       //
449 //       "chi20K0.="<<chi2K0<<
450 //       "chi2Gamma.="<<chi2Gamma<<
451 //       "chi2Lambda42.="<<chi2Lambda42<<
452 //       "chi2Lambda24.="<<chi2Lambda24<<
453 //       //
454 //       "chi20K0c.="<<chi2K0C<<
455 //       "chi2Gammac.="<<chi2GammaC<<
456 //       "chi2Lambda42c.="<<chi2Lambda42C<<
457 //       "chi2Lambda24c.="<<chi2Lambda24C<<
458 //       //
459 //       "v0K0C.="<<v0K0C<<
460 //       "v0GammaC.="<<v0GammaC<<
461 //       "v0Lambda42C.="<<v0Lambda42C<<
462 //       "v0Lambda24C.="<<v0Lambda24C<<
463 //       //
464 //       "massK0="<<massK0<<
465 //       "massGamma="<<massGamma<<
466 //       "massLambda42="<<massLambda42<<
467 //       "massLambda24="<<massLambda24<<
468 //       //
469 //       "timeK0="<<timeK0<<
470 //       "timeLambda42="<<timeLambda42<<
471 //       "timeLambda24="<<timeLambda24<<
472 //       "\n";
473 //     if (type==1) fGammas->AddLast(v0); 
474 //     //
475 //     //
476 //     //
477 //     delete v0K0;
478 //     delete v0Gamma;
479 //     delete v0Lambda42;
480 //     delete v0Lambda24;    
481 //     delete v0K0C;
482 //     delete v0GammaC;
483 //     delete v0Lambda42C;
484 //     delete v0Lambda24C;    
485 //   }
486 //   ProcessPI0(); 
487 // }
488
489
490
491
492 void AliTPCcalibV0::ProcessV0(Int_t ftype){
493   //
494   //
495   
496   if (! fGammas) fGammas = new TObjArray(10);
497   fGammas->Clear();
498   Int_t nV0s  = fV0s->GetEntries();
499   if (nV0s==0) return;
500   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
501   //
502   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
503     AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
504     AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track
505     AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track
506     
507     const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam();
508     const AliExternalTrackParam * paramInPos = trackP->GetInnerParam();
509   
510     if (!paramInPos) continue; // in case the inner paramters do not exist
511     if (!paramInNeg) continue;
512     // 
513     // 
514     //
515     AliKFParticle *v0K0       = Fit(primVtx,v0,-211,211);
516     AliKFParticle *v0Gamma    = Fit(primVtx,v0,11,-11);
517     AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211);
518     AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212);
519     //Set production vertex
520     v0K0->SetProductionVertex( primVtx );
521     v0Gamma->SetProductionVertex( primVtx );
522     v0Lambda42->SetProductionVertex( primVtx );
523     v0Lambda24->SetProductionVertex( primVtx );
524     Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
525     v0K0->GetMass( massK0,massSigma);
526     v0Gamma->GetMass( massGamma,massSigma);
527     v0Lambda42->GetMass( massLambda42,massSigma);
528     v0Lambda24->GetMass( massLambda24,massSigma);
529     Float_t chi2K0       = v0K0->GetChi2()/v0K0->GetNDF();
530     Float_t chi2Gamma    = v0Gamma->GetChi2()/v0Gamma->GetNDF();
531     Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
532     Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
533     //
534     // Mass Contrained params
535     //
536     AliKFParticle *v0K0C       = Fit(primVtx,v0,-211,211);
537     AliKFParticle *v0GammaC    = Fit(primVtx,v0,11,-11);
538     AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar
539     AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda
540     //   
541     v0K0C->SetProductionVertex( primVtx );
542     v0GammaC->SetProductionVertex( primVtx );
543     v0Lambda42C->SetProductionVertex( primVtx );
544     v0Lambda24C->SetProductionVertex( primVtx );
545
546     v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
547     v0GammaC->SetMassConstraint(0);
548     v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
549     v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
550     //    
551     Double_t timeK0, sigmaTimeK0;  
552     Double_t timeLambda42, sigmaTimeLambda42;  
553     Double_t timeLambda24, sigmaTimeLambda24;  
554     v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
555     //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
556     v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
557     v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
558     
559
560     //
561     Float_t chi2K0C       = v0K0C->GetChi2()/v0K0C->GetNDF();
562     if (chi2K0C<0) chi2K0C=100;
563     Float_t chi2GammaC    = v0GammaC->GetChi2()/v0GammaC->GetNDF();
564     if (chi2GammaC<0) chi2GammaC=100;
565     Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
566     if (chi2Lambda42C<0) chi2Lambda42C=100;
567     Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
568     if (chi2Lambda24C<0) chi2Lambda24C=100;
569     //
570     Float_t  minChi2C=99;
571     Int_t   type   =-1;
572     if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
573     if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
574     if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
575     if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
576     Float_t  minChi2=99;
577     Int_t   type0   =-1;
578     if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
579     if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
580     if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
581     if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
582     
583      // 0 is  negative particle; 1 is positive particle
584     Float_t betaGamma0 = 0;
585     Float_t betaGamma1 = 0;
586     
587     switch (type) {
588      case 0:
589       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
590       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
591       break;
592      case 1:
593       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass();
594       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass();
595       break;
596      case 2:
597       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass();
598       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
599       break;
600      case 3:
601       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
602       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass();
603       break;
604     }
605  
606     // cuts and histogram filling
607     Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2
608         
609     if (minChi2C < 2 && ftype == 1) {
610      //
611      if (chi2K0C < 10*minChi2C) numCand++;
612      if (chi2GammaC < 10*minChi2C) numCand++;
613      if (chi2Lambda42C < 10*minChi2C) numCand++;
614      if (chi2Lambda24C < 10*minChi2C) numCand++;
615      //
616      if (numCand < 2) {
617       if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal());
618       if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal());
619      }    
620     }
621     
622     //
623     //
624     // write output tree
625     if (minChi2>50) continue;
626     TTreeSRedirector *cstream = GetDebugStreamer();
627     if (cstream){
628       (*cstream)<<"V0"<<
629         "ftype="<<ftype<<
630         "v0.="<<v0<<
631         "trackN.="<<trackN<<
632         "trackP.="<<trackP<<
633         //
634         "betaGamma0="<<betaGamma0<<
635         "betaGamma1="<<betaGamma1<<
636         //
637         "type="<<type<<
638         "chi2C="<<minChi2C<<
639         "v0K0.="<<v0K0<<
640         "v0Gamma.="<<v0Gamma<<
641         "v0Lambda42.="<<v0Lambda42<<
642         "v0Lambda24.="<<v0Lambda24<<
643         //
644         "chi20K0.="<<chi2K0<<
645         "chi2Gamma.="<<chi2Gamma<<
646         "chi2Lambda42.="<<chi2Lambda42<<
647         "chi2Lambda24.="<<chi2Lambda24<<
648         //
649         "chi20K0c.="<<chi2K0C<<
650         "chi2Gammac.="<<chi2GammaC<<
651         "chi2Lambda42c.="<<chi2Lambda42C<<
652         "chi2Lambda24c.="<<chi2Lambda24C<<
653         //
654         "v0K0C.="<<v0K0C<<
655         "v0GammaC.="<<v0GammaC<<
656         "v0Lambda42C.="<<v0Lambda42C<<
657         "v0Lambda24C.="<<v0Lambda24C<<
658         //
659         "massK0="<<massK0<<
660         "massGamma="<<massGamma<<
661         "massLambda42="<<massLambda42<<
662         "massLambda24="<<massLambda24<<
663         //
664         "timeK0="<<timeK0<<
665         "timeLambda42="<<timeLambda42<<
666         "timeLambda24="<<timeLambda24<<
667         "\n";
668     }
669     if (type==1) fGammas->AddLast(v0); 
670     //
671     //
672     //
673     delete v0K0;
674     delete v0Gamma;
675     delete v0Lambda42;
676     delete v0Lambda24;    
677     delete v0K0C;
678     delete v0GammaC;
679     delete v0Lambda42C;
680     delete v0Lambda24C; 
681   }
682   ProcessPI0(); 
683 }
684
685
686
687 void AliTPCcalibV0::ProcessPI0(){
688   //
689   //
690   //
691   Int_t nentries = fGammas->GetEntries();
692   if (nentries<2) return;
693   // 
694   Double_t m0[3], m1[3];
695   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
696   for (Int_t i0=0; i0<nentries; i0++){
697     AliESDv0 *v00 = (AliESDv0*)fGammas->At(i0); 
698     v00->GetPxPyPz (m0[0], m0[1], m0[2]);
699     AliKFParticle *p00 = Fit(primVtx, v00, 11,-11);
700     p00->SetProductionVertex( primVtx );
701     p00->SetMassConstraint(0);
702     //
703     for (Int_t i1=i0; i1<nentries; i1++){
704       AliESDv0 *v01 = (AliESDv0*)fGammas->At(i1);
705       v01->GetPxPyPz (m1[0], m1[1], m1[2]);
706       AliKFParticle *p01 = Fit(primVtx, v01, 11,-11);
707       p01->SetProductionVertex( primVtx );
708       p01->SetMassConstraint(0);
709       if (v00->GetIndex(0) != v01->GetIndex(0) && 
710           v00->GetIndex(1) != v01->GetIndex(1)){
711         AliKFParticle pi0( *p00,*p01); 
712         pi0.SetProductionVertex(primVtx);
713         Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]);
714         Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]);
715         Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2])));
716         TTreeSRedirector *cstream = GetDebugStreamer();
717         if (cstream){
718           (*cstream)<<"PI0"<<
719             "v00.="<<v00<<
720             "v01.="<<v01<<
721             "mass="<<mass<<
722             "p00.="<<p00<<
723             "p01.="<<p01<<
724             "pi0.="<<&pi0<<
725             "\n";       
726         }
727       }
728       delete p01;
729     }
730     delete p00;
731   }
732 }
733
734
735
736
737
738 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
739   //
740   // Make KF Particle
741   //
742   AliKFParticle p1( *(v0->GetParamN()), PDG1 );
743   AliKFParticle p2( *(v0->GetParamP()), PDG2 );
744   AliKFParticle *V0 = new AliKFParticle;
745   Double_t x, y, z;
746   v0->GetXYZ(x,y,z );
747   V0->SetVtxGuess(x,y,z);
748   *(V0)+=p1;
749   *(V0)+=p2;
750   return V0;  
751 }
752
753 TH2F * AliTPCcalibV0::GetHistograms() {
754   //
755   //
756   //
757  return fTPCdEdx;
758 }
759
760
761
762
763 void AliTPCcalibV0::BinLogX(TH2F *h) {
764   //
765   //
766   //
767    TAxis *axis = h->GetXaxis();
768    int bins = axis->GetNbins();
769
770    Double_t from = axis->GetXmin();
771    Double_t to = axis->GetXmax();
772    Double_t *new_bins = new Double_t[bins + 1];
773    
774    new_bins[0] = from;
775    Double_t factor = pow(to/from, 1./bins);
776   
777    for (int i = 1; i <= bins; i++) {
778      new_bins[i] = factor * new_bins[i-1];
779    }
780    axis->Set(bins, new_bins);
781    delete [] new_bins;
782    
783 }
784
785
786
787