e4090b1a19655c71bc3b79cb496a202f407c9959
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibAlign.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 //     Class to make a internal alignemnt of TPC chambers                    //
20 //
21 //     Requierements - Warnings:
22 //     1. Before using this componenent the magnetic filed has to be set properly //
23 //     2. The systematic effects  - unlinearities has to be understood
24 //
25 //     If systematic and unlinearities are not under control
26 //     the alignment is just effective alignment. Not second order corrction
27 //     are calculated.
28 //    
29 //     The histograming of the edge effects and unlineratities integral part
30 //     of the component (currently only in debug stream)
31 //
32 //     3 general type of linear transformation investigated (see bellow)
33 //
34 //     By default only 6 parameter alignment to be used - other just for QA purposes
35
36 //     Different linear tranformation investigated
37 //     12 parameters - arbitrary linear transformation 
38 //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
39 //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
40 //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
41 //
42 //      9 parameters - scaling fixed to 1
43 //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
44 //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
45 //                     a20  a21  a22 a23     p[4]   p[5]  1      p[8] 
46 //
47 //      6 parameters - x-y rotation x-z, y-z tiliting
48 //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
49 //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
50 //                     a20  a21  a22 a23     p[1]   p[2]  1     p[5] 
51 //
52 //
53 //      Debug stream supported
54 //      0. Align    - The main output of the Alignment component
55 //                  - Used for visualization of the misalignment between sectors
56 //                  - Results of the missalignment fit and the mean and sigmas of histograms
57 //                   stored there
58 //      1. Tracklet - StreamLevel >1
59 //                  - Dump all information about tracklet match from sector1 to sector 2
60 //                  - Default histogram residulas created in parallel
61 //                  - Check this streamer in case of suspicious content of these histograms
62 //      2. Track    - StreamLevel>5  
63 //                  - For debugging of the edge effects
64 //                  - All information  - extrapolation inside of one sectors
65 //                  - Created in order to distinguish between unlinearities inside of o
66 //                    sector and  missalignment 
67    
68 //
69 //
70 /*
71   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
72   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
73   AliXRDPROOFtoolkit tool;
74   TChain * chain = tool.MakeChain("align.txt","Track",0,10200);
75   chain->Lookup();
76   TCut cutA("abs(tp1.fP[1]-tp2.fP[1])<0.3&&abs(tp1.fP[0]-tp2.fP[0])<0.15&&abs(tp1.fP[3]-tp2.fP[3])<0.01&&abs(tp1.fP[2]-tp2.fP[2])<0.01");
77   TCut cutS("s1%36==s2%36");
78   
79   .x ~/UliStyle.C
80   .x $ALICE_ROOT/macros/loadlibsREC.C
81
82   gSystem->Load("$ROOTSYS/lib/libXrdClient.so");
83   gSystem->Load("libProof");
84   gSystem->Load("libANALYSIS");
85   gSystem->Load("libSTAT");
86   gSystem->Load("libTPCcalib");
87   //
88   // compare reference
89   TFile fcalib("CalibObjects.root");
90   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
91
92   AliTPCcalibAlign * align = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
93   //
94   //
95   align->EvalFitters();
96   align->MakeTree("alignTree.root");
97   TFile falignTree("alignTree.root");
98   TTree * treeAlign = (TTree*)falignTree.Get("Align");
99    
100
101 */
102
103 ////
104 //// 
105
106 #include "TLinearFitter.h"
107 #include "AliTPCcalibAlign.h"
108 #include "AliExternalTrackParam.h"
109 #include "AliTPCTracklet.h"
110 #include "TH1D.h"
111 #include "TH2F.h"
112 #include "TVectorD.h"
113 #include "TTreeStream.h"
114 #include "TFile.h"
115 #include "TTree.h"
116 #include "TF1.h"
117 #include "TGraphErrors.h"
118 #include "AliTPCclusterMI.h"
119 #include "AliTPCseed.h"
120 #include "AliTracker.h"
121 #include "TClonesArray.h"
122
123
124 #include "TTreeStream.h"
125 #include <iostream>
126 #include <sstream>
127 using namespace std;
128
129 AliTPCcalibAlign* AliTPCcalibAlign::fgInstance = 0;
130 ClassImp(AliTPCcalibAlign)
131
132
133
134
135 AliTPCcalibAlign* AliTPCcalibAlign::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 AliTPCcalibAlign();
143   }
144   return fgInstance;
145 }
146
147
148
149
150 AliTPCcalibAlign::AliTPCcalibAlign()
151   :  AliTPCcalibBase(),
152      fDphiHistArray(72*72),
153      fDthetaHistArray(72*72),
154      fDyHistArray(72*72),
155      fDzHistArray(72*72),
156      //
157      fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
158      fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
159      fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
160      fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
161      fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
162      fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     
163      fFitterArray12(72*72),
164      fFitterArray9(72*72),
165      fFitterArray6(72*72),
166      fMatrixArray12(72*72),
167      fMatrixArray9(72*72),
168      fMatrixArray6(72*72)
169 {
170   //
171   // Constructor
172   //
173   for (Int_t i=0;i<72*72;++i) {
174     fPoints[i]=0;
175   }
176 }
177
178 AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
179   :AliTPCcalibBase(),  
180    fDphiHistArray(72*72),
181    fDthetaHistArray(72*72),
182    fDyHistArray(72*72),
183    fDzHistArray(72*72),
184    fDyPhiHistArray(72*72),      // array of residual histograms  y     -kYPhi
185    fDzThetaHistArray(72*72),    // array of residual histograms  z-z   -kZTheta
186    fDphiZHistArray(72*72),      // array of residual histograms  phi   -kPhiz
187    fDthetaZHistArray(72*72),    // array of residual histograms  theta -kThetaz
188    fDyZHistArray(72*72),        // array of residual histograms  y     -kYz
189    fDzZHistArray(72*72),        // array of residual histograms  z     -kZz     //
190    fFitterArray12(72*72),
191    fFitterArray9(72*72),
192    fFitterArray6(72*72),
193    fMatrixArray12(72*72),
194    fMatrixArray9(72*72),
195    fMatrixArray6(72*72)
196
197 {
198   //
199   // Constructor
200   //
201   SetName(name);
202   SetTitle(title);
203   for (Int_t i=0;i<72*72;++i) {
204     fPoints[i]=0;
205   }
206 }
207
208
209 AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
210   :AliTPCcalibBase(align),  
211    fDphiHistArray(align.fDphiHistArray),
212    fDthetaHistArray(align.fDthetaHistArray),
213    fDyHistArray(align.fDyHistArray),
214    fDzHistArray(align.fDzHistArray),
215    fDyPhiHistArray(align.fDyPhiHistArray),      // array of residual histograms  y     -kYPhi
216    fDzThetaHistArray(align.fDzThetaHistArray),    // array of residual histograms  z-z   -kZTheta
217    fDphiZHistArray(align.fDphiZHistArray),      // array of residual histograms  phi   -kPhiz
218    fDthetaZHistArray(align.fDthetaZHistArray),    // array of residual histograms  theta -kThetaz
219    fDyZHistArray(align.fDyZHistArray),        // array of residual histograms  y     -kYz
220    fDzZHistArray(align.fDzZHistArray),        // array of residual histograms  z     -kZz     
221    //
222    fFitterArray12(align.fFitterArray12),
223    fFitterArray9(align.fFitterArray9),
224    fFitterArray6(align.fFitterArray6),
225    
226    fMatrixArray12(align.fMatrixArray12),
227    fMatrixArray9(align.fMatrixArray9),
228    fMatrixArray6(align.fMatrixArray6)
229    
230 {
231   //
232   // copy constructor - copy also the content
233   //
234   TH1 * his = 0;
235   TObjArray * arr0=0;
236   const TObjArray *arr1=0;
237   for (Int_t index =0; index<72*72; index++){
238     for (Int_t iarray=0;iarray<10; iarray++){
239       if (iarray==kY){
240         arr0 = &fDyHistArray;
241         arr1 = &align.fDyHistArray;
242       }
243       if (iarray==kZ){
244         arr0 = &fDzHistArray;
245         arr1 = &align.fDzHistArray;
246       }
247       if (iarray==kPhi){
248         arr0 = &fDphiHistArray;
249         arr1 = &align.fDphiHistArray;
250       }
251       if (iarray==kTheta){
252         arr0 = &fDthetaHistArray;
253         arr1 = &align.fDthetaHistArray;
254       }
255       if (iarray==kYz){
256         arr0 = &fDyZHistArray;
257         arr1 = &align.fDyZHistArray;
258       }
259       if (iarray==kZz){
260         arr0 = &fDzZHistArray;
261         arr1 = &align.fDzZHistArray;
262       }
263       if (iarray==kPhiZ){
264         arr0 = &fDphiZHistArray;
265         arr1 = &align.fDphiZHistArray;
266       }
267       if (iarray==kThetaZ){
268         arr0 = &fDthetaZHistArray;
269         arr1 = &align.fDthetaZHistArray;
270       }
271
272       if (iarray==kYPhi){
273         arr0 = &fDyPhiHistArray;
274         arr1 = &align.fDyPhiHistArray;
275       }
276       if (iarray==kZTheta){
277         arr0 = &fDzThetaHistArray;
278         arr1 = &align.fDzThetaHistArray;
279       }
280
281       if (arr1->At(index)) {
282         his = (TH1*)arr1->At(index)->Clone();
283         his->SetDirectory(0);
284         arr0->AddAt(his,index);
285       }    
286     }
287   }
288 }
289
290
291 AliTPCcalibAlign::~AliTPCcalibAlign() {
292   //
293   // destructor
294   //
295 }
296
297 void AliTPCcalibAlign::Process(AliTPCseed *seed) {
298   //
299   // 
300   //
301   
302   TObjArray tracklets=
303     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
304                                     kFALSE,20,2);
305   //  TObjArray trackletsL=
306 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
307 //                                  kFALSE,20,2);
308 //   TObjArray trackletsQ=
309 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
310 //                                  kFALSE,20,2);
311 //   TObjArray trackletsR=
312 //     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
313 //                                  kFALSE,20,2);
314   tracklets.SetOwner();
315  //  trackletsL.SetOwner();
316 //   trackletsQ.SetOwner();
317 //   trackletsR.SetOwner();
318   if (tracklets.GetEntries()==2) {
319     AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
320     AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
321     if (t1->GetSector()>t2->GetSector()) {
322       AliTPCTracklet* tmp=t1;
323       t1=t2;
324       t2=tmp;
325     }
326     AliExternalTrackParam *common1=0,*common2=0;
327     if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
328       ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
329     delete common1;
330     delete common2;
331   }
332   
333 }
334
335 void AliTPCcalibAlign::Analyze(){
336   //
337   // Analyze function 
338   //
339   EvalFitters();
340 }
341
342
343 void AliTPCcalibAlign::Terminate(){
344   //
345   // Terminate function
346   // call base terminate + Eval of fitters
347   //
348   Info("AliTPCcalibAlign","Terminate");
349   EvalFitters();
350   AliTPCcalibBase::Terminate();
351 }
352
353
354
355
356 void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
357                                         const AliExternalTrackParam &tp2,
358                                         const AliTPCseed * seed,
359                                         Int_t s1,Int_t s2) {
360
361   //
362   //
363   //
364   //
365   //
366   // Process function to fill fitters
367   //
368   Double_t t1[5],t2[5];
369   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
370   Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
371   x1   =tp1.GetX();
372   y1   =tp1.GetY();
373   z1   =tp1.GetZ();
374   Double_t snp1=tp1.GetSnp();
375   dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
376   Double_t tgl1=tp1.GetTgl();
377   // dz/dx = 1/(cos(theta)*cos(phi))
378   dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1);
379   x2   =tp2.GetX();
380   y2   =tp2.GetY();
381   z2   =tp2.GetZ();
382   Double_t snp2=tp2.GetSnp();
383   dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
384   Double_t tgl2=tp2.GetTgl();
385   dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2);
386   Bool_t accept =   AcceptTracklet(tp1,tp2);  
387   //
388   //
389   //
390   if (fStreamLevel>1){
391     TTreeSRedirector *cstream = GetDebugStreamer();
392     if (cstream){
393       static TVectorD vec1(5);
394       static TVectorD vec2(5);
395       vec1.SetElements(t1);
396       vec2.SetElements(t2);
397       AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
398       AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
399       (*cstream)<<"Tracklet"<<
400         "run="<<fRun<<              //  run number
401         "event="<<fEvent<<          //  event number
402         "time="<<fTime<<            //  time stamp of event
403         "trigger="<<fTrigger<<      //  trigger
404         "mag="<<fMagF<<             //  magnetic field
405         "isOK="<<accept<<           //  flag - used for alignment
406         "tp1.="<<p1<<
407         "tp2.="<<p2<<
408         "v1.="<<&vec1<<
409         "v2.="<<&vec2<<
410         "s1="<<s1<<
411         "s2="<<s2<<
412         "\n";
413     }
414   }
415   if (!accept) return;
416
417   if (GetDebugLevel()>50) printf("Process track\n");
418   if (GetDebugLevel()>50) printf("Filling track\n");
419   //
420   // fill resolution histograms - previous cut included
421   ProcessDiff(tp1,tp2, seed,s1,s2);
422   FillHisto(tp1,tp2,s1,s2);  
423   ProcessAlign(t1,t2,s1,s2);
424 }
425
426 void AliTPCcalibAlign::ProcessAlign(Double_t * t1,
427                                     Double_t * t2,
428                                     Int_t s1,Int_t s2){
429   //
430   // Do intersector alignment
431   //
432   Process12(t1,t2,GetOrMakeFitter12(s1,s2));
433   Process9(t1,t2,GetOrMakeFitter9(s1,s2));
434   Process6(t1,t2,GetOrMakeFitter6(s1,s2));
435   ++fPoints[GetIndex(s1,s2)];
436 }
437
438 void AliTPCcalibAlign::ProcessTree(TTree * chainTracklet){
439   //
440   // Process the debug streamer tree
441   // Possible to modify selection criteria
442   // Used with entry list
443   //
444   TTreeSRedirector * cstream = new TTreeSRedirector("aligndump.root");
445
446   AliTPCcalibAlign *align = this;
447   //
448   TVectorD * vec1 = 0;
449   TVectorD * vec2 = 0;
450   AliExternalTrackParam * tp1 = new  AliExternalTrackParam;
451   AliExternalTrackParam * tp2 = new  AliExternalTrackParam;  
452   Int_t      s1 = 0;
453   Int_t      s2 = 0;
454   {
455     Int_t entries=chainTracklet->GetEntries();
456     for (Int_t i=0; i< entries; i++){
457       chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
458       chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
459       chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
460       chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
461       chainTracklet->GetBranch("s1")->SetAddress(&s1);
462       chainTracklet->GetBranch("s2")->SetAddress(&s2);      
463       chainTracklet->GetEntry(i);
464       if (s1==s2) continue;
465       if (i%100==0) printf("%d\t%d\t%d\t\n",i, s1,s2);
466       Bool_t accept =   AcceptTracklet(tp1,tp2);  
467       if (!accept) continue;
468       if (vec1->GetMatrixArray()){
469         align->FillHisto(*tp1,*tp2,s1,s2);
470         align->ProcessAlign(vec1->GetMatrixArray(),vec2->GetMatrixArray(),s1,s2);
471         (*cstream)<<"Tracklet"<<
472           "tp1.="<<tp1<<
473           "tp2.="<<tp2<<
474           "v1.="<<vec1<<
475           "v2.="<<vec2<<
476           "s1="<<s1<<
477           "s2="<<s2<<
478           "\n";
479       }
480     }
481   }
482   delete cstream;
483 }
484
485
486 Bool_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
487                                         const AliExternalTrackParam &p2){
488
489   //
490   // Accept pair of tracklets?
491   //
492   /*
493   // resolution cuts
494   TCut cutS0("sqrt(tp2.fC[0]+tp2.fC[0])<0.2");
495   TCut cutS1("sqrt(tp2.fC[2]+tp2.fC[2])<0.2");
496   TCut cutS2("sqrt(tp2.fC[5]+tp2.fC[5])<0.01");
497   TCut cutS3("sqrt(tp2.fC[9]+tp2.fC[9])<0.01");
498   TCut cutS4("sqrt(tp2.fC[14]+tp2.fC[14])<0.5");
499   TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
500   //
501   // parameters matching cuts
502   TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
503   TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
504   TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.03");
505   TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.03");
506   TCut cutP=cutP0+cutP1+cutP2+cutP3;    
507   */  
508   //
509   // resolution cuts
510   const Double_t *cp1 = p1.GetCovariance();
511   const Double_t *cp2 = p2.GetCovariance();
512   if (TMath::Abs(cp1[0]+cp2[0])>0.2) return kFALSE;
513   if (TMath::Abs(cp1[2]+cp2[2])>0.2) return kFALSE;
514   if (TMath::Abs(cp1[5]+cp2[5])>0.01) return kFALSE;
515   if (TMath::Abs(cp1[9]+cp2[9])>0.01) return kFALSE;
516   if (TMath::Abs(cp1[14]+cp2[14])>0.5) return kFALSE;
517
518   //parameters difference
519   const Double_t *tp1 = p1.GetParameter();
520   const Double_t *tp2 = p2.GetParameter();
521   if (TMath::Abs(tp1[0]-tp2[0])>0.6) return kFALSE;
522   if (TMath::Abs(tp1[1]-tp2[1])>0.6) return kFALSE;
523   if (TMath::Abs(tp1[2]-tp2[2])>0.03) return kFALSE;
524   if (TMath::Abs(tp1[3]-tp2[3])>0.03) return kFALSE;
525   //
526   if (TMath::Abs(tp2[1])>235) return kFALSE;
527   return kTRUE;
528 }
529
530
531
532 void  AliTPCcalibAlign::ProcessDiff(const AliExternalTrackParam &t1,
533                                     const AliExternalTrackParam &t2,
534                                     const AliTPCseed *seed,
535                                     Int_t s1,Int_t s2)
536 {
537   //
538   // Process local residuals function
539   // 
540   TVectorD vecX(160);
541   TVectorD vecY(160);
542   TVectorD vecZ(160);
543   TVectorD vecClY(160);
544   TVectorD vecClZ(160);
545   TClonesArray arrCl("AliTPCclusterMI",160);
546   arrCl.ExpandCreateFast(160);
547   Int_t count1=0, count2=0;
548   for (Int_t i=0;i<160;++i) {
549     AliTPCclusterMI *c=seed->GetClusterPointer(i);
550     vecX[i]=0;
551     vecY[i]=0;
552     vecZ[i]=0;
553     if (!c) continue;
554     AliTPCclusterMI & cl = (AliTPCclusterMI&) (*arrCl[i]);
555     if (c->GetDetector()!=s1 && c->GetDetector()!=s2) continue;
556     vecClY[i] = c->GetY();
557     vecClZ[i] = c->GetZ();
558     cl=*c;
559     const AliExternalTrackParam *par = (c->GetDetector()==s1)? &t1:&t2;
560     if (c->GetDetector()==s1) ++count1;
561     if (c->GetDetector()==s2) ++count2;
562     Double_t gxyz[3],xyz[3];
563     t1.GetXYZ(gxyz);
564     Float_t bz = AliTracker::GetBz(gxyz);
565     par->GetYAt(c->GetX(), bz, xyz[1]);
566     par->GetZAt(c->GetX(), bz, xyz[2]);
567     vecX[i] = c->GetX();
568     vecY[i]= xyz[1];
569     vecZ[i]= xyz[2];
570   }
571   //
572   //
573   if (fStreamLevel>5){
574     //
575     // huge output - cluster residuals to be investigated
576     //
577     TTreeSRedirector *cstream = GetDebugStreamer();
578     //AliTPCseed * t = (AliTPCseed*) seed;
579     //AliExternalTrackParam *p0 = &((AliExternalTrackParam&)seed);
580     AliExternalTrackParam *p1 = &((AliExternalTrackParam&)t1);
581     AliExternalTrackParam *p2 = &((AliExternalTrackParam&)t2);
582     /*
583       
584       Track->Draw("Cl[].fY-vtY.fElements:vtY.fElements-vtX.fElements*tan(pi/18.)>>his(100,-10,0)","Cl.fY!=0&&abs(Cl.fY-vtY.fElements)<1","prof");
585
586     */
587
588     if (cstream){
589       (*cstream)<<"Track"<<
590         "run="<<fRun<<              //  run number
591         "event="<<fEvent<<          //  event number
592         "time="<<fTime<<            //  time stamp of event
593         "trigger="<<fTrigger<<      //  trigger
594         "mag="<<fMagF<<             //  magnetic field
595         "Cl.="<<&arrCl<<
596         //"tp0.="<<p0<<
597         "tp1.="<<p1<<
598         "tp2.="<<p2<<
599         "vtX.="<<&vecX<<
600         "vtY.="<<&vecY<<
601         "vtZ.="<<&vecZ<<
602         "vcY.="<<&vecClY<<
603         "vcZ.="<<&vecClZ<<
604         "s1="<<s1<<
605         "s2="<<s2<<
606         "c1="<<count1<<
607         "c2="<<count2<<
608         "\n";
609     }
610   }
611 }
612
613
614
615
616 void AliTPCcalibAlign::Process12(const Double_t *t1,
617                                  const Double_t *t2,
618                                  TLinearFitter *fitter) {
619   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
620   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
621   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
622   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
623   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
624   //
625   //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
626   //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
627   //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
628
629
630
631   const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
632   const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
633
634   // TODO:
635   Double_t sy    = 0.1;
636   Double_t sz    = 0.1;
637   Double_t sdydx = 0.001;
638   Double_t sdzdx = 0.001;
639
640   Double_t p[12];
641   Double_t value;
642
643   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
644   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
645   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
646   for (Int_t i=0; i<12;i++) p[i]=0.;
647   p[3+0] = x1;          // a10
648   p[3+1] = y1;          // a11
649   p[3+2] = z1;          // a12
650   p[9+1] = 1.;          // a13
651   p[0+1] = y1*dydx2;    // a01
652   p[0+2] = z1*dydx2;    // a02
653   p[9+0] = dydx2;       // a03
654   value  = y2;
655   fitter->AddPoint(p,value,sy);
656
657   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
658   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
659   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
660   for (Int_t i=0; i<12;i++) p[i]=0.;
661   p[6+0] = x1;           // a20 
662   p[6+1] = y1;           // a21
663   p[6+2] = z1;           // a22
664   p[9+2] = 1.;           // a23
665   p[0+1] = y1*dzdx2;     // a01
666   p[0+2] = z1*dzdx2;     // a02
667   p[9+0] = dzdx2;        // a03
668   value  = z2;
669   fitter->AddPoint(p,value,sz);
670
671   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
672   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
673   for (Int_t i=0; i<12;i++) p[i]=0.;
674   p[3+0] = 1.;           // a10
675   p[3+1] = dydx1;        // a11
676   p[3+2] = dzdx1;        // a12
677   p[0+0] = -dydx2;       // a00
678   p[0+1] = -dydx1*dydx2; // a01
679   p[0+2] = -dzdx1*dydx2; // a02
680   value  = 0.;
681   fitter->AddPoint(p,value,sdydx);
682
683   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
684   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
685   for (Int_t i=0; i<12;i++) p[i]=0.;
686   p[6+0] = 1;            // a20
687   p[6+1] = dydx1;        // a21
688   p[6+2] = dzdx1;        // a22
689   p[0+0] = -dzdx2;       // a00
690   p[0+1] = -dydx1*dzdx2; // a01
691   p[0+2] = -dzdx1*dzdx2; // a02
692   value  = 0.;
693   fitter->AddPoint(p,value,sdzdx);
694 }
695
696 void AliTPCcalibAlign::Process9(Double_t *t1,
697                                 Double_t *t2,
698                                 TLinearFitter *fitter) {
699   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
700   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
701   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
702   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
703   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
704   //
705   //                     a00  a01  a02 a03     1      p[0]  p[1]   p[6]
706   //                     a10  a11  a12 a13 ==> p[2]   1     p[3]   p[7]
707   //                     a20  a21  a21 a23     p[4]   p[5]  1      p[8] 
708
709
710   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
711   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
712
713   // TODO:
714   Double_t sy    = 0.1;
715   Double_t sz    = 0.1;
716   Double_t sdydx = 0.001;
717   Double_t sdzdx = 0.001;
718   //
719   Double_t p[12];
720   Double_t value;
721
722   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
723   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
724   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
725   for (Int_t i=0; i<12;i++) p[i]=0.;
726   p[2]   += x1;           // a10
727   //p[]  +=1;             // a11
728   p[3]   += z1;           // a12    
729   p[7]   += 1;            // a13
730   p[0]   += y1*dydx2;     // a01
731   p[1]   += z1*dydx2;     // a02
732   p[6]   += dydx2;        // a03
733   value   = y2-y1;        //-a11
734   fitter->AddPoint(p,value,sy);
735   //
736   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
737   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
738   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
739   for (Int_t i=0; i<12;i++) p[i]=0.;
740   p[4]   += x1;           // a20 
741   p[5]   += y1;           // a21
742   //p[]  += z1;           // a22
743   p[8]   += 1.;           // a23
744   p[0]   += y1*dzdx2;     // a01
745   p[1]   += z1*dzdx2;     // a02
746   p[6]   += dzdx2;        // a03
747   value  = z2-z1;         //-a22
748   fitter->AddPoint(p,value,sz);
749
750   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
751   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
752   for (Int_t i=0; i<12;i++) p[i]=0.;
753   p[2]   += 1.;           // a10
754   //p[]  += dydx1;      // a11
755   p[3]   += dzdx1;        // a12
756   //p[]  += -dydx2;       // a00
757   p[0]   += -dydx1*dydx2; // a01
758   p[1]   += -dzdx1*dydx2; // a02
759   value  = -dydx1+dydx2;  // -a11 + a00
760   fitter->AddPoint(p,value,sdydx);
761
762   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
763   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
764   for (Int_t i=0; i<12;i++) p[i]=0.;
765   p[4]   += 1;            // a20
766   p[5]   += dydx1;        // a21
767   //p[]  += dzdx1;        // a22
768   //p[]  += -dzdx2;       // a00
769   p[0]   += -dydx1*dzdx2; // a01
770   p[1]   += -dzdx1*dzdx2; // a02
771   value  = -dzdx1+dzdx2;  // -a22 + a00
772   fitter->AddPoint(p,value,sdzdx);
773 }
774
775 void AliTPCcalibAlign::Process6(Double_t *t1,
776                                 Double_t *t2,
777                                 TLinearFitter *fitter) {
778   // x2    =  1  *x1 +-a01*y1 + 0      +a03
779   // y2    =  a01*x1 + 1  *y1 + 0      +a13
780   // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
781   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
782   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
783   //
784   //                     a00  a01  a02 a03     1     -p[0]  0     p[3]
785   //                     a10  a11  a12 a13 ==> p[0]   1     0     p[4]
786   //                     a20  a21  a21 a23     p[1]   p[2]  1     p[5] 
787
788   Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
789   Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
790
791   // TODO:
792   Double_t sy    = 0.1;
793   Double_t sz    = 0.1;
794   Double_t sdydx = 0.001;
795   Double_t sdzdx = 0.001;
796
797   Double_t p[12];
798   Double_t value;
799   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
800   // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
801   // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
802   for (Int_t i=0; i<12;i++) p[i]=0.;
803   p[0]   += x1;           // a10
804   //p[]  +=1;             // a11
805   //p[]  += z1;           // a12    
806   p[4]   += 1;            // a13
807   p[0]   += -y1*dydx2;    // a01
808   //p[]  += z1*dydx2;     // a02
809   p[3]   += dydx2;        // a03
810   value   = y2-y1;        //-a11
811   fitter->AddPoint(p,value,sy);
812   //
813   // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
814   // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
815   // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
816   for (Int_t i=0; i<12;i++) p[i]=0.;
817   p[1]   += x1;           // a20 
818   p[2]   += y1;           // a21
819   //p[]  += z1;           // a22
820   p[5]   += 1.;           // a23
821   p[0]   += -y1*dzdx2;    // a01
822   //p[]   += z1*dzdx2;     // a02
823   p[3]   += dzdx2;        // a03
824   value  = z2-z1;         //-a22
825   fitter->AddPoint(p,value,sz);
826
827   // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
828   // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
829   for (Int_t i=0; i<12;i++) p[i]=0.;
830   p[0]   += 1.;           // a10
831   //p[]  += dydx1;      // a11
832   //p[]   += dzdx1;        // a12
833   //p[]  += -dydx2;       // a00
834   p[0]   +=  dydx1*dydx2; // a01
835   //p[]   += -dzdx1*dydx2; // a02
836   value  = -dydx1+dydx2;  // -a11 + a00
837   fitter->AddPoint(p,value,sdydx);
838
839   // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
840   // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
841   for (Int_t i=0; i<12;i++) p[i]=0.;
842   p[1]   += 1;            // a20
843   p[2]   += dydx1;        // a21
844   //p[]  += dzdx1;        // a22
845   //p[]  += -dzdx2;       // a00
846   p[0]   +=  dydx1*dzdx2; // a01
847   //p[]  += -dzdx1*dzdx2; // a02
848   value  = -dzdx1+dzdx2;  // -a22 + a00
849   fitter->AddPoint(p,value,sdzdx);
850 }
851
852
853
854
855 void AliTPCcalibAlign::EvalFitters() {
856   //
857   // Analyze function 
858   // 
859   // Perform the fitting using linear fitters
860   //
861   Int_t kMinPoints =50;
862   TLinearFitter *f;
863   TFile fff("alignDebug.root","recreate");
864   for (Int_t s1=0;s1<72;++s1)
865     for (Int_t s2=0;s2<72;++s2){
866       if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
867         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
868         if (f->Eval()!=0) {
869           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
870           f->Write(Form("f12_%d_%d",s1,s2));
871         }else{
872           f->Write(Form("f12_%d_%d",s1,s2));
873         }
874       }
875       if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
876         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
877         if (f->Eval()!=0) {
878           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
879         }else{
880           f->Write(Form("f9_%d_%d",s1,s2));
881         }
882       }
883       if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
884         //      cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
885         if (f->Eval()!=0) {
886           cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
887         }else{
888           f->Write(Form("f6_%d_%d",s1,s2));
889         }
890       }
891     }
892   TMatrixD mat(4,4);
893   for (Int_t s1=0;s1<72;++s1)
894     for (Int_t s2=0;s2<72;++s2){
895       if (GetTransformation12(s1,s2,mat)){
896         fMatrixArray12.AddAt(mat.Clone(), GetIndex(s1,s2));
897       }
898       if (GetTransformation9(s1,s2,mat)){
899         fMatrixArray9.AddAt(mat.Clone(), GetIndex(s1,s2));
900       }
901       if (GetTransformation6(s1,s2,mat)){
902         fMatrixArray6.AddAt(mat.Clone(), GetIndex(s1,s2));
903       }
904     }
905   //this->Write("align");
906   
907 }
908
909 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
910   //
911   // get or make fitter - general linear transformation
912   //
913   static Int_t counter12=0;
914   static TF1 f12("f12","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]++x[9]++x[10]++x[11]");
915   TLinearFitter * fitter = GetFitter12(s1,s2);
916   if (fitter) return fitter;
917   //  fitter =new TLinearFitter(12,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]++x[9]++x[10]++x[11]");
918   fitter =new TLinearFitter(&f12,"");
919   fitter->StoreData(kFALSE);
920   fFitterArray12.AddAt(fitter,GetIndex(s1,s2)); 
921   counter12++;
922   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter12<<endl;
923   return fitter;
924 }
925
926 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
927   //
928   //get or make fitter - general linear transformation - no scaling
929   // 
930   static Int_t counter9=0;
931   static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
932   TLinearFitter * fitter = GetFitter9(s1,s2);
933   if (fitter) return fitter;
934   //  fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
935   fitter =new TLinearFitter(&f9,"");
936   fitter->StoreData(kFALSE);
937   fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
938   counter9++;
939   if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<"  :  "<<counter9<<endl;
940   return fitter;
941 }
942
943 TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
944   //
945   // get or make fitter  - 6 paramater linear tranformation
946   //                     - no scaling
947   //                     - rotation x-y
948   //                     - tilting x-z, y-z
949   static Int_t counter6=0;
950   static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
951   TLinearFitter * fitter = GetFitter6(s1,s2);
952   if (fitter) return fitter;
953   //  fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
954   fitter=new TLinearFitter(&f6,"");
955   fitter->StoreData(kTRUE);
956   fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
957   counter6++;
958   if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<"  :  "<<counter6<<endl;
959   return fitter;
960 }
961
962
963
964
965
966 Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
967   //
968   // GetTransformation matrix - 12 paramaters - generael linear transformation
969   //
970   if (!GetFitter12(s1,s2))
971     return false;
972   else {
973     TVectorD p(12);
974     GetFitter12(s1,s2)->GetParameters(p);
975     a.ResizeTo(4,4);
976     a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
977     a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
978     a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
979     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
980     return true;
981   } 
982 }
983
984 Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
985   //
986   // GetTransformation matrix - 9 paramaters - general linear transformation
987   //                            No scaling
988   //
989   if (!GetFitter9(s1,s2))
990     return false;
991   else {
992     TVectorD p(9);
993     GetFitter9(s1,s2)->GetParameters(p);
994     a.ResizeTo(4,4);
995     a[0][0]=1;    a[0][1]=p[0]; a[0][2]=p[1]; a[0][3]=p[6];
996     a[1][0]=p[2]; a[1][1]=1;    a[1][2]=p[3]; a[1][3]=p[7];
997     a[2][0]=p[4]; a[2][1]=p[5]; a[2][2]=1;    a[2][3]=p[8];
998     a[3][0]=0.;   a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
999     return true;
1000   } 
1001 }
1002
1003 Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
1004   //
1005   // GetTransformation matrix - 6  paramaters
1006   //                            3  translation
1007   //                            1  rotation -x-y  
1008   //                            2  tilting x-z y-z
1009   if (!GetFitter6(s1,s2))
1010     return false;
1011   else {
1012     TVectorD p(6);
1013     GetFitter6(s1,s2)->GetParameters(p);
1014     a.ResizeTo(4,4);
1015     a[0][0]=1;       a[0][1]=-p[0];a[0][2]=0;    a[0][3]=p[3];
1016     a[1][0]=p[0];    a[1][1]=1;    a[1][2]=0;    a[1][3]=p[4];
1017     a[2][0]=p[1];    a[2][1]=p[2]; a[2][2]=1;    a[2][3]=p[5];
1018     a[3][0]=0.;      a[3][1]=0.;   a[3][2]=0.;   a[3][3]=1.;
1019     return true;
1020   } 
1021 }
1022
1023 void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
1024                                         const AliExternalTrackParam &tp2,
1025                                         Int_t s1,Int_t s2) {
1026   //
1027   // Fill residual histograms
1028   // Innner-Outer
1029   // Left right - x-y
1030   // A-C side
1031   if (TMath::Abs(s2%36-s1%36)<2 || TMath::Abs(s2%18-s1%18)==0)  {  
1032     GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));    
1033     GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1034     GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
1035     GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
1036     //
1037     GetHisto(kPhiZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));    
1038     GetHisto(kThetaZ,s1,s2,kTRUE)->Fill(tp1.GetZ(),TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
1039     GetHisto(kYz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetY()-tp2.GetY());
1040     GetHisto(kZz,s1,s2,kTRUE)->Fill(tp1.GetZ(),tp1.GetZ()-tp2.GetZ());
1041     //
1042     GetHisto(kYPhi,s1,s2,kTRUE)->Fill(tp1.GetSnp(),tp1.GetY()-tp2.GetY());
1043     GetHisto(kZTheta,s1,s2,kTRUE)->Fill(tp1.GetTgl(),tp1.GetZ()-tp2.GetZ());
1044
1045
1046   }  
1047 }
1048
1049
1050
1051 TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
1052 {
1053   //
1054   // return specified residual histogram - it is only QA
1055   // if force specified the histogram and given histogram is not existing 
1056   //  new histogram is created
1057   //
1058   if (GetIndex(s1,s2)>=72*72) return 0;
1059   TObjArray *histoArray=0;
1060   switch (type) {
1061   case kY:
1062     histoArray = &fDyHistArray; break;
1063   case kZ:
1064     histoArray = &fDzHistArray; break;
1065   case kPhi:
1066     histoArray = &fDphiHistArray; break;
1067   case kTheta:
1068     histoArray = &fDthetaHistArray; break;
1069   case kYPhi:
1070     histoArray = &fDyPhiHistArray; break;
1071   case kZTheta:
1072     histoArray = &fDzThetaHistArray; break;
1073   case kYz:
1074     histoArray = &fDyZHistArray; break;
1075   case kZz:
1076     histoArray = &fDzZHistArray; break;
1077   case kPhiZ:
1078     histoArray = &fDphiZHistArray; break;
1079   case kThetaZ:
1080     histoArray = &fDthetaZHistArray; break;
1081   }
1082   TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
1083   if (histo) return histo;
1084   if (force==kFALSE) return 0; 
1085   //
1086   stringstream name;
1087   stringstream title;
1088   switch (type) {    
1089   case kY:
1090     name<<"hist_y_"<<s1<<"_"<<s2;
1091     title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
1092     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1093     break;
1094   case kZ:
1095     name<<"hist_z_"<<s1<<"_"<<s2;
1096     title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
1097     histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
1098     break;
1099   case kPhi:
1100     name<<"hist_phi_"<<s1<<"_"<<s2;
1101     title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
1102     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1103     break;
1104   case kTheta:
1105     name<<"hist_theta_"<<s1<<"_"<<s2;
1106     title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
1107     histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
1108     break;
1109     //
1110     //
1111   case kYPhi:
1112     name<<"hist_yphi_"<<s1<<"_"<<s2;
1113     title<<"Y Missalignment for sectors Phi"<<s1<<" and "<<s2;
1114     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-1,1,128,-0.3,0.3); // +/- 3 mm
1115     break;
1116   case kZTheta:
1117     name<<"hist_ztheta_"<<s1<<"_"<<s2;
1118     title<<"Z Missalignment for sectors Theta"<<s1<<" and "<<s2;
1119     histo = new TH2F(name.str().c_str(),title.str().c_str(),128,20,-1,1,-0.3,0.3); // +/- 3 mm
1120     break;
1121     //
1122     //
1123     //
1124   case kYz:
1125     name<<"hist_yz_"<<s1<<"_"<<s2;
1126     title<<"Y Missalignment for sectors Z"<<s1<<" and "<<s2;
1127     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1128     break;
1129   case kZz:
1130     name<<"hist_zz_"<<s1<<"_"<<s2;
1131     title<<"Z Missalignment for sectors Z"<<s1<<" and "<<s2;
1132     histo = new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.3,0.3); // +/- 3 mm
1133     break;
1134   case kPhiZ:
1135     name<<"hist_phiz_"<<s1<<"_"<<s2;
1136     title<<"Phi Missalignment for sectors Z"<<s1<<" and "<<s2;
1137     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1138     break;
1139   case kThetaZ:
1140     name<<"hist_thetaz_"<<s1<<"_"<<s2;
1141     title<<"Theta Missalignment for sectors Z"<<s1<<" and "<<s2;
1142     histo =new TH2F(name.str().c_str(),title.str().c_str(),20,-250,250,128,-0.01,0.01); // +/- 10 mrad
1143     break;
1144
1145
1146   }
1147   histo->SetDirectory(0);
1148   histoArray->AddAt(histo,GetIndex(s1,s2));
1149   return histo;
1150 }
1151
1152 TGraphErrors * AliTPCcalibAlign::MakeGraph(Int_t sec0, Int_t sec1, Int_t dsec, 
1153                                            Int_t i0, Int_t i1, FitType type) 
1154 {
1155   //
1156   //
1157   //
1158   TMatrixD mat;
1159   //TObjArray *fitArray=0;
1160   Double_t xsec[1000];
1161   Double_t ysec[1000];
1162   Int_t npoints=0;
1163   for (Int_t isec = sec0; isec<=sec1; isec++){
1164     Int_t isec2 = (isec+dsec)%72;    
1165     switch (type) {
1166     case k6:
1167       GetTransformation6(isec,isec2,mat);break;
1168     case k9:
1169       GetTransformation9(isec,isec2,mat);break;
1170     case k12:
1171       GetTransformation12(isec,isec2,mat);break;
1172     }
1173     xsec[npoints]=isec;
1174     ysec[npoints]=mat(i0,i1);
1175     ++npoints;
1176   }
1177   TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
1178   Char_t name[1000];
1179   sprintf(name,"Mat[%d,%d]  Type=%d",i0,i1,type);
1180   gr->SetName(name);
1181   return gr;
1182 }
1183
1184 void  AliTPCcalibAlign::MakeTree(const char *fname){
1185   //
1186   // make tree with alignment cosntant  -
1187   // For  QA visualization
1188   //
1189   /*
1190    TFile f("CalibObjects.root");
1191    TObjArray *array  = (TObjArray*)f.Get("TPCCalib");
1192    AliTPCcalibAlign   *alignTPC = (AliTPCcalibAlign   *)array->At(0);
1193    alignTPC->MakeTree("alignTree.root");
1194    TFile falign("alignTree.root");
1195    Align->Draw("dy")
1196    */
1197   const Int_t kMinPoints=50;
1198   TTreeSRedirector cstream(fname);
1199   for (Int_t s1=0;s1<72;++s1)
1200     for (Int_t s2=0;s2<72;++s2){
1201       TMatrixD m6;
1202       TMatrixD m9;
1203       TMatrixD m12;
1204       Double_t dy=0, dz=0, dphi=0,dtheta=0;
1205       Double_t sy=0, sz=0, sphi=0,stheta=0;
1206       Double_t ny=0, nz=0, nphi=0,ntheta=0;
1207       Double_t chi2v12=0, chi2v9=0, chi2v6=0;
1208       Int_t npoints=0;
1209       if (fPoints[GetIndex(s1,s2)]>kMinPoints){
1210         //
1211         //
1212         //
1213         TLinearFitter * fitter = 0;
1214         fitter = GetFitter12(s1,s2);
1215         npoints = fitter->GetNpoints();
1216         chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1217         //
1218         fitter = GetFitter9(s1,s2);
1219         npoints = fitter->GetNpoints();
1220         chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1221         //
1222         fitter = GetFitter6(s1,s2);
1223         npoints = fitter->GetNpoints();
1224         chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
1225         
1226         GetTransformation6(s1,s2,m6);
1227         GetTransformation9(s1,s2,m9);
1228         GetTransformation12(s1,s2,m12);
1229         TH1 * his=0;
1230         his = GetHisto(kY,s1,s2);
1231         if (his) { dy = his->GetMean(); sy = his->GetRMS(); ny = his->GetEntries();}
1232         his = GetHisto(kZ,s1,s2);
1233         if (his) { dz = his->GetMean(); sz = his->GetRMS(); nz = his->GetEntries();}
1234         his = GetHisto(kPhi,s1,s2);
1235         if (his) { dphi = his->GetMean(); sphi = his->GetRMS(); nphi = his->GetEntries();}
1236         his = GetHisto(kTheta,s1,s2);
1237         if (his) { dtheta = his->GetMean(); stheta = his->GetRMS(); ntheta = his->GetEntries();}
1238         //
1239       }
1240
1241       // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1242       // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1243       // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1244       // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1245       // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1246       //
1247       //                     a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
1248       //                     a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
1249       //                     a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
1250       
1251       //
1252       // 
1253       // dy:-(134*m6.fElements[4]+m6.fElements[7])
1254       // 
1255       // dphi:-(m6.fElements[4])
1256       //
1257       // dz:134*m6.fElements[8]+m6.fElements[11]
1258       //
1259       // dtheta:m6.fElements[8]
1260       //
1261       cstream<<"Align"<<
1262         "s1="<<s1<<     // reference sector
1263         "s2="<<s2<<     // sector to align
1264         "m6.="<<&m6<<   // tranformation matrix
1265         "m9.="<<&m9<<   // 
1266         "m12.="<<&m12<<
1267         "chi2v12="<<chi2v12<<
1268         "chi2v9="<<chi2v9<<
1269         "chi2v6="<<chi2v6<<
1270         //               histograms mean RMS and entries
1271         "dy="<<dy<<  
1272         "sy="<<sy<<
1273         "ny="<<ny<<
1274         "dz="<<dz<<
1275         "sz="<<sz<<
1276         "nz="<<nz<<
1277         "dphi="<<dphi<<
1278         "sphi="<<sphi<<
1279         "nphi="<<nphi<<
1280         "dtheta="<<dtheta<<
1281         "stheta="<<stheta<<
1282         "ntheta="<<ntheta<<
1283         "\n";
1284     }
1285
1286 }
1287
1288
1289 //_____________________________________________________________________
1290 Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
1291   //
1292   // merge function 
1293   //
1294   if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Merge");
1295   if (!list)
1296     return 0;  
1297   if (list->IsEmpty())
1298     return 1;
1299   
1300   TIterator* iter = list->MakeIterator();
1301   TObject* obj = 0;  
1302   iter->Reset();
1303   Int_t count=0;
1304   //
1305   TString str1(GetName());
1306   while((obj = iter->Next()) != 0)
1307     {
1308       AliTPCcalibAlign* entry = dynamic_cast<AliTPCcalibAlign*>(obj);
1309       if (entry == 0) continue; 
1310       if (str1.CompareTo(entry->GetName())!=0) continue;
1311       Add(entry);
1312       count++;
1313     } 
1314   return count;
1315 }
1316
1317
1318 void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
1319   //
1320   // Add entry - used for merging of compoents
1321   //
1322   for (Int_t i=0; i<72;i++){
1323     for (Int_t j=0; j<72;j++){
1324       if (align->fPoints[GetIndex(i,j)]<10) continue;
1325       fPoints[GetIndex(i,j)]+=align->fPoints[GetIndex(i,j)];
1326       //
1327       //
1328       //
1329       for (Int_t itype=0; itype<10; itype++){
1330         TH1 * his0=0, *his1=0;
1331         his0 = GetHisto((HistoType)itype,i,j);
1332         his1 = align->GetHisto((HistoType)itype,i,j);
1333         if (his1){
1334           if (his0) his0->Add(his1);
1335           else {
1336             his0 = GetHisto(kY,i,j,kTRUE);
1337             his0->Add(his1);
1338           }
1339         }       
1340       }      
1341     }
1342   }
1343   TLinearFitter *f0=0;
1344   TLinearFitter *f1=0;
1345   for (Int_t i=0; i<72;i++){
1346     for (Int_t j=0; j<72;j++){     
1347       if (align->fPoints[GetIndex(i,j)]<20) continue;
1348       //
1349       //
1350       // fitter12
1351       f0 =  GetFitter12(i,j);
1352       f1 =  align->GetFitter12(i,j);
1353       if (f1 &&f1->Eval()){
1354         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1355         else {
1356           fFitterArray12.AddAt(f1->Clone(),GetIndex(i,j));
1357         }
1358       }      
1359       //
1360       // fitter9
1361       f0 =  GetFitter9(i,j);
1362       f1 =  align->GetFitter9(i,j);
1363       if (f1&&f1->Eval()){
1364         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1365         else { 
1366           fFitterArray9.AddAt(f1->Clone(),GetIndex(i,j));
1367         }
1368       }      
1369       f0 =  GetFitter6(i,j);
1370       f1 =  align->GetFitter6(i,j);
1371       if (f1 &&f1->Eval()){
1372         if (f0&&f0->GetNpoints()>10) f0->Add(f1);
1373         else {
1374           fFitterArray6.AddAt(f1->Clone(),GetIndex(i,j));
1375         }
1376       }   
1377     }
1378   }
1379 }
1380
1381 Double_t AliTPCcalibAlign::Correct(Int_t type, Int_t value, Int_t s1, Int_t s2, Double_t x1, Double_t y1, Double_t z1, Double_t dydx1,Double_t dzdx1){
1382   //
1383   // GetTransformed value
1384   //
1385   //
1386   // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
1387   // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
1388   // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
1389   // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1390   // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1391   
1392   
1393   const TMatrixD * mat = GetTransformation(s1,s2,type);
1394   if (!mat) {
1395     if (value==0) return x1;
1396     if (value==1) return y1;
1397     if (value==2) return z1;
1398     if (value==3) return dydx1;
1399     if (value==4) return dzdx1;
1400     //
1401     if (value==5) return dydx1;
1402     if (value==6) return dzdx1;
1403     return 0;
1404   }
1405   Double_t valT=0;
1406
1407   if (value==0){
1408     valT = (*mat)(0,0)*x1+(*mat)(0,1)*y1+(*mat)(0,2)*z1+(*mat)(0,3);
1409   }
1410
1411   if (value==1){
1412     valT = (*mat)(1,0)*x1+(*mat)(1,1)*y1+(*mat)(1,2)*z1+(*mat)(1,3);
1413   }
1414   if (value==2){
1415     valT = (*mat)(2,0)*x1+(*mat)(2,1)*y1+(*mat)(2,2)*z1+(*mat)(2,3);
1416   }
1417   if (value==3){
1418     //    dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1419     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1 +(*mat)(1,2)*dzdx1;
1420     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1421   }
1422
1423   if (value==4){
1424     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)    
1425     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1 +(*mat)(2,2)*dzdx1;
1426     valT/= ((*mat)(0,0)    +(*mat)(0,1)*dydx1 +(*mat)(0,2)*dzdx1);
1427   }
1428   //
1429   if (value==5){
1430     // onlys shift in angle
1431     //    dydx2 =  (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
1432     valT =  (*mat)(1,0)    +(*mat)(1,1)*dydx1;
1433   }
1434
1435   if (value==6){
1436     // only shift in angle
1437     // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)    
1438     valT =  (*mat)(2,0)    +(*mat)(2,1)*dydx1;
1439   }
1440   //
1441   return valT;
1442 }
1443
1444
1445
1446 /*
1447   
1448
1449 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1450 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1451 AliXRDPROOFtoolkit tool;
1452 TChain * chainTr = tool.MakeChain("align.txt","Track",0,10200);
1453 chainTr->Lookup();
1454 TChain * chainTracklet = tool.MakeChain("align.txt","Tracklet",0,20);
1455 chainTracklet->Lookup();
1456
1457
1458 TCut cutS0("sqrt(tp2.fC[0]+tp2.fC[0])<0.6");
1459 TCut cutS1("sqrt(tp2.fC[2]+tp2.fC[2])<0.6");
1460 TCut cutS2("sqrt(tp2.fC[5]+tp2.fC[5])<0.04");
1461 TCut cutS3("sqrt(tp2.fC[9]+tp2.fC[9])<0.02");
1462 TCut cutS4("sqrt(tp2.fC[14]+tp2.fC[14])<0.5");
1463 // resolution cuts
1464 TCut cutS=cutS0+cutS1+cutS2+cutS3+cutS4;
1465
1466 TCut cutP0("abs(tp1.fP[0]-tp2.fP[0])<0.6");
1467 TCut cutP1("abs(tp1.fP[1]-tp2.fP[1])<0.6");
1468 TCut cutP2("abs(tp1.fP[2]-tp2.fP[2])<0.02");
1469 TCut cutP3("abs(tp1.fP[3]-tp2.fP[3])<0.01");
1470 TCut cutE("abs(tp2.fP[1])<235");
1471 TCut cutP=cutP0+cutP1+cutP2+cutP3+cutE;
1472
1473
1474 //
1475 //
1476 TCut cutA =  cutP+cutS;
1477 chainTr->Draw(">>listEL",cutA,"entryList");
1478 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1479 chainTr->SetEntryList(elist);
1480
1481
1482
1483
1484
1485
1486 TCut cutS("s1%36==s2%36");
1487
1488 TCut cutN("c1>32&&c2>60");
1489 TCut cutC0("sqrt(tp2.fC[0])<1");
1490
1491 TCut cutX("abs(tp2.fX-133.6)<2");
1492
1493 TCut cutA = cutP+cutN;
1494
1495
1496 TCut cutY("abs(vcY.fElements-vtY.fElements)<0.3&&vcY.fElements!=0")
1497 TCut cutZ("abs(vcZ.fElements-vtZ.fElements)<0.3&&vcZ.fElements!=0")
1498
1499
1500
1501
1502
1503
1504 TCut cutA =  cutP+cutS;
1505 chainTracklet->Draw(">>listEL",cutA,"entryList");
1506 TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1507 chainTracklet->SetEntryList(elist);
1508
1509 //
1510 TVectorD * vec1 = 0;
1511 TVectorD * vec2 = 0;
1512 AliExternalTrackParam * tp1 = 0;
1513 AliExternalTrackParam * tp2 = 0;
1514
1515 Int_t      s1 = 0;
1516 Int_t      s2 = 0;
1517 chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1518 chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1519 chainTracklet->GetBranch("s1")->SetAddress(&s1);
1520 chainTracklet->GetBranch("s2")->SetAddress(&s2);
1521
1522
1523 AliTPCcalibAlign align;
1524 {
1525 for (Int_t i=0; i< elist->GetN(); i++){
1526 //for (Int_t i=0; i<100000; i++){
1527 chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
1528 chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
1529 chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
1530 chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
1531 chainTracklet->GetBranch("s1")->SetAddress(&s1);
1532 chainTracklet->GetBranch("s2")->SetAddress(&s2);
1533
1534 chainTracklet->GetEntry(i);
1535 if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i,tentry, s1,s2);
1536 //vec1.Print();
1537 TLinearFitter * fitter = align.GetOrMakeFitter6(s1,s2);
1538 if (fitter) align.Process6(vec1->GetMatrixArray(),vec2->GetMatrixArray(),fitter);
1539 }
1540 }
1541
1542 */