.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / TPC / TPCcalib / AliTPCCorrectionFit.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17
18 /*
19   Responsible: marian.ivanov@cern.ch 
20   Tools for fitting of the space point distortion parameters.
21   Functionality  
22   
23  
24 1. Creation of the distortion maps from the residual histograms
25
26 2. Making fit trees
27
28 3. Space point distortion not directly observable. Instead a derived variables
29    like DCA at vertex, local y distortion in the TPC-*TOF,TRD,ITS) matching
30    in all 5 tracking parameters are obsereved.
31    In the AliTPCcorrection fir code we calculate the derivative of given variables
32    dO_{i}/dp_{i}
33
34 4. Global fit - later
35    d0 = sum{ki*dO_{i}/dp_{i}}  - linear fitting of the amplitudes ki
36
37  Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be 
38  extracted from the old macros 
39
40
41 */
42
43 #include "Riostream.h"
44 #include <fstream>
45 #include "TMap.h"
46 #include "TGraphErrors.h"
47 #include "AliExternalTrackParam.h"
48 #include "TFile.h"
49 #include "TGraph.h"
50 #include "TMultiGraph.h"
51 #include "TCanvas.h"
52 #include "THnSparse.h"
53 #include "THnBase.h"
54 #include "TProfile.h"
55 #include "TROOT.h"
56 #include "TLegend.h"
57 #include "TPad.h"
58 #include "TH2D.h"
59 #include "TH3D.h"
60 #include "AliTPCROC.h"
61 #include "AliTPCCalROC.h"
62 #include "AliESDfriend.h"
63 #include "AliTPCcalibTime.h"
64 #include "AliSplineFit.h"
65 #include "AliCDBMetaData.h"
66 #include "AliCDBId.h"
67 #include "AliCDBManager.h"
68 #include "AliCDBStorage.h"
69 #include "AliTPCcalibBase.h"
70 #include "AliTPCcalibDB.h"
71 #include "AliTPCcalibDButil.h"
72 #include "AliRelAlignerKalman.h"
73 #include "AliTPCParamSR.h"
74 #include "AliTPCcalibTimeGain.h"
75 #include "AliTPCcalibGainMult.h"
76 #include "AliSplineFit.h"
77 #include "AliTPCComposedCorrection.h"
78 #include "AliTPCExBTwist.h"
79 #include "AliTPCCalibGlobalMisalignment.h"
80 #include "TStatToolkit.h"
81 #include "TChain.h"
82 #include "TCut.h"
83 #include "AliTrackerBase.h"
84 #include "AliTPCCorrectionFit.h"
85 #include "AliTPCLaserTrack.h"
86 #include "TDatabasePDG.h"
87 #include "AliTPCcalibAlign.h"
88 #include "AliLog.h"
89 #include "AliRieman.h"
90
91 ClassImp(AliTPCCorrectionFit)
92
93 Double_t  AliTPCCorrectionFit::fgMaxChi2HelixAt=1;
94
95 AliTPCCorrectionFit::AliTPCCorrectionFit():
96   TNamed("TPCCorrectionFit","TPCCorrectionFit")
97 {
98   //
99   // default constructor
100   //
101 }
102
103 AliTPCCorrectionFit::~AliTPCCorrectionFit() {
104   //
105   // Destructor
106   //
107 }
108
109 /*
110 Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Float_t wt, Float_t t1, Float_t t2){
111   //
112   // Evaluation distortion at point using the linear approximation
113   //
114   //  refX  -  reference X where phi and snp is calculated
115   //  evalX  - ref X where the distortion is evaluated
116
117   if (wt<50){
118     AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);      
119     pcorr->SetOmegaTauT1T2(wt,t1,t2);
120   }
121   Double_t sector = 9*phi/TMath::Pi();
122   if (sector<0) sector+=18;
123   Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
124   Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
125   if (ptype==0) return y85+(y245-y85)*(evalX-85.)/(245.-85.);
126   if (ptype==2) return (y245-y85)/(245.-85.);
127   return 0;
128 }
129 */
130
131
132
133 Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps,Float_t wt, Float_t t1, Float_t t2, Bool_t *pstatus){
134   //
135   // Fit the distortion along the line with the parabolic model
136   // We assume that the track are primaries  - where the vertex is at (0,0,0)
137   //
138   // Parameters:
139   //  phi0  -  phi at the entrance of the TPC
140   //  snp   -  local inclination angle at the entrance of the TPC
141   //  refX  -  reference X where phi and snp is calculated
142   //  evalX  - ref X where the distortion is evaluated
143   //  theta -  inclination angle
144   //  corr  -  internal number of the correction and distortion 
145   //  ptype -  0 - local y distortion
146   //           //1 - local z distortion
147   //           2 - local derivative distortion
148   //           //3
149   //           4 - distortion in the curvature
150   //  nsteps - number of fit points
151   //
152   // return value -  distortion at point evalX with type ptype
153   //
154   static TLinearFitter fitter(3,"pol2"); 
155   fitter.ClearPoints();
156   if (nsteps<3) nsteps=3;
157   AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);        
158   if (pcorr==0) return 0;
159   if (wt<50){
160     if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
161   }
162   Double_t deltaX=(245-85)/(nsteps);
163   Double_t curv=2.*snp/refX;
164   Double_t dPhi0=TMath::ASin(snp);
165
166   for (Int_t istep=0; istep<(nsteps+1); istep++){
167     //
168     Double_t localX =85.+deltaX*istep;
169     Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
170     Double_t localPhi=phi0+dPhi;
171     Double_t sector = 9*localPhi/TMath::Pi();
172     if (sector<0) sector+=18;
173     if (sector>18) sector-=18;
174     Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
175     Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
176     Double_t x[1]={localX-dlocalX};
177     fitter.AddPoint(x,dy);
178   }
179   fitter.Eval();
180   Double_t par[3];
181   par[0]=fitter.GetParameter(0);
182   par[1]=fitter.GetParameter(1);
183   par[2]=fitter.GetParameter(2);
184   
185   Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
186   if (schi2> fgMaxChi2HelixAt){
187     ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("%s\tbad chi2:\t%2.2f",pcorr->GetName(),schi2).Data());
188     if (pstatus) (*pstatus)=kFALSE;    
189     return 0;
190   }
191   if (pstatus) (*pstatus)=kTRUE;
192   if (ptype==0) return par[0]+par[1]*evalX+par[2]*evalX*evalX;
193   if (ptype==2) return par[1]+2*par[2]*evalX;
194   if (ptype==4) return par[2];
195   if (ptype==5) return schi2;
196   return 0;
197 }
198
199
200 Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps, Float_t wt, Float_t t1, Float_t t2, Bool_t *pstatus){
201   //
202   // Fit the distortion along the line with the helix model
203   // FIXME - original trajectory to be changed - AliHelix to be used
204   // We assume that the track are primaries  - where the vertex is at (0,0,0)
205   //
206   // This procedure should emulate distortions as calibrated in the AliTPCcalibTime class
207   //     a.) AliTPCcalibTime::FillResHistoTPCTOF                     -  residuals, phi0 and snp  evaluated at reference refX  TOF clusters (refX=evalX)
208   //     b.) AliTPCcalibTime::FillResHistoTPCTRD                     -  residuals, phi0 and snp  evaluated at reference refX  TPCout radius (refX=evalX)
209   //     c.) AliTPCcalibTime::FillResHistoTPCITS                     -  residuals, phi0 and snp  evaluated at reference refX  TPCin radius (refX=evalX)
210   //     d.) AliTPCcalibTime::FillResHistoTPC  (vertex)              -  phi0 and snp  evaluated at reference refX  TPCin radius, residuals at vertex (refX!=evalX)
211   //         *         
212   // Parameters:
213   //  phi0  -  phi at the entrance of the TPC
214   //  snp   -  local inclination angle at the entrance of the TPC
215   //  refX  -  reference X where phi and snp is calculated
216   //  evalX  - ref X where the distortion is evaluated
217   //  theta -  inclination angle
218   //  corr  -  internal number of the correction and distortion 
219   //  ptype -  0 - local y distortion
220   //           //1 - local z distortion
221   //           2 - local derivative distortion
222   //           //3
223   //           4 - distortion in the curvature
224   //  nsteps - number of fit points
225   //
226   // return value -  distortion at point refX with type ptype
227   //
228   static TLinearFitter fitter(3,"pol2"); 
229   fitter.ClearPoints();
230   //
231   if (nsteps<4) nsteps=4;
232   Double_t deltaX=(245-85)/(nsteps);
233   AliRieman rieman(nsteps);
234   if (wt<50){
235     AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);          
236     if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
237   }
238   Double_t curv=2.*snp/refX;
239   Double_t dPhi0=TMath::ASin(snp);
240   for (Int_t istep=0; istep<(nsteps+1); istep++){
241     //
242     Double_t localX =85.+deltaX*istep;
243     Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
244     Double_t localPhi=phi0+dPhi;
245     Double_t sector = 9*localPhi/TMath::Pi();
246     if (sector<0) sector+=18;
247     if (sector>18) sector-=18;
248     Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
249     Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
250     Double_t x[1]={localX-dlocalX};
251     Double_t z=theta*x[0];
252     rieman.AddPoint(x[0],dy,z,0.1,0.1);
253     fitter.AddPoint(x,dy);
254   }
255   fitter.Eval();
256   rieman.Update();
257   //
258   Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
259   if (schi2>fgMaxChi2HelixAt){
260     ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("Bad chi2\t%2.2f",schi2).Data());
261     if (pstatus) (*pstatus)=kFALSE;    
262     return 0;
263   }
264   if (pstatus) (*pstatus)=kTRUE;    
265   if (ptype==0) return rieman.GetYat(evalX);
266   if (ptype==2) return rieman.GetDYat(evalX);
267   if (ptype==4) return rieman.GetC();
268   return 0;
269 }
270
271
272
273
274 void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){
275   //
276   // Create cluster distortion map
277   //
278   TFile *falign = TFile::Open("CalibObjects.root");
279   TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign");
280   if (!arrayAlign) {
281     AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
282     return;
283   }
284   AliTPCcalibAlign * align =  (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC");
285   if (!align) {
286       AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
287     return;
288   }
289   TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root");
290
291   THnBase * hdY = (THnBase*)align->GetClusterDelta(0);
292   //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1);
293   AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz);
294   //  AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz);
295
296   const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
297   for (Int_t ihis=0; ihis<4; ihis++){
298     THnSparse * hisAlign =align->GetTrackletDelta(ihis);
299     AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz);
300   }
301   delete pcstream;
302   delete falign;
303 }
304
305
306
307
308
309 void  AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){
310   //
311   // Make cluster residual map from the n-dimensional histogram
312   // hisInput supposed to have given format:
313   //     4 Dim:
314   //      delta, 
315   //      sector 
316   //      localX
317   //      kZ
318   // Vertex position assumed to be at (0,0,0)          
319   //
320   //TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
321   //
322   Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
323   Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
324   Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
325   TF1 *fgaus=0;
326   TH3F * hisResMap3D = 
327     new TH3F("his3D","his3D",
328              nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
329              nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
330              nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
331   hisResMap3D->GetXaxis()->SetTitle("sector");
332   hisResMap3D->GetYaxis()->SetTitle("localX");
333   hisResMap3D->GetZaxis()->SetTitle("kZ");
334
335   TH2F * hisResMap2D[4] ={0,0,0,0};
336   for (Int_t i=0; i<4; i++){
337     hisResMap2D[i]=
338       new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
339                nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
340                nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
341     hisResMap2D[i]->GetXaxis()->SetTitle("sector");
342     hisResMap2D[i]->GetYaxis()->SetTitle("localX");
343   }
344   //
345   //
346   //
347   TF1 * f1= 0;
348   Int_t axis0[4]={0,1,2,3};
349   Int_t axis1[4]={0,1,2,3};
350   Int_t counter=0;
351   for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
352     // phi- sector  range
353     hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
354     THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0); 
355     Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
356     //
357     for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
358       // local x range
359       // kz fits
360       his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
361       THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1); 
362       Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
363       //      
364       //A side
365       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
366       TH1 * hisA = his2->Projection(0);
367       Double_t meanA= hisA->GetMean();
368       Double_t rmsA= hisA->GetRMS();
369       Double_t entriesA= hisA->GetEntries();
370       delete hisA;
371       //C side
372       his2->GetAxis(3)->SetRangeUser(0.01,0.3);
373       TH1 * hisC = his2->Projection(0);
374       Double_t meanC= hisC->GetMean();
375       Double_t rmsC= hisC->GetRMS();
376       Double_t entriesC= hisC->GetEntries();
377       delete hisC;
378       his2->GetAxis(3)->SetRangeUser(-1.2,1.2);      
379       TH2 * hisAC       = his2->Projection(0,3);
380       TProfile *profAC  = hisAC->ProfileX(); 
381       delete hisAC;
382       profAC->Fit("pol1","QNR","QNR",0.05,1);
383       if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
384       Double_t offsetA=f1->GetParameter(0);
385       Double_t slopeA=f1->GetParameter(1);
386       Double_t offsetAE=f1->GetParError(0);
387       Double_t slopeAE=f1->GetParError(1); 
388       Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
389       profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
390       f1=(TF1*)gROOT->FindObject("pol1");
391       Double_t offsetC=f1->GetParameter(0);
392       Double_t slopeC=f1->GetParameter(1); 
393       Double_t offsetCE=f1->GetParError(0);
394       Double_t slopeCE=f1->GetParError(1); 
395       Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
396       if (counter%50==0) printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
397       counter++;
398       (*pcstream)<<"deltaFit"<<
399         "sector="<<sector<<
400         "localX="<<localX<<
401         "meanA="<<meanA<<
402         "rmsA="<<rmsA<<
403         "entriesA="<<entriesA<<
404         "meanC="<<meanC<<
405         "rmsC="<<rmsC<<
406         "entriesC="<<entriesC<<
407         "offsetA="<<offsetA<<
408         "slopeA="<<slopeA<<
409         "offsetAE="<<offsetAE<<
410         "slopeAE="<<slopeAE<<
411         "chi2A="<<chi2A<<
412         "offsetC="<<offsetC<<
413         "slopeC="<<slopeC<<
414         "offsetCE="<<offsetCE<<
415         "slopeCE="<<slopeCE<<
416         "chi2C="<<chi2C<<
417         "\n";
418       //
419       hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
420       hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
421       hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
422       hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
423       
424       for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
425         Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
426         if (TMath::Abs(kZ)<0.05) continue;  // crossing 
427         his2->GetAxis(3)->SetRange(ibin3,ibin3);
428         if (TMath::Abs(kZ)>0.15){
429           his2->GetAxis(3)->SetRange(ibin3,ibin3);
430         }
431         TH1 * his = his2->Projection(0);
432         Double_t mean= his->GetMean();
433         Double_t rms= his->GetRMS();
434         Double_t entries= his->GetEntries();
435         //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
436         hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
437         Double_t phi=TMath::Pi()*sector/9;
438         if (phi>TMath::Pi()) phi+=TMath::Pi();
439         Double_t meanG=0;
440         Double_t rmsG=0;
441         if (entries>50){
442           if (!fgaus) {     
443             his->Fit("gaus","Q","goff");
444             fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
445           }
446           if (fgaus) {
447             his->Fit(fgaus,"Q","goff");
448             meanG=fgaus->GetParameter(1);
449             rmsG=fgaus->GetParameter(2);
450           }
451         }
452         Double_t dsec=sector-Int_t(sector)-0.5;
453         Double_t snp=dsec*TMath::Pi()/9.;
454         (*pcstream)<<"delta"<<
455           "ptype="<<ptype<<
456           "dtype="<<dtype<<
457           "sector="<<sector<<
458           "dsec="<<dsec<<
459           "snp="<<snp<<
460           "phi="<<phi<<
461           "localX="<<localX<<
462           "kZ="<<kZ<<
463           "theta="<<kZ<<
464           "mean="<<mean<<
465           "rms="<<rms<<
466           "meanG="<<meanG<<
467           "rmsG="<<rmsG<<
468           "entries="<<entries<<
469           "meanA="<<meanA<<
470           "rmsA="<<rmsA<<
471           "entriesA="<<entriesA<<
472           "meanC="<<meanC<<
473           "rmsC="<<rmsC<<
474           "entriesC="<<entriesC<<
475           "offsetA="<<offsetA<<
476           "slopeA="<<slopeA<<
477           "chi2A="<<chi2A<<
478           "offsetC="<<offsetC<<
479           "slopeC="<<slopeC<<
480           "chi2C="<<chi2C<<
481           "\n";
482         delete his;
483       }
484       delete his2;
485     }
486     delete his1;
487   }
488   hisResMap3D->Write();
489   hisResMap2D[0]->Write();
490   hisResMap2D[1]->Write();
491   hisResMap2D[2]->Write();
492   hisResMap2D[3]->Write();
493   //  delete pcstream;
494 }
495
496
497
498 void   AliTPCCorrectionFit::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ, Double_t bz){
499   //
500   // make a distortion map out ou fthe residual histogram
501   // Results are written to the debug streamer - pcstream
502   // Parameters:
503   //   his0       - input (4D) residual histogram
504   //   pcstream   - file to write the tree
505   //   run        - run number
506   //   refX       - track matching reference X
507   //   type       - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
508   // THnSparse axes:
509   // OBJ: TAxis     #Delta  #Delta
510   // OBJ: TAxis     tanTheta        tan(#Theta)
511   // OBJ: TAxis     phi     #phi
512   // OBJ: TAxis     snp     snp
513
514   // marian.ivanov@cern.ch
515   const Int_t kMinEntries=10;
516   Int_t idim[4]={0,1,2,3};
517   //
518   //
519   //
520   Int_t nbins3=his0->GetAxis(3)->GetNbins();
521   Int_t first3=his0->GetAxis(3)->GetFirst();
522   Int_t last3 =his0->GetAxis(3)->GetLast();
523   //
524   for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){   // axis 3 - local angle
525     his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
526     Double_t      x3= his0->GetAxis(3)->GetBinCenter(ibin3);
527     THnSparse * his3= his0->Projection(3,idim);         //projected histogram according selection 3
528     //
529     Int_t nbins2    = his3->GetAxis(2)->GetNbins();
530     Int_t first2    = his3->GetAxis(2)->GetFirst();
531     Int_t last2     = his3->GetAxis(2)->GetLast();
532     //
533     for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){   // axis 2 - phi
534       his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
535       Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
536       THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
537       Int_t nbins1     = his2->GetAxis(1)->GetNbins();
538       Int_t first1     = his2->GetAxis(1)->GetFirst();
539       Int_t last1      = his2->GetAxis(1)->GetLast();
540       for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - theta
541         //
542         Double_t       x1= his2->GetAxis(1)->GetBinCenter(ibin1);
543         his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
544         if (TMath::Abs(x1)<0.1){
545           if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
546           if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
547         }
548         if (TMath::Abs(x1)<0.06){
549           his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
550         }
551         TH1 * hisDelta = his2->Projection(0);
552         //
553         Double_t entries = hisDelta->GetEntries();
554         Double_t mean=0, rms=0;
555         if (entries>kMinEntries){
556           mean    = hisDelta->GetMean(); 
557           rms = hisDelta->GetRMS(); 
558         }
559         Double_t sector = 9.*x2/TMath::Pi();
560         if (sector<0) sector+=18;
561         Double_t dsec = sector-Int_t(sector)-0.5;
562         Double_t z=refX*x1;
563         (*pcstream)<<hname<<
564           "run="<<run<<
565           "bz="<<bz<<
566           "theta="<<x1<<
567           "phi="<<x2<<
568           "z="<<z<<            // dummy z
569           "snp="<<x3<<
570           "entries="<<entries<<
571           "mean="<<mean<<
572           "rms="<<rms<<
573           "refX="<<refX<<   // track matching refernce plane
574           "type="<<type<<   //
575           "sector="<<sector<<
576           "dsec="<<dsec<<
577           "\n";
578         delete hisDelta;
579         //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
580       }
581       delete his2;
582     }
583     delete his3;
584   }
585 }
586
587
588
589 void   AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
590   //
591   // make a distortion map out ou fthe residual histogram
592   // Results are written to the debug streamer - pcstream
593   // Parameters:
594   //   his0       - input (4D) residual histogram
595   //   pcstream   - file to write the tree
596   //   run        - run number
597   //   refX       - track matching reference X
598   //   type       - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
599   // marian.ivanov@cern.ch
600   //
601   //  Histo axeses
602   //   Collection name='TObjArray', class='TObjArray', size=16
603   //  0. OBJ: TAxis     #Delta  #Delta
604   //  1. OBJ: TAxis     N_{cl}  N_{cl}
605   //  2. OBJ: TAxis     dca_{r} (cm)    dca_{r} (cm)
606   //  3. OBJ: TAxis     z (cm)  z (cm)
607   //  4. OBJ: TAxis     sin(#phi)       sin(#phi)
608   //  5. OBJ: TAxis     tan(#theta)     tan(#theta)
609   //  6. OBJ: TAxis     1/pt (1/GeV)    1/pt (1/GeV)
610   //  7. OBJ: TAxis     pt (GeV)        pt (GeV)
611   //  8. OBJ: TAxis     alpha   alpha
612   const Int_t kMinEntries=10;
613   //
614   //  1. make default selections
615   //
616   TH1 * hisDelta=0;
617   Int_t idim0[4]={0 , 5, 8,  3};   // delta, theta, alpha, z
618   hisInput->GetAxis(1)->SetRangeUser(110,190);   //long tracks
619   hisInput->GetAxis(2)->SetRangeUser(-10,35);    //tracks close to beam pipe
620   hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
621   hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
622   hisDelta= hisInput->Projection(0);
623   hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
624   delete hisDelta;
625   THnSparse *his0=  hisInput->Projection(4,idim0);
626   //
627   // 2. Get mean in diferent bins
628   //
629   Int_t nbins1=his0->GetAxis(1)->GetNbins();
630   Int_t first1=his0->GetAxis(1)->GetFirst();
631   Int_t last1 =his0->GetAxis(1)->GetLast();
632   //
633   Double_t bz=AliTrackerBase::GetBz();
634   Int_t idim[4]={0,1, 2,  3};  // delta, theta,alpha,z
635   //
636   for (Int_t ibin1=first1; ibin1<=last1; ibin1++){   //axis 1 - theta
637     //
638     Double_t       x1= his0->GetAxis(1)->GetBinCenter(ibin1);  
639     his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
640     //
641     THnSparse * his1 = his0->Projection(4,idim);  // projected histogram according range1
642     Int_t nbins3     = his1->GetAxis(3)->GetNbins();
643     Int_t first3     = his1->GetAxis(3)->GetFirst();
644     Int_t last3      = his1->GetAxis(3)->GetLast();
645     //
646     for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){   // axis 3 - z at "vertex"
647       his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
648       Double_t      x3= his1->GetAxis(3)->GetBinCenter(ibin3);
649       if (ibin3<first3) {
650         his1->GetAxis(3)->SetRangeUser(-1,1);
651         x3=0;
652       }
653       THnSparse * his3= his1->Projection(4,idim);         //projected histogram according selection 3
654       Int_t nbins2    = his3->GetAxis(2)->GetNbins();
655       Int_t first2    = his3->GetAxis(2)->GetFirst();
656       Int_t last2     = his3->GetAxis(2)->GetLast();
657       //
658       for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
659         his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
660         Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
661         hisDelta = his3->Projection(0);
662         //
663         Double_t entries = hisDelta->GetEntries();
664         Double_t mean=0, rms=0;
665         if (entries>kMinEntries){
666           mean    = hisDelta->GetMean(); 
667           rms = hisDelta->GetRMS(); 
668         }
669         Double_t sector = 9.*x2/TMath::Pi();
670         if (sector<0) sector+=18;
671         Double_t dsec = sector-Int_t(sector)-0.5;
672         Double_t snp=0;  // dummy snp - equal 0
673         (*pcstream)<<hname<<
674           "run="<<run<<
675           "bz="<<bz<<            // magnetic field
676           "theta="<<x1<<         // theta
677           "phi="<<x2<<           // phi (alpha)
678           "z="<<x3<<             // z at "vertex"
679           "snp="<<snp<<          // dummy snp
680           "entries="<<entries<<  // entries in bin
681           "mean="<<mean<<        // mean
682           "rms="<<rms<<
683           "refX="<<refX<<        // track matching refernce plane
684           "type="<<type<<        // parameter type
685           "sector="<<sector<<    // sector
686           "dsec="<<dsec<<        // dummy delta sector
687           "\n";
688         delete hisDelta;
689         printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
690       }
691       delete his3;
692     }
693     delete his1;
694   }
695   delete his0;
696 }
697
698
699
700 void   AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){
701   //
702   // make a distortion map out of the residual histogram
703   // Results are written to the debug streamer - pcstream
704   // Parameters:
705   //   his0       - input (4D) residual histogram
706   //   pcstream   - file to write the tree
707   //   run        - run number
708   //   type       - 0- y 1-z,2 -snp, 3-theta
709   //   bz         - magnetic field
710   // marian.ivanov@cern.ch
711
712   //Collection name='TObjArray', class='TObjArray', size=16
713   //0  OBJ: TAxis     delta   delta
714   //1  OBJ: TAxis     phi     phi
715   //2  OBJ: TAxis     localX  localX
716   //3  OBJ: TAxis     kY      kY
717   //4  OBJ: TAxis     kZ      kZ
718   //5  OBJ: TAxis     is1     is1
719   //6  OBJ: TAxis     is0     is0
720   //7. OBJ: TAxis     z       z
721   //8. OBJ: TAxis     IsPrimary       IsPrimary
722
723   const Int_t kMinEntries=10;
724   THnSparse * hisSector0=0;
725   TH1 * htemp=0;    // histogram to calculate mean value of parameter
726   //  Double_t bz=AliTrackerBase::GetBz();
727
728   //
729   // Loop over pair of sector:
730   // isPrim         - 8  ==> 8
731   // isec0          - 6  ==> 7
732   //   isec1        - 5  ==> 6
733   //     refX       - 2  ==> 5
734   //
735   //     phi        - 1  ==> 4
736   //       z        - 7  ==> 3
737   //         snp    - 3  ==> 2
738   //           theta- 4  ==> 1
739   //                  0  ==> 0;           
740   for (Int_t isec0=0; isec0<72; isec0++){
741     Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
742     //
743     //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4);  // select secondaries only ? - get out later ?
744     hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
745     hisSector0=hisInput->Projection(7,index0);
746     //
747     //
748     for (Int_t isec1=isec0+1; isec1<72; isec1++){    
749       //if (isec1!=isec0+36) continue;
750       if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
751       printf("Sectors %d\t%d\n",isec1,isec0);
752       hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);      
753       TH1 * hisX=hisSector0->Projection(5);
754       Double_t refX= hisX->GetMean();
755       delete hisX;
756       TH1 *hisDelta=hisSector0->Projection(0);
757       Double_t dmean = hisDelta->GetMean();
758       Double_t drms = hisDelta->GetRMS();
759       hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
760       delete hisDelta;
761       //
762       //  1. make default selections
763       //
764       Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
765       THnBase *hisSector1=  hisSector0->ProjectionND(5,idim0);
766       //
767       // 2. Get mean in diferent bins
768       //
769       Int_t idim[5]={0, 1, 2,  3, 4};  // {delta, theta-1,snp-2 ,z-3, phi-4}
770       //
771       //      Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
772       Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
773       Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
774       //
775       for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){   //axis 4 - phi
776         //
777         // Phi loop
778         //
779         Double_t       xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);         
780         Double_t psec    = (9*xPhi/TMath::Pi());
781         if (psec<0) psec+=18;
782         Bool_t isOK0=kFALSE;
783         Bool_t isOK1=kFALSE;
784         if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.)  isOK0=kTRUE;
785         if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.)  isOK1=kTRUE;
786         if (!isOK0) continue;
787         if (!isOK1) continue;
788         //
789         hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
790         if (isec1!=isec0+36) {
791           hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
792         }
793         //
794         htemp = hisSector1->Projection(4);
795         xPhi=htemp->GetMean();
796         delete htemp;
797         THnBase * hisPhi = hisSector1->ProjectionND(4,idim);
798         //Int_t nbinsZ     = hisPhi->GetAxis(3)->GetNbins();
799         Int_t firstZ     = hisPhi->GetAxis(3)->GetFirst();
800         Int_t lastZ      = hisPhi->GetAxis(3)->GetLast();
801         //
802         for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){   // axis 3 - z
803           //
804           // Z loop
805           //
806           hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
807           if (isec1!=isec0+36) {
808             hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));     
809           }
810           htemp = hisPhi->Projection(3);
811           Double_t      xZ= htemp->GetMean();
812           delete htemp;
813           THnBase * hisZ= hisPhi->ProjectionND(3,idim);         
814           //projected histogram according selection 3 -z
815           //
816           //
817           //Int_t nbinsSnp    = hisZ->GetAxis(2)->GetNbins();
818           Int_t firstSnp    = hisZ->GetAxis(2)->GetFirst();
819           Int_t lastSnp     = hisZ->GetAxis(2)->GetLast();
820           for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){   // axis 2 - snp
821             //
822             // Snp loop
823             //
824             hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
825             if (isec1!=isec0+36) {
826               hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
827             }
828             htemp = hisZ->Projection(2);
829             Double_t      xSnp= htemp->GetMean();
830             delete htemp;
831             THnBase * hisSnp= hisZ->ProjectionND(2,idim);         
832             //projected histogram according selection 2 - snp
833             
834             //Int_t nbinsTheta    = hisSnp->GetAxis(1)->GetNbins();
835             Int_t firstTheta    = hisSnp->GetAxis(1)->GetFirst();
836             Int_t lastTheta     = hisSnp->GetAxis(1)->GetLast();
837             //
838             for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){  // axis1 theta
839               
840               
841               hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
842               if (isec1!=isec0+36) {
843                  hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));             
844               }
845               htemp = hisSnp->Projection(1);          
846               Double_t xTheta=htemp->GetMean();
847               delete htemp;
848               hisDelta = hisSnp->Projection(0);
849               //
850               Double_t entries = hisDelta->GetEntries();
851               Double_t mean=0, rms=0;
852               if (entries>kMinEntries){
853                 mean    = hisDelta->GetMean(); 
854                 rms = hisDelta->GetRMS(); 
855               }
856               Double_t sector = 9.*xPhi/TMath::Pi();
857               if (sector<0) sector+=18;
858               Double_t dsec = sector-Int_t(sector)-0.5;
859               Int_t dtype=1;  // TPC alignment type
860               (*pcstream)<<hname<<
861                 "run="<<run<<
862                 "bz="<<bz<<             // magnetic field
863                 "ptype="<<type<<         // parameter type
864                 "dtype="<<dtype<<         // parameter type
865                 "isec0="<<isec0<<       // sector 0 
866                 "isec1="<<isec1<<       // sector 1             
867                 "sector="<<sector<<     // sector as float
868                 "dsec="<<dsec<<         // delta sector
869                 //
870                 "theta="<<xTheta<<      // theta
871                 "phi="<<xPhi<<          // phi (alpha)        
872                 "z="<<xZ<<              // z
873                 "snp="<<xSnp<<          // snp
874                 
875                 //
876                 "entries="<<entries<<  // entries in bin
877                 "mean="<<mean<<        // mean
878                 "rms="<<rms<<          // rms 
879                 "refX="<<refX<<        // track matching reference plane
880                 "\n";
881               delete hisDelta;
882               //printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean);
883               //
884             }//ibinTheta
885             delete hisSnp;
886           } //ibinSnp
887           delete hisZ;
888         }//ibinZ
889         delete hisPhi;
890       }//ibinPhi
891       delete hisSector1;      
892     }//isec1
893     delete hisSector0;
894   }//isec0
895 }
896
897
898
899
900
901 // void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
902 //   //
903 //   // Make a laser fit tree for global minimization
904 //   //  
905 //   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();  
906 //   AliTPCCorrection * correction = calib->GetTPCComposedCorrection();  
907 //   if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());  
908 //   correction->AddVisualCorrection(correction,0);  //register correction
909
910 //   //  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
911 //   //AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
912 //   //
913 //   const Double_t cutErrY=0.05;
914 //   const Double_t kSigmaCut=4;
915 //   //  const Double_t cutErrZ=0.03;
916 //   const Double_t kEpsilon=0.00000001;
917 //   //  const Double_t kMaxDist=1.;  // max distance - space correction
918 //   TVectorD *vecdY=0;
919 //   TVectorD *vecdZ=0;
920 //   TVectorD *veceY=0;
921 //   TVectorD *veceZ=0;
922 //   AliTPCLaserTrack *ltr=0;
923 //   AliTPCLaserTrack::LoadTracks();
924 //   tree->SetBranchAddress("dY.",&vecdY);
925 //   tree->SetBranchAddress("dZ.",&vecdZ);
926 //   tree->SetBranchAddress("eY.",&veceY);
927 //   tree->SetBranchAddress("eZ.",&veceZ);
928 //   tree->SetBranchAddress("LTr.",&ltr);
929 //   Int_t entries= tree->GetEntries();
930 //   TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
931 //   Double_t bz=AliTrackerBase::GetBz();
932 //   // 
933 //   //  Double_t globalXYZ[3];
934 //   //Double_t globalXYZCorr[3];
935 //   for (Int_t ientry=0; ientry<entries; ientry++){
936 //     tree->GetEntry(ientry);
937 //     if (!ltr->GetVecGX()){
938 //       ltr->UpdatePoints();
939 //     }
940 //     //
941 //     TVectorD fit10(5);
942 //     TVectorD fit5(5);
943 //     printf("Entry\t%d\n",ientry);
944 //     for (Int_t irow0=0; irow0<158; irow0+=1){
945 //       //       
946 //       TLinearFitter fitter10(4,"hyp3");
947 //       TLinearFitter fitter5(2,"hyp1");
948 //       Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
949 //       if (sector<0) continue;
950 //       //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
951
952 //       Double_t refX= (*ltr->GetVecLX())[irow0];
953 //       Int_t firstRow1 = TMath::Max(irow0-10,0);
954 //       Int_t lastRow1  = TMath::Min(irow0+10,158);
955 //       Double_t padWidth=(irow0<64)?0.4:0.6;
956 //       // make long range fit
957 //       for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
958 //      if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
959 //      if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
960 //      if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
961 //      Double_t idealX= (*ltr->GetVecLX())[irow1];
962 //      Double_t idealY= (*ltr->GetVecLY())[irow1];
963 //      //      Double_t idealZ= (*ltr->GetVecLZ())[irow1];
964 //      Double_t gx= (*ltr->GetVecGX())[irow1];
965 //      Double_t gy= (*ltr->GetVecGY())[irow1];
966 //      Double_t gz= (*ltr->GetVecGZ())[irow1];
967 //      Double_t measY=(*vecdY)[irow1]+idealY;
968 //      Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
969 //      // deltaR = R distorted -R ideal
970 //      Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
971 //      fitter10.AddPoint(xxx,measY,1);
972 //       }
973 //       Bool_t isOK=kTRUE;
974 //       Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
975 //       Double_t mean10  =0;//   fitter10.GetParameter(0);
976 //       Double_t slope10  =0;//   fitter10.GetParameter(0);
977 //       Double_t cosPart10  = 0;//  fitter10.GetParameter(2);
978 //       Double_t sinPart10   =0;//  fitter10.GetParameter(3); 
979
980 //       if (fitter10.GetNpoints()>10){
981 //      fitter10.Eval();
982 //      rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
983 //      mean10      =   fitter10.GetParameter(0);
984 //      slope10     =   fitter10.GetParameter(1);
985 //      cosPart10   =   fitter10.GetParameter(2);
986 //      sinPart10   =  fitter10.GetParameter(3); 
987 //      //
988 //      // make short range fit
989 //      //
990 //      for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
991 //        if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
992 //        if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
993 //        if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
994 //        Double_t idealX= (*ltr->GetVecLX())[irow1];
995 //        Double_t idealY= (*ltr->GetVecLY())[irow1];
996 //        //      Double_t idealZ= (*ltr->GetVecLZ())[irow1];
997 //        Double_t gx= (*ltr->GetVecGX())[irow1];
998 //        Double_t gy= (*ltr->GetVecGY())[irow1];
999 //        Double_t gz= (*ltr->GetVecGZ())[irow1];
1000 //        Double_t measY=(*vecdY)[irow1]+idealY;
1001 //        Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
1002 //        // deltaR = R distorted -R ideal 
1003 //        Double_t expY= mean10+slope10*(idealX+deltaR-refX);
1004 //        if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
1005 //        //
1006 //        Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
1007 //        Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
1008 //        fitter5.AddPoint(xxx,measY-corr,1);
1009 //      }     
1010 //       }else{
1011 //      isOK=kFALSE;
1012 //       }
1013 //       if (fitter5.GetNpoints()<8) isOK=kFALSE;
1014
1015 //       Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
1016 //       Double_t offset5  =0;//  fitter5.GetParameter(0);
1017 //       Double_t slope5   =0;//  fitter5.GetParameter(0); 
1018 //       if (isOK){
1019 //      fitter5.Eval();
1020 //      rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
1021 //      offset5  =  fitter5.GetParameter(0);
1022 //      slope5   =  fitter5.GetParameter(0); 
1023 //       }
1024 //       //
1025 //       Double_t dtype=5;
1026 //       Double_t ptype=0;
1027 //       Double_t phi   =(*ltr->GetVecPhi())[irow0];
1028 //       Double_t theta =ltr->GetTgl();
1029 //       Double_t mean=(vecdY)->GetMatrixArray()[irow0];
1030 //       Double_t gx=0,gy=0,gz=0;
1031 //       Double_t snp = (*ltr->GetVecP2())[irow0];
1032 //       Int_t bundle= ltr->GetBundle();
1033 //       Int_t id= ltr->GetId();
1034 //       //      Double_t rms = err->GetMatrixArray()[irow];
1035 //       //
1036 //       gx = (*ltr->GetVecGX())[irow0];
1037 //       gy = (*ltr->GetVecGY())[irow0];
1038 //       gz = (*ltr->GetVecGZ())[irow0];
1039 //       Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
1040 //       fitter10.GetParameters(fit10);
1041 //       fitter5.GetParameters(fit5);      
1042 //       Double_t idealY= (*ltr->GetVecLY())[irow0];
1043 //       Double_t measY=(*vecdY)[irow0]+idealY;
1044 //       Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
1045 //       if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
1046 //       //
1047 //       (*pcstream)<<"fitFull"<<  // dumpe also intermediate results
1048 //      "bz="<<bz<<         // magnetic filed used
1049 //      "dtype="<<dtype<<   // detector match type
1050 //      "ptype="<<ptype<<   // parameter type
1051 //      "theta="<<theta<<   // theta
1052 //      "phi="<<phi<<       // phi 
1053 //      "snp="<<snp<<       // snp
1054 //      "sector="<<sector<<
1055 //      "bundle="<<bundle<<
1056 // //   //      "dsec="<<dsec<<
1057 //      "refX="<<refX<<      // reference radius
1058 //      "gx="<<gx<<         // global position
1059 //      "gy="<<gy<<         // global position
1060 //      "gz="<<gz<<         // global position
1061 //      "dRrec="<<dRrec<<      // delta Radius in reconstruction
1062 //      "id="<<id<<     //bundle
1063 //      "rms10="<<rms10<<
1064 //      "rms5="<<rms5<<
1065 //      "fit10.="<<&fit10<<
1066 //      "fit5.="<<&fit5<<
1067 //      "measY="<<measY<<
1068 //      "mean="<<mean<<
1069 //      "idealY="<<idealY<<
1070 //      "corr="<<corr<<
1071 //      "isOK="<<isOK<<
1072 //      "\n";
1073 //     }
1074 //   }
1075 //   delete pcstream;
1076 // }
1077
1078
1079 void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1080   //
1081   // Make a fit tree:
1082   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1083   // calculates partial distortions
1084   // Partial distortion is stored in the resulting tree
1085   // Output is storred in the file distortion_<dettype>_<partype>.root
1086   // Partial  distortion is stored with the name given by correction name
1087   //
1088   //
1089   // Parameters of function:
1090   // input     - input tree
1091   // dtype     - distortion type 0 - ITSTPC,  1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF,  4 - TPCTPC track crossing 
1092   // ppype     - parameter type
1093   // corrArray - array with partial corrections
1094   // step      - skipe entries  - if 1 all entries processed - it is slow
1095   // debug     0 if debug on also space points dumped - it is slow
1096
1097   const Double_t kMaxSnp = 0.85;  
1098   const Double_t kcutSnp=0.25;
1099   const Double_t kcutTheta=1.;
1100   const Double_t kRadiusTPC=85;
1101   //  AliTPCROC *tpcRoc =AliTPCROC::Instance();  
1102   //
1103   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1104   //  const Double_t kB2C=-0.299792458e-3;
1105   const Int_t kMinEntries=20; 
1106   Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
1107   Float_t refX;  
1108   Int_t run;
1109   tinput->SetBranchAddress("run",&run);
1110   tinput->SetBranchAddress("theta",&theta);
1111   tinput->SetBranchAddress("phi", &phi);
1112   tinput->SetBranchAddress("snp",&snp);
1113   tinput->SetBranchAddress("mean",&mean);
1114   tinput->SetBranchAddress("rms",&rms);
1115   tinput->SetBranchAddress("entries",&entries);
1116   tinput->SetBranchAddress("sector",&sector);
1117   tinput->SetBranchAddress("dsec",&dsec);
1118   tinput->SetBranchAddress("refX",&refX);
1119   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
1120   //
1121   Int_t nentries=tinput->GetEntries();
1122   Int_t ncorr=corrArray->GetEntries();
1123   Double_t corrections[100]={0}; //
1124   Double_t tPar[5];
1125   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1126   Int_t dir=0;
1127   if (dtype==5 || dtype==6) dtype=4;
1128   if (dtype==0) { dir=-1;}
1129   if (dtype==1) { dir=1;}
1130   if (dtype==2) { dir=-1;}
1131   if (dtype==3) { dir=1;}
1132   if (dtype==4) { dir=-1;}
1133   //
1134   for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1135     tinput->GetEntry(ientry);
1136     if (TMath::Abs(snp)>kMaxSnp) continue;
1137     tPar[0]=0;
1138     tPar[1]=theta*refX;
1139     if (dtype==2)  tPar[1]=theta*kRadiusTPC;
1140     tPar[2]=snp;
1141     tPar[3]=theta;
1142     tPar[4]=(gRandom->Rndm()-0.5)*0.02;  // should be calculated - non equal to 0
1143     if (dtype==4){
1144       // tracks crossing CE
1145       tPar[1]=0;   // track at the CE
1146       //if (TMath::Abs(theta) <0.05) continue;  // deep cross
1147     }
1148
1149     if (TMath::Abs(snp) >kcutSnp) continue;
1150     if (TMath::Abs(theta) >kcutTheta) continue;
1151     printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1152     Double_t bz=AliTrackerBase::GetBz();
1153     if (dtype !=4) { //exclude TPC  - for TPC mainly non primary tracks
1154       if (dtype!=2 && TMath::Abs(bz)>0.1 )  tPar[4]=snp/(refX*bz*kB2C*2);
1155       
1156       if (dtype==2 && TMath::Abs(bz)>0.1 )  {
1157         tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1158         // snp at the TPC inner radius in case the vertex match used
1159       }
1160     }
1161     //
1162     tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1163     AliExternalTrackParam track(refX,phi,tPar,cov);
1164     Double_t xyz[3];
1165     track.GetXYZ(xyz);
1166     Int_t id=0;
1167     Double_t pt=1./tPar[4];
1168     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1169     //if (ptype==4 &&bz<0) mean*=-1;  // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
1170     Double_t refXD=refX;
1171     (*pcstream)<<"fit"<<
1172       "run="<<run<<       // run number
1173       "bz="<<bz<<         // magnetic filed used
1174       "dtype="<<dtype<<   // detector match type
1175       "ptype="<<ptype<<   // parameter type
1176       "theta="<<theta<<   // theta
1177       "phi="<<phi<<       // phi 
1178       "snp="<<snp<<       // snp
1179       "mean="<<mean<<     // mean dist value
1180       "rms="<<rms<<       // rms
1181       "sector="<<sector<<
1182       "dsec="<<dsec<<
1183       "refX="<<refXD<<         // referece X as double
1184       "gx="<<xyz[0]<<         // global position at reference
1185       "gy="<<xyz[1]<<         // global position at reference
1186       "gz="<<xyz[2]<<         // global position at reference   
1187       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1188       "pt="<<pt<<            // pt
1189       "id="<<id<<             // track id
1190       "entries="<<entries;// number of entries in bin
1191     //
1192     Bool_t isOK=kTRUE;
1193     if (entries<kMinEntries) isOK=kFALSE;
1194     //
1195     if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1196       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1197       corrections[icorr]=0;
1198       if (entries>kMinEntries){
1199         AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1200         AliExternalTrackParam *trackOut = 0;
1201         if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1202         if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1203         if (dtype==0) {dir= -1;}
1204         if (dtype==1) {dir=  1;}
1205         if (dtype==2) {dir= -1;}
1206         if (dtype==3) {dir=  1;}
1207         //
1208         if (trackOut){
1209           if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1210           if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
1211           if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1212           //      trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1213           //      
1214           corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1215           delete trackOut;      
1216         }else{
1217           corrections[icorr]=0;
1218           isOK=kFALSE;
1219         }
1220         //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out
1221       }      
1222       (*pcstream)<<"fit"<<
1223         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
1224     }
1225   
1226     if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1227       //
1228       // special case of the TPC tracks crossing the CE
1229       //
1230       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1231       corrections[icorr]=0;
1232       if (entries>kMinEntries){
1233         AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
1234         AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet 
1235         AliExternalTrackParam *trackOut0 = 0;
1236         AliExternalTrackParam *trackOut1 = 0;
1237         //
1238         if (debug)  trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1239         if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1240         if (debug)  trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1241         if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1242         //
1243         if (trackOut0 && trackOut1){
1244           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp))  isOK=kFALSE;
1245           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1246           if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1247           if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1248           //
1249           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1250           if (!trackIn1.Rotate(trackIn0.GetAlpha()))  isOK=kFALSE;
1251           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1252           if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;       
1253           if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1254           //
1255           corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1256           corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1257           if (isOK)
1258             if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
1259                 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
1260                 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
1261                 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
1262                 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
1263                 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
1264                 ){
1265               isOK=kFALSE;
1266             }             
1267           delete trackOut0;      
1268           delete trackOut1;       
1269         }else{
1270           corrections[icorr]=0;
1271           isOK=kFALSE;
1272         }
1273         //
1274         //if (ptype==4 &&bz<0) corrections[icorr]*=-1;  // interpret as curvature - commented out no in lookup
1275       }      
1276       (*pcstream)<<"fit"<<
1277         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
1278     }
1279     //
1280     (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1281   }
1282
1283
1284   delete pcstream;
1285 }
1286
1287
1288
1289 void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1290   //
1291   // Make a fit tree:
1292   // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1293   // calculates partial distortions
1294   // Partial distortion is stored in the resulting tree
1295   // Output is storred in the file distortion_<dettype>_<partype>.root
1296   // Partial  distortion is stored with the name given by correction name
1297   //
1298   //
1299   // Parameters of function:
1300   // input     - input tree
1301   // dtype     - distortion type 10 - IROC-OROC 
1302   // ppype     - parameter type
1303   // corrArray - array with partial corrections
1304   // step      - skipe entries  - if 1 all entries processed - it is slow
1305   // debug     0 if debug on also space points dumped - it is slow
1306
1307   const Double_t kMaxSnp = 0.8;  
1308   const Int_t kMinEntries=200; 
1309   //  AliTPCROC *tpcRoc =AliTPCROC::Instance();  
1310   //
1311   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1312   //  const Double_t kB2C=-0.299792458e-3;
1313   Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
1314   Int_t isec1, isec0;
1315   Double_t refXD;
1316   Float_t refX;
1317   Int_t run;
1318   tinput->SetBranchAddress("run",&run);
1319   tinput->SetBranchAddress("theta",&theta);
1320   tinput->SetBranchAddress("phi", &phi);
1321   tinput->SetBranchAddress("snp",&snp);
1322   tinput->SetBranchAddress("mean",&mean);
1323   tinput->SetBranchAddress("rms",&rms);
1324   tinput->SetBranchAddress("entries",&entries);
1325   tinput->SetBranchAddress("sector",&sector);
1326   tinput->SetBranchAddress("dsec",&dsec);
1327   tinput->SetBranchAddress("refX",&refXD);
1328   tinput->SetBranchAddress("z",&globalZ);
1329   tinput->SetBranchAddress("isec0",&isec0);
1330   tinput->SetBranchAddress("isec1",&isec1);
1331   TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
1332   //
1333   Int_t nentries=tinput->GetEntries();
1334   Int_t ncorr=corrArray->GetEntries();
1335   Double_t corrections[100]={0}; //
1336   Double_t tPar[5];
1337   Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1338   Int_t dir=0;
1339   //
1340   for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1341     tinput->GetEntry(ientry);
1342     refX=refXD;
1343     Int_t id=-1;
1344     if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1;  // IROC-OROC - opposite side
1345     if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2;  // IROC-OROC - same side
1346     if (dtype==10  && id==-1) continue;
1347     //
1348     dir=-1;
1349     tPar[0]=0;
1350     tPar[1]=globalZ;
1351     tPar[2]=snp;
1352     tPar[3]=theta;
1353     tPar[4]=(gRandom->Rndm()-0.1)*0.2;  //
1354     Double_t pt=1./tPar[4];
1355     //
1356     printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1357     Double_t bz=AliTrackerBase::GetBz();
1358     AliExternalTrackParam track(refX,phi,tPar,cov);    
1359     Double_t xyz[3],xyzIn[3],xyzOut[3];
1360     track.GetXYZ(xyz);
1361     track.GetXYZAt(85,bz,xyzIn);    
1362     track.GetXYZAt(245,bz,xyzOut);    
1363     Double_t phiIn  = TMath::ATan2(xyzIn[1],xyzIn[0]);
1364     Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
1365     Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
1366     Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
1367     Int_t sectorIn  = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
1368     Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
1369     //
1370     Bool_t isOK=kTRUE; 
1371     if (sectorIn!=sectorOut) isOK=kFALSE;  // requironment - cluster in the same sector
1372     if (sectorIn!=sectorRef) isOK=kFALSE;  // requironment - cluster in the same sector
1373     if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE;  // requironment - minimal amount of tracks in bin
1374     // Do downscale
1375     if (TMath::Abs(theta)>1) isOK=kFALSE;
1376     //
1377     Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1378     //
1379     (*pcstream)<<"fit"<<
1380       "run="<<run<<       //run
1381       "bz="<<bz<<         // magnetic filed used
1382       "dtype="<<dtype<<   // detector match type
1383       "ptype="<<ptype<<   // parameter type
1384       "theta="<<theta<<   // theta
1385       "phi="<<phi<<       // phi 
1386       "snp="<<snp<<       // snp
1387       "mean="<<mean<<     // mean dist value
1388       "rms="<<rms<<       // rms
1389       "sector="<<sector<<
1390       "dsec="<<dsec<<
1391       "refX="<<refXD<<         // referece X
1392       "gx="<<xyz[0]<<         // global position at reference
1393       "gy="<<xyz[1]<<         // global position at reference
1394       "gz="<<xyz[2]<<         // global position at reference   
1395       "dRrec="<<dRrec<<      // delta Radius in reconstruction
1396       "pt="<<pt<<      //pt
1397       "id="<<id<<             // track id
1398       "entries="<<entries;// number of entries in bin
1399     //
1400     AliExternalTrackParam *trackOut0 = 0;
1401     AliExternalTrackParam *trackOut1 = 0;
1402     AliExternalTrackParam *ptrackIn0 = 0;
1403     AliExternalTrackParam *ptrackIn1 = 0;
1404
1405     for (Int_t icorr=0; icorr<ncorr; icorr++) {
1406       //
1407       // special case of the TPC tracks crossing the CE
1408       //
1409       AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1410       corrections[icorr]=0;
1411       if (entries>kMinEntries &&isOK){
1412         AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
1413         AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
1414         ptrackIn1=&trackIn0;
1415         ptrackIn0=&trackIn1;
1416         //
1417         if (debug)  trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1418         if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1419         if (debug)  trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1420         if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1421         //
1422         if (trackOut0 && trackOut1){
1423           //
1424           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp))  isOK=kFALSE;
1425           if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1426           // rotate all tracks to the same frame
1427           if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1428           if (!trackIn1.Rotate(trackIn0.GetAlpha()))  isOK=kFALSE;
1429           if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;       
1430           //
1431           if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1432           if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1433           if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1434           //
1435           corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1436           corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1437           (*pcstream)<<"fitDebug"<< // just to debug the correction
1438             "mean="<<mean<<
1439             "pIn0.="<<ptrackIn0<<
1440             "pIn1.="<<ptrackIn1<<
1441             "pOut0.="<<trackOut0<<
1442             "pOut1.="<<trackOut1<<
1443             "refX="<<refXD<<
1444             "\n";
1445           delete trackOut0;      
1446           delete trackOut1;      
1447         }else{
1448           corrections[icorr]=0;
1449           isOK=kFALSE;
1450         }
1451       }      
1452       (*pcstream)<<"fit"<<
1453         Form("%s=",corr->GetName())<<corrections[icorr];   // dump correction value
1454     }
1455     //
1456     (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1457   }
1458   delete pcstream;
1459 }
1460
1461