.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / TPC / TPCcalib / 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 "AliTPCcalibV0.h"
39 #include "AliV0.h"
40 #include "TRandom.h"
41 #include "TTreeStream.h"
42 #include "AliTPCcalibDB.h"
43 #include "AliTPCCorrection.h"
44 #include "AliGRPObject.h"
45 #include "AliTPCTransform.h"
46 #include "AliAnalysisManager.h"
47
48
49
50 ClassImp(AliTPCcalibV0)
51
52
53 AliTPCcalibV0::AliTPCcalibV0() : 
54    AliTPCcalibBase(),
55    fV0Tree(0),
56    fHPTTree(0),
57    fStack(0),
58    fESD(0),
59    fPdg(0),
60    fParticles(0),
61    fV0s(0),
62    fGammas(0)
63 {
64   
65 }   
66 AliTPCcalibV0::AliTPCcalibV0(const Text_t *name, const Text_t *title):
67    AliTPCcalibBase(),
68    fV0Tree(0),
69    fHPTTree(0),
70    fStack(0),
71    fESD(0),
72    fPdg(0),
73    fParticles(0),
74    fV0s(0),
75    fGammas(0)
76 {
77   fPdg = new TDatabasePDG;       
78   // create output histograms
79   SetName(name);
80   SetTitle(title);
81 }   
82
83 AliTPCcalibV0::~AliTPCcalibV0(){
84   //
85   //
86   //
87   delete fV0Tree;
88   delete fHPTTree;
89 }
90
91
92
93
94
95 void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd){
96   //
97   //
98   //
99   fESD = esd;
100   AliKFParticle::SetField(esd->GetMagneticField());
101   if (TMath::Abs(AliTracker::GetBz())<1) return;  
102   DumpToTree(esd);
103   DumpToTreeHPT(esd);
104 }
105
106 void  AliTPCcalibV0::DumpToTreeHPT(AliESDEvent *esd){
107   //
108   // Dump V0s fith full firend information to the 
109   // 
110   if (TMath::Abs(AliTracker::GetBz())<1) return;
111   const Int_t kMinCluster=110;
112   const Float_t kMinPt   =4.;
113   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
114 //   if (!esdFriend) {
115 //     Printf("ERROR: esdFriend not available");
116 //     return;
117 //   }
118   //
119   Int_t ntracks=esd->GetNumberOfTracks();
120   for (Int_t i=0;i<ntracks;++i) {
121     Bool_t isOK=kFALSE;
122     AliESDtrack *track = esd->GetTrack(i);
123     if (track->GetTPCncls()<kMinCluster) continue;
124     if (TMath::Abs(AliTracker::GetBz())>1){ // cut on momenta if measured
125       if (track->Pt()>kMinPt) isOK=kTRUE;
126     }
127     if (TMath::Abs(AliTracker::GetBz())<1){  // require primary track for the B field OFF data
128       Bool_t isAccepted=kTRUE;
129       if (!track->IsOn(AliESDtrack::kITSrefit)) isAccepted=kFALSE;
130       if (!track->IsOn(AliESDtrack::kTPCrefit)) isAccepted=kFALSE;
131       if (!track->IsOn(AliESDtrack::kTOFout))   isAccepted=kFALSE;
132       Float_t dvertex[2],cvertex[3]; 
133       track->GetImpactParametersTPC(dvertex,cvertex);
134       if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>20) isAccepted=kFALSE;
135       if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>20) isAccepted=kFALSE;
136       track->GetImpactParameters(dvertex,cvertex);
137       if (TMath::Abs(dvertex[0]/TMath::Sqrt(cvertex[0]+0.01))>10) isAccepted=kFALSE;
138       if (TMath::Abs(dvertex[1]/TMath::Sqrt(TMath::Abs(cvertex[2]+0.01)))>10) isAccepted=kFALSE;
139       if (!isAccepted) isOK=kFALSE;
140     } 
141     if ( track->GetTPCsignal()>100 && track->GetInnerParam()->Pt()>1 ){
142       if (track->IsOn(AliESDtrack::kITSin)||track->IsOn(AliESDtrack::kTRDout)||track->IsOn(AliESDtrack::kTOFin))
143         isOK=kTRUE;
144       if (isOK){
145         TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
146         Int_t eventNumber = esd->GetEventNumberInFile();
147         Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0; 
148         Bool_t hasITS=(track->GetNcls(0)>2);
149         printf("DUMPIONTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->GetInnerParam()->Pt()*track->GetTPCsignal()/50., eventNumber,hasFriend,hasITS);
150       }
151     }
152     if (!isOK) continue;
153     TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
154     Int_t eventNumber = esd->GetEventNumberInFile();    
155     Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(i)!=0):0;
156     Bool_t hasITS=(track->GetNcls(0)>2);    
157     printf("DUMPHPTTrack:%s|%f|%d|%d|%d\n",filename.Data(),track->Pt(), eventNumber,hasFriend,hasITS);
158     //
159     if (!esdFriend) continue;
160     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
161     if (!friendTrack) continue;
162
163     if (!isOK) continue;
164     //
165
166     TObject *calibObject;
167     AliTPCseed *seed = 0;
168     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
169       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
170     }
171     if (!seed) continue;
172       if (!fHPTTree) {
173       fHPTTree = new TTree("HPT","HPT");
174       fHPTTree->SetDirectory(0);
175     }
176     if (fHPTTree->GetEntries()==0){
177       //
178       fHPTTree->SetDirectory(0);
179       fHPTTree->Branch("t.",&track);
180       fHPTTree->Branch("ft.",&friendTrack);
181       fHPTTree->Branch("s.",&seed);
182     }else{
183       fHPTTree->SetBranchAddress("t.",&track);
184       fHPTTree->SetBranchAddress("ft.",&friendTrack);
185       fHPTTree->SetBranchAddress("s.",&seed);
186     }
187     fHPTTree->Fill();
188     //
189   }
190 }
191
192
193
194 void  AliTPCcalibV0::DumpToTree(AliESDEvent *esd){
195   //
196   // Dump V0s fith full firend information to the 
197   // 
198   Int_t nV0s  = fESD->GetNumberOfV0s();
199   const Int_t kMinCluster=110;
200   const Double_t kDownscale=0.01;
201   const Float_t kMinPt   =1.0;
202   const Float_t kMinMinPt   =0.7;
203   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(esd->FindListObject("AliESDfriend"));
204   //
205   
206   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
207     Bool_t isOK=kFALSE;
208     AliESDv0 * v0 = (AliESDv0*) esd->GetV0(ivertex);
209     AliESDtrack * track0 = fESD->GetTrack(v0->GetIndex(0)); // negative track
210     AliESDtrack * track1 = fESD->GetTrack(v0->GetIndex(1)); // positive track 
211     if (track0->GetTPCNcls()<kMinCluster) continue;
212     if (track0->GetKinkIndex(0)>0) continue;    
213     if (track1->GetTPCNcls()<kMinCluster) continue;
214     if (track1->GetKinkIndex(0)>0) continue;
215     if (v0->GetOnFlyStatus()==kFALSE) continue;
216     //
217     if (TMath::Min(track0->Pt(),track1->Pt())<kMinMinPt) continue;
218     //
219     //
220     if (TMath::Max(track0->Pt(),track1->Pt())>kMinPt) isOK=kTRUE;
221     if (gRandom->Rndm()<kDownscale) isOK=kTRUE;  
222     if (!isOK) continue;
223     //
224     TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
225     Int_t eventNumber = esd->GetEventNumberInFile();
226     Bool_t hasITS=(track0->GetNcls(0)+ track1->GetNcls(0)>4);
227     printf("DUMPHPTV0:%s|%f|%d|%d|%d\n",filename.Data(), (TMath::Min(track0->Pt(),track1->Pt())), eventNumber,(esdFriend!=0), hasITS);
228     //
229     if (!esdFriend) continue;
230     //
231     
232     //
233     AliESDfriendTrack *ftrack0 = esdFriend->GetTrack(v0->GetIndex(0));
234     if (!ftrack0) continue;
235     AliESDfriendTrack *ftrack1 = esdFriend->GetTrack(v0->GetIndex(1));
236     if (!ftrack1) continue;
237     //
238     TObject *calibObject;
239     AliTPCseed *seed0 = 0;
240     AliTPCseed *seed1 = 0;
241     for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
242       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
243     }
244     for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
245       if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
246     }
247     if (!seed0) continue;
248     if (!seed1) continue;
249     AliExternalTrackParam * paramIn0 = (AliExternalTrackParam *)track0->GetInnerParam();
250     AliExternalTrackParam * paramIn1 = (AliExternalTrackParam *)track1->GetInnerParam();
251     if (!paramIn0) continue;
252     if (!paramIn1) continue;
253     //
254     //
255     if (!fV0Tree) {
256       fV0Tree = new TTree("V0s","V0s");
257       fV0Tree->SetDirectory(0);
258     }
259     if (fV0Tree->GetEntries()==0){
260       //
261       fV0Tree->SetDirectory(0);
262       fV0Tree->Branch("v0.",&v0);
263       fV0Tree->Branch("t0.",&track0);
264       fV0Tree->Branch("t1.",&track1);
265       fV0Tree->Branch("ft0.",&ftrack0);
266       fV0Tree->Branch("ft1.",&ftrack1);
267       fV0Tree->Branch("s0.",&seed0);
268       fV0Tree->Branch("s1.",&seed1);
269     }else{
270       fV0Tree->SetBranchAddress("v0.",&v0);
271       fV0Tree->SetBranchAddress("t0.",&track0);
272       fV0Tree->SetBranchAddress("t1.",&track1);
273       fV0Tree->SetBranchAddress("ft0.",&ftrack0);
274       fV0Tree->SetBranchAddress("ft1.",&ftrack1);
275       fV0Tree->SetBranchAddress("s0.",&seed0);
276       fV0Tree->SetBranchAddress("s1.",&seed1);
277     }
278     fV0Tree->Fill();
279   }
280 }
281
282
283 Long64_t AliTPCcalibV0::Merge(TCollection *const li) {
284
285   TIterator* iter = li->MakeIterator();
286   AliTPCcalibV0* cal = 0;
287
288   while ((cal = (AliTPCcalibV0*)iter->Next())) {
289     if (cal->fV0Tree){
290       if (!fV0Tree) {
291         fV0Tree = new TTree("V0s","V0s");
292         fV0Tree->SetDirectory(0);
293       }
294       if (cal->fV0Tree->GetEntries()>0) AliTPCcalibV0::AddTree(cal->fV0Tree);
295       if (cal->fHPTTree->GetEntries()>0) AliTPCcalibV0::AddTreeHPT(cal->fHPTTree);
296     }    
297   }
298   return 0;
299 }
300
301
302 void AliTPCcalibV0::AddTree(TTree * treeInput){
303   //
304   // Add the content of tree: 
305   // Notice automatic copy of tree in ROOT does not work for such complicated tree
306   //  
307   return ;
308   AliESDv0 * v0 = new AliESDv0;
309   Double_t kMinPt=0.8;
310   AliESDtrack * track0 = 0; // negative track
311   AliESDtrack * track1 = 0; // positive track 
312   AliESDfriendTrack *ftrack0 = 0;
313   AliESDfriendTrack *ftrack1 = 0;
314   AliTPCseed *seed0 = 0;
315   AliTPCseed *seed1 = 0;
316   treeInput->SetBranchStatus("ft0.",kFALSE);
317   treeInput->SetBranchStatus("ft1.",kFALSE);
318   TDatabasePDG pdg;
319   Double_t massK0= pdg.GetParticle("K0")->Mass();
320   Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
321
322   Int_t entries= treeInput->GetEntries();
323   for (Int_t i=0; i<entries; i++){
324     treeInput->SetBranchAddress("v0.",&v0);
325     treeInput->SetBranchAddress("t0.",&track0);
326     treeInput->SetBranchAddress("t1.",&track1);
327     treeInput->SetBranchAddress("ft0.",&ftrack0);
328     treeInput->SetBranchAddress("ft1.",&ftrack1);
329     treeInput->SetBranchAddress("s0.",&seed0);
330     treeInput->SetBranchAddress("s1.",&seed1);
331     if (fV0Tree->GetEntries()==0){
332       fV0Tree->SetDirectory(0);
333       fV0Tree->Branch("v0.",&v0);
334       fV0Tree->Branch("t0.",&track0);
335       fV0Tree->Branch("t1.",&track1);
336       fV0Tree->Branch("ft0.",&ftrack0);
337       fV0Tree->Branch("ft1.",&ftrack1);
338       fV0Tree->Branch("s0.",&seed0);
339       fV0Tree->Branch("s1.",&seed1);
340     }else{
341       fV0Tree->SetBranchAddress("v0.",&v0);
342       fV0Tree->SetBranchAddress("t0.",&track0);
343       fV0Tree->SetBranchAddress("t1.",&track1);
344       fV0Tree->SetBranchAddress("ft0.",&ftrack0);
345       fV0Tree->SetBranchAddress("ft1.",&ftrack1);
346       fV0Tree->SetBranchAddress("s0.",&seed0);
347       fV0Tree->SetBranchAddress("s1.",&seed1);
348     }
349     //
350     treeInput->GetEntry(i);
351     //ftrack0->GetCalibContainer()->SetOwner(kTRUE);
352     //ftrack1->GetCalibContainer()->SetOwner(kTRUE);
353     Bool_t isOK=kTRUE;
354     if (v0->GetOnFlyStatus()==kFALSE) isOK=kFALSE;
355     if (track0->GetTPCncls()<100) isOK=kFALSE;
356     if (track1->GetTPCncls()<100) isOK=kFALSE;    
357     if (TMath::Min(seed0->Pt(),seed1->Pt())<kMinPt) isOK=kFALSE;
358     if (TMath::Min(track0->Pt(),track1->Pt())<kMinPt) isOK=kFALSE;
359     Bool_t isV0=kFALSE;    
360     if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.05)     isV0=kTRUE;
361     if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.05) isV0=kTRUE; 
362     if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.05) isV0=kTRUE;
363     if (TMath::Abs(v0->GetEffMass(0,0))<0.02) isV0=kFALSE; //reject electrons
364     if (!isV0) isOK=kFALSE;
365     if (isOK) fV0Tree->Fill();
366     delete v0;
367     delete track0;
368     delete track1;
369     delete ftrack0;
370     delete ftrack1;
371     delete seed0;
372     delete seed1;
373     v0=0;
374     track0=0;
375     track1=0;
376     ftrack0=0;
377     ftrack1=0;
378     seed0=0;
379     seed1=0;
380   }
381 }
382
383 void AliTPCcalibV0::AddTreeHPT(TTree * treeInput){
384   //
385   // Add the content of tree: 
386   // Notice automatic copy of tree in ROOT does not work for such complicated tree
387   //  
388   return ;
389   AliESDtrack *track = 0;
390   AliESDfriendTrack *friendTrack = 0;
391   AliTPCseed *seed = 0;
392   if (!treeInput) return;
393   if (treeInput->GetEntries()==0) return;
394   //
395   Int_t entries= treeInput->GetEntries();  
396   //
397   for (Int_t i=0; i<entries; i++){
398     track=0;
399     friendTrack=0;
400     seed=0;
401     //
402     treeInput->SetBranchAddress("t.",&track);
403     treeInput->SetBranchAddress("ft.",&friendTrack);
404     treeInput->SetBranchAddress("s.",&seed);
405     treeInput->GetEntry(i);
406     //
407     if (fHPTTree->GetEntries()==0){
408       fHPTTree->SetDirectory(0);
409       fHPTTree->Branch("t.",&track);
410       fHPTTree->Branch("ft.",&friendTrack);
411       fHPTTree->Branch("s.",&seed);
412     }else{
413       fHPTTree->SetBranchAddress("t.",&track);
414       fHPTTree->SetBranchAddress("ft.",&friendTrack);
415       fHPTTree->SetBranchAddress("s.",&seed);
416     }    
417     Bool_t isOK=kTRUE;
418     if (!track->IsOn(AliESDtrack::kITSrefit)) isOK=kFALSE;
419     if (!track->IsOn(AliESDtrack::kTOFout)) isOK=kFALSE;
420     if (isOK) fHPTTree->Fill();
421     //
422     delete track;
423     delete friendTrack;
424     delete seed;
425   }
426 }
427
428
429 void AliTPCcalibV0::MakeFitTreeTrack(const TObjArray * corrArray, Double_t ptCut, Int_t /*run*/){
430   //
431   // Make a fit tree
432   //
433   // 0. Loop over selected tracks
434   // 1. Loop over all transformation - refit the track with and without the
435   //    transformtation
436   // 2. Dump the matching paramaeters to the debugStremer
437   //
438   
439   //Connect input
440   const Int_t kMinNcl=120;
441   TFile f("TPCV0Objects.root");
442   AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
443   TTree * treeInput = v0TPC->GetHPTTree();
444   TTreeSRedirector *pcstream = new TTreeSRedirector("fitHPT.root");
445   AliESDtrack *track = 0;
446   AliESDfriendTrack *friendTrack = 0;
447   AliTPCseed *seed = 0;
448   if (!treeInput) return;
449   if (treeInput->GetEntries()==0) return;
450   //
451   treeInput->SetBranchAddress("t.",&track);
452   treeInput->SetBranchAddress("ft.",&friendTrack);
453   treeInput->SetBranchAddress("s.",&seed);
454   //
455   Int_t ncorr=0;
456   if (corrArray) ncorr = corrArray->GetEntries();
457   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
458  //  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
459 //   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
460 //   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd());
461   //
462   //
463   //  
464   Int_t ntracks= treeInput->GetEntries();
465   for (Int_t itrack=0; itrack<ntracks; itrack++){
466     treeInput->GetEntry(itrack);
467     if (!track) continue;
468     if (seed->Pt()<ptCut) continue;
469     if (track->Pt()<ptCut) continue;
470     if (track->GetTPCncls()<kMinNcl) continue;
471     //
472     // Reapply transformation
473     //
474     for (Int_t irow=0; irow<159; irow++){
475       AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
476       if (cluster &&cluster->GetX()>10){
477         Double_t x0[3]={ static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
478         Int_t index0[1]={cluster->GetDetector()};
479         transform->Transform(x0,index0,0,1);
480         cluster->SetX(x0[0]);
481         cluster->SetY(x0[1]);
482         cluster->SetZ(x0[2]);
483         //
484       }
485     }    
486     //
487     AliExternalTrackParam* paramInner=0;
488     AliExternalTrackParam* paramOuter=0;
489     AliExternalTrackParam* paramIO=0;
490     Bool_t isOK=kTRUE;
491     for (Int_t icorr=-1; icorr<ncorr; icorr++){
492       //
493       AliTPCCorrection *corr = 0;
494       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
495       AliExternalTrackParam * trackInner = RefitTrack(seed, corr,85,134,0.1);      
496       AliExternalTrackParam * trackIO = RefitTrack(seed, corr,245,85,0.1);      
497       AliExternalTrackParam * trackOuter = RefitTrack(seed, corr,245,134,0.1 ); 
498       if (trackInner&&trackOuter&&trackIO){
499         trackOuter->Rotate(trackInner->GetAlpha());
500         trackOuter->PropagateTo(trackInner->GetX(),AliTracker::GetBz());
501         if (icorr<0) {
502           paramInner=trackInner;
503           paramOuter=trackOuter;
504           paramIO=trackIO;
505           paramIO->Rotate(seed->GetAlpha());
506           paramIO->PropagateTo(seed->GetX(),AliTracker::GetBz());
507         }
508       }else{
509         isOK=kFALSE;
510       }
511       
512     }
513     if (paramOuter&& paramInner) {
514       //      Bool_t isOK=kTRUE;
515       if (paramInner->GetSigmaY2()>0.01) isOK&=kFALSE;
516       if (paramOuter->GetSigmaY2()>0.01) isOK&=kFALSE;
517       if (paramInner->GetSigmaZ2()>0.01) isOK&=kFALSE;
518       if (paramOuter->GetSigmaZ2()>0.01) isOK&=kFALSE;      
519       (*pcstream)<<"fit"<<
520         "s.="<<seed<<
521         "io.="<<paramIO<<
522         "pIn.="<<paramInner<<
523         "pOut.="<<paramOuter;      
524     }
525     //
526   }
527   delete pcstream;
528   /*
529     .x ~/rootlogon.C
530     Int_t run=117112;
531     .x ../ConfigCalibTrain.C(run)
532     .L ../AddTaskTPCCalib.C
533     ConfigOCDB(run)
534     TFile fexb("../../RegisterCorrectionExB.root");
535     AliTPCComposedCorrection *cc=  (AliTPCComposedCorrection*) fexb.Get("ComposedExB");
536     cc->Init();
537     cc->Print("DA"); // Print used correction classes
538     TObjArray *array = cc->GetCorrections()
539     AliTPCcalibV0::MakeFitTreeTrack(array,2,run);
540    
541    */
542 }
543
544 void AliTPCcalibV0::MakeFitTreeV0(const TObjArray * corrArray, Double_t ptCut, Int_t run){
545   //
546   // Make a fit tree
547   //
548   // 0. Loop over selected tracks
549   // 1. Loop over all transformation - refit the track with and without the
550   //    transformtation
551   // 2. Dump the matching paramaeters to the debugStremer
552   //
553   
554   //Connect input
555   TFile f("TPCV0Objects.root");
556   AliTPCcalibV0 *v0TPC = (AliTPCcalibV0*) f.Get("v0TPC");
557   TTree * treeInput = v0TPC->GetV0Tree();
558   TTreeSRedirector *pcstream = new TTreeSRedirector("fitV0.root");
559   AliESDv0 *v0 = 0;
560   AliESDtrack *track0 = 0;
561   AliESDfriendTrack *friendTrack0 = 0;
562   AliTPCseed *seed0 = 0;
563   AliTPCseed *s0 = 0;
564   AliESDtrack *track1 = 0;
565   AliESDfriendTrack *friendTrack1 = 0;
566   AliTPCseed *s1 = 0;
567   AliTPCseed *seed1 = 0;
568   if (!treeInput) return;
569   if (treeInput->GetEntries()==0) return;
570   //
571   treeInput->SetBranchAddress("v0.",&v0);
572   treeInput->SetBranchAddress("t0.",&track0);
573   treeInput->SetBranchAddress("ft0.",&friendTrack0);
574   treeInput->SetBranchAddress("s0.",&s0);
575   treeInput->SetBranchAddress("t1.",&track1);
576   treeInput->SetBranchAddress("ft1.",&friendTrack1);
577   treeInput->SetBranchAddress("s1.",&s1);
578   //
579   TDatabasePDG pdg;
580   Int_t ncorr=0;
581   if (corrArray) ncorr = corrArray->GetEntries();
582   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
583   Double_t massK0= pdg.GetParticle("K0")->Mass();
584   Double_t massLambda= pdg.GetParticle("Lambda0")->Mass();
585   Double_t massPion=pdg.GetParticle("pi+")->Mass();
586   Double_t massProton=pdg.GetParticle("proton")->Mass();
587   Int_t pdgPion=pdg.GetParticle("pi+")->PdgCode();
588   Int_t pdgProton=pdg.GetParticle("proton")->PdgCode();
589   Double_t rmass0=0;
590   Double_t rmass1=0;
591   Double_t massV0=0;
592   Int_t    pdg0=0;
593   Int_t    pdg1=0;
594   //
595   //
596   //  
597   Int_t nv0s= treeInput->GetEntries();
598   for (Int_t iv0=0; iv0<nv0s; iv0++){
599     Int_t  v0Type=0;
600     Int_t isK0=0;
601     Int_t isLambda=0;
602     Int_t isAntiLambda=0;
603     treeInput->GetEntry(iv0);
604     if (TMath::Abs(v0->GetEffMass(2,2)-massK0)<0.03) {isK0=1; v0Type=1;} //select K0s    
605     if (TMath::Abs(v0->GetEffMass(4,2)-massLambda)<0.01) {isLambda=1; v0Type=2;} //select Lambda   
606     if (TMath::Abs(v0->GetEffMass(2,4)-massLambda)<0.01) {isAntiLambda=1;v0Type=3;} //select Anti Lambda
607     if (isK0+isLambda+isAntiLambda!=1) continue;
608     rmass0=massPion;
609     rmass1=massPion;
610     pdg0=pdgPion;
611     pdg1=pdgPion;
612     if (isLambda) {rmass0=massProton; pdg0=pdgProton;}
613     if (isAntiLambda) {rmass1=massProton; pdg1=pdgProton;}
614     massV0=massK0;
615     if (isK0==0) massV0=massLambda;
616     //
617     if (!s0) continue;
618     seed0=(s0->GetSign()>0)?s0:s1;
619     seed1=(s0->GetSign()>0)?s1:s0;
620     if (seed0->GetZ()*seed1->GetZ()<0) continue; //remove membrane crossed tracks
621     if (seed0->Pt()<ptCut) continue;
622     if (seed1->Pt()<ptCut) continue;
623     //
624     // Reapply transformation
625     //
626     for  (Int_t itype=0; itype<2; itype++){
627       AliTPCseed * seed = (itype==0) ? seed0: seed1;      
628       for (Int_t irow=0; irow<159; irow++){
629         AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
630         if (cluster &&cluster->GetX()>10){
631           Double_t x0[3]={ static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
632           Int_t index0[1]={cluster->GetDetector()};
633           transform->Transform(x0,index0,0,1);
634           cluster->SetX(x0[0]);
635           cluster->SetY(x0[1]);
636           cluster->SetZ(x0[2]);
637           //
638         }
639       }
640     }   
641     Bool_t isOK=kTRUE;
642     Double_t radius = v0->GetRr();
643     Double_t xyz[3];
644     v0->GetXYZ(xyz[0],xyz[1],xyz[2]);
645     Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
646     TObjArray arrayV0in(ncorr+1);
647     TObjArray arrayV0io(ncorr+1);
648     TObjArray arrayT0(ncorr+1);
649     TObjArray arrayT1(ncorr+1);
650     arrayV0in.SetOwner(kTRUE);
651     arrayV0io.SetOwner(kTRUE);
652     //
653     for (Int_t icorr=-1; icorr<ncorr; icorr++){
654       AliTPCCorrection *corr =0;
655       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
656       //
657       AliExternalTrackParam * trackInner0 = RefitTrack(seed0, corr,160,85,rmass0);      
658       AliExternalTrackParam * trackIO0    = RefitTrack(seed0, corr,245,85,rmass0);      
659       AliExternalTrackParam * trackInner1 = RefitTrack(seed1, corr,160,85,rmass1);      
660       AliExternalTrackParam * trackIO1    = RefitTrack(seed1, corr,245,85,rmass1);      
661       if (!trackInner0) isOK=kFALSE;
662       if (!trackInner1) isOK=kFALSE;
663       if (!trackIO0)    isOK=kFALSE;
664       if (!trackIO1)    isOK=kFALSE;
665       if (isOK){
666         if (!trackInner0->Rotate(alpha)) isOK=kFALSE;
667         if (!trackInner1->Rotate(alpha)) isOK=kFALSE;
668         if (!trackIO0->Rotate(alpha)) isOK=kFALSE;
669         if (!trackIO1->Rotate(alpha)) isOK=kFALSE;
670         //
671         if (!AliTracker::PropagateTrackToBxByBz(trackInner0, radius, rmass0, 1, kFALSE)) isOK=kFALSE; 
672         if (!AliTracker::PropagateTrackToBxByBz(trackInner1, radius, rmass1, 1, kFALSE)) isOK=kFALSE; 
673         if (!AliTracker::PropagateTrackToBxByBz(trackIO0, radius, rmass0, 1, kFALSE)) isOK=kFALSE; 
674         if (!AliTracker::PropagateTrackToBxByBz(trackIO1, radius, rmass1, 1, kFALSE)) isOK=kFALSE; 
675         if (!isOK) continue;
676         arrayT0.AddAt(trackIO0->Clone(),icorr+1);
677         arrayT1.AddAt(trackIO1->Clone(),icorr+1);
678         Int_t charge=TMath::Nint(trackIO0->GetSign());
679         AliKFParticle pin0( *trackInner0,  pdg0*charge);
680         AliKFParticle pin1( *trackInner1, -pdg1*charge);
681         AliKFParticle pio0( *trackIO0,  pdg0*charge);
682         AliKFParticle pio1( *trackIO1, -pdg1*charge);
683         AliKFParticle v0in;
684         AliKFParticle v0io;
685         v0in+=pin0;
686         v0in+=pin1;
687         v0io+=pio0;
688         v0io+=pio1;
689         arrayV0in.AddAt(v0in.Clone(),icorr+1);
690         arrayV0io.AddAt(v0io.Clone(),icorr+1);
691       }
692     }
693     if (!isOK) continue;
694     //
695     //AliKFParticle* pin0= (AliKFParticle*)arrayV0in.At(0);
696     AliKFParticle* pio0= (AliKFParticle*)arrayV0io.At(0);
697     AliExternalTrackParam *param0=(AliExternalTrackParam *)arrayT0.At(0);
698     AliExternalTrackParam *param1=(AliExternalTrackParam *)arrayT1.At(0);
699     Double_t mass0=0, mass0E=0; 
700     pio0->GetMass( mass0,mass0E);
701     //
702     Double_t mean=mass0-massV0;
703     if (TMath::Abs(mean)>0.05) continue;
704     Double_t mass[10000];
705     //
706     Int_t dtype=30;  // id for V0
707     Int_t ptype=5;   // id for invariant mass
708     //    Int_t id=TMath::Nint(100.*(param0->Pt()-param1->Pt())/(param0->Pt()+param1->Pt()));      // K0s V0 asymetry
709     Int_t id=Int_t(1000.*(param0->Pt()-param1->Pt()));      // K0s V0 asymetry
710     Double_t gx,gy,gz, px,py,pz;
711     Double_t pt = v0->Pt();
712     v0->GetXYZ(gx,gy,gz);
713     v0->GetPxPyPz(px,py,pz);
714     Double_t theta=pz/TMath::Sqrt(px*px+py*py);
715     Double_t phi=TMath::ATan2(py,px);
716     Double_t snp=0.5*(seed0->GetSnp()+seed1->GetSnp());
717     Double_t sector=9*phi/TMath::Pi();
718     Double_t dsec=sector-TMath::Nint(sector);
719     Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
720     //Int_t nentries=v0Type;
721     Double_t bz=AliTracker::GetBz();
722     Double_t dRrec=0;
723     (*pcstream)<<"fitDebug"<< 
724       "id="<<id<<
725       "v0.="<<v0<<
726       "mean="<<mean<<
727       "rms="<<mass0E<<
728       "pio0.="<<pio0<<
729       "p0.="<<param0<<
730       "p1.="<<param1;
731     (*pcstream)<<"fit"<<  // dump valus for fit
732       "run="<<run<<       //run number
733       "bz="<<bz<<         // magnetic filed used
734       "dtype="<<dtype<<   // detector match type 30
735       "ptype="<<ptype<<   // parameter type
736       "theta="<<theta<<   // theta
737       "phi="<<phi<<       // phi 
738       "snp="<<snp<<       // snp
739       "mean="<<mean<<     // mean dist value
740       "rms="<<mass0E<<       // rms
741       "sector="<<sector<<
742       "dsec="<<dsec<<
743       //
744       "refX="<<refX<<      // reference radius
745       "gx="<<gx<<         // global position
746       "gy="<<gy<<         // global position
747       "gz="<<gz<<         // global position
748       "dRrec="<<dRrec<<      // delta Radius in reconstruction
749       "pt="<<pt<<         // pt of the particle
750       "id="<<id<<     //delta of the momenta      
751       "entries="<<v0Type;//  type of the V0
752     for (Int_t icorr=0; icorr<ncorr; icorr++){
753       AliTPCCorrection *corr =0;
754       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
755       //      AliKFParticle* pin= (AliKFParticle*)arrayV0in.At(icorr+1);
756       AliKFParticle* pio= (AliKFParticle*)arrayV0io.At(icorr+1);
757       AliExternalTrackParam *par0=(AliExternalTrackParam *)arrayT0.At(icorr+1);
758       AliExternalTrackParam *par1=(AliExternalTrackParam *)arrayT1.At(icorr+1);
759       Double_t massE=0; 
760       pio->GetMass( mass[icorr],massE);
761       mass[icorr]-=mass0;
762       (*pcstream)<<"fit"<< 
763         Form("%s=",corr->GetName())<<mass[icorr];
764       (*pcstream)<<"fitDebug"<< 
765         Form("%s=",corr->GetName())<<mass[icorr]<<
766         Form("%sp0.=",corr->GetName())<<par0<<
767         Form("%sp1=",corr->GetName())<<par1;
768     }
769     (*pcstream)<<"fit"<< "isOK="<<isOK<<"\n";        
770     (*pcstream)<<"fitDebug"<< "isOK="<<isOK<<"\n";        
771   }
772   delete pcstream;
773 }
774
775 AliExternalTrackParam * AliTPCcalibV0::RefitTrack(AliTPCseed *seed, AliTPCCorrection * corr, Double_t xstart, Double_t xstop, Double_t mass){
776   //
777   // Refit the track:
778   //    seed   - tpc track with cluster
779   //    corr   - distrotion/correction class  - apllied to the points
780   //    xstart - radius to start propagate/update
781   //    xstop  - radius to stop propagate/update
782   // 
783   const Double_t kResetCov=20.;
784   const Double_t kSigma=5.;
785   Double_t covar[15];
786   for (Int_t i=0;i<15;i++) covar[i]=0;
787   covar[0]=kSigma*kSigma;
788   covar[2]=kSigma*kSigma;
789   covar[5]=kSigma*kSigma/Float_t(150*150);
790   covar[9]=kSigma*kSigma/Float_t(150*150);
791   covar[14]=1*1;
792   // 
793   AliExternalTrackParam *refit  = new AliExternalTrackParam(*seed);
794   refit->PropagateTo(xstart, AliTracker::GetBz());
795   refit->AddCovariance(covar);
796   refit->ResetCovariance(kResetCov);
797   Double_t xmin = TMath::Min(xstart,xstop);
798   Double_t xmax = TMath::Max(xstart,xstop);
799   Int_t ncl=0;
800   //
801   Bool_t isOK=kTRUE;
802   for (Int_t index0=0; index0<159; index0++){
803     Int_t irow= (xstart<xstop)? index0:159-index0;
804     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);  //cluster in local system
805     if (!cluster) continue;
806     if (cluster->GetX()<xmin) continue;
807     if (cluster->GetX()>xmax) continue;
808     Double_t alpha = TMath::Pi()*(cluster->GetDetector()+0.5)/9.;
809     if (!refit->Rotate(alpha)) isOK=kFALSE;
810     Double_t x     = cluster->GetX();
811     Double_t y     = cluster->GetY();
812     Double_t z     = cluster->GetZ();
813     if (corr){      
814       Float_t xyz[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};  // original position
815       corr->DistortPointLocal(xyz,cluster->GetDetector());
816       x=xyz[0];
817       y=xyz[1];
818       z=xyz[2];
819     }
820     if (!AliTracker::PropagateTrackToBxByBz(refit, x,mass,1.,kFALSE)) isOK=kFALSE;
821     if (!isOK) continue;
822     Double_t cov[3]={0.01,0.,0.01};
823     Double_t yz[2]={y,z};
824     if (!refit->Update(yz,cov)) isOK=kFALSE;    
825     ncl++;
826   }
827   if (!AliTracker::PropagateTrackToBxByBz(refit, xstop, mass,1.,kTRUE)) isOK=kFALSE;
828   //  
829   if (!isOK) {
830     delete refit;
831     return 0;
832   }  
833   return refit;
834 }
835
836
837
838
839
840 //
841 // Obsolete part
842 //
843
844
845
846
847 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
848   //
849   // Make KF Particle
850   //
851   AliKFParticle p1( *(v0->GetParamN()), PDG1 );
852   AliKFParticle p2( *(v0->GetParamP()), PDG2 );
853   AliKFParticle *V0 = new AliKFParticle;
854   Double_t x, y, z;
855   v0->GetXYZ(x,y,z );
856   V0->SetVtxGuess(x,y,z);
857   *(V0)+=p1;
858   *(V0)+=p2;
859   return V0;  
860 }
861
862
863
864
865 void AliTPCcalibV0::BinLogX(TH2F *h) {
866   //
867   //
868   //
869    TAxis *axis = h->GetXaxis();
870    int bins = axis->GetNbins();
871
872    Double_t from = axis->GetXmin();
873    Double_t to = axis->GetXmax();
874    Double_t *new_bins = new Double_t[bins + 1];
875    
876    new_bins[0] = from;
877    Double_t factor = pow(to/from, 1./bins);
878   
879    for (int i = 1; i <= bins; i++) {
880      new_bins[i] = factor * new_bins[i-1];
881    }
882    axis->Set(bins, new_bins);
883    delete [] new_bins;   
884 }
885
886
887
888 void AliTPCcalibV0::FilterV0s(AliESDEvent* event){
889   //
890   // 
891   TDatabasePDG pdg;  
892   const Double_t kChi2Cut=20;
893   const Double_t kMinR=2;
894   const Double_t ptCut=0.2;
895   const Int_t kMinNcl=110;
896   //
897   Int_t nv0 = event->GetNumberOfV0s(); 
898   AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex();
899   AliKFVertex kfvertex=*vertex;
900   //
901   for (Int_t iv0=0;iv0<nv0;iv0++){
902     AliESDv0 *v0 = event->GetV0(iv0);
903     if (!v0) continue;
904     if (v0->GetPindex()<0) continue;
905     if (v0->GetNindex()<0) continue;
906     if (TMath::Max(v0->GetPindex(), v0->GetNindex())>event->GetNumberOfTracks()) continue;
907     //
908     //   
909     AliExternalTrackParam pp=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamP())):(*(v0->GetParamN()));
910     AliExternalTrackParam pn=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamN())):(*(v0->GetParamP()));
911     AliKFParticle kfp1( pp, 211 );
912     AliKFParticle kfp2( pn, -211 );
913     AliKFParticle *v0KFK0 = new AliKFParticle(kfp1,kfp2);
914     AliKFParticle *v0KFK0CV = new AliKFParticle(*v0KFK0);
915     v0KFK0CV->SetProductionVertex(kfvertex);
916     v0KFK0CV->TransportToProductionVertex();
917     AliKFParticle *v0KFK0CVM = new AliKFParticle(*v0KFK0CV);
918     v0KFK0CVM->SetMassConstraint(pdg.GetParticle("K_S0")->Mass());
919     Double_t chi2K0 = v0KFK0CV->GetChi2();
920     //    Double_t chi2K0M= v0KFK0CVM->GetChi2();    
921     //Double_t massK0=v0KFK0CV->GetMass();
922     if (chi2K0>kChi2Cut) continue;
923     if (v0->GetRr()>kMinR) continue;
924     //
925     Double_t maxPt = TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
926     Double_t effMass22=v0->GetEffMass(2,2);
927     Double_t effMass42=v0->GetEffMass(4,2);
928     Double_t effMass24=v0->GetEffMass(2,4);
929     Bool_t isV0= kFALSE;
930     isV0|=TMath::Abs(effMass22-pdg.GetParticle("K_S0")->Mass())<0.1;
931     isV0|=TMath::Abs(effMass42-pdg.GetParticle("Lambda0")->Mass())<0.1;
932     isV0|=TMath::Abs(effMass24-pdg.GetParticle("Lambda0")->Mass())<0.1;
933
934     Double_t sign= v0->GetParamP()->GetSign()* v0->GetParamN()->GetSign();
935     if (sign<0&&v0->GetOnFlyStatus()>0.5&&maxPt>ptCut&&isV0){
936       AliESDtrack * trackP = event->GetTrack(v0->GetPindex());
937       AliESDtrack * trackN = event->GetTrack(v0->GetNindex());
938       if (!trackN) continue;
939       if (!trackP) continue;
940       Int_t nclP= (Int_t)trackP->GetTPCClusterInfo(2,1);
941       Int_t nclN= (Int_t)trackN->GetTPCClusterInfo(2,1);
942       if (TMath::Min(nclP,nclN)<kMinNcl) continue;
943       Double_t eta = TMath::Max(TMath::Abs(trackP->Eta()), TMath::Abs(trackN->Eta()));
944       Double_t ncls = TMath::Min(TMath::Abs(trackP->GetNcls(0)), TMath::Abs(trackN->GetNcls(0)));
945       if (eta<0.8&&ncls>2){
946         //      printf("%d\t%f\t%f\t%d\t%d\t%f\t%f\n",i, v0->Pt(), maxPt, v0->GetNindex(),v0->GetPindex(),v0->GetRr(),effMass22);       
947         (*fDebugStreamer)<<"v0tree"<<
948           "v0.="<<v0<<
949           "tp.="<<trackP<<
950           "tm.="<<trackN<<
951           //
952           "v.="<<vertex<<
953           "ncls="<<ncls<<
954           "maxPt="<<maxPt<<
955           "\n";        
956       }
957     }
958   }
959 }