AliTPCcalibTimeGain.cxx - Adding the Gamma conversion selected electorns
[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 #include "TRandom.h"
42 #include "TTreeStream.h"
43 #include "AliTPCcalibDB.h"
44 #include "AliTPCCorrection.h"
45 #include "AliGRPObject.h"
46 #include "AliTPCTransform.h"
47
48
49
50
51
52 ClassImp(AliTPCcalibV0)
53
54
55 AliTPCcalibV0::AliTPCcalibV0() : 
56    AliTPCcalibBase(),
57    fV0Tree(0),
58    fHPTTree(0),
59    fStack(0),
60    fESD(0),
61    fPdg(0),
62    fParticles(0),
63    fV0s(0),
64    fGammas(0)
65 {
66   
67 }   
68 AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
69    AliTPCcalibBase(),
70    fV0Tree(0),
71    fHPTTree(0),
72    fStack(0),
73    fESD(0),
74    fPdg(0),
75    fParticles(0),
76    fV0s(0),
77    fGammas(0)
78 {
79   fPdg = new TDatabasePDG;       
80   // create output histograms
81   SetName(name);
82   SetTitle(title);
83 }   
84
85 AliTPCcalibV0::~AliTPCcalibV0(){
86   //
87   //
88   //
89   delete fV0Tree;
90   delete fHPTTree;
91 }
92
93
94
95
96
97 void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){
98   //
99   //
100   //
101   fESD = esd;
102   AliKFParticle::SetField(esd->GetMagneticField());
103   if (TMath::Abs(AliTracker::GetBz())<1) return;  
104   DumpToTree(esd);
105   DumpToTreeHPT(esd);
106   return;
107   //
108   MakeV0s();
109   if (stack) {
110     fStack = stack;
111     MakeMC();
112   }else{
113     fStack =0;
114   }
115 }
116
117 void  AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){
118   //
119   // Dump V0s fith full firend information to the 
120   // 
121   if (TMath::Abs(AliTracker::GetBz())<1) return;
122   const Int_t kMinCluster=110;
123   const Float_t kMinPt   =3.;
124   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
125   if (!esdFriend) {
126     Printf("ERROR: esdFriend not available");
127     return;
128   }
129   //
130   Int_t ntracks=esd->GetNumberOfTracks();
131   for (Int_t i=0;i<ntracks;++i) {
132     Bool_t isOK=kFALSE;
133     AliESDtrack *track = esd->GetTrack(i);
134     if (track->GetTPCncls()<kMinCluster) continue;
135     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
136     if (!friendTrack) continue;
137     if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured
138       if (track->Pt()>kMinPt) isOK=kTRUE;
139     }
140     if (TMath::Abs(AliTracker::GetBz())<1){  // require primary track for the B field OFF data
141       Bool_t isAccepted=kTRUE;
142       if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE;
143       if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE;
144       if (!track->IsOn(AliESDtrack::kTOFout))   isAccepted=kFALSE;
145       Float_t dvertex[2],cvertex[3]; 
146       track->GetImpactParametersTPC(dvertex,cvertex);
147       if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE;
148       if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE;
149       track->GetImpactParameters(dvertex,cvertex);
150       if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE;
151       if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE;
152       if (!isAccepted) isOK=kFALSE;
153     }
154     
155     if (!isOK) continue;
156     //
157
158     TObject *calibObject;
159     AliTPCseed *seed = 0;
160     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
161       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
162     }
163     if (!seed) continue;
164       if (!fHPTTree) {
165       fHPTTree = new TTree("HPT","HPT");
166       fHPTTree->SetDirectory(0);
167     }
168     if (fHPTTree->GetEntries()==0){
169       //
170       fHPTTree->SetDirectory(0);
171       fHPTTree->Branch("t.",&track);
172       fHPTTree->Branch("ft.",&friendTrack);
173       fHPTTree->Branch("s.",&seed);
174     }else{
175       fHPTTree->SetBranchAddress("t.",&track);
176       fHPTTree->SetBranchAddress("ft.",&friendTrack);
177       fHPTTree->SetBranchAddress("s.",&seed);
178     }
179     fHPTTree->Fill();
180     //
181   }
182 }
183
184
185
186 void  AliTPCcalibV0::DumpToTree(AliESDEvent *esd){
187   //
188   // Dump V0s fith full firend information to the 
189   // 
190   Int_t nV0s  = fESD->GetNumberOfV0s();
191   const Int_t kMinCluster=110;
192   const Double_t kDownscale=0.01;
193   const Float_t kMinR    =0;
194   const Float_t kMinPt   =1.0;
195   const Float_t kMinMinPt   =0.7;
196   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
197   if (!esdFriend) {
198     Printf("ERROR: esdFriend not available");
199     return;
200   }
201   //
202   
203   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
204     Bool_t isOK=kFALSE;
205     AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
206     AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
207     AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track 
208     if (track0->GetTPCNcls()<kMinCluster) continue;
209     if (track0->GetKinkIndex(0)>0) continue;    
210     if (track1->GetTPCNcls()<kMinCluster) continue;
211     if (track1->GetKinkIndex(0)>0) continue;
212     if (v0->GetOnFlyStatus()==kFALSE) continue;
213     //
214     if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
215     //
216     //
217     if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
218     if (gRandom->Rndm()<kDownscale) isOK=kTRUE;  
219     if (!isOK) continue;
220     //
221     AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0));
222     if (!ftrack0) continue;
223     AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1));
224     if (!ftrack1) continue;
225     //
226     TObject *calibObject;
227     AliTPCseed *seed0 = 0;
228     AliTPCseed *seed1 = 0;
229     for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
230       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
231     }
232     for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
233       if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
234     }
235     if (!seed0) continue;
236     if (!seed1) continue;
237     AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam();
238     AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam();
239     if (!paramIn0) continue;
240     if (!paramIn1) continue;
241     //
242     //
243     if (!fV0Tree) {
244       fV0Tree = new TTree("V0s","V0s");
245       fV0Tree->SetDirectory(0);
246     }
247     if (fV0Tree->GetEntries()==0){
248       //
249       fV0Tree->SetDirectory(0);
250       fV0Tree->Branch("v0.",&v0);
251       fV0Tree->Branch("t0.",&track0);
252       fV0Tree->Branch("t1.",&track1);
253       fV0Tree->Branch("ft0.",&ftrack0);
254       fV0Tree->Branch("ft1.",&ftrack1);
255       fV0Tree->Branch("s0.",&seed0);
256       fV0Tree->Branch("s1.",&seed1);
257     }else{
258       fV0Tree->SetBranchAddress("v0.",&v0);
259       fV0Tree->SetBranchAddress("t0.",&track0);
260       fV0Tree->SetBranchAddress("t1.",&track1);
261       fV0Tree->SetBranchAddress("ft0.",&ftrack0);
262       fV0Tree->SetBranchAddress("ft1.",&ftrack1);
263       fV0Tree->SetBranchAddress("s0.",&seed0);
264       fV0Tree->SetBranchAddress("s1.",&seed1);
265     }
266     fV0Tree->Fill();
267   }
268 }
269
270
271 Long64_t AliTPCcalibV0::Merge(TCollection *const li) {
272
273   TIterator* iter = li->MakeIterator();
274   AliTPCcalibV0* cal = 0;
275
276   while ((cal = (AliTPCcalibV0*)iter->Next())) {
277     if (cal->fV0Tree){
278       if (!fV0Tree) {
279         fV0Tree = new TTree("V0s","V0s");
280         fV0Tree->SetDirectory(0);
281       }
282       if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree);
283       if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree);
284     }    
285   }
286   return 0;
287 }
288
289
290 void AliTPCcalibV0::AddTree(TTree * treeInput){
291   //
292   // Add the content of tree: 
293   // Notice automatic copy of tree in ROOT does not work for such complicated tree
294   //  
295   AliESDv0 * v0 = new AliESDv0;
296   Double_t kMinPt=0.8;
297   AliESDtrack * track0 = 0; // negative track
298   AliESDtrack * track1 = 0; // positive track 
299   AliESDfriendTrack *ftrack0 = 0;
300   AliESDfriendTrack *ftrack1 = 0;
301   AliTPCseed *seed0 = 0;
302   AliTPCseed *seed1 = 0;
303   treeInput->SetBranchStatus("ft0.",kFALSE);
304   treeInput->SetBranchStatus("ft1.",kFALSE);
305   TDatabasePDG pdg;
306   Double_t massK0= pdg.GetParticle("K0")->Mass();
307   Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
308
309   Int_t entries= treeInput->GetEntries();
310   for (Int_t i=0; i<entries; i++){
311     treeInput->SetBranchAddress("v0.",&v0);
312     treeInput->SetBranchAddress("t0.",&track0);
313     treeInput->SetBranchAddress("t1.",&track1);
314     treeInput->SetBranchAddress("ft0.",&ftrack0);
315     treeInput->SetBranchAddress("ft1.",&ftrack1);
316     treeInput->SetBranchAddress("s0.",&seed0);
317     treeInput->SetBranchAddress("s1.",&seed1);
318     if (fV0Tree->GetEntries()==0){
319       fV0Tree->SetDirectory(0);
320       fV0Tree->Branch("v0.",&v0);
321       fV0Tree->Branch("t0.",&track0);
322       fV0Tree->Branch("t1.",&track1);
323       fV0Tree->Branch("ft0.",&ftrack0);
324       fV0Tree->Branch("ft1.",&ftrack1);
325       fV0Tree->Branch("s0.",&seed0);
326       fV0Tree->Branch("s1.",&seed1);
327     }else{
328       fV0Tree->SetBranchAddress("v0.",&v0);
329       fV0Tree->SetBranchAddress("t0.",&track0);
330       fV0Tree->SetBranchAddress("t1.",&track1);
331       fV0Tree->SetBranchAddress("ft0.",&ftrack0);
332       fV0Tree->SetBranchAddress("ft1.",&ftrack1);
333       fV0Tree->SetBranchAddress("s0.",&seed0);
334       fV0Tree->SetBranchAddress("s1.",&seed1);
335     }
336     //
337     treeInput->GetEntry(i);
338     ftrack0->GetCalibContainer()->SetOwner(kTRUE);
339     ftrack1->GetCalibContainer()->SetOwner(kTRUE);
340     Bool_t isOK=kTRUE;
341     if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE;
342     if (track0->GetTPCncls()<100) isOK=kFALSE;
343     if (track1->GetTPCncls()<100) isOK=kFALSE;    
344     if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE;
345     if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
346     Bool_t isV0=kFALSE;    
347     if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05)     isV0=kTRUE;
348     if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE; 
349     if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE;
350     if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons
351     if (!isV0) isOK=kFALSE;
352     if (isOK) fV0Tree->Fill();
353     delete v0;
354     delete track0;
355     delete track1;
356     delete ftrack0;
357     delete ftrack1;
358     delete seed0;
359     delete seed1;
360     v0=0;
361     track0=0;
362     track1=0;
363     ftrack0=0;
364     ftrack1=0;
365     seed0=0;
366     seed1=0;
367   }
368 }
369
370 void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){
371   //
372   // Add the content of tree: 
373   // Notice automatic copy of tree in ROOT does not work for such complicated tree
374   //  
375   AliESDtrack *track = 0;
376   AliESDfriendTrack *friendTrack = 0;
377   AliTPCseed *seed = 0;
378   if (!treeInput) return;
379   if (treeInput->GetEntries()==0) return;
380   //
381   Int_t entries= treeInput->GetEntries();  
382   //
383   for (Int_t i=0; i<entries; i++){
384     track=0;
385     friendTrack=0;
386     seed=0;
387     //
388     treeInput->SetBranchAddress("t.",&track);
389     treeInput->SetBranchAddress("ft.",&friendTrack);
390     treeInput->SetBranchAddress("s.",&seed);
391     treeInput->GetEntry(i);
392     //
393     if (fHPTTree->GetEntries()==0){
394       fHPTTree->SetDirectory(0);
395       fHPTTree->Branch("t.",&track);
396       fHPTTree->Branch("ft.",&friendTrack);
397       fHPTTree->Branch("s.",&seed);
398     }else{
399       fHPTTree->SetBranchAddress("t.",&track);
400       fHPTTree->SetBranchAddress("ft.",&friendTrack);
401       fHPTTree->SetBranchAddress("s.",&seed);
402     }    
403     Bool_t isOK=kTRUE;
404     if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE;
405     if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE;
406     if (isOK) fHPTTree->Fill();
407     //
408     delete track;
409     delete friendTrack;
410     delete seed;
411   }
412 }
413
414
415 void AliTPCcalibV0::MakeMC(){
416   //
417   // MC comparison 
418   //    1. Select interesting particles
419   //    2. Assign the recosntructed particles 
420   //
421   //1. Select interesting particles
422   const Float_t kMinP   = 0.2;
423   const Float_t kMinPt  = 0.1;
424   const Float_t kMaxR   = 0.5;
425   const Float_t kMaxTan = 1.2;
426   const Float_t kMaxRad = 150;
427   //
428   if (!fParticles) fParticles = new TObjArray;
429   TParticle *part=0;
430   //  
431   Int_t entries = fStack->GetNtrack();
432   for (Int_t ipart=0; ipart<entries; ipart++){
433     part = fStack->Particle(ipart);
434     if (!part) continue;
435     if (part->P()<kMinP) continue;
436     if (part->R()>kMaxR) continue;
437     if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue;
438     Bool_t isInteresting = kFALSE;
439     if (part->GetPdgCode()==22) isInteresting =kTRUE;
440     if (part->GetPdgCode()==310) isInteresting =kTRUE;
441     if (part->GetPdgCode()==111) isInteresting =kTRUE;
442     if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE;
443
444     //
445     if (!isInteresting) continue;    
446     fParticles->AddLast(new TParticle(*part));
447   }
448   if (fParticles->GetEntries()<1) {
449     return;
450   }
451   //
452   //
453   //
454   Int_t sentries=fParticles->GetEntries();;
455   for (Int_t ipart=0; ipart<sentries; ipart++){
456     part = (TParticle*)fParticles->At(ipart);
457     TParticle *p0 = 0;
458     TParticle *p1 = 0;
459
460     Int_t nold =0;
461     Int_t nnew =0;
462     Int_t id0  = part->GetDaughter(0);
463     Int_t id1  = part->GetDaughter(1);    
464     if (id0>=fStack->GetNtrack() ) id0*=-1;
465     if (id1>=fStack->GetNtrack() ) id1*=-1;
466     Bool_t findable = kTRUE;
467     if (id0<0 || id1<0) findable = kFALSE;
468     Int_t charge =0; 
469     if (findable){
470       p0 = fStack->Particle(id0);
471       if (p0->R()>kMaxRad) findable = kFALSE;
472       if (p0->Pt()<kMinPt) findable = kFALSE;
473       if (p0->Vz()>250) findable= kFALSE;
474       if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
475       if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE;
476       else
477         if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++;
478           
479       p1 = fStack->Particle(id1);
480       if (p1->R()>kMaxRad) findable = kFALSE;
481       if (p1->Pt()<kMinPt) findable = kFALSE;
482       if (TMath::Abs(p1->Vz())>250) findable= kFALSE;
483       if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
484       if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE;
485       else
486         if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++;
487                           
488     }
489     //   (*fDebugStream)<<"MC0"<<
490     //       "P.="<<part<<
491     //       "findable="<<findable<<
492     //       "id0="<<id0<<
493     //       "id1="<<id1<<
494     //       "\n";
495     if (!findable) continue;
496     Float_t minpt = TMath::Min(p0->Pt(), p1->Pt());
497     Int_t type=-1;
498     
499     //
500     // 
501     AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
502     for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
503       AliESDv0 * v0 = fESD->GetV0(ivertex);
504       // select coresponding track
505       AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
506       if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue;
507       AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
508       if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue;
509       TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel()));
510       TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel()));
511       //
512       //
513       if ( v0->GetOnFlyStatus()) nnew++;
514       if (!v0->GetOnFlyStatus()) nold++;
515       if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 &&  TMath::Abs(pp->GetPdgCode())==11) 
516         type =1;
517       if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 &&  TMath::Abs(pp->GetPdgCode())==211) 
518         type =0;
519       if (part->GetPdgCode()==3122){
520         if (TMath::Abs(pn->GetPdgCode())==210 ) type=2;
521         else type=3;
522       }
523       AliKFParticle *v0kf       = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
524       v0kf->SetProductionVertex( primVtx );
525       AliKFParticle *v0kfc       = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
526       v0kfc->SetProductionVertex( primVtx );
527       v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass());
528       Float_t chi2 = v0kf->GetChi2();
529       Float_t chi2C = v0kf->GetChi2();
530       //
531       //
532       TTreeSRedirector *cstream = GetDebugStreamer();
533       if (cstream){
534         (*cstream)<<"MCRC"<<
535           "P.="<<part<<
536           "type="<<type<<
537           "chi2="<<chi2<<
538           "chi2C="<<chi2C<<
539           "minpt="<<minpt<<
540           "id0="<<id0<<
541           "id1="<<id1<<
542           "Pn.="<<pn<<
543           "Pp.="<<pp<<
544           "tn.="<<trackN<<
545           "tp.="<<trackP<<
546           "nold.="<<nold<<
547           "nnew.="<<nnew<<
548           "v0.="<<v0<<
549           "v0kf.="<<v0kf<<
550           "v0kfc.="<<v0kfc<<
551           "\n";     
552         delete v0kf;
553         delete v0kfc; 
554         //
555       }
556       
557       if (cstream){
558         (*cstream)<<"MC"<<
559           "P.="<<part<<
560           "charge="<<charge<<
561           "type="<<type<<
562           "minpt="<<minpt<<
563           "id0="<<id0<<
564           "id1="<<id1<<
565           "P0.="<<p0<<
566           "P1.="<<p1<<
567           "nold="<<nold<<
568           "nnew="<<nnew<<
569           "\n";
570       }
571     }
572     fParticles->Delete(); 
573   }
574 }
575
576 void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t run){
577   //
578   // Make a fit tree
579   //
580   // 0. Loop over selected tracks
581   // 1. Loop over all transformation - refit the track with and without the
582   //    transformtation
583   // 2. Dump the matching paramaeters to the debugStremer
584   //
585   
586   //Connect input
587   const Int_t kMinNcl=120;
588   TFile f("TPCV0Objects.root");
589   AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
590   TTree * treeInput = v0TPC->GetHPTTree();
591   TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root");
592   AliESDtrack *track = 0;
593   AliESDfriendTrack *friendTrack = 0;
594   AliTPCseed *seed = 0;
595   if (!treeInput) return;
596   if (treeInput->GetEntries()==0) return;
597   //
598   treeInput->SetBranchAddress("t.",&track);
599   treeInput->SetBranchAddress("ft.",&friendTrack);
600   treeInput->SetBranchAddress("s.",&seed);
601   //
602   Int_t ncorr=0;
603   if (corrArray) ncorr = corrArray->GetEntries();
604   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
605   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
606   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
607   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
608   //
609   //
610   //  
611   Int_t ntracks= treeInput->GetEntries();
612   for (Int_t itrack=0; itrack<ntracks; itrack++){
613     treeInput->GetEntry(itrack);
614     if (!track) continue;
615     if (seed->Pt()<ptCut) continue;
616     if (track->Pt()<ptCut) continue;
617     if (track->GetTPCncls()<kMinNcl) continue;
618     //
619     // Reapply transformation
620     //
621     for (Int_t irow=0; irow<159; irow++){
622       AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
623       if (cluster &&cluster->GetX()>10){
624         Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
625         Int_t index0[1]={cluster->GetDetector()};
626         transform->Transform(x0,index0,0,1);
627         cluster->SetX(x0[0]);
628         cluster->SetY(x0[1]);
629         cluster->SetZ(x0[2]);
630         //
631       }
632     }    
633     //
634     Double_t xref=134;
635     AliExternalTrackParam* paramInner=0;
636     AliExternalTrackParam* paramOuter=0;
637     AliExternalTrackParam* paramIO=0;
638     Bool_t isOK=kTRUE;
639     for (Int_t icorr=-1; icorr<ncorr; icorr++){
640       //
641       AliTPCCorrection *corr = 0;
642       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
643       AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1);      
644       AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1);      
645       AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 ); 
646       if (trackInner&&trackOuter&&trackIO){
647         trackOuter->Rotate(trackInner->GetAlpha());
648         trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz());
649         if (icorr<0) {
650           paramInner=trackInner;
651           paramOuter=trackOuter;
652           paramIO=trackIO;
653           paramIO->Rotate(seed->GetAlpha());
654           paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz());
655         }
656       }else{
657         isOK=kFALSE;
658       }
659       
660     }
661     if (paramOuter&& paramInner) {
662       //      Bool_t isOK=kTRUE;
663       if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE;
664       if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE;
665       if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE;
666       if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE;      
667       (*pcstream)<<"fit"<<
668         "s.="<<seed<<
669         "io.="<<paramIO<<
670         "pIn.="<<paramInner<<
671         "pOut.="<<paramOuter;      
672     }
673     //
674   }
675   delete pcstream;
676   /*
677     .x ~/rootlogon.C
678     Int_t run=117112;
679     .x ../ConfigCalibTrain.C(run)
680     .L ../AddTaskTPCCalib.C
681     ConfigOCDB(run)
682     TFile fexb("../../RegisterCorrectionExB.root");
683     AliTPCComposedCorrection *cc=  (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
684     cc->Init();
685     cc->Print("DA"); // Print used correction classes
686     TObjArray *array = cc->GetCorrections()
687     AliTPCcalibV0::MakeFitTreeTrack(array,2,run);
688    
689    */
690 }
691
692 void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){
693   //
694   // Make a fit tree
695   //
696   // 0. Loop over selected tracks
697   // 1. Loop over all transformation - refit the track with and without the
698   //    transformtation
699   // 2. Dump the matching paramaeters to the debugStremer
700   //
701   
702   //Connect input
703   const Int_t kMinNcl=120;
704   TFile f("TPCV0Objects.root");
705   AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
706   TTree * treeInput = v0TPC->GetV0Tree();
707   TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root");
708   AliESDv0 *v0 = 0;
709   AliESDtrack *track0 = 0;
710   AliESDfriendTrack *friendTrack0 = 0;
711   AliTPCseed *seed0 = 0;
712   AliTPCseed *s0 = 0;
713   AliESDtrack *track1 = 0;
714   AliESDfriendTrack *friendTrack1 = 0;
715   AliTPCseed *s1 = 0;
716   AliTPCseed *seed1 = 0;
717   if (!treeInput) return;
718   if (treeInput->GetEntries()==0) return;
719   //
720   treeInput->SetBranchAddress("v0.",&v0);
721   treeInput->SetBranchAddress("t0.",&track0);
722   treeInput->SetBranchAddress("ft0.",&friendTrack0);
723   treeInput->SetBranchAddress("s0.",&s0);
724   treeInput->SetBranchAddress("t1.",&track1);
725   treeInput->SetBranchAddress("ft1.",&friendTrack1);
726   treeInput->SetBranchAddress("s1.",&s1);
727   //
728   TDatabasePDG pdg;
729   Int_t ncorr=0;
730   if (corrArray) ncorr = corrArray->GetEntries();
731   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
732   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
733   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
734   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
735   Double_t massK0= pdg.GetParticle("K0")->Mass();
736   Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
737   Double_t massPion=pdg.GetParticle("pi+")->Mass();
738   Double_t massProton=pdg.GetParticle("proton")->Mass();
739   Double_t pdgPion=pdg.GetParticle("pi+")->PdgCode();
740   Double_t pdgProton=pdg.GetParticle("proton")->PdgCode();
741   Double_t mass0=0;
742   Double_t mass1=0;
743   Double_t massV0=0;
744   Int_t    pdg0=0;
745   Int_t    pdg1=0;
746   //
747   //
748   //  
749   Int_t nv0s= treeInput->GetEntries();
750   for (Int_t iv0=0; iv0<nv0s; iv0++){
751     Int_t  v0Type=0;
752     Int_t isK0=0;
753     Int_t isLambda=0;
754     Int_t isAntiLambda=0;
755     treeInput->GetEntry(iv0);
756     if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s    
757     if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda   
758     if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda
759     if (isK0+isLambda+isAntiLambda!=1) continue;
760     mass0=massPion;
761     mass1=massPion;
762     pdg0=pdgPion;
763     pdg1=pdgPion;
764     if (isLambda) {mass0=massProton; pdg0=pdgProton;}
765     if (isAntiLambda) {mass1=massProton; pdg1=pdgProton;}
766     massV0=massK0;
767     if (isK0==0) massV0=massLambda;
768     //
769     if (!s0) continue;
770     seed0=(s0->GetSign()>0)?s0:s1;
771     seed1=(s0->GetSign()>0)?s1:s0;
772     if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks
773     if (seed0->Pt()<ptCut) continue;
774     if (seed1->Pt()<ptCut) continue;
775     //
776     // Reapply transformation
777     //
778     for  (Int_t itype=0; itype<2; itype++){
779       AliTPCseed * seed = (itype==0) ? seed0: seed1;      
780       for (Int_t irow=0; irow<159; irow++){
781         AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
782         if (cluster &&cluster->GetX()>10){
783           Double_t x0[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
784           Int_t index0[1]={cluster->GetDetector()};
785           transform->Transform(x0,index0,0,1);
786           cluster->SetX(x0[0]);
787           cluster->SetY(x0[1]);
788           cluster->SetZ(x0[2]);
789           //
790         }
791       }
792     }   
793     Bool_t isOK=kTRUE;
794     Double_t radius = v0->GetRr();
795     Double_t xyz[3];
796     v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
797     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
798     TObjArray arrayV0in(ncorr+1);
799     TObjArray arrayV0io(ncorr+1);
800     TObjArray arrayT0(ncorr+1);
801     TObjArray arrayT1(ncorr+1);
802     arrayV0in.SetOwner(kTRUE);
803     arrayV0io.SetOwner(kTRUE);
804     //
805     for (Int_t icorr=-1; icorr<ncorr; icorr++){
806       AliTPCCorrection *corr =0;
807       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
808       //
809       AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,mass0);      
810       AliExternalTrackParam * trackIO0    = RefitTrack(seed0, corr,245,85,mass0);      
811       AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,mass1);      
812       AliExternalTrackParam * trackIO1    = RefitTrack(seed1, corr,245,85,mass1);      
813       if (!trackInner0) isOK=kFALSE;
814       if (!trackInner1) isOK=kFALSE;
815       if (!trackIO0)    isOK=kFALSE;
816       if (!trackIO1)    isOK=kFALSE;
817       if (isOK){
818         if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
819         if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
820         if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
821         if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
822         //
823         if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, mass0, 1, kFALSE)) isOK=kFALSE; 
824         if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, mass1, 1, kFALSE)) isOK=kFALSE; 
825         if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, mass0, 1, kFALSE)) isOK=kFALSE; 
826         if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, mass1, 1, kFALSE)) isOK=kFALSE; 
827         if (!isOK) continue;
828         arrayT0.AddAt(trackIO0->Clone(),icorr+1);
829         arrayT1.AddAt(trackIO1->Clone(),icorr+1);
830         Int_t charge=(trackIO0->GetSign());
831         AliKFParticle pin0( *trackInner0,  pdg0*charge);
832         AliKFParticle pin1( *trackInner1, -pdg1*charge);
833         AliKFParticle pio0( *trackIO0,  pdg0*charge);
834         AliKFParticle pio1( *trackIO1, -pdg1*charge);
835         AliKFParticle v0in;
836         AliKFParticle v0io;
837         v0in+=pin0;
838         v0in+=pin1;
839         v0io+=pio0;
840         v0io+=pio1;
841         arrayV0in.AddAt(v0in.Clone(),icorr+1);
842         arrayV0io.AddAt(v0io.Clone(),icorr+1);
843       }
844     }
845     if (!isOK) continue;
846     //
847     //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0);
848     AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0);
849     AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0);
850     AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0);
851     Double_t mass0=0, mass0E=0; 
852     pio0->GetMass( mass0,mass0E);
853     //
854     Double_t mean=mass0-massV0;
855     if (TMath::Abs(mean)>0.05) continue;
856     Double_t mass[10000];
857     //
858     Int_t dtype=30;  // id for V0
859     Int_t ptype=5;   // id for invariant mass
860     //    Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt()));      // K0s V0 asymetry
861     Int_t id=1000.*(param0->Pt()-param1->Pt());      // K0s V0 asymetry
862     Double_t gx,gy,gz, px,py,pz;
863     Double_t pt = v0->Pt();
864     v0->GetXYZ(gx,gy,gz);
865     v0->GetPxPyPz(px,py,pz);
866     Double_t theta=pz/TMath::Sqrt(px*px+py*py);
867     Double_t phi=TMath::ATan2(py,px);
868     Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp());
869     Double_t sector=9*phi/TMath::Pi();
870     Double_t dsec=sector-TMath::Nint(sector);
871     Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
872     //Int_t nentries=v0Type;
873     Double_t bz=AliTracker::GetBz();
874     Double_t dRrec=0;
875     (*pcstream)<<"fitDebug"<< 
876       "id="<<id<<
877       "v0.="<<v0<<
878       "mean="<<mean<<
879       "rms="<<mass0E<<
880       "pio0.="<<pio0<<
881       "p0.="<<param0<<
882       "p1.="<<param1;
883     (*pcstream)<<"fit"<<  // dump valus for fit
884       "run="<<run<<       //run number
885       "bz="<<bz<<         // magnetic filed used
886       "dtype="<<dtype<<   // detector match type 30
887       "ptype="<<ptype<<   // parameter type
888       "theta="<<theta<<   // theta
889       "phi="<<phi<<       // phi 
890       "snp="<<snp<<       // snp
891       "mean="<<mean<<     // mean dist value
892       "rms="<<mass0E<<       // rms
893       "sector="<<sector<<
894       "dsec="<<dsec<<
895       //
896       "refX="<<refX<<      // reference radius
897       "gx="<<gx<<         // global position
898       "gy="<<gy<<         // global position
899       "gz="<<gz<<         // global position
900       "dRrec="<<dRrec<<      // delta Radius in reconstruction
901       "pt="<<pt<<         // pt of the particle
902       "id="<<id<<     //delta of the momenta      
903       "entries="<<v0Type;//  type of the V0
904     for (Int_t icorr=0; icorr<ncorr; icorr++){
905       AliTPCCorrection *corr =0;
906       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
907       AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
908       AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
909       AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
910       AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
911       Double_t massE=0; 
912       pio->GetMass( mass[icorr],massE);
913       mass[icorr]-=mass0;
914       (*pcstream)<<"fit"<< 
915         Form("%s=",corr->GetName())<<mass[icorr];
916       (*pcstream)<<"fitDebug"<< 
917         Form("%s=",corr->GetName())<<mass[icorr]<<
918         Form("%sp0.=",corr->GetName())<<par0<<
919         Form("%sp1=",corr->GetName())<<par1;
920     }
921     (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";        
922     (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";        
923   }
924   delete pcstream;
925 }
926
927 AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
928   //
929   // Refit the track:
930   //    seed   - tpc track with cluster
931   //    corr   - distrotion/correction class  - apllied to the points
932   //    xstart - radius to start propagate/update
933   //    xstop  - radius to stop propagate/update
934   // 
935   const Double_t kResetCov=20.;
936   const Double_t kSigma=5.;
937   Double_t covar[15];
938   for (Int_t i=0;i<15;i++) covar[i]=0;
939   covar[0]=kSigma*kSigma;
940   covar[2]=kSigma*kSigma;
941   covar[5]=kSigma*kSigma/Float_t(150*150);
942   covar[9]=kSigma*kSigma/Float_t(150*150);
943   covar[14]=1*1;
944   // 
945   AliExternalTrackParam *refit  = new AliExternalTrackParam(*seed);
946   refit->PropagateTo(xstart, AliTracker::GetBz());
947   refit->AddCovariance(covar);
948   refit->ResetCovariance(kResetCov);
949   Double_t xmin = TMath::Min(xstart,xstop);
950   Double_t xmax = TMath::Max(xstart,xstop);
951   Int_t ncl=0;
952   //
953   Bool_t isOK=kTRUE;
954   for (Int_t index0=0; index0<159; index0++){
955     Int_t irow= (xstart<xstop)? index0:159-index0;
956     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);  //cluster in local system
957     if (!cluster) continue;
958     if (cluster->GetX()<xmin) continue;
959     if (cluster->GetX()>xmax) continue;
960     Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.;
961     if (!refit->Rotate(alpha)) isOK=kFALSE;
962     Double_t x     = cluster->GetX();
963     Double_t y     = cluster->GetY();
964     Double_t z     = cluster->GetZ();
965     if (corr){      
966       Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};  // original position
967       corr->DistortPointLocal(xyz,cluster->GetDetector());
968       x=xyz[0];
969       y=xyz[1];
970       z=xyz[2];
971     }
972     if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE;
973     if (!isOK) continue;
974     Double_t cov[3]={0.01,0.,0.01};
975     Double_t yz[2]={y,z};
976     if (!refit->Update(yz,cov)) isOK=kFALSE;    
977     ncl++;
978   }
979   if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE;
980   //  
981   if (!isOK) {
982     delete refit;
983     return 0;
984   }  
985   return refit;
986 }
987
988
989
990
991
992 //
993 // Obsolete part
994 //
995
996
997 void AliTPCcalibV0::MakeV0s(){
998   //
999   // 
1000   //
1001   const Int_t kMinCluster=40;
1002   const Float_t kMinR    =0;
1003   if (! fV0s) fV0s = new TObjArray(10);
1004   fV0s->Clear();
1005   //
1006   // Old V0 finder
1007   //
1008   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
1009     AliESDv0 * v0 = fESD->GetV0(ivertex);
1010     if (v0->GetOnFlyStatus()) continue;
1011     fV0s->AddLast(v0);
1012   }
1013   ProcessV0(0);
1014   fV0s->Clear(0);
1015   //
1016   // MI V0 finder
1017   //
1018   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
1019     AliESDv0 * v0 = fESD->GetV0(ivertex);
1020     if (!v0->GetOnFlyStatus()) continue;
1021     fV0s->AddLast(v0);
1022   }
1023   ProcessV0(1);
1024   fV0s->Clear();
1025   return;
1026   //
1027   // combinatorial
1028   //
1029   Int_t ntracks = fESD->GetNumberOfTracks();
1030   for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
1031     AliESDtrack * track0 = fESD->GetTrack(itrack0);
1032     if (track0->GetSign()>0) continue;
1033     if ( track0->GetTPCNcls()<kMinCluster) continue;
1034     if (track0->GetKinkIndex(0)>0) continue;    
1035     //
1036     for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
1037       AliESDtrack * track1 = fESD->GetTrack(itrack1);
1038       if (track1->GetSign()<0) continue;
1039       if ( track1->GetTPCNcls()<kMinCluster) continue;
1040       if (track1->GetKinkIndex(0)>0) continue;
1041       //
1042       //      AliExternalTrackParam param0(*track0);
1043       // AliExternalTrackParam param1(*track1);
1044       AliV0 vertex;
1045       vertex.SetParamN(*track0);
1046       vertex.SetParamP(*track1);
1047       Float_t xyz[3];
1048       xyz[0] = fESD->GetPrimaryVertex()->GetXv();
1049       xyz[1] = fESD->GetPrimaryVertex()->GetYv();
1050       xyz[2] = fESD->GetPrimaryVertex()->GetZv();
1051       vertex.Update(xyz);
1052       if (vertex.GetRr()<kMinR) continue;
1053       if (vertex.GetDcaV0Daughters()>1.) continue;
1054       if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue;
1055       // if (vertex.GetPointAngle()<0.9) continue;
1056       vertex.SetIndex(0,itrack0);
1057       vertex.SetIndex(1,itrack1);      
1058       fV0s->AddLast(new AliV0(vertex));
1059     }
1060   }
1061   ProcessV0(2);
1062   for (Int_t i=0;i<fV0s->GetEntries(); i++) delete fV0s->At(i);
1063   fV0s->Clear();
1064 }
1065
1066
1067
1068
1069
1070
1071
1072
1073 void AliTPCcalibV0::ProcessV0(Int_t ftype){
1074   //
1075   // Obsolete
1076   //  
1077   if (! fGammas) fGammas = new TObjArray(10);
1078   fGammas->Clear();
1079   Int_t nV0s  = fV0s->GetEntries();
1080   if (nV0s==0) return;
1081   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
1082   //
1083   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
1084     AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
1085     AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track
1086     AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track
1087     
1088     const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam();
1089     const AliExternalTrackParam * paramInPos = trackP->GetInnerParam();
1090   
1091     if (!paramInPos) continue; // in case the inner paramters do not exist
1092     if (!paramInNeg) continue;
1093     // 
1094     // 
1095     //
1096     AliKFParticle *v0K0       = Fit(primVtx,v0,-211,211);
1097     AliKFParticle *v0Gamma    = Fit(primVtx,v0,11,-11);
1098     AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211);
1099     AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212);
1100     //Set production vertex
1101     v0K0->SetProductionVertex( primVtx );
1102     v0Gamma->SetProductionVertex( primVtx );
1103     v0Lambda42->SetProductionVertex( primVtx );
1104     v0Lambda24->SetProductionVertex( primVtx );
1105     Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
1106     v0K0->GetMass( massK0,massSigma);
1107     v0Gamma->GetMass( massGamma,massSigma);
1108     v0Lambda42->GetMass( massLambda42,massSigma);
1109     v0Lambda24->GetMass( massLambda24,massSigma);
1110     Float_t chi2K0       = v0K0->GetChi2()/v0K0->GetNDF();
1111     Float_t chi2Gamma    = v0Gamma->GetChi2()/v0Gamma->GetNDF();
1112     Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
1113     Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
1114     //
1115     // Mass Contrained params
1116     //
1117     AliKFParticle *v0K0C       = Fit(primVtx,v0,-211,211);
1118     AliKFParticle *v0GammaC    = Fit(primVtx,v0,11,-11);
1119     AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar
1120     AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda
1121     //   
1122     v0K0C->SetProductionVertex( primVtx );
1123     v0GammaC->SetProductionVertex( primVtx );
1124     v0Lambda42C->SetProductionVertex( primVtx );
1125     v0Lambda24C->SetProductionVertex( primVtx );
1126
1127     v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
1128     v0GammaC->SetMassConstraint(0);
1129     v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
1130     v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
1131     //    
1132     Double_t timeK0, sigmaTimeK0;  
1133     Double_t timeLambda42, sigmaTimeLambda42;  
1134     Double_t timeLambda24, sigmaTimeLambda24;  
1135     v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
1136     //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
1137     v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
1138     v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
1139     
1140
1141     //
1142     Float_t chi2K0C       = v0K0C->GetChi2()/v0K0C->GetNDF();
1143     if (chi2K0C<0) chi2K0C=100;
1144     Float_t chi2GammaC    = v0GammaC->GetChi2()/v0GammaC->GetNDF();
1145     if (chi2GammaC<0) chi2GammaC=100;
1146     Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
1147     if (chi2Lambda42C<0) chi2Lambda42C=100;
1148     Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
1149     if (chi2Lambda24C<0) chi2Lambda24C=100;
1150     //
1151     Float_t  minChi2C=99;
1152     Int_t   type   =-1;
1153     if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
1154     if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
1155     if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
1156     if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
1157     Float_t  minChi2=99;
1158     Int_t   type0   =-1;
1159     if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
1160     if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
1161     if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
1162     if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
1163     
1164      // 0 is  negative particle; 1 is positive particle
1165     Float_t betaGamma0 = 0;
1166     Float_t betaGamma1 = 0;
1167     
1168     switch (type) {
1169      case 0:
1170       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
1171       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
1172       break;
1173      case 1:
1174       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass();
1175       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass();
1176       break;
1177      case 2:
1178       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass();
1179       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
1180       break;
1181      case 3:
1182       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
1183       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass();
1184       break;
1185     }
1186  
1187     // cuts and histogram filling
1188     Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2
1189         
1190     if (minChi2C < 2 && ftype == 1) {
1191      //
1192      if (chi2K0C < 10*minChi2C) numCand++;
1193      if (chi2GammaC < 10*minChi2C) numCand++;
1194      if (chi2Lambda42C < 10*minChi2C) numCand++;
1195      if (chi2Lambda24C < 10*minChi2C) numCand++;
1196      //
1197     }
1198     
1199     //
1200     //
1201     // write output tree
1202     if (minChi2>50) continue;
1203     TTreeSRedirector *cstream = GetDebugStreamer();
1204     if (cstream){
1205       (*cstream)<<"V0"<<
1206         "ftype="<<ftype<<
1207         "v0.="<<v0<<
1208         "trackN.="<<trackN<<
1209         "trackP.="<<trackP<<
1210         //
1211         "betaGamma0="<<betaGamma0<<
1212         "betaGamma1="<<betaGamma1<<
1213         //
1214         "type="<<type<<
1215         "chi2C="<<minChi2C<<
1216         "v0K0.="<<v0K0<<
1217         "v0Gamma.="<<v0Gamma<<
1218         "v0Lambda42.="<<v0Lambda42<<
1219         "v0Lambda24.="<<v0Lambda24<<
1220         //
1221         "chi20K0.="<<chi2K0<<
1222         "chi2Gamma.="<<chi2Gamma<<
1223         "chi2Lambda42.="<<chi2Lambda42<<
1224         "chi2Lambda24.="<<chi2Lambda24<<
1225         //
1226         "chi20K0c.="<<chi2K0C<<
1227         "chi2Gammac.="<<chi2GammaC<<
1228         "chi2Lambda42c.="<<chi2Lambda42C<<
1229         "chi2Lambda24c.="<<chi2Lambda24C<<
1230         //
1231         "v0K0C.="<<v0K0C<<
1232         "v0GammaC.="<<v0GammaC<<
1233         "v0Lambda42C.="<<v0Lambda42C<<
1234         "v0Lambda24C.="<<v0Lambda24C<<
1235         //
1236         "massK0="<<massK0<<
1237         "massGamma="<<massGamma<<
1238         "massLambda42="<<massLambda42<<
1239         "massLambda24="<<massLambda24<<
1240         //
1241         "timeK0="<<timeK0<<
1242         "timeLambda42="<<timeLambda42<<
1243         "timeLambda24="<<timeLambda24<<
1244         "\n";
1245     }
1246     if (type==1) fGammas->AddLast(v0); 
1247     //
1248     //
1249     //
1250     delete v0K0;
1251     delete v0Gamma;
1252     delete v0Lambda42;
1253     delete v0Lambda24;    
1254     delete v0K0C;
1255     delete v0GammaC;
1256     delete v0Lambda42C;
1257     delete v0Lambda24C; 
1258   }
1259   ProcessPI0(); 
1260 }
1261
1262
1263
1264 void AliTPCcalibV0::ProcessPI0(){
1265   //
1266   //
1267   //
1268   Int_t nentries = fGammas->GetEntries();
1269   if (nentries<2) return;
1270   // 
1271   Double_t m0[3], m1[3];
1272   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
1273   for (Int_t i0=0; i0<nentries; i0++){
1274     AliESDv0 *v00 = (AliESDv0*)fGammas->At(i0); 
1275     v00->GetPxPyPz (m0[0], m0[1], m0[2]);
1276     AliKFParticle *p00 = Fit(primVtx, v00, 11,-11);
1277     p00->SetProductionVertex( primVtx );
1278     p00->SetMassConstraint(0);
1279     //
1280     for (Int_t i1=i0; i1<nentries; i1++){
1281       AliESDv0 *v01 = (AliESDv0*)fGammas->At(i1);
1282       v01->GetPxPyPz (m1[0], m1[1], m1[2]);
1283       AliKFParticle *p01 = Fit(primVtx, v01, 11,-11);
1284       p01->SetProductionVertex( primVtx );
1285       p01->SetMassConstraint(0);
1286       if (v00->GetIndex(0) != v01->GetIndex(0) && 
1287           v00->GetIndex(1) != v01->GetIndex(1)){
1288         AliKFParticle pi0( *p00,*p01); 
1289         pi0.SetProductionVertex(primVtx);
1290         Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]);
1291         Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]);
1292         Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2])));
1293         TTreeSRedirector *cstream = GetDebugStreamer();
1294         if (cstream){
1295           (*cstream)<<"PI0"<<
1296             "v00.="<<v00<<
1297             "v01.="<<v01<<
1298             "mass="<<mass<<
1299             "p00.="<<p00<<
1300             "p01.="<<p01<<
1301             "pi0.="<<&pi0<<
1302             "\n";       
1303         }
1304       }
1305       delete p01;
1306     }
1307     delete p00;
1308   }
1309 }
1310
1311
1312
1313
1314
1315 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
1316   //
1317   // Make KF Particle
1318   //
1319   AliKFParticle p1( *(v0->GetParamN()), PDG1 );
1320   AliKFParticle p2( *(v0->GetParamP()), PDG2 );
1321   AliKFParticle *V0 = new AliKFParticle;
1322   Double_t x, y, z;
1323   v0->GetXYZ(x,y,z );
1324   V0->SetVtxGuess(x,y,z);
1325   *(V0)+=p1;
1326   *(V0)+=p2;
1327   return V0;  
1328 }
1329
1330
1331
1332
1333 void AliTPCcalibV0::BinLogX(TH2F *h) {
1334   //
1335   //
1336   //
1337    TAxis *axis = h->GetXaxis();
1338    int bins = axis->GetNbins();
1339
1340    Double_t from = axis->GetXmin();
1341    Double_t to = axis->GetXmax();
1342    Double_t *new_bins = new Double_t[bins + 1];
1343    
1344    new_bins[0] = from;
1345    Double_t factor = pow(to/from, 1./bins);
1346   
1347    for (int i = 1; i <= bins; i++) {
1348      new_bins[i] = factor * new_bins[i-1];
1349    }
1350    axis->Set(bins, new_bins);
1351    delete [] new_bins;
1352    
1353 }
1354
1355
1356
1357