Coding Conventions
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibUnlinearity.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   Class for histogramming of cluster residuals
18   and fitting of the unlinearities. To be used only for data without
19   magnetic field. The laser tracks should be also rejected.
20   //
21
22   Track fitting:
23   The track is fitted using linear and parabolic model 
24   The edge clusters are removed from the fit.
25   Edge pad-row - first and last 15 padrows removed from the fit
26
27   Unlinearities fitting:
28   Unlinearities at the edge aproximated using two exponential decays.
29
30   Model:
31   dz = dz0(r,z) +dr(r,z)*tan(theta) 
32   dy =          +dr(r,z)*tan(phi)
33
34    
35   .x ~/NimStyle.C
36   gSystem->Load("libANALYSIS");
37   gSystem->Load("libTPCcalib");
38   TFile fcalib("CalibObjects.root");
39   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
40   AliTPCcalibUnlinearity * calibUnlin = ( AliTPCcalibUnlinearity *)array->FindObject("calibUnlinearity");
41   //
42   
43 */
44
45
46 #include "TLinearFitter.h"
47
48 #include "Riostream.h"
49 #include "TChain.h"
50 #include "TTree.h"
51 #include "TH1F.h"
52 #include "TH2F.h"
53 #include "TH3F.h"
54 #include "THnSparse.h"
55 #include "TList.h"
56 #include "TMath.h"
57 #include "TCanvas.h"
58 #include "TFile.h"
59 #include "TF1.h"
60 #include "TVectorD.h"
61 #include "TProfile.h"
62 #include "TGraphErrors.h"
63 #include "TCanvas.h"
64 #include "TCut.h"
65
66
67 #include "AliTPCclusterMI.h"
68 #include "AliTPCseed.h"
69 #include "AliESDVertex.h"
70 #include "AliESDEvent.h"
71 #include "AliESDfriend.h"
72 #include "AliESDInputHandler.h"
73 #include "AliAnalysisManager.h"
74 #include "AliTracker.h"
75 #include "AliMagF.h"
76 #include "AliTPCCalROC.h"
77
78 #include "AliLog.h"
79
80
81 #include "TTreeStream.h"
82 #include "AliTPCTracklet.h"
83 #include "TTimeStamp.h"
84 #include "AliTPCcalibDB.h"
85 #include "AliDCSSensorArray.h"
86 #include "AliDCSSensor.h"
87 #include "TLinearFitter.h"
88 #include "AliTPCClusterParam.h"
89 #include "TStatToolkit.h"
90
91
92 #include "AliTPCcalibUnlinearity.h"
93 #include "AliTPCPointCorrection.h"
94
95
96 ClassImp(AliTPCcalibUnlinearity)
97
98 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
99   AliTPCcalibBase(),
100   fInit(kFALSE),
101   fDiffHistoLineY(0),
102   fDiffHistoParY(0),
103   fDiffHistoLineZ(0),
104   fDiffHistoParZ(0),
105   fFittersOutR(38),
106   fFittersOutZ(38),
107   fParamsOutR(38),
108   fParamsOutZ(38),
109   fErrorsOutR(38),
110   fErrorsOutZ(38),
111   fDistRPHIPlus(74),
112   fDistRPHIMinus(74),
113   fFitterQuadrantY(16*38),
114   fFitterQuadrantPhi(16*38)
115 {
116   //
117   // Defualt constructor
118   //
119 }
120
121 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
122   AliTPCcalibBase(name,title),
123   fInit(kFALSE),
124   fDiffHistoLineY(0),
125   fDiffHistoParY(0),  
126   fDiffHistoLineZ(0),
127   fDiffHistoParZ(0),  
128   fFittersOutR(38),
129   fFittersOutZ(38),
130   fParamsOutR(38),
131   fParamsOutZ(38),
132   fErrorsOutR(38),
133   fErrorsOutZ(38),
134   fDistRPHIPlus(74),
135   fDistRPHIMinus(74),
136   fFitterQuadrantY(16*38),
137   fFitterQuadrantPhi(16*38)
138 {
139   //
140   // Non default constructor
141   //
142   MakeHisto();
143 }
144
145 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
146   //
147   //
148   //
149   if (fDiffHistoLineZ) delete fDiffHistoLineY;
150   if (fDiffHistoParZ)  delete fDiffHistoParY;
151   if (fDiffHistoLineZ) delete fDiffHistoLineZ;
152   if (fDiffHistoParZ)  delete fDiffHistoParZ;
153 }
154
155
156
157
158
159 void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
160   //
161   // 
162   //  
163   if (!fInit) {
164     Init();   // work around for PROOF - initialize firs time used
165   }
166   const Int_t  kdrow=10;
167   const Int_t kMinCluster=35;
168   if (!seed) return;
169   if (TMath::Abs(fMagF)>0.01) return;   // use only data without field
170   Int_t counter[72];
171   for (Int_t i=0;i<72;i++) counter[i]=0;
172   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
173     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
174     if (!cluster) continue;
175     if (cluster->GetType()<0) continue;
176     counter[cluster->GetDetector()]++;
177   }
178   //
179   for (Int_t is=0; is<72;is++){
180     if (counter[is]>kMinCluster) ProcessDiff(seed, is);
181     if (counter[is]>kMinCluster) {
182       AlignOROC(seed, is);
183     }
184   }
185 }
186
187
188 void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
189   //
190   // process differnce of the cluster in respect with linear and parabolic fit
191   //
192   const  Double_t kXmean=134;
193   const Int_t     kdrow=10;
194   const Float_t kSagitaCut = 1;
195   static TLinearFitter fy1(2,"pol1");
196   static TLinearFitter fy2(3,"pol2");
197   static TLinearFitter fz1(2,"pol1");
198   static TLinearFitter fz2(3,"pol2");
199   //
200   static TVectorD py1(2);
201   static TVectorD py2(3);
202   static TVectorD pz1(2);
203   static TVectorD pz2(3);  
204   //  
205   fy1.ClearPoints();
206   fy2.ClearPoints();
207   fz1.ClearPoints();
208   fz2.ClearPoints();
209   //
210   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
211     AliTPCclusterMI *c=track->GetClusterPointer(irow);
212     if (!c) continue;
213     if (c->GetDetector()!=isec) continue;
214     if (c->GetType()<0) continue;
215     Double_t dx = c->GetX()-kXmean;
216     Double_t x[2]={dx, dx*dx};
217     fy1.AddPoint(x,c->GetY());
218     fy2.AddPoint(x,c->GetY());
219     fz1.AddPoint(x,c->GetZ());
220     fz2.AddPoint(x,c->GetZ());
221   }
222   fy1.Eval(); fy1.GetParameters(py1);
223   fy2.Eval(); fy2.GetParameters(py2);
224   fz1.Eval(); fz1.GetParameters(pz1);
225   fz2.Eval(); fz2.GetParameters(pz2);
226   //
227   
228   for (Int_t irow=0;irow<159;irow++) {
229     AliTPCclusterMI *c=track->GetClusterPointer(irow);
230     if (!c) continue;
231     if (c->GetDetector()!=isec) continue;
232     Double_t dx = c->GetX()-kXmean;
233     Double_t y1 = py1[0]+py1[1]*dx;
234     Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
235     //
236     Double_t z1 = pz1[0]+pz1[1]*dx;
237     Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
238     //
239     Double_t edgeMinus= -y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
240     Double_t edgePlus =  y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
241     Int_t    npads = AliTPCROC::Instance()->GetNPads(isec,c->GetRow());
242     Float_t    dpad  = TMath::Min(c->GetPad(), npads-1- c->GetPad());
243     Float_t dy =c->GetY()-y1;
244     //
245     // Corrections 
246     //    
247     AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
248     Double_t corrclY=0, corrtrY=0, corrR=0;
249     corrclY =  corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
250                                        c->GetY(),c->GetY(), c->GetZ(), py1[1],c->GetMax(),2.5); 
251     corrtrY =  corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
252                                        c->GetY(),y1, c->GetZ(), py1[1], c->GetMax(),2.5);     
253     corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,c->GetX(),c->GetY(),c->GetZ(),isec);
254     //
255     //
256     //
257     if (fStreamLevel>1){
258       TTreeSRedirector *cstream = GetDebugStreamer();
259       if (cstream){
260         (*cstream)<<"Diff"<<
261           "run="<<fRun<<              //  run number
262           "event="<<fEvent<<          //  event number
263           "time="<<fTime<<            //  time stamp of event
264           "trigger="<<fTrigger<<      //  trigger
265           "mag="<<fMagF<<             //  magnetic field          
266           "isec="<<isec<<             // current sector
267           "npads="<<npads<<           // number of pads at given sector
268           "dpad="<<dpad<<             // distance to edge pad
269           //
270           "crR="<<corrR<<              // radial correction
271           "cclY="<<corrclY<<            // edge effect correction cl
272           "ctrY="<<corrtrY<<            // edge effect correction using track
273           //
274           "Cl.="<<c<<
275           "y1="<<y1<<
276           "y2="<<y2<<
277           "z1="<<z1<<
278           "z2="<<z2<<
279           "py1.="<<&py1<<
280           "py2.="<<&py2<<
281           "pz1.="<<&pz1<<
282           "pz2.="<<&pz2<<
283           "eP="<<edgePlus<<
284           "eM="<<edgeMinus<<
285           "\n";
286       }
287     }
288     if (TMath::Min(edgeMinus,edgePlus)<6){
289       AddPointRPHI(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
290       TTreeSRedirector *cstream = GetDebugStreamer();
291       if (TMath::Abs(dy)<2 && cstream){
292         (*cstream)<<"rphi"<<
293           "isec="<<isec<<             // current sector
294           "npads="<<npads<<           // number of pads at given sector
295           "dpad="<<dpad<<             // distance to edge pad
296           "eP="<<edgePlus<<           // distance to edge 
297           "eM="<<edgeMinus<<          // distance to edge minus
298           //
299           "crR="<<corrR<<              // radial correction
300           "cclY="<<corrclY<<            // edge effect correction cl
301           "ctrY="<<corrtrY<<            // edge effect correction using track
302           //
303           "dy="<<dy<<                 // dy
304           "Cl.="<<c<<
305           "y1="<<y1<<
306           "y2="<<y2<<
307           "z1="<<z1<<
308           "z2="<<z2<<
309           "py1.="<<&py1<<
310           "pz1.="<<&pz1<<
311           "\n";
312       }
313     }
314     if (c->GetType()<0) continue;  // don't use edge rphi cluster
315     //
316     //
317     Double_t x[10];
318     x[0]=c->GetDetector();
319     x[1]=c->GetRow();
320     x[2]=c->GetY()/c->GetX();
321     x[3]=c->GetZ();
322     //
323     // apply sagita cut
324     //
325     if (TMath::Abs(py2[2]*150*150/4)<kSagitaCut && 
326         TMath::Abs(pz2[2]*150*150/4)<kSagitaCut){
327       //
328       x[4]=py1[1];
329       x[5]=c->GetY()-y1;
330       fDiffHistoLineY->Fill(x);
331       x[5]=c->GetY()-y2;
332       //fDiffHistoParY->Fill(x);
333       //
334       x[4]=pz1[1];
335       x[5]=c->GetY()-z1;
336       fDiffHistoLineZ->Fill(x);
337       x[5]=c->GetY()-z2;
338       //fDiffHistoParZ->Fill(x);      
339       // Apply sagita cut
340       AddPoint(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
341     }
342              
343   }
344 }
345
346
347 void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
348   //
349   // fit the tracklet in OROC sepatately in Quadrant
350   //
351   //
352   if (isec<36) return; 
353   const Int_t kMinClusterF=35;
354   const Int_t kMinClusterQ=10;
355   //
356   const  Int_t     kdrow1 =8;       // rows to skip at the end      
357   const  Int_t     kdrow0 =2;        // rows to skip at beginning  
358   const  Float_t   kedgey=3;
359   const  Float_t   kMaxDist=0.5;
360   const  Float_t   kMaxCorrY=0.1;
361   const  Float_t   kPRFWidth = 0.6;   //cut 2 sigma of PRF
362   //
363   //
364   AliTPCROC * roc = AliTPCROC::Instance();
365   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
366   const  Double_t kXmean=roc->GetPadRowRadii(isec,53);
367   //
368   // full track fit parameters
369   // 
370   TLinearFitter fyf(2,"pol1");
371   TLinearFitter fzf(2,"pol1");
372   TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
373   Int_t nf=0;
374   //
375   // make full fit as reference
376   //
377   for (Int_t iter=0; iter<2; iter++){
378     fyf.ClearPoints();
379     for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
380       AliTPCclusterMI *c=track->GetClusterPointer(irow);
381       if (!c) continue;
382       if (c->GetDetector()!=isec) continue;
383       if (c->GetRow()<kdrow0) continue;
384       Double_t dx = c->GetX()-kXmean;
385       Double_t x[2]={dx, dx*dx};
386       if (iter==0 &&c->GetType()<0) continue;
387       if (iter==1){     
388         Double_t yfit  =  pyf[0]+pyf[1]*dx;
389         Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
390         if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
391         if (dedge<kedgey) continue;
392         Double_t corrtrY =  
393           corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
394                                   c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
395         if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
396       }
397       fyf.AddPoint(x,c->GetY(),0.1);
398       fzf.AddPoint(x,c->GetZ(),0.1);
399     }
400     nf = fyf.GetNpoints();
401     if (nf<kMinClusterF) return;   // not enough points - skip 
402     fyf.Eval(); 
403     fyf.GetParameters(pyf); 
404     fyf.GetErrors(peyf);
405     fzf.Eval(); 
406     fzf.GetParameters(pzf); 
407     fzf.GetErrors(pezf);
408   }
409   //
410   // Make Fitters and params for 3 fitters
411   //
412   TLinearFitter *fitters[4];
413   Int_t npoints[4];
414   TVectorD params[4];
415   TVectorD errors[4];
416   Double_t chi2C[4];
417   for (Int_t i=0;i<4;i++) {
418     fitters[i] = new TLinearFitter(2,"pol1");
419     npoints[i]=0;
420     params[i].ResizeTo(2);
421     errors[i].ResizeTo(2);
422   }
423   //
424   // Update fitters
425   //
426   for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
427     AliTPCclusterMI *c=track->GetClusterPointer(irow);
428     if (!c) continue;
429     if (c->GetDetector()!=isec) continue;
430     if (c->GetRow()<kdrow0) continue;      
431     Double_t dx = c->GetX()-kXmean;
432     Double_t x[2]={dx, dx*dx};
433     Double_t yfit  =  pyf[0]+pyf[1]*dx;
434     Double_t dedge =  c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
435     if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
436     if (dedge<kedgey) continue;
437     Double_t corrtrY =  
438       corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
439                               c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
440     if (TMath::Abs(corrtrY)>kMaxCorrY) continue;  
441     if (dx<0){
442       if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
443       if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
444     }
445     if (dx>0){
446       if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
447       if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
448     }
449   }
450
451   for (Int_t i=0;i<4;i++) {
452     npoints[i] = fitters[i]->GetNpoints();
453     if (npoints[i]>=kMinClusterQ){
454       fitters[i]->Eval();
455       Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
456       chi2C[i]=chi2Fac;
457       fitters[i]->GetParameters(params[i]);
458       fitters[i]->GetErrors(errors[i]);
459       // renormalize errors
460       errors[i][0]*=chi2Fac;
461       errors[i][1]*=chi2Fac;
462     }
463   }
464   //
465   // Fill fitters
466   //
467   for (Int_t i0=0;i0<4;i0++){
468     for (Int_t i1=0;i1<4;i1++){
469       if (i0==i1) continue;
470       if(npoints[i0]<kMinClusterQ) continue;
471       if(npoints[i1]<kMinClusterQ) continue;
472       Int_t index0=i0*4+i1;
473       Double_t xy[1];
474       Double_t xfi[1];
475       xy[0] = pyf[1];   
476       xfi[0] = pyf[1];
477       //
478       Double_t dy = params[i1][0]-params[i0][0];
479       Double_t sy = TMath::Sqrt(errors[i1][0]*errors[i1][0]+errors[i0][0]*errors[i0][0]);
480       Double_t dphi = params[i1][1]-params[i0][1];
481       Double_t sphi = TMath::Sqrt(errors[i1][1]*errors[i1][1]+errors[i0][1]*errors[i0][1]);
482       //
483       Int_t sector = isec-36;
484       //
485       if (TMath::Abs(dy/(sy+0.1))>5.) continue;          // 5 sigma cut
486       if (TMath::Abs(dphi/(sphi+0.004))>5.) continue;    // 5 sigma cut 
487       TLinearFitter * fitterY,*fitterPhi;
488       fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*sector+index0);
489       if (TMath::Abs(dy/(sy+0.1))<5.)  fitterY->AddPoint(xy,dy,sy);
490       fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*36+index0);
491       if (TMath::Abs(dy/(sy+0.1))<5.)  fitterY->AddPoint(xy,dy,sy);
492       //
493       fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*sector+index0);
494       if (TMath::Abs(dphi/(sphi+0.001))<5.)  fitterPhi->AddPoint(xfi,dphi,sphi);
495       fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*36+index0);
496       if (TMath::Abs(dphi/(sphi+0.001))<5.)  fitterPhi->AddPoint(xfi,dphi,sphi);
497     }
498   }
499   //
500   // Dump debug information
501   //
502   if (fStreamLevel>0){    
503     TTreeSRedirector *cstream = GetDebugStreamer();      
504     Int_t sector = isec-36;
505     if (cstream){      
506       for (Int_t i0=0;i0<4;i0++){
507         for (Int_t i1=i0;i1<4;i1++){
508           if (i0==i1) continue;
509           if(npoints[i0]<kMinClusterQ) continue;
510           if(npoints[i1]<kMinClusterQ) continue;
511           (*cstream)<<"fitQ"<<
512             "run="<<fRun<<              //  run number
513             "event="<<fEvent<<          //  event number
514             "time="<<fTime<<            //  time stamp of event
515             "trigger="<<fTrigger<<      //  trigger
516             "mag="<<fMagF<<             //  magnetic field        
517             "sec="<<sector<<            // current sector from 0 to 36
518             "isec="<<isec<<             //  current sector
519             // full fit
520             "nf="<<nf<<                 //  total number of points
521             "pyf.="<<&pyf<<             //  full OROC fit y
522             "pzf.="<<&pzf<<             //  full OROC fit z
523             // quadrant fit
524             "i0="<<i0<<                 // quadrant number
525             "i1="<<i1<<                 
526             "n0="<<npoints[i0]<<        // number of points
527             "n1="<<npoints[i1]<<
528             "p0.="<<&params[i0]<<       // parameters
529             "p1.="<<&params[i1]<< 
530             "e0.="<<&errors[i0]<<       // errors
531             "e1.="<<&errors[i1]<<
532             "chi0="<<chi2C[i0]<<       // chi2s
533             "chi1="<<chi2C[i1]<<
534
535             "\n";
536         }    
537       }
538     }
539   }
540 }
541
542
543
544
545
546 void AliTPCcalibUnlinearity::MakeHisto(){
547   //
548   //
549   //
550   TString  axisName[10];
551   Double_t xmin[10],  xmax[10];
552   Int_t    nbins[10];
553   //
554   //
555   axisName[0]="sector";
556   xmin[0]=0; xmax[0]=72; nbins[0]=72;
557   //          
558   axisName[1]="row";
559   xmin[1]=0; xmax[1]=96; nbins[1]=96;
560   //
561   axisName[2]="tphi";            // tan phi - ly/lx 
562   xmin[2]=-0.17; xmax[2]=0.17; nbins[2]=5;
563   //  
564   axisName[3]="lz";              // global z        
565   xmin[3]=-250; xmax[3]=250; nbins[3]=10;
566   //
567   axisName[4]="k";              // tangent - in angle
568   xmin[4]=-0.5; xmax[4]=0.5; nbins[4]=10;
569   //
570   //
571   axisName[5]="delta";          // delta
572   xmin[5]=-0.5; xmax[5]=0.5; nbins[5]=100;
573   //
574   //
575   fDiffHistoLineY = new THnSparseF("hnDistHistoLineY","hnDistHistoLineY",6,nbins,xmin,xmax);
576   fDiffHistoParY = new THnSparseF("hnDistHistoParY","hnDistHistoParY",6,nbins,xmin,xmax);
577   fDiffHistoLineZ = new THnSparseF("hnDistHistoLineZ","hnDistHistoLineZ",6,nbins,xmin,xmax);
578   fDiffHistoParZ = new THnSparseF("hnDistHistoParZ","hnDistHistoParZ",6,nbins,xmin,xmax);
579  //  for (Int_t i=0; i<8;i++){
580 //     fDiffHistoLineY->GetAxis(i)->SetName(axisName[i].Data());
581 //     fDiffHistoParY->GetAxis(i)->SetName(axisName[i].Data());
582 //     fDiffHistoLineY->GetAxis(i)->SetTitle(axisName[i].Data());
583 //     fDiffHistoParY->GetAxis(i)->SetTitle(axisName[i].Data());
584 //     fDiffHistoLineZ->GetAxis(i)->SetName(axisName[i].Data());
585 //     fDiffHistoParZ->GetAxis(i)->SetName(axisName[i].Data());
586 //     fDiffHistoLineZ->GetAxis(i)->SetTitle(axisName[i].Data());
587 //     fDiffHistoParZ->GetAxis(i)->SetTitle(axisName[i].Data());
588 //   }
589   
590   //
591   // 
592   //
593   char hname[1000];
594   TH2F * his=0;
595   for (Int_t isec=0;isec<74;isec++){
596     sprintf(hname,"DeltaRPhiPlus_Sector%d",isec);
597     his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
598     his->SetDirectory(0);
599     fDistRPHIPlus.AddAt(his,isec);
600     sprintf(hname,"DeltaRPhiMinus_Sector%d",isec);
601     his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
602     his->SetDirectory(0);
603     fDistRPHIMinus.AddAt(his,isec);
604   }
605 }
606
607
608 void AliTPCcalibUnlinearity::Terminate(){
609   //
610   // Terminate function
611   // call base terminate + Eval of fitters
612   //
613   Info("AliTPCcalibUnlinearities","Terminate");
614   EvalFitters();
615   DumpTree();
616   AliTPCcalibBase::Terminate();
617 }
618
619 Long64_t AliTPCcalibUnlinearity::Merge(TCollection *list) {
620   //
621   // merge results
622   //
623   TIterator* iter = list->MakeIterator();
624   AliTPCcalibUnlinearity* cal = 0;
625   //
626   while ((cal = (AliTPCcalibUnlinearity*)iter->Next())) {
627     if (!cal->InheritsFrom(AliTPCcalibUnlinearity::Class())) {
628       return -1;
629     }    
630     Add(cal);
631   }
632   return 0;  
633 }
634
635 void    AliTPCcalibUnlinearity::Add(AliTPCcalibUnlinearity * calib){
636   //
637   //
638   //
639   if (!fInit) Init();
640   if (fDiffHistoLineY && calib->fDiffHistoLineY){
641     fDiffHistoLineY->Add(calib->fDiffHistoLineY);
642   }
643   if (fDiffHistoParY && calib->fDiffHistoParY){
644     fDiffHistoParY->Add(calib->fDiffHistoParY);
645   }
646   if (fDiffHistoLineZ && calib->fDiffHistoLineZ){
647     fDiffHistoLineZ->Add(calib->fDiffHistoLineZ);
648   }
649   if (fDiffHistoParZ && calib->fDiffHistoParZ){
650     fDiffHistoParZ->Add(calib->fDiffHistoParZ);
651   }
652   for (Int_t isec=0;isec<38;isec++){
653     TLinearFitter *f0r = (TLinearFitter*)fFittersOutR.At(isec);
654     TLinearFitter *f1r = (TLinearFitter*)(calib->fFittersOutR.At(isec));
655     if (f0r&&f1r) f0r->Add(f1r);
656     TLinearFitter *f0z = (TLinearFitter*)fFittersOutZ.At(isec);
657     TLinearFitter *f1z = (TLinearFitter*)(calib->fFittersOutZ.At(isec));
658     if (f0z&&f1z) f0z->Add(f1z);
659   }
660
661   for (Int_t isec=0;isec<16*38;isec++){
662     TLinearFitter *f0y = (TLinearFitter*)fFitterQuadrantY.At(isec);
663     TLinearFitter *f1y = (TLinearFitter*)(calib->fFitterQuadrantY.At(isec));
664     if (f0y&&f1y) f0y->Add(f1y);
665     TLinearFitter *f0phi = (TLinearFitter*)fFitterQuadrantPhi.At(isec);
666     TLinearFitter *f1phi = (TLinearFitter*)(calib->fFitterQuadrantPhi.At(isec));
667     if (f0phi&&f1phi) f0phi->Add(f1phi);
668   }
669   
670   
671   for (Int_t isec=0;isec<74;isec++){
672     TH2F * h0p= (TH2F*)(fDistRPHIPlus.At(isec));
673     TH2F * h1p= (TH2F*)(calib->fDistRPHIPlus.At(isec));
674     if (h0p&&h1p) h0p->Add(h1p);
675     TH2F * h0m= (TH2F*)(fDistRPHIMinus.At(isec));
676     TH2F * h1m= (TH2F*)(calib->fDistRPHIMinus.At(isec));
677     if (h0m&&h1m) h0m->Add(h1m);
678   }
679
680  
681 }
682
683
684
685 void AliTPCcalibUnlinearity::DumpTree(const char *fname){
686   //
687   //
688   // 
689   TTreeSRedirector *cstream =new TTreeSRedirector(fname);
690   if (!cstream) return;
691   //  
692   THnSparse *his=0;  
693   Double_t position[10];
694   Double_t value; 
695   Int_t *bins = new Int_t[10];
696   //
697   for (Int_t ihis=0; ihis<2; ihis++){
698     his =  (ihis==0)? fDiffHistoLineY:fDiffHistoParY;
699     if (!his) continue;
700     //
701     Int_t ndim = his->GetNdimensions();
702     //
703     for (Long64_t i = 0; i < his->GetNbins(); ++i) {
704       value = his->GetBinContent(i, bins);
705       for (Int_t idim = 0; idim < ndim; idim++) {
706         position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
707       }      
708       (*cstream)<<"Resol"<<
709         "ihis="<<ihis<<
710         "bincont="<<value;
711       for (Int_t idim=0;idim<ndim;idim++){
712         (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
713       }
714       (*cstream)<<"Resol"<<
715         "\n";      
716     }
717   }
718   delete cstream;
719 }
720
721
722 void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t ky, Double_t kz, Int_t npoints){
723   //
724   //
725   // Simple distortion fit in outer sectors
726   //
727   // Model - 2 exponential decay toward the center of chamber
728   //       - Decay length 10 and 5 cm
729   //       - Decay amplitude - Z dependent 
730   //
731  
732   Double_t side   = (-1+((sector%36)<18)*2); // A side +1   C side -1
733   Double_t dzt    = (cz-tz)*side;
734   Double_t radius = TMath::Sqrt(cx*cx+cy*cy);  
735   AliTPCROC * roc = AliTPCROC::Instance();
736   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
737   Double_t dout = xout-radius;
738   if (dout>30) return;
739   //drift
740   Double_t dr   = 0.5 - TMath::Abs(cz/250.);
741   Double_t dy = cy-ty;
742   Double_t eout10 = TMath::Exp(-dout/10.);
743   Double_t eout5 = TMath::Exp(-dout/5.);
744   //
745   Double_t eoutp  = (eout10+eout5)*0.5;
746   Double_t eoutm  = (eout10-eout5)*0.5;
747
748   //
749   Double_t xxxR[6], xxxZ[6], xxxRZ[6];
750   //TString fstring="";
751   xxxZ[0]=eoutp;                //fstring+="eoutp++";  
752   xxxZ[1]=eoutm;               //fstring+="eoutm++";  
753   xxxZ[2]=dr*eoutp;            //fstring+="dr*eoutp++";  
754   xxxZ[3]=dr*eoutm;            //fstring+="dr*eoutm++";  
755   xxxZ[4]=dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
756   xxxZ[5]=dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
757   //
758   xxxR[0]=ky*eoutp;                //fstring+="eoutp++";  
759   xxxR[1]=ky*eoutm;               //fstring+="eoutm++";  
760   xxxR[2]=ky*dr*eoutp;            //fstring+="dr*eoutp++";  
761   xxxR[3]=ky*dr*eoutm;            //fstring+="dr*eoutm++";  
762   xxxR[4]=ky*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
763   xxxR[5]=ky*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
764   //
765   xxxRZ[0]=kz*eoutp;                //fstring+="eoutp++";  
766   xxxRZ[1]=kz*eoutm;               //fstring+="eoutm++";  
767   xxxRZ[2]=kz*dr*eoutp;            //fstring+="dr*eoutp++";  
768   xxxRZ[3]=kz*dr*eoutm;            //fstring+="dr*eoutm++";  
769   xxxRZ[4]=kz*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
770   xxxRZ[5]=kz*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
771   //
772   TLinearFitter * fitter=0;
773   Double_t err=0.1;
774   for (Int_t ip=0; ip<npoints; ip++){
775     //
776     fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
777     fitter->AddPoint(xxxR,dy,err);
778     //fitter->AddPoint(xxxRZ,dz);
779     fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
780     fitter->AddPoint(xxxZ,dzt,err);
781     //
782     if (side>0){
783       fitter = (TLinearFitter*)fFittersOutR.At(36);
784       fitter->AddPoint(xxxR,dy,err);
785       //fitter->AddPoint(xxxRZ,dz);
786       fitter = (TLinearFitter*)fFittersOutZ.At(36);
787       fitter->AddPoint(xxxZ,dzt,err);
788     }
789     if (side<0){
790       fitter = (TLinearFitter*)fFittersOutR.At(37);
791       fitter->AddPoint(xxxR,dy,err);
792       //fitter->AddPoint(xxxRZ,dz);
793       fitter = (TLinearFitter*)fFittersOutZ.At(37);
794       fitter->AddPoint(xxxZ,dzt,err);
795     }
796   }
797 }
798 void AliTPCcalibUnlinearity::AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz,  Double_t /*ky*/, Double_t /*kz*/, Int_t npoints){
799   //
800   //
801   //
802   Float_t kMaxDz = 0.5;  // cut on z in cm
803   Double_t dy = cy-ty;
804   Double_t dz = cz-tz;
805   if (TMath::Abs(dz)>kMaxDz) return;
806   Double_t edgeMinus= -ty+cx*TMath::Tan(TMath::Pi()/18.);
807   Double_t edgePlus =  ty+cx*TMath::Tan(TMath::Pi()/18.);
808   //
809   TH2F * his =0;
810   his = (TH2F*)fDistRPHIPlus.At(sector);
811   his->Fill(edgePlus,dy,npoints);
812   his = (TH2F*)((sector<36)? fDistRPHIPlus.At(72):fDistRPHIPlus.At(73));
813   his->Fill(edgePlus,dy,npoints);
814   his = (TH2F*)fDistRPHIMinus.At(sector);
815   his->Fill(edgeMinus,dy,npoints);
816   his = (TH2F*)((sector<36)? fDistRPHIMinus.At(72):fDistRPHIMinus.At(73));
817   his->Fill(edgeMinus,dy,npoints);
818 }
819
820
821
822 void AliTPCcalibUnlinearity::Init(){
823   //
824   //
825   //
826   //
827   MakeHisto();
828   // Make outer fitters
829   TLinearFitter * fitter = 0;
830   for (Int_t ifit=0; ifit<38; ifit++){
831     fitter = new TLinearFitter(7,"hyp6");
832     fitter->StoreData(kFALSE);
833     fFittersOutR.AddAt(fitter,ifit);
834     fitter = new TLinearFitter(7,"hyp6");
835     fitter->StoreData(kFALSE);
836     fFittersOutZ.AddAt(fitter,ifit);
837   }
838   for (Int_t ifit=0; ifit<16*38;ifit++){
839     fitter = new TLinearFitter(2,"hyp1");
840     fitter->StoreData(kFALSE);
841     fFitterQuadrantY.AddAt(fitter,ifit);
842     fitter = new TLinearFitter(2,"hyp1");
843     fitter->StoreData(kFALSE);
844     fFitterQuadrantPhi.AddAt(fitter,ifit);    
845   }  
846   fInit= kTRUE;
847 }
848
849 void AliTPCcalibUnlinearity::EvalFitters(){
850   //
851   // 
852   // Evaluate fitters
853   // 
854   TLinearFitter * fitter = 0;
855   TVectorD vec;
856   for (Int_t ifit=0; ifit<38; ifit++){
857     fitter=(TLinearFitter*)fFittersOutR.At(ifit);
858     if (fitter) {
859       fitter->Eval();
860       fitter->GetParameters(vec);
861       fParamsOutR.AddAt(vec.Clone(),ifit);
862       fitter->GetErrors(vec);
863       fErrorsOutR.AddAt(vec.Clone(),ifit);
864     }
865     fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
866     if (fitter) {
867       fitter->Eval();
868       fitter->GetParameters(vec);
869       fParamsOutZ.AddAt(vec.Clone(),ifit);
870       fitter->GetErrors(vec);
871       fErrorsOutZ.AddAt(vec.Clone(),ifit);
872     }
873   }
874 }
875
876
877
878
879
880 void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
881   //
882   //
883   //
884   //  TTree * tree = chainUnlinP;
885   AliTPCcalibUnlinearity *calib = this;
886   //
887   AliTPCclusterMI * cl=0;
888   Double_t ty,tz;
889   TVectorD *vy=0, *vz=0;
890   TVectorD *vy2=0, *vz2=0;
891   Long64_t nentries = tree->GetEntries();
892   if (nmaxPoints>0) nentries = TMath::Min(nentries,nmaxPoints);
893   //
894   {
895   for (Long64_t i=0; i<nentries; i++){
896     tree->SetBranchAddress("Cl.",&cl);
897     tree->SetBranchAddress("y1",&ty);
898     tree->SetBranchAddress("z1",&tz);
899     tree->SetBranchAddress("py1.",&vy);
900     tree->SetBranchAddress("pz1.",&vz);
901     tree->SetBranchAddress("py2.",&vy2);
902     tree->SetBranchAddress("pz2.",&vz2);
903     //
904     tree->GetEntry(i);
905     if (!cl) continue;
906     if (i%10000==0) printf("%d\n",(Int_t)i);
907     Int_t row= cl->GetRow();
908     if (cl->GetDetector()>36) row+=64;
909     calib->AddPointRPHI(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
910                       ty,tz,(*vy)[1],(*vz)[1],1);
911
912     if (cl->GetType()<0) continue; 
913     Double_t dy = cl->GetY()-ty;
914     Double_t dz = cl->GetZ()-tz;
915     //
916     //
917     if (TMath::Abs(dy)>0.25) continue;
918     if (TMath::Abs(dz)>0.25) continue;
919     
920     if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
921     if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
922       // Apply sagita cut
923     if (cl->GetType()>=0) {
924       calib->AddPoint(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
925                       ty,tz,(*vy)[1],(*vz)[1],1);
926     }
927   }
928   }
929 }
930
931
932 void  AliTPCcalibUnlinearity::MakeQPosNormAll(TTree * chainUnlinD, AliTPCClusterParam * param, Int_t maxPoints){
933   //
934   // Make position corrections
935   // for the moment Only using debug streamer 
936   // chainUnlinD  - debug tree
937   // param     - parameters to be updated
938   // maxPoints - maximal number of points using for fit
939   // verbose   - print info flag
940   //
941   // Current implementation - only using debug streamers
942   // 
943   
944   /*    
945   //Defaults
946   Int_t maxPoints=100000;
947   */
948   /*
949     Usage: 
950     //0. Load libraries
951     gSystem->Load("libANALYSIS");
952     gSystem->Load("libSTAT");
953     gSystem->Load("libTPCcalib");
954     
955
956     //1. Load Parameters to be modified:
957     //e.g:
958     
959     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
960     AliCDBManager::Instance()->SetRun(0) 
961     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
962     //
963     //2. Load the debug streamer
964     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
965     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
966     AliXRDPROOFtoolkit tool;  
967     TChain * chainUnlin = tool.MakeChain("unlin.txt","Resol",0,10200);
968     chainUnlin->Lookup();
969     TChain * chainUnlinD = tool.MakeChain("unlin.txt","Diff",0,10200);
970     chainUnlinD->Lookup();
971
972     //3. Do fits and store results
973     // 
974     AliTPCcalibUnlinearity::MakeQPosNormAll(chainUnlinD,param,200000,0) ;
975     TFile f("paramout.root","recreate");
976     param->Write("clusterParam");
977     f.Close();
978     //4. Verification
979     TFile f2("paramout.root");
980     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
981     param2->SetInstance(param2);
982     chainUnlinD->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line 
983     
984    */
985
986
987   TStatToolkit toolkit;
988   Double_t chi2;
989   TVectorD fitParamY0;
990   TVectorD fitParamY1;
991   TVectorD fitParamZ0;
992   TVectorD fitParamZ1;
993   TMatrixD covMatrix;
994   Int_t npoints;
995   
996   chainUnlinD->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
997   chainUnlinD->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
998   chainUnlinD->SetAlias("sp","(sin(dp*pi)-dp*pi)");
999   chainUnlinD->SetAlias("st","(sin(dt)-dt)");
1000   chainUnlinD->SetAlias("dq","(-1+(Cl.fZ>0)*2)/Cl.fMax");
1001   //
1002   chainUnlinD->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1003   //
1004   //
1005   //
1006   TCut cutE("Cl.fType>=0");
1007   TCut cutDy("abs(Cl.fY-y1)<0.4");
1008   TCut cutDz("abs(Cl.fZ-z1)<0.4");
1009   TCut cutY2("abs(py2.fElements[2])*150^2/4<1");
1010   TCut cutZ2("abs(pz2.fElements[2])*150^2/4<1");
1011   TCut cutAll = cutE+cutDy+cutDz+cutY2+cutZ2;
1012   //
1013   TCut cutA("Cl.fZ>0");
1014   TCut cutC("Cl.fZ<0");
1015
1016   TString fstringY="";  
1017   //
1018   fstringY+="(dp)++";            //1
1019   fstringY+="(dp)*di++";         //2
1020   fstringY+="(sp)++";            //3
1021   fstringY+="(sp)*di++";         //4
1022   fstringY+="(dq)++";            //5
1023   TString fstringZ="";  
1024   fstringZ+="(dt)++";            //1
1025   fstringZ+="(dt)*di++";         //2
1026   fstringZ+="(st)++";            //3
1027   fstringZ+="(st)*di++";         //4
1028   fstringZ+="(dq)++";            //5
1029   //
1030   // Z corrections
1031   //
1032   TString *strZ0 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
1033   printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1034   strZ0->Tokenize("++")->Print();
1035   param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone();
1036   //
1037   TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
1038   //
1039   printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1040   strZ1->Tokenize("++")->Print();
1041   param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone();
1042   param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone();
1043   //
1044   // Y corrections
1045   //   
1046   TString *strY0 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
1047   printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); 
1048   strY0->Tokenize("++")->Print();
1049   param->PosYcor(0) = (TVectorD*) fitParamY0.Clone();
1050   
1051
1052   TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
1053   //
1054   printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));  
1055   strY1->Tokenize("++")->Print();
1056   param->PosYcor(1) = (TVectorD*) fitParamY1.Clone();
1057   param->PosYcor(2) = (TVectorD*) fitParamY1.Clone();
1058   //
1059   //
1060   //
1061   chainUnlinD->SetAlias("fitZ0",strZ0->Data());
1062   chainUnlinD->SetAlias("fitZ1",strZ1->Data());
1063   chainUnlinD->SetAlias("fitY0",strY0->Data());
1064   chainUnlinD->SetAlias("fitY1",strY1->Data());
1065 }
1066
1067
1068
1069
1070
1071
1072