M AliTPCcalibBase.h - SetRun function
[u/mrichter/AliRoot.git] / TPC / AliTPCkalmanAlign.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   TPC Kalman filter implementation for TPC sector alignment.
18   Output of the AliTPCcalibAlign is used as a input for TPC global alignment.
19   In AliTPCcalibAlign histograms - track parameter matching at sector boundaries are created.
20   Each sector is aligned with 5 neighborhoud (sectors)
21     1. Up-down    - 1
22     2. Left-right - 4
23
24   Sector alignment parameters are obtained finding the alignment parameters, minimizing
25   misalignmnet for all piars fo sectors.
26
27   Global minimization-  MakeGlobalAlign
28   
29
30   Example usage:
31   gSystem->Load("libANALYSIS");
32   gSystem->Load("libTPCcalib");
33   //
34   Int_t run=117220;
35   .x ConfigCalibTrain.C(run)
36   calibDB = AliTPCcalibDB::Instance()
37
38   AliTPCkalmanAlign kalmanAlign("TPC align", "TPC align");  // create the object
39   kalmanAlign.ReadAlign("CalibObjects.root");               // read the calibration form file         
40   kalmanAlign.MakeGlobalAlign();                            // do kalman alignment
41   kalmanAlign.DrawDeltaAlign();                             // make QA plot
42   //
43   
44
45 */
46 #include "TMath.h"
47 #include "TTreeStream.h"
48 #include "TGraph.h"
49 #include "TGraphErrors.h"
50 #include "TVectorD.h"
51 #include "TClonesArray.h"
52 #include "TCut.h"
53 #include "TChain.h"
54 #include "AliXRDPROOFtoolkit.h"
55 #include "TLegend.h"
56 #include "TCanvas.h"
57
58 #include "AliTPCcalibDB.h"
59 #include "AliTPCCalROC.h"
60 #include "AliCDBMetaData.h"
61 #include "AliCDBId.h"
62 #include "AliCDBManager.h"
63 #include "AliCDBStorage.h"
64 #include "AliCDBEntry.h"
65 #include "AliAlignObjParams.h"
66 #include "AliTPCROC.h"
67 #include "AliTracker.h"
68 #include "TFile.h"
69 #include "TLinearFitter.h"
70 #include "AliTPCcalibAlign.h"
71 #include "TH1.h"
72 #include "AliTPCCalPad.h"
73 #include "AliTPCkalmanAlign.h"
74 #include "TStatToolkit.h"
75 #include "AliTPCPreprocessorOnline.h"
76 #include "TPostScript.h"
77 #include "AliGRPObject.h"
78
79 AliTPCkalmanAlign::AliTPCkalmanAlign():
80   TNamed(),
81   fCalibAlign(0),     // kalman alignnmnt
82   fOriginalAlign(0),   // original alignment 0 read for the OCDB
83   fNewAlign(0),
84   fPadTime0(0),
85   fFitCEGlobal(0),
86   fFitCELocal(0)
87 {
88   //
89   // Default constructor
90   //
91   for (Int_t i=0; i<4; i++){
92     fDelta1D[i]=0;
93     fCovar1D[i]=0;
94   }
95 }
96
97 AliTPCkalmanAlign::AliTPCkalmanAlign(const char* name, const char* title): 
98   TNamed(name,title),
99   fCalibAlign(0),     // kalman alignnmnt
100   fOriginalAlign(0),   // original alignment 0 read for the OCDB
101   fNewAlign(0),
102   fPadTime0(0),
103   fFitCEGlobal(0),
104   fFitCELocal(0) 
105 {
106   //
107   // Default constructor
108   //
109   for (Int_t i=0; i<4; i++){
110     fDelta1D[i]=0;
111     fCovar1D[i]=0;
112   }
113   fFitCEGlobal = new TObjArray(6); 
114   fFitCELocal  = new TObjArray(6); 
115   for (Int_t ipar=0; ipar<6;ipar++){
116     fFitCEGlobal->AddAt(new TVectorD(36),ipar);
117     fFitCELocal->AddAt(new TVectorD(36),ipar);
118   }
119 }
120
121 void AliTPCkalmanAlign::ReadAlign(const char *fname){
122   //
123   // Read old alignment used in the reco
124   // and the residual histograms
125   // WE ASSUME that the OCDB path is set in the same way as done in the calibration
126   //
127   TFile fcalib(fname);
128   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
129   fCalibAlign=0;
130   if (array) fCalibAlign=( AliTPCcalibAlign *)array->FindObject("alignTPC");
131   fCalibAlign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
132   //
133   // old alignment used
134   AliCDBEntry * cdbEntry= AliCDBManager::Instance()->Get("TPC/Align/Data");
135   fOriginalAlign =0;
136   if (cdbEntry){
137     fOriginalAlign = (TClonesArray*)cdbEntry->GetObject();
138   }
139   
140 }
141
142 void AliTPCkalmanAlign::BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t mean, Double_t sigma){
143   //
144   // Book Align 1D parameters and covariance
145   //
146   param.ResizeTo(72,1);
147   covar.ResizeTo(72,72);
148   for (Int_t i=0;i<72;i++){
149     param(i,0)=mean;
150     for (Int_t j=0;j<72;j++) covar(i,j)=0;      
151     covar(i,i)=sigma*sigma;
152   }
153 }
154
155
156 void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &vecXk, TMatrixD &covXk){
157   //
158   // Update 1D kalman filter
159   //
160   const Int_t knMeas=1;
161   const Int_t knElem=72;
162   TMatrixD mat1(72,72);            // update covariance matrix
163   TMatrixD matHk(1,knElem);        // vector to mesurement
164   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
165   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
166   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
167   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
168   TMatrixD covXk2(knElem,knElem);  // helper matrix
169   TMatrixD covXk3(knElem,knElem);  // helper matrix
170   TMatrixD vecZk(1,1);
171   TMatrixD measR(1,1);
172   vecZk(0,0)=delta;
173   measR(0,0)=sigma*sigma;
174   //
175   // reset matHk
176   for (Int_t iel=0;iel<knElem;iel++) 
177     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
178   //mat1
179   for (Int_t iel=0;iel<knElem;iel++) {
180     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
181     mat1(iel,iel)=1;
182   }
183   //
184   matHk(0, s1)=1;
185   matHk(0, s2)=-1;
186   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
187   matHkT=matHk.T(); matHk.T();
188   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
189   matSk.Invert();
190   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
191   vecXk += matKk*vecYk;                    //  updated vector 
192   covXk2= (mat1-(matKk*matHk));
193   covXk3 =  covXk2*covXk;          
194   covXk = covXk3;  
195 }
196
197
198
199
200 void AliTPCkalmanAlign::MakeGlobalAlign(){
201   //
202   // Combine all pairs of fitters and make global alignemnt
203   //
204   
205   AliTPCkalmanAlign &kalmanAlign=*this;
206   TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
207   Int_t run = AliCDBManager::Instance()->GetRun();
208   AliGRPObject * grp = AliTPCcalibDB::Instance()->GetGRP(run);
209   Float_t bz = AliTracker::GetBz();
210   UInt_t timeS = grp->GetTimeStart();
211   UInt_t timeE = grp->GetTimeEnd();
212   UInt_t  time = (timeS+timeE)/2;
213   
214   //
215   // get ce info
216   //
217   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
218   TVectorD paramCE[72];
219   TMatrixD covarCE[72];
220   Int_t  statUpDown=0;   // statistic up down
221   Int_t  statLeftRight=0;   // statistic left-right
222   Float_t chi2;
223   for (Int_t isec=0; isec<72; isec++){
224     AliTPCCalROC * roc0  = padTime0->GetCalROC(isec);
225     roc0->GlobalFit(0,kFALSE,paramCE[isec],covarCE[isec],chi2,0);
226     (*pcstream)<<"ceFit"<<
227       "isec="<<isec<<
228       "p0.="<<&paramCE[isec]<<
229       "\n";
230   }
231
232   DumpOldAlignment(pcstream);
233   const Int_t kMinEntries=400;
234   TMatrixD vec[5];
235   TMatrixD cov[5];
236   kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
237   kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
238   kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
239   kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
240   //
241   TVectorD delta(5);
242   TVectorD rms(5);
243   TVectorD stat(5);
244   TH1 * his=0;
245   for (Int_t is0=0;is0<72;is0++)
246     for (Int_t is1=0;is1<72;is1++){
247       //
248       //
249       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
250       if (!his) continue;
251       if (his->GetEntries()<kMinEntries) continue;
252       delta[0]=his->GetMean();
253       rms[0]=his->GetRMS();
254       stat[0]=his->GetEntries();
255       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
256       //     
257       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
258       if (!his) continue;
259       delta[1]=his->GetMean();
260       rms[1]=his->GetRMS();
261       stat[1]=his->GetEntries();
262       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
263       //
264       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
265       if (!his) continue;
266       delta[2] = his->GetMean();
267       rms[2]=his->GetRMS();
268       stat[2]=his->GetEntries();
269       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
270       //
271       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
272       if (!his) continue;
273       delta[3] = his->GetMean();       
274       rms[3]=his->GetRMS();
275       stat[3]=his->GetEntries();
276       kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
277       if (is1==is0+36) statUpDown+=Int_t(stat[0]);
278       if (is1==is0+35) statLeftRight+=Int_t(stat[0]);
279     }  
280
281   fDelta1D[0] = (TMatrixD*)vec[0].Clone();
282   fDelta1D[1] = (TMatrixD*)vec[1].Clone();
283   fDelta1D[2] = (TMatrixD*)vec[2].Clone();
284   fDelta1D[3] = (TMatrixD*)vec[3].Clone();
285   //
286   fCovar1D[0] = (TMatrixD*)cov[0].Clone();
287   fCovar1D[1] = (TMatrixD*)cov[1].Clone();
288   fCovar1D[2] = (TMatrixD*)cov[2].Clone();
289   fCovar1D[3] = (TMatrixD*)cov[3].Clone();
290   statUpDown/=36;
291   statLeftRight/=36;
292   MakeNewAlignment(kTRUE);
293   FitCE();
294   for (Int_t is0=0;is0<72;is0++)
295     for (Int_t is1=0;is1<72;is1++){
296       Bool_t isPair=kFALSE;
297       if (TMath::Abs(is0%18-is1%18)<2) isPair=kTRUE;
298       if (TMath::Abs(is0%18-is1%18)==17) isPair=kTRUE;
299       if (!isPair) continue;
300       stat[0]=0; stat[1]=0; stat[2]=0; stat[3]=0; 
301       //
302       //
303       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
304       if (his){
305         delta[0]=his->GetMean();
306         rms[0]=his->GetRMS();
307         stat[0]=his->GetEntries();
308       }
309       //     
310       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
311       if (his) {
312         delta[1]=his->GetMean();
313         rms[1]=his->GetRMS();
314         stat[1]=his->GetEntries();
315       }
316       //
317       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
318       if (his){
319         delta[2] = his->GetMean();
320         rms[2]=his->GetRMS();
321         stat[2]=his->GetEntries();
322       }
323       //
324       his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
325       if (his){
326         delta[3] = his->GetMean();       
327         rms[3]=his->GetRMS();
328         stat[3]=his->GetEntries();
329       }
330       TVectorD fceG[8],fceL[6];
331       if (fFitCEGlobal)
332         for (Int_t ipar=0; ipar<8;ipar++){
333           fceG[ipar].ResizeTo(36);
334           if (ipar<6) fceL[ipar].ResizeTo(36);
335           if (fFitCEGlobal->At(ipar)){
336             fceG[ipar]=*((TVectorD*)fFitCEGlobal->At(ipar));
337             if (ipar<6){
338               fceL[ipar]=*((TVectorD*)fFitCELocal->At(ipar));
339             }
340           }           
341         }
342       
343       (*pcstream)<<"kalmanAlignDebug"<<
344         "run="<<run<<            // runs
345         "time="<<time<<          // time
346         "timeE="<<timeE<<        // sart of tun time
347         "timeS="<<timeS<<        // end od run time
348         "bz="<<bz<<
349         "is0="<<is0<<
350         "is1="<<is1<<
351         "delta.="<<&delta<<      // alignment deltas
352         "rms.="<<&rms<<          // rms
353         "stat.="<<&stat<<
354         "vec0.="<<&vec[0]<<
355         "vec1.="<<&vec[1]<<
356         "vec2.="<<&vec[2]<<
357         "vec3.="<<&vec[3]<<
358         "pceIn0.="<<&paramCE[is0%36]<<        // default CE parameters
359         "pceOut0.="<<&paramCE[is0%36+36]<<
360         "pceIn1.="<<&paramCE[is1%36]<<
361         "pceOut1.="<<&paramCE[is1%36+36]<<
362         //                     // current CE parameters form last calibration - not used in Reco
363         "fceG0.="<<&fceG[0]<<  // global fit of CE  - offset
364         "fceG1.="<<&fceG[1]<<  // global fit of CE  - Gy gradient 
365         "fceG2.="<<&fceG[2]<<  // global fit of CE  - Gx gradient
366         "fceG3.="<<&fceG[3]<<  // global fit of CE  - IROC-OROC offset
367         "fceG4.="<<&fceG[4]<<  // global fit of CE  - commont slope LX
368         "fceG5.="<<&fceG[5]<<  // global fit of CE  - delta slope LX
369         "fceG6.="<<&fceG[6]<<  // global fit of CE  - common slope LY
370         "fceG7.="<<&fceG[7]<<  // global fit of CE  - delta slope LY
371         //
372         "fceL0.="<<&fceL[0]<<  // local fit of CE  - offset to mean
373         "fceL1.="<<&fceL[1]<<  // local fit of CE  - IROC-OROC offset
374         "fceL2.="<<&fceL[2]<<  // local fit of CE  - common slope LX
375         "fceL3.="<<&fceL[3]<<  // local fit of CE  - delta slope  LX
376         "fceL4.="<<&fceL[4]<<  // local fit of CE  - common slope LY
377         "fceL5.="<<&fceL[5]<<  // local fit of CE  - delta slope  LY
378         "\n";
379     }
380   
381   (*pcstream)<<"runSummary"<<
382     "run="<<run<<                      // run number 
383     "bz="<<bz<<                        // bz field
384     "statUpDown="<<statUpDown<<        // stat up-down
385     "statLeftRight="<<statLeftRight<<  // stat left-right
386     "\n";
387
388   delete pcstream;
389 }
390
391
392
393
394
395
396 void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad  *pad, Int_t ustartRun, Int_t uendRun,  const char* storagePath ){
397   //
398   // Update OCDB
399   // .x ConfigCalibTrain.C(117116)
400   // AliTPCcalibDB::Instance()->GetPulserTmean()
401   // pad->Add(-pad->GetMean())
402   AliCDBMetaData *metaData= new AliCDBMetaData();
403   metaData->SetObjectClassName("TObjArray");
404   metaData->SetResponsible("Marian Ivanov");
405   metaData->SetBeamPeriod(1);
406   metaData->SetAliRootVersion("05-25-01"); //root version
407   metaData->SetComment("Calibration of the z - time misalignment");
408   AliCDBId* id1=NULL;
409   id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
410   AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
411   gStorage->Put(pad, (*id1), metaData);
412 }
413
414
415
416 void AliTPCkalmanAlign::DrawDeltaAlign(){
417   //
418   // Draw the reuslts of the alignment
419   // Residual misalignment in respect with previous alignment shown
420   //
421   //
422   TFile f("AliTPCkalmanAlign.root","update");
423   TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
424   TH1::AddDirectory(0);
425   // 
426   treeDelta->SetAlias("sector","is0");
427   treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
428   treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
429   treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
430   treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
431   //
432   treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
433   treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
434   treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
435   treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
436   //
437   treeDelta->SetMarkerStyle(25);
438   treeDelta->SetMarkerColor(4);
439   treeDelta->SetLineColor(4);
440   const char *type[3]={"up-down","left-right","right-left"};
441   const char *gropt[3]={"alp","lp","lp"};
442   const char *varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10*vec0.fElements[is0]:sector"};
443   const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
444   const char *varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10*vec1.fElements[is0]:sector"};
445   const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
446   const char *varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000*vec3.fElements[is0]:sector"};
447   const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
448   const char *varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000*vec2.fElements[is0]:sector"};
449   const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
450   TLegend *legend = 0;
451   Int_t grcol[3]={2,1,4};
452   Int_t entries=0;
453   TGraph *grDelta[3]={0,0,0};
454   TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
455   TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
456   TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
457   TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
458   canvasDy->Divide(2,2);
459   canvasDz->Divide(2,2);
460   canvasDtheta->Divide(2,2);
461   canvasDphi->Divide(2,2);
462
463   //
464   // Dy
465   //
466   canvasDy->cd(1);
467   treeDelta->Draw("dYmeas:dYfit");
468   for (Int_t itype=0; itype<3; itype++){
469     canvasDy->cd(itype+2);
470     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
471     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
472     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
473     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
474     entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
475     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
476     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
477     for (Int_t i=0; i<3; i++) {
478       grDelta[i]->SetMaximum(1.5);
479       grDelta[i]->SetMinimum(-1);
480       grDelta[i]->SetTitle(type[i]);
481       grDelta[i]->SetMarkerColor(grcol[i]); 
482       grDelta[i]->SetLineColor(grcol[i]); 
483       grDelta[i]->SetMarkerStyle(25+i); 
484       grDelta[i]->GetXaxis()->SetTitle("sector"); 
485       grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)"); 
486       if (itype==2 && i>0) continue;
487       grDelta[i]->Draw(gropt[i]); 
488       legend->AddEntry(grDelta[i]);
489     }
490     legend->Draw();
491   }
492   //
493   // Dz
494   //
495   canvasDz->cd();
496   canvasDz->cd(1);
497   treeDelta->Draw("dZmeas:dZfit");
498   for (Int_t itype=0; itype<3; itype++){
499     canvasDz->cd(itype+2);
500     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
501     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
502     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
503     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
504     entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
505     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
506     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
507     for (Int_t i=0; i<3; i++) {
508       grDelta[i]->SetMaximum(1.5);
509       grDelta[i]->SetMinimum(-1);
510       grDelta[i]->SetTitle(type[i]);
511       grDelta[i]->SetMarkerColor(grcol[i]); 
512       grDelta[i]->SetLineColor(grcol[i]); 
513       grDelta[i]->SetMarkerStyle(25+i); 
514       grDelta[i]->GetXaxis()->SetTitle("sector"); 
515       grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)"); 
516       if (itype==2 && i>0) continue;
517       grDelta[i]->Draw(gropt[i]); 
518       legend->AddEntry(grDelta[i]);
519     }
520     legend->Draw();
521   }
522
523   //
524   // Dtheta
525   //
526   canvasDtheta->cd(1);
527   treeDelta->Draw("dThetameas:dThetafit");
528   for (Int_t itype=0; itype<3; itype++){
529     canvasDtheta->cd(itype+2);
530     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
531     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
532     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
533     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
534     entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
535     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
536     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
537     for (Int_t i=0; i<3; i++) {
538       grDelta[i]->SetMaximum(4.);
539       grDelta[i]->SetMinimum(-3.);
540       grDelta[i]->SetTitle(type[i]);
541       grDelta[i]->SetMarkerColor(grcol[i]); 
542       grDelta[i]->SetLineColor(grcol[i]); 
543       grDelta[i]->SetMarkerStyle(25+i); 
544       grDelta[i]->GetXaxis()->SetTitle("sector"); 
545       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)"); 
546       if (itype==2 && i>0) continue;
547       grDelta[i]->Draw(gropt[i]); 
548       legend->AddEntry(grDelta[i]);
549     }
550     legend->Draw();
551   }
552
553   //
554   // Dphi
555   //
556   canvasDphi->cd();
557   canvasDphi->cd(1);
558   treeDelta->Draw("dPhimeas:dPhifit");
559   for (Int_t itype=0; itype<3; itype++){
560     canvasDphi->cd(itype+2);
561     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
562     grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
563     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
564     grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
565     entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
566     grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
567     legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
568     for (Int_t i=0; i<3; i++) {
569       grDelta[i]->SetMaximum(4.);
570       grDelta[i]->SetMinimum(-3.);
571       grDelta[i]->SetTitle(type[i]);
572       grDelta[i]->SetMarkerColor(grcol[i]); 
573       grDelta[i]->SetLineColor(grcol[i]); 
574       grDelta[i]->SetMarkerStyle(25+i); 
575       grDelta[i]->GetXaxis()->SetTitle("sector"); 
576       grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)"); 
577       if (itype==2 && i>0) continue;
578       grDelta[i]->Draw(gropt[i]); 
579       legend->AddEntry(grDelta[i]);
580     }
581     legend->Draw();
582   }
583   //
584   //
585   f.cd();
586   canvasDy->Write();
587   canvasDz->Write();
588   canvasDtheta->Write();
589   canvasDphi->Write();
590   //  
591   //
592   //
593   TPostScript *ps = new TPostScript("alignReport.ps", 112); 
594   ps->NewPage();
595   canvasDy->Draw();
596   ps->NewPage();
597   canvasDz->Draw();
598   ps->NewPage();
599   canvasDtheta->Draw();
600   ps->NewPage();
601   canvasDphi->Draw();
602   ps->Close();
603   delete ps;
604 }
605
606
607
608 void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
609   // Dump the content of old alignemnt
610   // Expected that the old alignmnet is loaded
611   //
612   if (!fOriginalAlign) return;
613   //
614   TVectorD localTrans(3);
615   TVectorD globalTrans(3);
616   TVectorD localRot(3);
617   TVectorD globalRot(3);
618   AliGeomManager::ELayerID idLayer;
619   Int_t idModule=0;
620   //
621   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
622     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
623     params->GetVolUID(idLayer,idModule);
624     params->GetLocalTranslation(localTrans.GetMatrixArray());
625     params->GetLocalAngles(localRot.GetMatrixArray());
626     params->GetTranslation(globalTrans.GetMatrixArray());
627     params->GetAngles(globalRot.GetMatrixArray());
628     Int_t sector=idModule;
629     if (idLayer>7) sector+=36;
630     (*pcstream)<<"oldAlign"<<
631       //"idLayer="<<idLayer<<
632       "idModule="<<idModule<<
633       "sector="<<sector<<
634       "lT.="<<&localTrans<<
635       "gT.="<<&localTrans<<
636       "lR.="<<&localRot<<
637       "gR.="<<&globalRot<<
638       "\n";
639   }
640 }
641
642
643 void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
644   //
645   // make a new Alignment entry
646   //
647   if (!fOriginalAlign) return;
648   //
649   TVectorD localTrans(3);
650   TVectorD globalTrans(3);
651   TVectorD localRot(3);
652   TVectorD globalRot(3);
653   //
654   TVectorD localTransNew(3);   // new entries
655   TVectorD globalTransNew(3);
656   TVectorD localRotNew(3);
657   TVectorD globalRotNew(3);
658   //
659   AliGeomManager::ELayerID idLayer;
660   Int_t idModule=0;
661   //
662   fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
663   for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
664     AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
665     //AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
666     params->GetVolUID(idLayer,idModule);
667     Int_t sector=(Int_t)idModule;
668     if (idLayer>7) sector+=36;
669     params->GetLocalTranslation(localTrans.GetMatrixArray());
670     params->GetLocalAngles(localRot.GetMatrixArray());
671     params->GetTranslation(globalTrans.GetMatrixArray());
672     params->GetAngles(globalRot.GetMatrixArray());
673     //
674     //
675     //
676     if (badd){ // addition if 
677       localTransNew=localTrans;
678       localRotNew=localRot;
679     }
680     localTransNew[1]=localTransNew[1]-((*fDelta1D[0])(sector,0));
681     localRot[0]  =localRot[0]-(*fDelta1D[2])(sector,0);
682     //
683     if (pcstream) (*pcstream)<<"alignParams"<<
684       //"idLayer="<<idLayer<<
685       "idModule="<<idModule<<
686       "sector="<<sector<<
687       "olT.="<<&localTrans<<
688       "olR.="<<&localRot<<
689       "ogT.="<<&localTrans<<
690       "ogR.="<<&globalRot<<
691       "nlT.="<<&localTransNew<<
692       "nlR.="<<&localRotNew<<
693       "ngT.="<<&localTransNew<<
694       "ngR.="<<&globalRotNew<<
695       "\n";
696   }
697 }
698
699
700
701 void AliTPCkalmanAlign::DrawAlignmentTrends(){
702   //
703   // Draw trends of alingment variables
704   //
705   /*
706     //1.  Create a list of  align data
707     //
708     //2. Filter list
709     AliXRDPROOFtoolkit::FilterList("align.list","AliTPCkalmanAlign.root kalmanAlignDebug",0);
710
711   */
712   AliXRDPROOFtoolkit toolkit;
713   TChain * chain = toolkit.MakeChainRandom("align.list.Good","kalmanAlignDebug",0,2000);
714   TChain * chainRef = toolkit.MakeChainRandom("alignRef.list","kalmanAlignDebug",0,2000);
715   chain->AddFriend(chainRef,"R");
716   chainRef->AddFriend(chainRef,"T");
717   //cuts
718   TCut cutS="stat.fElements[0]>200&&stat.fElements[1]>200&&stat.fElements[3]>200&&stat.fElements[3]>200";   //statistic in the bin
719   TCut cutST="T.stat.fElements[0]>200&&T.stat.fElements[1]>200&&T.stat.fElements[3]>200&&T.stat.fElements[3]>200";   //statistic in the bin
720   //  TTree *tree  = chain->CopyTree(cutS);
721   //TTree *treeR = chainRef->CopyTree(cutST);
722
723   TCanvas * canvasDy= new TCanvas("canvasDy","canvasDy");
724   TH1 *his=0;
725   TLegend *legend = 0;
726   //  Int_t grcol[3]={2,1,4};
727   
728   legend = new TLegend(0.7,0.6,0.9,0.9, "Alignment #Delta_{y}- Up-Down");
729   for (Int_t isec=0; isec<18; isec+=2){
730     chain->SetMarkerColor(1+(isec%5));
731     chain->SetMarkerStyle(isec+20);
732     chain->Draw("10*(delta.fElements[0]-R.delta.fElements[0]):run",cutS+Form("is1==is0+36&&is0==%d",isec),"profgoff");
733     his = (TH1*)(chain->GetHistogram()->Clone());
734     his->SetName(Form("#Delta_{Y} sector %d",isec));
735     his->SetTitle(Form("#Delta_{Y} sector %d",isec));
736     his->SetMaximum(1.);
737     his->SetMinimum(-1.);
738     his->GetYaxis()->SetTitle("#Delta_{y} (mm)");
739     his->GetXaxis()->SetTitle("run Number");
740     if (isec==0) his->Draw("");
741     if (isec>0) his->Draw("same");
742     legend->AddEntry(his);
743   }
744   legend->Draw();
745   canvasDy->Draw();
746 }
747
748
749
750
751
752
753 void AliTPCkalmanAlign::FitCE(){
754   //
755   // fit CE 
756   // 1. Global fit - gy and gx
757   // 2. Local X fit common 
758   // 3. Sector fit 
759   //
760   AliTPCPreprocessorOnline * preprocesor = new AliTPCPreprocessorOnline;
761   //
762   AliTPCCalPad *padTime0 = AliTPCcalibDB::Instance()->GetPadTime0();
763   AliTPCCalPad *padNoise = AliTPCcalibDB::Instance()->GetPadNoise();
764   AliTPCCalPad * ceTmean = AliTPCcalibDB::Instance()->GetCETmean();  // CE information
765   AliTPCCalPad * ceTrms  = AliTPCcalibDB::Instance()->GetCETrms();
766   AliTPCCalPad * ceQmean = AliTPCcalibDB::Instance()->GetCEQmean();  
767   AliTPCCalPad * pulserTmean = AliTPCcalibDB::Instance()->GetPulserTmean(); //
768   AliTPCCalPad * pulserTrms = AliTPCcalibDB::Instance()->GetPulserTrms();
769   AliTPCCalPad * pulserQmean = AliTPCcalibDB::Instance()->GetPulserQmean();
770   AliTPCCalPad * dmap0  = AliTPCcalibDB::Instance()->GetDistortionMap(0);   // distortion maps
771   AliTPCCalPad * dmap1  = AliTPCcalibDB::Instance()->GetDistortionMap(1);
772   AliTPCCalPad * dmap2  = AliTPCcalibDB::Instance()->GetDistortionMap(2);
773   pulserTmean->Add(-pulserTmean->GetMean());
774   //
775   preprocesor->AddComponent(padTime0->Clone());
776   preprocesor->AddComponent(padNoise->Clone());
777   preprocesor->AddComponent(pulserTmean->Clone());
778   preprocesor->AddComponent(pulserQmean->Clone());
779   preprocesor->AddComponent(pulserTrms->Clone());
780   preprocesor->AddComponent(ceTmean->Clone());
781   preprocesor->AddComponent(ceQmean->Clone());
782   preprocesor->AddComponent(ceTrms->Clone());
783   preprocesor->AddComponent(dmap0->Clone());
784   preprocesor->AddComponent(dmap1->Clone());
785   preprocesor->AddComponent(dmap2->Clone());
786   preprocesor->DumpToFile("cetmean.root");
787
788   TCut cutNoise="abs(PadNoise.fElements/PadNoise_Median-1)<0.3";
789   TCut cutPulserT="abs(PulserTrms.fElements/PulserTrms_Median-1)<0.2";
790   TCut cutPulserQ="abs(PulserQmean.fElements/PulserQmean_Median-1)<0.2";
791   TCut cutCEQ="CEQmean.fElements>50";
792   TCut cutCET="abs(CETmean.fElements)<2";
793   TCut cutAll=cutNoise+cutPulserT+cutPulserQ+cutCEQ+cutCET;
794   //
795   //
796   TFile * f = new TFile("cetmean.root");
797   TTree * chain = (TTree*) f->Get("calPads");
798   Int_t entries = chain->Draw("1",cutAll,"goff");
799   if (entries<200000) return;  // no calibration available - pulser or CE or noise
800
801   TStatToolkit toolkit;
802   Double_t chi2=0;
803   Int_t    npoints=0;
804   TVectorD param;
805   TMatrixD covar;
806   //
807   // make a aliases
808   AliTPCkalmanAlign::MakeAliasCE(chain);
809   TString  fstringG="";              // global part
810   //
811   fstringG+="Gy++";                  // 1 - global y
812   fstringG+="Gx++";                  // 2 - global x
813   // 
814   fstringG+="isin++";                // 3 -delta IROC-OROC offset
815   fstringG+="Lx++";                  // 4 -common slope 
816   fstringG+="Lx*isin++";             // 5 -delta slope 
817   fstringG+="Ly++";                  // 6 -common slope 
818   fstringG+="Ly*isin++";             // 7 -delta slope 
819   TVectorD vecG[2];
820   TString * strFitG=0;
821   TString * strFitLX=0;
822   //
823   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideA"+cutAll, chi2,npoints,vecG[0],covar,-1,0, 10000000, kFALSE);
824   chain->SetAlias("tfitGA",strFitG->Data());
825   strFitG->Tokenize("++")->Print();
826   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
827   //
828   strFitG = TStatToolkit::FitPlane(chain,"deltaT", fstringG.Data(),"sideC"+cutAll, chi2,npoints,vecG[1],covar,-1,0, 10000000, kFALSE);
829   chain->SetAlias("tfitGC",strFitG->Data());
830   strFitG->Tokenize("++")->Print();
831   printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
832   //
833   AliTPCCalPad *padFitG =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[0],vecG[1]);
834   AliTPCCalPad *padFitLX=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[0],vecG[1]);
835   // swap a side and c side
836   AliTPCCalPad *padFitGSwap =AliTPCCalPad::CreateCalPadFit("1++gy/500.++gx/500.++0+++0++0++0++0",vecG[1],vecG[0]);
837   AliTPCCalPad *padFitLXSwap=AliTPCCalPad::CreateCalPadFit("0++0++0++(sector<36)++(lx-133)/100++(sector<36)*(lx-133)/100.++(ly)/100++(sector<36)*(ly)/100.",vecG[1],vecG[0]);
838   padFitG->SetName("CEG");
839   padFitLX->SetName("CELX");
840   padFitGSwap->SetName("CEGS");
841   padFitLXSwap->SetName("CELXS");
842   preprocesor->AddComponent(padFitG->Clone());
843   preprocesor->AddComponent(padFitLX->Clone());
844   preprocesor->AddComponent(padFitGSwap->Clone());
845   preprocesor->AddComponent(padFitLXSwap->Clone());
846   preprocesor->DumpToFile("cetmean.root");   // add it to the file
847   //
848   // make local fits
849   //
850   f = new TFile("cetmean.root");
851   chain = (TTree*) f->Get("calPads");
852   AliTPCkalmanAlign::MakeAliasCE(chain);
853   TString  fstringL="";              // local fit
854   //                                 // 0. delta common
855   fstringL+="isin++";                // 1. delta IROC-OROC offset
856   fstringL+="Lx++";                  // 2. common slope 
857   fstringL+="Lx*isin++";             // 3. delta slope 
858   fstringL+="Ly++";                  // 4. common slope 
859   fstringL+="Ly*isin++";             // 5. delta slope 
860   TVectorD vecL[36];
861   TVectorD dummy(6);
862   AliTPCCalPad *padFitLCE = new AliTPCCalPad("LocalCE","LocalCE");
863   AliTPCCalPad *padFitTmpCE;
864   for (Int_t isec=0; isec<36; isec++){
865     TCut cutSector=Form("(sector%36)==%d",isec);
866     strFitLX = TStatToolkit::FitPlane(chain,"deltaT-CEG.fElements-CELX.fElements", fstringL.Data(),cutSector+cutAll+"abs(deltaT-CEG.fElements-CELX.fElements)<0.4", chi2,npoints,vecL[isec],covar,-1,0, 10000000, kFALSE);
867     printf("sec=%d\tchi2=%f\n",isec,TMath::Sqrt(chi2/npoints));
868     //
869     TString fitL=Form("((sector%36)==%d)++((sector%36)==%d)*(sector<36)++((sector%36)==%d)*(lx-133)/100.++((sector%36)==%d)*(sector<36)*(lx-133)/100.++((sector%36)==%d)*(ly)/100.++((sector%36)==%d)*(sector<36)*(ly)/100.",isec,isec,isec,isec);
870     if (isec<18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),vecL[isec],dummy);
871     if (isec>=18) padFitTmpCE=AliTPCCalPad::CreateCalPadFit(fitL.Data(),dummy,vecL[isec]);
872     padFitLCE->Add(padFitTmpCE);
873   }
874   //
875   padFitLCE->SetName("CELocal");
876   preprocesor->AddComponent(padFitLCE->Clone());
877   preprocesor->DumpToFile("cetmean.root");   // add it to the file
878   //
879   // write data to array
880   //
881   fFitCELocal  = new TObjArray(6); 
882   fFitCEGlobal = new TObjArray(8); 
883   for (Int_t ipar=0; ipar<8;ipar++){
884     //
885     fFitCEGlobal->AddAt(new TVectorD(36),ipar);
886     TVectorD &fvecG = *((TVectorD*)fFitCEGlobal->At(ipar));
887     for (Int_t isec=0; isec<36;isec++){      
888       if (isec<18)  fvecG[isec]=vecG[0][ipar];
889       if (isec>=18) fvecG[isec]=vecG[1][ipar];    
890     }
891     if (ipar<6){
892       fFitCELocal->AddAt(new TVectorD(36),ipar);
893       TVectorD &fvecL = *((TVectorD*)fFitCELocal->At(ipar));
894       for (Int_t isec=0; isec<36;isec++){      
895         fvecL[isec]=vecL[isec][ipar];
896       }
897     }
898   }
899 }
900
901 void AliTPCkalmanAlign::MakeAliasCE(TTree * chain){
902   //
903   // make a aliases of pad variables
904   //
905   chain->SetAlias("side","(-1+(sector%36<18)*2)");
906   chain->SetAlias("sideA","(sector%36<18)");
907   chain->SetAlias("sideC","(sector%36>=18)");
908   chain->SetAlias("isin","(sector<36)");
909   chain->SetAlias("deltaT","CETmean.fElements-PulserTmean.fElements");
910   chain->SetAlias("timeP","PulserTmean.fElements");
911   chain->SetAlias("Gy","(gy.fElements/500.)");
912   chain->SetAlias("Gx","(gx.fElements/500.)");
913   chain->SetAlias("Lx","(lx.fElements-133)/100.");   // lx in meters
914   chain->SetAlias("Ly","(ly.fElements)/100.");
915   chain->SetAlias("La","(ly.fElements/lx.fElements/0.155)");
916   chain->SetAlias("deltaT","(CETmean.fElements-PulserTmean.fElements)");
917 }