Extacting the OCDB in a separate module. The detectors have write permission in the...
[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
36   //  
37 */
38
39
40 #include "TLinearFitter.h"
41
42 #include "Riostream.h"
43 #include "TChain.h"
44 #include "TTree.h"
45 #include "TH1F.h"
46 #include "TH2F.h"
47 #include "TH3F.h"
48 #include "THnSparse.h"
49 #include "TList.h"
50 #include "TMath.h"
51 #include "TCanvas.h"
52 #include "TFile.h"
53 #include "TF1.h"
54 #include "TVectorD.h"
55 #include "TProfile.h"
56 #include "TGraphErrors.h"
57 #include "TCanvas.h"
58 #include "TCut.h"
59
60
61 #include "AliTPCclusterMI.h"
62 #include "AliTPCseed.h"
63 #include "AliESDVertex.h"
64 #include "AliESDEvent.h"
65 #include "AliESDfriend.h"
66 #include "AliESDInputHandler.h"
67 #include "AliAnalysisManager.h"
68 #include "AliTracker.h"
69 #include "AliMagF.h"
70 #include "AliTPCCalROC.h"
71
72 #include "AliLog.h"
73
74
75 #include "TTreeStream.h"
76 #include "AliTPCTracklet.h"
77 #include "TTimeStamp.h"
78 #include "AliTPCcalibDB.h"
79 #include "AliTPCcalibLaser.h"
80 #include "AliDCSSensorArray.h"
81 #include "AliDCSSensor.h"
82 #include "TLinearFitter.h"
83 #include "AliTPCClusterParam.h"
84 #include "TStatToolkit.h"
85
86
87 #include "AliTPCcalibUnlinearity.h"
88
89 AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::fgInstance = 0;
90
91 ClassImp(AliTPCcalibUnlinearity)
92
93 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
94   AliTPCcalibBase(),
95   fDiffHistoLine(0),
96   fDiffHistoPar(0),
97   fFittersOutR(38),
98   fFittersOutZ(38),
99   fParamsOutR(38),
100   fParamsOutZ(38),
101   fErrorsOutR(38),
102   fErrorsOutZ(38)
103 {
104   //
105   // Defualt constructor
106   //
107 }
108
109 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
110   AliTPCcalibBase(name,title),
111   fDiffHistoLine(0),
112   fDiffHistoPar(0),  
113   fFittersOutR(38),
114   fFittersOutZ(38),
115   fParamsOutR(38),
116   fParamsOutZ(38),
117   fErrorsOutR(38),
118   fErrorsOutZ(38)
119 {
120   //
121   // Non default constructor
122   //
123   MakeFitters();
124 }
125
126 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
127   //
128   //
129   //
130   if (fDiffHistoLine) delete fDiffHistoLine;
131   if (fDiffHistoPar)  delete fDiffHistoPar;
132 }
133
134
135 AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::Instance()
136 {
137   //
138   // Singleton implementation
139   // Returns an instance of this class, it is created if neccessary
140   //
141   if (fgInstance == 0){
142     fgInstance = new AliTPCcalibUnlinearity();
143   }
144   return fgInstance;
145 }
146
147
148
149
150 void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
151   //
152   // 
153   //  
154   const Int_t  kdrow=10;
155   const Int_t kMinCluster=35;
156   if (!fDiffHistoLine) MakeHisto();    
157   if (!seed) return;
158   if (TMath::Abs(fMagF)>0.01) return;   // use only data without field
159   Int_t counter[72];
160   for (Int_t i=0;i<72;i++) counter[i]=0;
161   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
162     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
163     if (!cluster) continue;
164     if (cluster->GetType()<0) continue;
165     counter[cluster->GetDetector()]++;
166   }
167   //
168   for (Int_t is=0; is<72;is++){
169     if (counter[is]>kMinCluster) ProcessDiff(seed, is);
170   }
171 }
172
173
174 void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
175   //
176   // process differnce of the cluster in respect with linear and parabolic fit
177   //
178   const  Double_t kXmean=134;
179   const Int_t     kdrow=10;
180   static TLinearFitter fy1(2,"pol1");
181   static TLinearFitter fy2(3,"pol2");
182   static TLinearFitter fz1(2,"pol1");
183   static TLinearFitter fz2(3,"pol2");
184   //
185   static TVectorD py1(2);
186   static TVectorD py2(3);
187   static TVectorD pz1(2);
188   static TVectorD pz2(3);  
189   //  
190   fy1.ClearPoints();
191   fy2.ClearPoints();
192   fz1.ClearPoints();
193   fz2.ClearPoints();
194   //
195   for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
196     AliTPCclusterMI *c=track->GetClusterPointer(irow);
197     if (!c) continue;
198     if (c->GetDetector()!=isec) continue;
199     if (c->GetType()<0) continue;
200     Double_t dx = c->GetX()-kXmean;
201     Double_t x[2]={dx, dx*dx};
202     fy1.AddPoint(x,c->GetY());
203     fy2.AddPoint(x,c->GetY());
204     fz1.AddPoint(x,c->GetZ());
205     fz2.AddPoint(x,c->GetZ());
206   }
207   fy1.Eval(); fy1.GetParameters(py1);
208   fy2.Eval(); fy2.GetParameters(py2);
209   fz1.Eval(); fz1.GetParameters(pz1);
210   fz2.Eval(); fz2.GetParameters(pz2);
211   //
212   
213   for (Int_t irow=0;irow<159;irow++) {
214     AliTPCclusterMI *c=track->GetClusterPointer(irow);
215     if (!c) continue;
216     if (c->GetDetector()!=isec) continue;
217     if (c->GetType()<0) continue;
218     Double_t dx = c->GetX()-kXmean;
219     Double_t y1 = py1[0]+py1[1]*dx;
220     Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
221     //
222     Double_t z1 = pz1[0]+pz1[1]*dx;
223     Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
224     //
225     //
226     Double_t x[10];
227     x[0]=isec;
228     x[1]=irow;
229     x[2]=c->GetY();
230     x[3]=c->GetZ();
231     x[4]=py1[1];
232     x[5]=pz1[1];
233     x[6]=py2[2]*150*150/4; // sagita 150 cm
234     //
235     x[7]=c->GetY()-y1;
236     x[8]=c->GetZ()-z1;
237     fDiffHistoLine->Fill(x);
238     x[7]=c->GetY()-y2;
239     x[8]=c->GetZ()-z2;
240     fDiffHistoPar->Fill(x);   
241
242     if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(pz2[2]*150*150/4)<1){
243       // Apply sagita cut
244       AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
245                py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
246     }
247              
248     if (fStreamLevel>1){
249       TTreeSRedirector *cstream = GetDebugStreamer();
250       if (cstream){
251         (*cstream)<<"Diff"<<
252           "run="<<fRun<<              //  run number
253           "event="<<fEvent<<          //  event number
254           "time="<<fTime<<            //  time stamp of event
255           "trigger="<<fTrigger<<      //  trigger
256           "mag="<<fMagF<<             //  magnetic field          
257           "isec="<<isec<<
258           "Cl.="<<c<<
259           "y1="<<y1<<
260           "y2="<<y2<<
261           "z1="<<z1<<
262           "z2="<<z2<<
263           "py1.="<<&py1<<
264           "py2.="<<&py2<<
265           "pz1.="<<&pz1<<
266           "pz2.="<<&pz2<<
267           "\n";
268       }
269     }
270   }
271 }
272
273
274 void AliTPCcalibUnlinearity::MakeHisto(){
275   //
276   //
277   //
278   TString  axisName[10];
279   Double_t xmin[10],  xmax[10];
280   Int_t    nbins[10];
281   //
282   //
283   axisName[0]="sector";
284   xmin[0]=0; xmax[0]=72; nbins[0]=72;
285   //
286   axisName[1]="row";
287   xmin[1]=0; xmax[1]=159; nbins[1]=159;
288   //
289   axisName[2]="ly";
290   xmin[2]=-50; xmax[2]=50; nbins[2]=10;
291   //
292   axisName[3]="lz";
293   xmin[3]=-250; xmax[3]=250; nbins[3]=50;
294   //
295   axisName[4]="p2";
296   xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
297   //
298   axisName[5]="p3";
299   xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
300   //
301   axisName[6]="p4";
302   xmin[6]=-2; xmax[6]=2; nbins[6]=20;
303   //
304   axisName[7]="dy";
305   xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
306   //
307   axisName[8]="dz";
308   xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
309   //
310   fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
311   fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
312   for (Int_t i=0; i<9;i++){
313     fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
314     fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
315     fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
316     fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
317   }
318 }
319
320
321 void AliTPCcalibUnlinearity::Terminate(){
322   //
323   // Terminate function
324   // call base terminate + Eval of fitters
325   //
326   Info("AliTPCcalibUnlinearities","Terminate");
327   EvalFitters();
328   DumpTree();
329   AliTPCcalibBase::Terminate();
330 }
331
332
333 void AliTPCcalibUnlinearity::DumpTree(){
334   //
335   //
336   // 
337   if (fStreamLevel==0) return;
338   TTreeSRedirector *cstream = GetDebugStreamer();
339   if (!cstream) return;
340   //  
341   THnSparse *his=0;  
342   Double_t position[10];
343   Double_t value; 
344   Int_t *bins = new Int_t[10];
345   //
346   for (Int_t ihis=0; ihis<2; ihis++){
347     his =  (ihis==0)? fDiffHistoLine:fDiffHistoPar;
348     if (!his) continue;
349     //
350     Int_t ndim = his->GetNdimensions();
351     //
352     for (Long64_t i = 0; i < his->GetNbins(); ++i) {
353       value = his->GetBinContent(i, bins);
354       for (Int_t idim = 0; idim < ndim; idim++) {
355         position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
356       }      
357       (*cstream)<<"Resol"<<
358         "ihis="<<ihis<<
359         "bincont="<<value;
360       for (Int_t idim=0;idim<ndim;idim++){
361         (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
362       }
363       (*cstream)<<"Resol"<<
364         "\n";      
365     }
366   }
367 }
368
369
370 void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints){
371   //
372   //
373   // Simple distortion fit in outer sectors
374   //
375   // Model - 2 exponential decay toward the center of chamber
376   //       - Decay length 10 and 5 cm
377   //       - Decay amplitude - Z dependent 
378   //
379   /*
380   chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
381   chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
382   chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
383   chainUnlin->SetAlias("dzt","dz*side");
384   chainUnlin->SetAlias("dr","(1-abs(lz/250))");
385   chainUnlin->SetAlias("dout","(159-row)*1.5");
386   chainUnlin->SetAlias("din","row*0.75");
387   chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
388   chainUnlin->SetAlias("eout5","exp(-(dout)/5.)");  
389   */
390
391   Double_t side   = (-1+((sector%36)<18)*2); // A side +1   C side -1
392   Double_t dzt    = dz*side;
393   Double_t dout   = (159-row)*1.5;  // distance to the last pad row - (valid only for long pads)
394   if (dout>30) return; // use only edge pad rows
395   dr-=0.5;             // dr shifted to the middle - reduce numerical instabilities
396
397   Double_t eout10 = TMath::Exp(-dout/10.);
398   Double_t eout5 = TMath::Exp(-dout/5.);
399   //
400   Double_t eoutp  = (eout10+eout5)*0.5;
401   Double_t eoutm  = (eout10-eout5)*0.5;
402
403   //
404   Double_t xxxR[6], xxxZ[6], xxxRZ[6];
405   //TString fstring="";
406   xxxZ[0]=eoutp;                //fstring+="eoutp++";  
407   xxxZ[1]=eoutm;               //fstring+="eoutm++";  
408   xxxZ[2]=dr*eoutp;            //fstring+="dr*eoutp++";  
409   xxxZ[3]=dr*eoutm;            //fstring+="dr*eoutm++";  
410   xxxZ[4]=dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
411   xxxZ[5]=dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
412   //
413   xxxR[0]=p2*eoutp;                //fstring+="eoutp++";  
414   xxxR[1]=p2*eoutm;               //fstring+="eoutm++";  
415   xxxR[2]=p2*dr*eoutp;            //fstring+="dr*eoutp++";  
416   xxxR[3]=p2*dr*eoutm;            //fstring+="dr*eoutm++";  
417   xxxR[4]=p2*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
418   xxxR[5]=p2*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
419   //
420   xxxRZ[0]=p3*eoutp;                //fstring+="eoutp++";  
421   xxxRZ[1]=p3*eoutm;               //fstring+="eoutm++";  
422   xxxRZ[2]=p3*dr*eoutp;            //fstring+="dr*eoutp++";  
423   xxxRZ[3]=p3*dr*eoutm;            //fstring+="dr*eoutm++";  
424   xxxRZ[4]=p3*dr*dr*eoutp;         //fstring+="dr*dr*eoutp++";  
425   xxxRZ[5]=p3*dr*dr*eoutm;         //fstring+="dr*dr*eoutm++";    
426   //
427   TLinearFitter * fitter=0;
428   Double_t err=0.1;
429   for (Int_t ip=0; ip<npoints; ip++){
430     //
431     fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
432     fitter->AddPoint(xxxR,dy,err);
433     //fitter->AddPoint(xxxRZ,dz);
434     fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
435     fitter->AddPoint(xxxZ,dzt,err);
436     //
437     if (side>0){
438       fitter = (TLinearFitter*)fFittersOutR.At(36);
439       fitter->AddPoint(xxxR,dy,err);
440       //fitter->AddPoint(xxxRZ,dz);
441       fitter = (TLinearFitter*)fFittersOutZ.At(36);
442       fitter->AddPoint(xxxZ,dzt,err);
443     }
444     if (side<0){
445       fitter = (TLinearFitter*)fFittersOutR.At(37);
446       fitter->AddPoint(xxxR,dy,err);
447       //fitter->AddPoint(xxxRZ,dz);
448       fitter = (TLinearFitter*)fFittersOutZ.At(37);
449       fitter->AddPoint(xxxZ,dzt,err);
450     }
451   }
452 }
453
454
455 void AliTPCcalibUnlinearity::MakeFitters(){
456   //
457   //
458   //
459   // Make outer fitters
460   TLinearFitter * fitter = 0;
461   for (Int_t ifit=0; ifit<38; ifit++){
462     fitter = new TLinearFitter(7,"hyp6");
463     fitter->StoreData(kFALSE);
464     fFittersOutR.AddAt(fitter,ifit);
465     fitter = new TLinearFitter(7,"hyp6");
466     fitter->StoreData(kFALSE);
467     fFittersOutZ.AddAt(fitter,ifit);
468   }
469 }
470
471 void AliTPCcalibUnlinearity::EvalFitters(){
472   //
473   // 
474   // Evaluate fitters
475   // 
476   TLinearFitter * fitter = 0;
477   TVectorD vec;
478   for (Int_t ifit=0; ifit<38; ifit++){
479     fitter=(TLinearFitter*)fFittersOutR.At(ifit);
480     if (fitter) {
481       fitter->Eval();
482       fitter->GetParameters(vec);
483       fParamsOutR.AddAt(vec.Clone(),ifit);
484       fitter->GetErrors(vec);
485       fErrorsOutR.AddAt(vec.Clone(),ifit);
486     }
487     fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
488     if (fitter) {
489       fitter->Eval();
490       fitter->GetParameters(vec);
491       fParamsOutZ.AddAt(vec.Clone(),ifit);
492       fitter->GetErrors(vec);
493       fErrorsOutZ.AddAt(vec.Clone(),ifit);
494     }
495   }
496 }
497
498 Double_t AliTPCcalibUnlinearity::GetDr(Int_t sector, Double_t dout, Double_t dr){
499   //
500   //
501   //
502   TVectorD * vec = GetParamOutR(sector);
503   if (!vec) return 0;
504   dr-=0.5;        // dr=0 at the middle drift length
505   Double_t eout10 = TMath::Exp(-dout/10.);
506   Double_t eout5 = TMath::Exp(-dout/5.);                    
507   Double_t eoutp  = (eout10+eout5)*0.5;
508   Double_t eoutm  = (eout10-eout5)*0.5;
509
510   Double_t result=0;
511   result+=(*vec)[1]*eoutp;
512   result+=(*vec)[2]*eoutm;
513   result+=(*vec)[3]*eoutp*dr;
514   result+=(*vec)[4]*eoutm*dr;
515   result+=(*vec)[5]*eoutp*dr*dr;
516   result+=(*vec)[6]*eoutm*dr*dr;
517   return result;
518 }
519
520
521 Double_t AliTPCcalibUnlinearity::GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
522   //
523   //
524   //
525   Double_t alpha    = TMath::ATan2(gy,gx);
526   if (alpha<0)  alpha+=TMath::Pi()*2;
527   Int_t lsector = Int_t(18*alpha/(TMath::Pi()*2));
528   alpha = (lsector+0.5)*TMath::Pi()/9.;
529   //
530   Double_t r[3];
531   r[0]=gx; r[1]=gy;
532   Double_t cs=TMath::Cos(-alpha), sn=TMath::Sin(-alpha), x=r[0];
533   r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
534   //
535   AliTPCROC * roc = AliTPCROC::Instance();
536   Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
537   Double_t dout = xout-r[0];
538   if (dout<0) return 0;
539   Double_t dr   = 1-TMath::Abs(gz/250.);
540   if (gz<0) lsector+=18;
541   if (stype==0) lsector = (gz>0) ? 36:37;
542   if (stype<0) return lsector;  // 
543   return GetDr(lsector,dout,dr);
544 }
545
546
547 Double_t AliTPCcalibUnlinearity::GetDz(Int_t sector, Double_t dout, Double_t dr){
548   //
549   //
550   //
551   TVectorD * vec = GetParamOutZ(sector);
552   if (!vec) return 0;
553   dr-=0.5; // 0 at the middle
554   Double_t eout10 = TMath::Exp(-dout/10.);
555   Double_t eout5 = TMath::Exp(-dout/5.);
556   Double_t eoutp  = (eout10+eout5)*0.5;
557   Double_t eoutm  = (eout10-eout5)*0.5;
558
559   
560   Double_t result=0;
561   result+=(*vec)[1]*eoutp;
562   result+=(*vec)[2]*eoutm;
563   result+=(*vec)[3]*eoutp*dr;
564   result+=(*vec)[4]*eoutm*dr;
565   result+=(*vec)[5]*eoutp*dr*dr;
566   result+=(*vec)[6]*eoutm*dr*dr;
567   return result;
568 }
569
570
571 Double_t      AliTPCcalibUnlinearity::SGetDr(Int_t sector, Double_t dout, Double_t dr){
572   //
573   //
574   //
575   return fgInstance->GetDr(sector,dout,dr);
576 }
577 Double_t      AliTPCcalibUnlinearity::SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
578   //
579   //
580   //
581   return fgInstance->GetGDr(stype,gx,gy,gz);
582 }
583 Double_t      AliTPCcalibUnlinearity::SGetDz(Int_t sector, Double_t dout, Double_t dr){
584   //
585   //
586   //
587   return fgInstance->GetDz(sector,dout,dr);
588 }
589
590
591
592
593 void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
594   //
595   //
596   //
597   //  TTree * tree = chainUnlinP;
598   AliTPCcalibUnlinearity *calib = this;
599   //
600   AliTPCclusterMI * cl=0;
601   Double_t ty,tz;
602   TVectorD *vy=0, *vz=0;
603   TVectorD *vy2=0, *vz2=0;
604   Long64_t nentries = tree->GetEntries();
605   if (nmaxPoints>0) nentries = TMath::Min(nentries,nmaxPoints);
606   //
607   {
608   for (Long64_t i=0; i<nentries; i++){
609     tree->SetBranchAddress("Cl.",&cl);
610     tree->SetBranchAddress("y1",&ty);
611     tree->SetBranchAddress("z1",&tz);
612     tree->SetBranchAddress("py1.",&vy);
613     tree->SetBranchAddress("pz1.",&vz);
614     tree->SetBranchAddress("py2.",&vy2);
615     tree->SetBranchAddress("pz2.",&vz2);
616     //
617     tree->GetEntry(i);
618     if (!cl) continue;
619     if (i%10000==0) printf("%d\n",(Int_t)i);
620     Int_t row= cl->GetRow();
621     if (cl->GetDetector()>36) row+=64;
622     if (cl->GetType()<0) continue; 
623     Double_t dy = cl->GetY()-ty;
624     Double_t dz = cl->GetZ()-tz;
625     Double_t dr = 1-TMath::Abs(cl->GetZ()/250);
626     //
627     //
628     if (TMath::Abs(dy)>0.25) continue;
629     if (TMath::Abs(dz)>0.25) continue;
630     
631     if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
632     if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
633       // Apply sagita cut
634     calib->AddPoint(cl->GetDetector(), row, dz, dy,
635                     (*vy)[1],(*vz)[1], dr, 1);
636   }
637   }
638 }
639
640
641 void  AliTPCcalibUnlinearity::MakeQPosNormAll(TTree * chainUnlinD, AliTPCClusterParam * param, Int_t maxPoints){
642   //
643   // Make position corrections
644   // for the moment Only using debug streamer 
645   // chainUnlinD  - debug tree
646   // param     - parameters to be updated
647   // maxPoints - maximal number of points using for fit
648   // verbose   - print info flag
649   //
650   // Current implementation - only using debug streamers
651   // 
652   
653   /*    
654   //Defaults
655   Int_t maxPoints=100000;
656   */
657   /*
658     Usage: 
659     //0. Load libraries
660     gSystem->Load("libANALYSIS");
661     gSystem->Load("libSTAT");
662     gSystem->Load("libTPCcalib");
663     
664
665     //1. Load Parameters to be modified:
666     //e.g:
667     
668     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
669     AliCDBManager::Instance()->SetRun(0) 
670     AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
671     //
672     //2. Load the debug streamer
673     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
674     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
675     AliXRDPROOFtoolkit tool;  
676     TChain * chainUnlin = tool.MakeChain("unlin.txt","Resol",0,10200);
677     chainUnlin->Lookup();
678     TChain * chainUnlinD = tool.MakeChain("unlin.txt","Diff",0,10200);
679     chainUnlinD->Lookup();
680
681     //3. Do fits and store results
682     // 
683     AliTPCcalibUnlinearity::MakeQPosNormAll(chainUnlinD,param,200000,0) ;
684     TFile f("paramout.root","recreate");
685     param->Write("clusterParam");
686     f.Close();
687     //4. Verification
688     TFile f2("paramout.root");
689     AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
690     param2->SetInstance(param2);
691     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 
692     
693    */
694
695
696   TStatToolkit toolkit;
697   Double_t chi2;
698   TVectorD fitParamY0;
699   TVectorD fitParamY1;
700   TVectorD fitParamZ0;
701   TVectorD fitParamZ1;
702   TMatrixD covMatrix;
703   Int_t npoints;
704   
705   chainUnlinD->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
706   chainUnlinD->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
707   chainUnlinD->SetAlias("sp","(sin(dp*pi)-dp*pi)");
708   chainUnlinD->SetAlias("st","(sin(dt)-dt)");
709   chainUnlinD->SetAlias("dq","(-1+(Cl.fZ>0)*2)/Cl.fMax");
710   //
711   chainUnlinD->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
712   //
713   //
714   //
715   TCut cutE("Cl.fType>=0");
716   TCut cutDy("abs(Cl.fY-y1)<0.4");
717   TCut cutDz("abs(Cl.fZ-z1)<0.4");
718   TCut cutY2("abs(py2.fElements[2])*150^2/4<1");
719   TCut cutZ2("abs(pz2.fElements[2])*150^2/4<1");
720   TCut cutAll = cutE+cutDy+cutDz+cutY2+cutZ2;
721   //
722   TCut cutA("Cl.fZ>0");
723   TCut cutC("Cl.fZ<0");
724
725   TString fstringY="";  
726   //
727   fstringY+="(dp)++";            //1
728   fstringY+="(dp)*di++";         //2
729   fstringY+="(sp)++";            //3
730   fstringY+="(sp)*di++";         //4
731   fstringY+="(dq)++";            //5
732   TString fstringZ="";  
733   fstringZ+="(dt)++";            //1
734   fstringZ+="(dt)*di++";         //2
735   fstringZ+="(st)++";            //3
736   fstringZ+="(st)*di++";         //4
737   fstringZ+="(dq)++";            //5
738   //
739   // Z corrections
740   //
741   TString *strZ0 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
742   printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
743   strZ0->Tokenize("++")->Print();
744   param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
745   //
746   TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
747   //
748   printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
749   strZ1->Tokenize("++")->Print();
750   param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
751   param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
752   //
753   // Y corrections
754   //   
755   TString *strY0 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
756   printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); 
757   strY0->Tokenize("++")->Print();
758   param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
759   
760
761   TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
762   //
763   printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));  
764   strY1->Tokenize("++")->Print();
765   param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
766   param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
767   //
768   //
769   //
770   chainUnlinD->SetAlias("fitZ0",strZ0->Data());
771   chainUnlinD->SetAlias("fitZ1",strZ1->Data());
772   chainUnlinD->SetAlias("fitY0",strY0->Data());
773   chainUnlinD->SetAlias("fitY1",strY1->Data());
774 }
775
776
777
778
779
780
781