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