AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FitAlignCombined.C
Go to the documentation of this file.
1 
37 #if !defined(__CINT__) || defined(__MAKECINT__)
38 #include "TH1D.h"
39 #include "TH2F.h"
40 #include "THnSparse.h"
41 #include "TVectorD.h"
42 #include "TTreeStream.h"
43 #include "TFile.h"
44 #include "TChain.h"
45 #include "AliTPCcalibAlign.h"
46 #include "AliTPCROC.h"
47 #include "AliTPCCalROC.h"
48 #include "AliTPCCalPad.h"
49 #include "TF1.h"
50 #include "TH2.h"
51 #include "TH3.h"
52 #include "TROOT.h"
53 #include "TProfile.h"
54 #include "AliTPCPreprocessorOnline.h"
55 #include "AliTPCcalibDB.h"
56 #include "AliTPCkalmanAlign.h"
57 #include "TPostScript.h"
58 #include "TLegend.h"
59 #include "TCut.h"
60 #include "TCanvas.h"
61 #include "TStyle.h"
62 #include "AliLog.h"
63 #include "AliTPCExBEffectiveSector.h"
64 #include "AliTPCComposedCorrection.h"
65 #include "AliTPCCalibGlobalMisalignment.h"
66 //
67 #include "AliCDBMetaData.h"
68 #include "AliCDBId.h"
69 #include "AliCDBManager.h"
70 #include "AliCDBStorage.h"
71 #include "AliCDBEntry.h"
72 #include "TStatToolkit.h"
73 #include "AliAlignObjParams.h"
74 #include "AliTPCParam.h"
75 #endif
76 
77 
78 void DrawFitQA();
79 
80 AliTPCROC *proc = AliTPCROC::Instance();
81 AliTPCCalibGlobalMisalignment *combAlignOCDBOld=0; // old misalignment from OCDB - used for reconstruction
82 AliTPCCalibGlobalMisalignment *combAlignOCDBNew=0; // new misalignment
83 AliTPCCalibGlobalMisalignment *combAlignGlobal=0; // delta misalignment - global part
84 AliTPCCalibGlobalMisalignment *combAlignLocal=0; // delta misalignment - sector/local part
85 //
86 
87 TTreeSRedirector *pcstream= 0; // Workspace to store the fit parameters and QA histograms
88 //
89 TChain * chain=0;
90 TChain * chainZ=0; // trees with z measurement
91 TTree * tree=0;
92 TTree *itsdy=0;
93 TTree *itsdyP=0;
94 TTree *itsdyM=0;
95 TTree *itsdp=0;
96 TTree *itsdphiP=0;
97 TTree *itsdphiM=0;
98 
99 TTree *itsdz=0;
100 TTree *itsdt=0;
101 TTree *tofdy=0;
102 TTree *trddy=0;
103 //
104 TTree *vdy=0;
105 TTree *vdyP=0;
106 TTree *vdyM=0;
107 TTree *vds=0;
108 TTree *vdz=0;
109 TTree *vdt=0;
110 
111 TCut cutS="entries>1000&&abs(snp)<0.2&&abs(theta)<1.";
112 
113 
114 void RegisterAlignFunction(){
125 
126  TGeoHMatrix matrixX;
127  TGeoHMatrix matrixY;
128  TGeoHMatrix matrixZ;
129  TGeoRotation rot0;
130  TGeoRotation rot1;
131  TGeoRotation rot2; //transformation matrices
132  TGeoRotation rot90; //transformation matrices
133  matrixX.SetDx(0.1); matrixY.SetDy(0.1); matrixZ.SetDz(0.1); //1 mm translation
134  rot0.SetAngles(0.001*TMath::RadToDeg(),0,0);
135  rot1.SetAngles(0,0.001*TMath::RadToDeg(),0);
136  rot2.SetAngles(0,0,0.001*TMath::RadToDeg());
137  //how to get rot02 ?
138  rot90.SetAngles(0,90,0);
139  rot2.MultiplyBy(&rot90,kTRUE);
140  rot90.SetAngles(0,-90,0);
141  rot2.MultiplyBy(&rot90,kFALSE);
142  AliTPCCalibGlobalMisalignment *alignRot0 =new AliTPCCalibGlobalMisalignment;
143  alignRot0->SetAlignGlobal(&rot0);
144  AliTPCCalibGlobalMisalignment *alignRot1=new AliTPCCalibGlobalMisalignment;
145  alignRot1->SetAlignGlobal(&rot1);
146  AliTPCCalibGlobalMisalignment *alignRot2=new AliTPCCalibGlobalMisalignment;
147  alignRot2->SetAlignGlobal(&rot2);
148  //
149  AliTPCCalibGlobalMisalignment *alignTrans0 =new AliTPCCalibGlobalMisalignment;
150  alignTrans0->SetAlignGlobal(&matrixX);
151  AliTPCCalibGlobalMisalignment *alignTrans1=new AliTPCCalibGlobalMisalignment;
152  alignTrans1->SetAlignGlobal(&matrixY);
153  AliTPCCalibGlobalMisalignment *alignTrans2=new AliTPCCalibGlobalMisalignment;
154  alignTrans2->SetAlignGlobal(&matrixZ);
155  AliTPCCorrection::AddVisualCorrection(alignTrans0 ,0);
156  AliTPCCorrection::AddVisualCorrection(alignTrans1 ,1);
157  AliTPCCorrection::AddVisualCorrection(alignTrans2 ,2);
158  AliTPCCorrection::AddVisualCorrection(alignRot0 ,3);
159  AliTPCCorrection::AddVisualCorrection(alignRot1 ,4);
160  AliTPCCorrection::AddVisualCorrection(alignRot2 ,5);
161  TObjArray arrAlign(6);
162  arrAlign.AddAt(alignTrans0->Clone(),0);
163  arrAlign.AddAt(alignTrans1->Clone(),1);
164  arrAlign.AddAt(alignTrans2->Clone(),2);
165  arrAlign.AddAt(alignRot0->Clone(),3);
166  arrAlign.AddAt(alignRot1->Clone(),4);
167  arrAlign.AddAt(alignRot2->Clone(),5);
168  //combAlign.SetCorrections((TObjArray*)arrAlign.Clone());
169 }
170 
171 AliTPCCalibGlobalMisalignment * MakeAlignFunctionGlobal(TVectorD paramYGlobal){
176 
177  AliTPCCalibGlobalMisalignment *alignGlobal =new AliTPCCalibGlobalMisalignment;
178  TGeoHMatrix matGlobal; // global parameters
179  TGeoHMatrix matDelta; // delta A side - C side
180  //
181  TGeoHMatrix matrixX;
182  TGeoHMatrix matrixY;
183  TGeoHMatrix matrixZ;
184  TGeoRotation rot0;
185  TGeoRotation rot1;
186  TGeoRotation rot2; //transformation matrices
187  TGeoRotation rot90; //transformation matrices
188  //
189  //
190  matrixX.SetDx(0.1*paramYGlobal[1]);
191  matrixY.SetDy(0.1*paramYGlobal[2]);
192  rot0.SetAngles(0.001*TMath::RadToDeg()*paramYGlobal[3],0,0);
193  rot1.SetAngles(0,0.001*TMath::RadToDeg()*paramYGlobal[4],0);
194  rot2.SetAngles(0,0,0.001*TMath::RadToDeg()*paramYGlobal[5]);
195  rot90.SetAngles(0,90,0);
196  rot2.MultiplyBy(&rot90,kTRUE);
197  rot90.SetAngles(0,-90,0);
198  rot2.MultiplyBy(&rot90,kFALSE);
199  matGlobal.Multiply(&matrixX);
200  matGlobal.Multiply(&matrixY);
201  matGlobal.Multiply(&rot0);
202  matGlobal.Multiply(&rot1);
203  matGlobal.Multiply(&rot2);
204  alignGlobal->SetAlignGlobal((TGeoMatrix*)matGlobal.Clone());
205  AliTPCCorrection::AddVisualCorrection(alignGlobal ,100);
206  //
207  /*
208  chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame");
209  chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame");
210 
211  chain->Draw("deltaG:fitYGlobal",cutS+"type==0");
212  chain->Draw("deltaG:fitYGlobal",cutS+"type==2");
213 
214  */
215 
216  return alignGlobal;
217 }
218 
219 
220 AliTPCCalibGlobalMisalignment * MakeAlignFunctionSector(TVectorD paramYLocal){
227 
228  AliTPCCalibGlobalMisalignment *alignLocal =new AliTPCCalibGlobalMisalignment;
229  TGeoHMatrix matrixX;
230  TGeoHMatrix matrixY;
231  TGeoHMatrix matrixZ;
232  TGeoRotation rot0;
233  TGeoRotation rot1;
234  TGeoRotation rot2; //transformation matrices
235  TGeoRotation rot90; //transformation matrices
236  TObjArray * array = new TObjArray(72);
237  TVectorD vecLY(72);
238  TVectorD vecGY(72);
239  TVectorD vecGX(72);
240  TVectorD vecPhi(72);
241  TVectorD vecSec(72);
242  //
243  //
244  {for (Int_t isec=0; isec<72; isec++){
245  Int_t sector=isec%18;
246  Int_t offset=sector*4+1;
247  if ((isec%36)>=18) offset+=2;
248  Double_t phi= (Double_t(sector)+0.5)*TMath::Pi()/9.;
249  //
250  Double_t dly = paramYLocal[offset];
251  Double_t dphi= paramYLocal[offset+1];
252  Double_t dgx = TMath::Cos(phi)*0+TMath::Sin(phi)*dly;
253  Double_t dgy = TMath::Sin(phi)*0-TMath::Cos(phi)*dly;
254  vecSec[isec]=isec;
255  vecLY[isec]=dly;
256  vecGX[isec]=dgx;
257  vecGY[isec]=dgy;
258  vecPhi[isec]=dphi;
259  //
260  matrixX.SetDx(dgx);
261  matrixY.SetDy(dgy);
262  rot0.SetAngles(0.001*TMath::RadToDeg()*dphi,0,0);
263  TGeoHMatrix matrixSector; // global parameters
264  matrixSector.Multiply(&matrixX);
265  matrixSector.Multiply(&matrixY);
266  matrixSector.Multiply(&rot0);
267  array->AddAt(matrixSector.Clone(),isec);
268  }}
269  alignLocal->SetAlignSectors(array);
270  AliTPCCorrection::AddVisualCorrection(alignLocal,101);
271  //
272  //
273  TGraph * graphLY = new TGraph(72,vecSec.GetMatrixArray(), vecLY.GetMatrixArray());
274  TGraph * graphGX = new TGraph(72,vecSec.GetMatrixArray(), vecGX.GetMatrixArray());
275  TGraph * graphGY = new TGraph(72,vecSec.GetMatrixArray(), vecGY.GetMatrixArray());
276  TGraph * graphPhi = new TGraph(72,vecSec.GetMatrixArray(), vecPhi.GetMatrixArray());
277  graphLY->SetMarkerStyle(25); graphGX->SetMarkerStyle(25); graphGY->SetMarkerStyle(25); graphPhi->SetMarkerStyle(25);
278  graphLY->SetName("DeltaLY"); graphLY->SetTitle("#Delta_{ly} (mm)");
279  graphGX->SetName("DeltaGX"); graphGX->SetTitle("#Delta_{gx} (mm)");
280  graphGY->SetName("DeltaGY"); graphGY->SetTitle("#Delta_{gy} (mm)");
281  graphPhi->SetName("DeltaPhi"); graphPhi->SetTitle("#Delta_{phi} (mrad)");
282  //
283  graphLY->Write("grDeltaLY");
284  graphGX->Write("grDeltaGX");
285  graphGY->Write("grDeltaGY");
286  graphPhi->Write("grDeltaPhi");
287  // Check:
288  /*
289  graphLY.Draw("alp")
290  chain->Draw("mean-deltaG:int(sector)",cutS+"type==0&&refX==0&&theta>0","profsame");
291  chain->Draw("mean-deltaG:int(sector+18)",cutS+"type==0&&refX==0&&theta<0","profsame");
292 
293  graphPhi.Draw("alp")
294  chain->Draw("1000*(mean-deltaG):int(sector)",cutS+"type==2&&refX==0&&theta>0","profsame");
295  chain->Draw("1000*(mean-deltaG):int(sector+18)",cutS+"type==2&&refX==0&&theta<0","profsame");
296 
297  //
298  chain->Draw("deltaL:fitYLocal",cutS+"type==2&&refX==0",""); // phi shift - OK
299  chain->Draw("deltaL:fitYLocal",cutS+"type==0&&refX==0",""); // OK
300  chain->Draw("deltaL:fitYLocal",cutS+"type==0",""); // OK
301 
302  */
303 
304  return alignLocal;
305 }
306 
307 
308 
309 
310 
311 void LoadTrees(){
317 
318  TFile *f0 = new TFile("../mergeField0/mean.root");
319  TFile *fP= new TFile("../mergePlus/mean.root");
320  TFile *fM= new TFile("../mergeMinus/mean.root");
321  //
322  itsdy=(TTree*)f0->Get("ITSdy");
323  itsdyP=(TTree*)fP->Get("ITSdy");
324  itsdyM=(TTree*)fM->Get("ITSdy");
325  itsdp=(TTree*)f0->Get("ITSdsnp");
326  itsdphiP=(TTree*)fP->Get("ITSdsnp");
327  itsdphiM=(TTree*)fM->Get("ITSdsnp");
328 
329  itsdz=(TTree*)f0->Get("ITSdz");
330  itsdt=(TTree*)f0->Get("ITSdtheta");
331  tofdy=(TTree*)f0->Get("TOFdy");
332  trddy=(TTree*)f0->Get("TRDdy");
333  //
334  vdy=(TTree*)f0->Get("Vertexdy");
335  vdyP=(TTree*)fP->Get("Vertexdy");
336  vdyM=(TTree*)fM->Get("Vertexdy");
337  vds=(TTree*)f0->Get("Vertexdsnp");
338  vdz=(TTree*)f0->Get("Vertexdz");
339  vdt=(TTree*)f0->Get("Vertexdtheta");
340  chain = new TChain("mean","mean");
341  chainZ = new TChain("mean","mean");
342  chain->AddFile("../mergeField0/mean.root",10000000,"ITSdy");
343  chain->AddFile("../mergeField0/mean.root",10000000,"ITSdsnp");
344  chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdy");
345  chain->AddFile("../mergeField0/mean.root",10000000,"Vertexdsnp");
346  chain->AddFile("../mergeField0/mean.root",10000000,"TOFdy");
347  chain->AddFile("../mergeField0/mean.root",10000000,"TRDdy");
348  //
349  chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdz");
350  chainZ->AddFile("../mergeField0/mean.root",10000000,"Vertexdtheta");
351  chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdz");
352  chainZ->AddFile("../mergeField0/mean.root",10000000,"ITSdtheta");
353 
354  //
355  itsdy->AddFriend(itsdp,"Phi");
356  itsdy->AddFriend(itsdphiP,"PhiP");
357  itsdy->AddFriend(itsdphiM,"PhiM");
358  itsdy->AddFriend(itsdz,"Z");
359  itsdy->AddFriend(itsdt,"T");
360  itsdy->AddFriend(itsdyP,"YP");
361  itsdy->AddFriend(itsdyM,"YM");
362  //
363  itsdy->AddFriend(vdy,"V");
364  itsdy->AddFriend(vdyP,"VP");
365  itsdy->AddFriend(vdyP,"VM");
366  itsdy->AddFriend(tofdy,"TOF");
367  itsdy->AddFriend(trddy,"TRD");
368  itsdy->AddFriend(vds,"VPhi");
369  itsdy->AddFriend(vdz,"VZ");
370  itsdy->AddFriend(vdt,"VT");
371  itsdy->AddFriend(tofdy,"TOF.");
372  itsdy->SetMarkerStyle(25);
373  itsdy->SetMarkerSize(0.4);
374  tree=itsdy;
375  tree->SetAlias("side","(-1+2*(theta>0))");
376  chain->SetAlias("side","(-1+2*(theta>0))");
377  chain->SetMarkerStyle(25);
378  chain->SetMarkerSize(0.5);
379 }
380 
381 
382 void FitAlignCombinedCorr(){
384 
385  RegisterAlignFunction();
386  LoadTrees();
387  combAlignOCDBOld = AliTPCCalibGlobalMisalignment::CreateOCDBAlign();
388  //
389 
390  pcstream= new TTreeSRedirector("fitAlignCombined.root");
391 
392  TStatToolkit toolkit;
393  Double_t chi2=0;
394  Int_t npoints=0;
395  // Fit vectors
396  // global fits
397  TVectorD vecSec(37);
398  TVectorD paramYGlobal;
399  TVectorD paramYLocal;
400  TMatrixD covar;
401  // make Aliases
402  {for (Int_t isec=0; isec<18; isec++){ // sectors
403  tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
404  chain->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
405  }}
406  for (Int_t isec=-1;isec<36; isec++){
407  vecSec[isec+1]=isec;
408  }
409  chain->SetAlias("err","((type==0)*0.1+(type==2)*0.001)");
410  // delta phi
411  chain->SetAlias("dpgx","(0+0)");
412  chain->SetAlias("dpgy","(0+0)");
413  chain->SetAlias("dpgr0","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,3)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,3))/160)");
414  chain->SetAlias("dpgr1","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,4)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,4))/160)");
415  chain->SetAlias("dpgr2","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,5)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,5))/160)");
416  chain->SetAlias("deltaPG","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,100)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,100))/160)");
417  chain->SetAlias("deltaPL","((AliTPCCorrection::GetCorrSector(sector,245,theta,1,101)-AliTPCCorrection::GetCorrSector(sector,85,theta,1,101))/160)");
418  // delta y at 85 cm
419  chain->SetAlias("dygxT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,0)+0)");
420  chain->SetAlias("dygyT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,1)+0)");
421  chain->SetAlias("dyr0T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,3)+0)");
422  chain->SetAlias("dyr1T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,4)+0)");
423  chain->SetAlias("dyr2T","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,5)+0)");
424  chain->SetAlias("deltaYGT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,100)+0)");
425  chain->SetAlias("deltaYLT","(AliTPCCorrection::GetCorrSector(sector,85,theta,1,101)+0)");
426  //
427  // delta y at reference X
428  chain->SetAlias("dygx","(dygxT+0)"); // due global x shift
429  chain->SetAlias("dygy","(dygyT+0)"); // due global y shift
430  chain->SetAlias("dyr0","(dyr0T+dpgr0*(refX-85.))"); // due rotation 0
431  chain->SetAlias("dyr1","(dyr1T+dpgr1*(refX-85.))"); // due rotation 1
432  chain->SetAlias("dyr2","(dyr2T+dpgr2*(refX-85.))"); // due rotation 2
433  chain->SetAlias("deltaYG","(deltaYGT+deltaPG*(refX-85.))"); // due global distortion
434  chain->SetAlias("deltaYL","(deltaYLT+deltaPL*(refX-85.))"); // due local distortion
435 
436  chain->SetAlias("deltaG","(type==0)*deltaYG+(type==2)*deltaPG"); //alias to global function
437  chain->SetAlias("deltaL","(type==0)*deltaYL+(type==2)*deltaPL"); //alias to local function
438  //
439  TString fstringGlobal="";
440  // global part
441  fstringGlobal+="((type==0)*dygx+0)++";
442  fstringGlobal+="((type==0)*dygy+0)++";
443  fstringGlobal+="((type==0)*dyr0+(type==2)*dpgr0)++";
444  fstringGlobal+="((type==0)*dyr1+(type==2)*dpgr1)++";
445  fstringGlobal+="((type==0)*dyr2+(type==2)*dpgr2)++";
446  //
447  // Make global fits
448  //
449  TString *strFitYGlobal = TStatToolkit::FitPlane(chain,"mean:err", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYGlobal,covar,-1,0, 10000000, kTRUE);
450  combAlignGlobal =MakeAlignFunctionGlobal(paramYGlobal);
451  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
452  chain->SetAlias("fitYGlobal",strFitYGlobal->Data());
453  paramYGlobal.Print();
454  paramYGlobal.Write("paramYGlobal");
455  //
456  //
457  //
458  //
459  {for (Int_t isec=0; isec<18; isec++){
460  //A side
461  chain->SetAlias(Form("glyASec%d",isec),Form("((type==0)*sec%d*(theta>0))",isec)); // ly shift
462  chain->SetAlias(Form("gxASec%d",isec),Form("((type==0)*dygx*sec%d*(theta>0))",isec)); // gx shift
463  chain->SetAlias(Form("gyASec%d",isec),Form("((type==0)*dygy*sec%d*(theta>0))",isec)); // gy shift
464  chain->SetAlias(Form("gpASec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta>0)",isec)); // phi rotation
465  //C side
466  chain->SetAlias(Form("glyCSec%d",isec),Form("((type==0)*sec%d*(theta<0))",isec)); // ly shift
467  chain->SetAlias(Form("gxCSec%d",isec),Form("((type==0)*dygx*sec%d*(theta<0))",isec));
468  chain->SetAlias(Form("gyCSec%d",isec),Form("((type==0)*dygy*sec%d*(theta<0))",isec));
469  chain->SetAlias(Form("gpCSec%d",isec),Form("(((type==0)*dyr0+(type==2)*dpgr0)*sec%d)*(theta<0)",isec)); // phi rotation
470  }}
471 
472  TString fstringLocal="";
473  {for (Int_t isec=0; isec<18; isec++){
474  //fstringLocal+=Form("((type==0)*dygx*sec%d*(theta>0))++",isec);
475  // fstringLocal+=Form("(gxASec%d)++",isec);
476  // fstringLocal+=Form("(gyASec%d)++",isec);
477  fstringLocal+=Form("(glyASec%d)++",isec);
478  fstringLocal+=Form("(gpASec%d)++",isec);
479  fstringLocal+=Form("(glyCSec%d)++",isec);
480  //fstringLocal+=Form("(gxCSec%d)++",isec);
481  //fstringLocal+=Form("(gyCSec%d)++",isec);
482  fstringLocal+=Form("(gpCSec%d)++",isec);
483  }}
484 
485  TString *strFitYLocal = TStatToolkit::FitPlane(chain,"(mean-deltaG):err", fstringLocal.Data(),cutS+"", chi2,npoints,paramYLocal,covar,-1,0, 10000000, kTRUE);
486  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
487  chain->SetAlias("fitYLocal",strFitYLocal->Data());
488  combAlignLocal =MakeAlignFunctionSector(paramYLocal);
489  //paramYLocal.Print();
490  paramYLocal.Write("paramYLocal");
491  //
492  // scaling - why do we need it?
493  TString fstringScale="";
494  fstringScale+="((type==1)*refX+(type==2))*dsec++";
495  fstringScale+="((type==1)*refX+(type==2))*dsec*cos(phi)++";
496  fstringScale+="((type==1)*refX+(type==2))*dsec*sin(phi)++";
497  fstringScale+="((type==1)*refX+(type==2))*dsec*side++";
498  fstringScale+="((type==1)*refX+(type==2))*dsec*side*cos(phi)++";
499  fstringScale+="((type==1)*refX+(type==2))*dsec*side*sin(phi)++";
500  //
501  //
502  //
503  pcstream->GetFile()->cd();
504  DrawFitQA();
505  //
506  AliTPCCalibGlobalMisalignment * combAlignOCDBNew0 =( AliTPCCalibGlobalMisalignment *)combAlignOCDBOld->Clone();
507  combAlignOCDBNew0->AddAlign(*combAlignGlobal);
508  combAlignOCDBNew0->AddAlign(*combAlignLocal);
509  combAlignOCDBNew= AliTPCCalibGlobalMisalignment::CreateMeanAlign(combAlignOCDBNew0);
510  //
511  combAlignGlobal->SetName("fcombAlignGlobal");
512  combAlignLocal->SetName("fcombAlignLocal");
513  combAlignOCDBOld->SetName("fcombAlignOCDBOld");
514  combAlignOCDBNew->SetName("fcombAlignOCDBNew");
515  combAlignOCDBNew0->SetName("fcombAlignOCDBNew0");
516  //
517  combAlignGlobal->Write("fcombAlignGlobal");
518  combAlignLocal->Write("fcombAlignLocal");
519  combAlignOCDBOld->Write("fcombAlignOCDBOld");
520  combAlignOCDBNew->Write("fcombAlignOCDBNew");
521  combAlignOCDBNew0->Write("fcombAlignOCDBNew0");
522  combAlignOCDBNew->Print();
523  delete pcstream;
524 }
525 
526 void DrawFitQA(){
528 
529  TCanvas c;
530  c.SetLeftMargin(0.15);
531  chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&refX==0","");
532  chain->GetHistogram()->SetName("TPC_VertexDphiGlobal1D");
533  chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)");
534  chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)");
535  chain->GetHistogram()->Fit("gaus","qnr");
536  chain->GetHistogram()->Write();
537  //
538  chain->Draw("1000*(mean-deltaG)>>his(100,-1.5,1.5)",cutS+"type==2&&abs(refX-85)<2","");
539  chain->GetHistogram()->SetName("TPC_ITSDphiGlobal1D");
540  chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)");
541  chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{#phi} (mrad)");
542  chain->GetHistogram()->Fit("gaus","qnr");
543  chain->GetHistogram()->Write();
544  //
545  chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-85)<2","");
546  chain->GetHistogram()->SetName("TPC_ITSDRPhiGlobal1D");
547  chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)");
548  chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)");
549  chain->GetHistogram()->Fit("gaus","qnr");
550  chain->GetHistogram()->Write();
551  //
552  chain->Draw("10*(mean-deltaG)>>his(100,-1,1)",cutS+"type==0&&abs(refX-0)<2","");
553  chain->GetHistogram()->SetName("TPC-VertexDRPhiGlobal1D");
554  chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)");
555  chain->GetHistogram()->GetXaxis()->SetTitle("#Delta_{r#phi} (mm)");
556  chain->GetHistogram()->Fit("gaus","qnr");
557  chain->GetHistogram()->Write();
558  //
559  // Make QA plot 3D
560  //
561  chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz");
562  chain->GetHistogram()->SetName("TPC_VertexDphi3DGlobal3D");
563  chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Global Fit)");
564  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
565  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
566  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
567  chain->GetHistogram()->Draw("colz");
568  chain->GetHistogram()->Write();
569  //
570  chain->Draw("1000*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz");
571  chain->GetHistogram()->SetName("TPC_ITSDphi3DGlobal3D");
572  chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Global Fit)");
573  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
574  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
575  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
576  chain->GetHistogram()->Draw("colz");
577  chain->GetHistogram()->Write();
578  //
579  chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz");
580  chain->GetHistogram()->SetName("TPC_ITSDRphi3DGlobal3D");
581  chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Global Fit)");
582  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
583  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
584  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
585  chain->GetHistogram()->Draw("colz");
586  chain->GetHistogram()->Write();
587  //
588  chain->Draw("10*(mean-deltaG):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz");
589  chain->GetHistogram()->SetName("TPC_VertexDRphi3DGlobal3D");
590  chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Global Fit)");
591  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
592  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
593  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
594  chain->GetHistogram()->Draw("colz");
595  chain->GetHistogram()->Write();
596  //
597  //
598  //
599  chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&refX==0","colz");
600  chain->GetHistogram()->SetName("TPC_VertexDphi3DLocal3D");
601  chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{#phi} - (Local Fit)");
602  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
603  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
604  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
605  chain->GetHistogram()->Draw("colz");
606  chain->GetHistogram()->Write();
607  //
608  chain->Draw("1000*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==2&&abs(refX-85)<2","colz");
609  chain->GetHistogram()->SetName("TPC_ITSDphi3DLocal3D");
610  chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{#phi} - (Local Fit)");
611  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
612  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
613  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
614  chain->GetHistogram()->Draw("colz");
615  chain->GetHistogram()->Write();
616  //
617  chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-85)<2","colz");
618  chain->GetHistogram()->SetName("TPC_ITSDRphi3DLocal3D");
619  chain->GetHistogram()->SetTitle("TPC-ITS #Delta_{r#phi} - (Local Fit)");
620  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
621  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
622  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
623  chain->GetHistogram()->Draw("colz");
624  chain->GetHistogram()->Write();
625  //
626  chain->Draw("10*(mean-deltaG-deltaL):sector:abs(theta)",cutS+"theta>0&&type==0&&abs(refX-0)<2","colz");
627  chain->GetHistogram()->SetName("TPC_VertexDRphiLocal3D");
628  chain->GetHistogram()->SetTitle("TPC-Vertex #Delta_{r#phi} - (Local Fit)");
629  chain->GetHistogram()->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
630  chain->GetHistogram()->GetXaxis()->SetTitle("sector");
631  chain->GetHistogram()->GetZaxis()->SetTitle("tan(#Theta)");
632  chain->GetHistogram()->Draw("colz");
633  chain->GetHistogram()->Write();
634 }
635 
636 
637 
638 
639 void FitAlignCombined0(){
641 
642  TTreeSRedirector *pcstream= new TTreeSRedirector("fitAlignCombined.root");
643 
644  TStatToolkit toolkit;
645  Double_t chi2=0;
646  Int_t npoints=0;
647  // Fit vectors
648  // global fits
649  TVectorD vecSec(37);
650  TVectorD paramYVGlobal;
651  TVectorD paramYITSGlobal;
652  TVectorD paramPhiVGlobal;
653  TVectorD paramPhiITSGlobal;
654  // local fits
655  TVectorD paramYVLocal;
656  TVectorD paramPhiVLocal;
657  TVectorD paramYITSLocal;
658  TVectorD paramPhiITSLocal;
659  TMatrixD covar;
660  // make Aliases
661  {for (Int_t isec=0; isec<18; isec++){ // sectors
662  tree->SetAlias(Form("sec%d",isec), Form("(abs(sector-%f)<0.5)",isec+0.5));
663  }}
664  for (Int_t isec=-1;isec<36; isec++){
665  vecSec[isec+1]=isec;
666  }
667  //
668  TString fstringGlobal="";
669  fstringGlobal+="side++";
670  fstringGlobal+="theta++";
671  fstringGlobal+="cos(phi)*(theta>0)++"; // GX - A side
672  fstringGlobal+="sin(phi)*(theta>0)++"; // GY - A side
673  fstringGlobal+="cos(phi)*(theta<0)++"; // GX - C side
674  fstringGlobal+="sin(phi)*(theta<0)++"; // GY - C side
675  fstringGlobal+="cos(phi)*(theta)++"; // theta - get rid of rotation
676  fstringGlobal+="sin(phi)*(theta)++"; // theta - get rid of rotation
677  //
678  // Local Fit
679  //
680  TString fstringLocal="";
681  {for (Int_t sector=0; sector<18; sector++){
682  fstringLocal+=Form("((sec%d)*theta>0)++",sector);
683  fstringLocal+=Form("((sec%d)*theta<0)++",sector);
684  }}
685  //
686  // Make global fits
687  //
688  TString *strFitYVGlobal = TStatToolkit::FitPlane(tree,"V.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYVGlobal,covar,-1,0, 10000000, kFALSE);
689  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
690  tree->SetAlias("fitYVGlobal",strFitYVGlobal->Data());
691  tree->Draw("V.mean:fitYVGlobal",cutS);
692  //
693  TString *strFitPhiVGlobal = TStatToolkit::FitPlane(tree,"VPhi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiVGlobal,covar,-1,0, 10000000, kFALSE);
694  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
695  tree->SetAlias("fitPhiVGlobal",strFitPhiVGlobal->Data());
696  tree->Draw("VPhi.mean:fitPhiVGlobal",cutS);
697  //
698  TString *strFitYITSGlobal = TStatToolkit::FitPlane(tree,"mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramYITSGlobal,covar,-1,0, 10000000, kFALSE);
699  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
700  tree->SetAlias("fitYITSGlobal",strFitYITSGlobal->Data());
701  tree->Draw("mean:fitYITSGlobal",cutS);
702  TString *strFitPhiITSGlobal = TStatToolkit::FitPlane(tree,"Phi.mean", fstringGlobal.Data(),cutS+"", chi2,npoints,paramPhiITSGlobal,covar,-1,0, 10000000, kFALSE);
703  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
704  tree->SetAlias("fitPhiITSGlobal",strFitPhiITSGlobal->Data());
705  tree->Draw("Phi.mean:fitPhiITSGlobal",cutS);
706  //
707  // Residuals to the global fit
708  //
709  tree->SetAlias("dyVG", "V.mean-fitYVGlobal");
710  tree->SetAlias("dphiVG","(VPhi.mean-fitPhiVGlobal)");
711  tree->SetAlias("dyITSG","mean-fitYITSGlobal");
712  tree->SetAlias("dphiITSG","(Phi.mean-fitPhiITSGlobal)");
713 
714  TString *strFitYVLocal = TStatToolkit::FitPlane(tree,"dyVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYVLocal,covar,-1,0, 10000000, kFALSE);
715  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
716  tree->SetAlias("fitYVLocal",strFitYVLocal->Data());
717  tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz");
718  //
719  TString *strFitPhiVLocal = TStatToolkit::FitPlane(tree,"dphiVG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiVLocal,covar,-1,0, 10000000, kFALSE);
720  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
721  tree->SetAlias("fitPhiVLocal",strFitPhiVLocal->Data());
722  tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz");
723  //
724  //
725  //
726  TString *strFitYITSLocal = TStatToolkit::FitPlane(tree,"dyITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramYITSLocal,covar,-1,0, 10000000, kFALSE);
727  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
728  tree->SetAlias("fitYITSLocal",strFitYITSLocal->Data());
729  tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
730  //
731  TString *strFitPhiITSLocal = TStatToolkit::FitPlane(tree,"dphiITSG", fstringLocal.Data(),cutS+"", chi2,npoints,paramPhiITSLocal,covar,-1,0, 10000000, kFALSE);
732  printf("chi2=%f\n",TMath::Sqrt(chi2/npoints));
733  tree->SetAlias("fitPhiITSLocal",strFitPhiITSLocal->Data());
734  tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
735 
736  //
737  TH1D *hisA = 0;
738  TH1D *hisC = 0;
739  TVectorD vSec(36);
740  TVectorD vecDPhi(36);
741  TVectorD vecDLY(36);
742  TVectorD vecDGX(36);
743  TVectorD vecDGY(36);
744 
745  //
746  tree->SetAlias("phiMean","(fitPhiVLocal+fitPhiITSLocal+fitPhiVGlobal+fitPhiITSGlobal)*0.5");
747  tree->SetAlias("yMeanV","((fitYITSLocal+fitYITSGlobal-85*phiMean)+(fitYVLocal+fitYVGlobal))*0.5");
748  tree->SetAlias("yMeanITS","((fitYITSLocal+fitYITSGlobal)+(fitYVLocal+fitYVGlobal+85*phiMean))*0.5");
749 
750 
751  tree->Draw("phiMean:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof");
752  hisA = (TH1D*) tree->GetHistogram()->Clone();
753  tree->Draw("phiMean:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof");
754  hisC = (TH1D*) tree->GetHistogram()->Clone();
755 
756  for (Int_t isec=0; isec<36; isec++){
757  vSec[isec]=isec;
758  if (isec<18) vecDPhi[isec]=hisA->GetBinContent(isec%18+1);
759  if (isec>=18) vecDPhi[isec]=hisC->GetBinContent(isec%18+1);
760  }
761  tree->Draw("yMeanV:int(sector)>>hhhA(18,0.18)",cutS+"abs(theta)<0.15&&theta>0","prof");
762  hisA = (TH1D*) tree->GetHistogram()->Clone();
763  tree->Draw("yMeanV:int(sector)>>hhhC(18,0.18)",cutS+"abs(theta)<0.15&&theta<0","prof");
764  hisC = (TH1D*) tree->GetHistogram()->Clone();
765  //
766  for (Int_t isec=0; isec<36; isec++){
767  if (isec<18) vecDLY[isec]=hisA->GetBinContent(isec%18+1);
768  if (isec>=18) vecDLY[isec]=hisC->GetBinContent(isec%18+1);
769  vecDGX[isec]=-vecDLY[isec]*TMath::Sin(TMath::Pi()*(isec+0.5)/9.);
770  vecDGY[isec]= vecDLY[isec]*TMath::Cos(TMath::Pi()*(isec+0.5)/9.);
771  }
772 
773 
774 
775 
776  //
777  // Store results
778  //
779  {(*pcstream)<<"fitLocal"<<
780  //global fits
781  "sec.="<<&vSec<<
782  "pYVGlobal.="<<&paramYVGlobal<<
783  "pYITSGlobal.="<<&paramYITSGlobal<<
784  "pPhiVGlobal.="<<&paramPhiVGlobal<<
785  "pPhiITSGlobal.="<<&paramPhiITSGlobal<<
786  // local fits
787  "pYVLocal.="<<&paramYVLocal<<
788  "pPhiVLocal.="<<&paramPhiVLocal<<
789  "pYITSLocal.="<<&paramYITSLocal<<
790  "pPhiITSLocal.="<<&paramPhiITSLocal<<
791  //
792  // combined
793  "vDPhi.="<<&vecDPhi<<
794  "vDLY.="<<&vecDLY<<
795  "vDGX.="<<&vecDGX<<
796  "vDGY.="<<&vecDGY<<
797  "\n";
798  }
799  paramYVGlobal.Write("paramYVGlobal");
800  paramYITSGlobal.Write("paramYITSGlobal");
801  paramPhiVGlobal.Write("paramPhiVGlobal");
802  paramPhiITSGlobal.Write("paramPhiITSGlobal");
803  // local fits
804  paramYVLocal.Write("paramYVLocal");
805  paramPhiVLocal.Write("paramPhiVLocal");
806  paramYITSLocal.Write("paramYITSLocal");
807  paramPhiITSLocal.Write("paramPhiITSLocal");
808  vecDPhi.Write("vecDPhi");
809  vecDLY.Write("vecDLY");
810  vecDGX.Write("vecDGX");
811  vecDGY.Write("vecDGY");
812 
813 
814 
815  tree->Draw("dyVG:sector:abs(theta)",cutS+"theta>0","colz");
816  tree->GetHistogram()->Write("fitYVGlobal");
817  tree->Draw("dphiVG:sector:abs(theta)",cutS+"theta>0","colz");
818  tree->GetHistogram()->Write("fitPhiVGlobal");
819  tree->Draw("dyITSG:sector:abs(theta)",cutS+"theta>0","colz");
820  tree->GetHistogram()->Write("fitYITSGlobal");
821  tree->Draw("dphiITSG:sector:abs(theta)",cutS+"theta>0","colz");
822  tree->GetHistogram()->Write("fitPhiITSGlobal");
823  //
824  tree->Draw("dyVG-fitYVLocal:sector:abs(theta)",cutS+"theta>0","colz");
825  tree->GetHistogram()->Write("fitYVLocal");
826  tree->Draw("dphiVG-fitPhiVLocal:sector:abs(theta)",cutS+"theta>0","colz");
827  tree->GetHistogram()->Write("fitPhiVLocal");
828  tree->Draw("dyITSG-fitYITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
829  tree->GetHistogram()->Write("fitYITSLocal");
830  tree->Draw("dphiITSG-fitPhiITSLocal:sector:abs(theta)",cutS+"theta>0","colz");
831  tree->GetHistogram()->Write("fitPhiITSLocal");
832  delete pcstream;
833 }
834 
835 void FitAlignCombined(){
842 
843  TFile fcalib("../mergeField0/TPCAlignObjects.root");
844  AliTPCcalibAlign * align = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
845 
846  TCut cutQ="entries>1000&&abs(theta)<0.5&&abs(snp)<0.2&&YP.entries>50&&YM.entries>50";
847  TH1F his1("hdeltaY1","hdeltaY1",100,-0.5,0.5);
848  TMatrixD vecAlign(72,1);
849  TMatrixD covAlign(72,72);
850  TMatrixD vecAlignY(72,1);
851  TMatrixD covAlignY(72,72);
852  TMatrixD vecAlignTheta(72,1);
853  TMatrixD covAlignTheta(72,72);
854  TMatrixD vecAlignZ(72,1);
855  TMatrixD covAlignZ(72,72);
856  AliTPCkalmanAlign::BookAlign1D(vecAlign,covAlign,0,0.002);
857  AliTPCkalmanAlign::BookAlign1D(vecAlignY,covAlignY,0,0.1);
858  AliTPCkalmanAlign::BookAlign1D(vecAlignTheta,covAlignTheta,0,0.002);
859  AliTPCkalmanAlign::BookAlign1D(vecAlignZ,covAlignZ,0,0.1);
860  TVectorD vecITSY(72);
861  TVectorD vecITSYPM(72);
862  TVectorD vecITSPhi(72);
863  TVectorD vecITSPhiPM(72);
864  TVectorD vecVY(72);
865  TVectorD vecVS(72);
866  TVectorD vecITSTheta(72);
867  TVectorD vecVTheta(72);
868  {for (Int_t isec0=0; isec0<36; isec0++){
869  Double_t phi0=(isec0%18+0.5)*TMath::Pi()/9.;
870  if (phi0>TMath::Pi()) phi0-=TMath::TwoPi();
871  Int_t iside0=(isec0%36<18)? 0:1;
872  TCut cutSector=Form("abs(%f-phi)<0.14",phi0);
873  TCut cutSide = (iside0==0)? "theta>0":"theta<0";
874  //
875  itsdy->Draw("mean",cutQ+cutSector+cutSide);
876  Double_t meanITSY=itsdy->GetHistogram()->GetMean();
877  vecITSY[isec0]=meanITSY;
878  vecITSY[isec0+36]=meanITSY;
879  //
880  itsdy->Draw("(YP.mean+YM.mean)*0.5",cutQ+cutSector+cutSide);
881  Double_t meanITSYPM=itsdy->GetHistogram()->GetMean();
882  vecITSYPM[isec0]=meanITSYPM;
883  vecITSYPM[isec0+36]=meanITSYPM;
884  //
885  itsdy->Draw("Phi.mean",cutQ+cutSector+cutSide);
886  Double_t meanITSPhi=itsdy->GetHistogram()->GetMean();
887  vecITSPhi[isec0]=meanITSPhi;
888  vecITSPhi[isec0+36]=meanITSPhi;
889  //
890  itsdy->Draw("(PhiP.mean+PhiM.mean)*0.5",cutQ+cutSector+cutSide);
891  Double_t meanITSPhiPM=itsdy->GetHistogram()->GetMean();
892  vecITSPhiPM[isec0]=meanITSPhiPM;
893  vecITSPhiPM[isec0+36]=meanITSPhiPM;
894  //
895  //
896  itsdy->Draw("VPhi.mean",cutQ+cutSector+cutSide);
897  Double_t meanVS=itsdy->GetHistogram()->GetMean();
898  vecVS[isec0]=meanVS;
899  vecVS[isec0+36]=meanVS;
900  //
901  itsdy->Draw("V.mean",cutQ+cutSector+cutSide);
902  Double_t meanVY=itsdy->GetHistogram()->GetMean();
903  vecVY[isec0]=meanVY;
904  vecVY[isec0+36]=meanVY;
905  //
906  itsdy->Draw("T.mean",cutQ+cutSector+cutSide);
907  Double_t meanITST=itsdy->GetHistogram()->GetMean();
908  vecITSTheta[isec0]=meanITST;
909  vecITSTheta[isec0+36]=meanITST;
910  //
911  itsdy->Draw("VT.mean",cutQ+cutSector+cutSide);
912  Double_t meanVT=itsdy->GetHistogram()->GetMean();
913  vecVTheta[isec0]=meanVT;
914  vecVTheta[isec0+36]=meanVT;
915  }
916  }
917  AliTPCROC * roc = AliTPCROC::Instance();
918  Double_t fXIROC = (roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(0,roc->GetNRows(0)-1))*0.5;
919  Double_t fXOROC = (roc->GetPadRowRadii(36,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
920  Double_t fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
921 
922  TTreeSRedirector *pcstream=new TTreeSRedirector("combAlign.root");
923 
924  {
925  for (Int_t iter=0; iter<5; iter++){
926  for (Int_t isec0=0; isec0<72; isec0++){
927  for (Int_t isec1=0; isec1<72; isec1++){
928  //if (isec0%36!=isec0%36) continue;
929  if (iter==0 && isec0%36!=isec1%36) continue;
930  if (iter==1 && isec0%18!=isec1%18) continue;
931  TH1 * his = align->GetHisto(AliTPCcalibAlign::kY,isec0,isec1);
932  TH1 * hisPhi = align->GetHisto(AliTPCcalibAlign::kPhi,isec0,isec1);
933  TH1 * hisTheta = align->GetHisto(AliTPCcalibAlign::kTheta,isec0,isec1);
934  TH1 * hisZ = align->GetHisto(AliTPCcalibAlign::kZ,isec0,isec1);
935  Int_t side0=(isec0%36)<18 ? 1:-1;
936  Int_t side1=(isec1%36)<18 ? 1:-1;
937  Double_t weightTPC= 0.002;
938  if (isec0%18==isec1%18) weightTPC=0.0005;
939  if (isec0%36==isec1%36) weightTPC=0.00005;
940  if (side0!=side1) weightTPC*=2.;
941  Double_t weightTPCT= 0.001;
942  if (isec0%18==isec1%18) weightTPC=0.0005;
943  if (isec0%36==isec1%36) weightTPC=0.00005;
944  if (TMath::Abs(isec0%18-isec1%18)==1) weightTPC = 0.0005;
945  if (TMath::Abs(isec0%36-isec1%36)==1) weightTPC = 0.0001;
946 
947  //
948  Double_t meanITS0=vecITSY[isec0];
949  Double_t meanITSPM0=vecITSYPM[isec0];
950  Double_t meanITSPhi0=vecITSPhi[isec0];
951  Double_t meanITSPhiPM0=vecITSPhiPM[isec0];
952  Double_t meanVY0=vecVY[isec0];
953  Double_t meanVPhi0=vecVS[isec0];
954  Double_t meanITSTheta0=vecITSTheta[isec0];
955  Double_t meanVTheta0=vecVTheta[isec0];
956  //
957  Double_t meanITS1=vecITSY[isec1];
958  Double_t meanITSPM1=vecITSYPM[isec1];
959  Double_t meanITSPhi1=vecITSPhi[isec1];
960  Double_t meanITSPhiPM1=vecITSPhiPM[isec1];
961  Double_t meanVY1=vecVY[isec1];
962  Double_t meanVPhi1=vecVS[isec1];
963  Double_t meanITSTheta1=vecITSTheta[isec1];
964  Double_t meanVTheta1=vecVTheta[isec1];
965  //
966  Double_t meanPhi0 = (meanITSPhi0+meanVPhi0+2*meanITS0/83.)/4.;
967  Double_t meanPhi1 = (meanITSPhi1+meanVPhi1+2*meanITS1/83.)/4.;
968  //
969  //
970  if (iter>2 &&isec0==isec1){
971  AliTPCkalmanAlign::UpdateAlign1D(-meanITS0/83, 0.001, isec0, vecAlign,covAlign);
972  AliTPCkalmanAlign::UpdateAlign1D(-meanITSPhi0, 0.001, isec0, vecAlign,covAlign);
973  AliTPCkalmanAlign::UpdateAlign1D(-meanVPhi0, 0.001, isec0, vecAlign,covAlign);
974  //
975  AliTPCkalmanAlign::UpdateAlign1D(-meanPhi0, 0.001, isec0, vecAlign,covAlign);
976  AliTPCkalmanAlign::UpdateAlign1D(-meanVTheta0, 0.001, isec0, vecAlignTheta,covAlignTheta);
977  AliTPCkalmanAlign::UpdateAlign1D(-meanITSTheta0, 0.001, isec0, vecAlignTheta,covAlignTheta);
978  }
979  if (iter>2){
980  AliTPCkalmanAlign::UpdateAlign1D(meanPhi1-meanPhi0, 0.001, isec0,isec1, vecAlign,covAlign);
981  AliTPCkalmanAlign::UpdateAlign1D(meanITSPhiPM1-meanITSPhiPM0, 0.001, isec0,isec1, vecAlign,covAlign);
982  AliTPCkalmanAlign::UpdateAlign1D((meanITSPM1-meanITSPM0)/83., 0.001, isec0,isec1, vecAlign,covAlign);
983  AliTPCkalmanAlign::UpdateAlign1D(meanVTheta1-meanVTheta0, 0.001, isec0,isec1, vecAlignTheta,covAlignTheta);
984  AliTPCkalmanAlign::UpdateAlign1D(meanITSTheta1-meanITSTheta0, 0.001, isec0,isec1, vecAlignTheta,covAlignTheta);
985  }
986  //
987  if (!his) continue;
988  if (!hisPhi) continue;
989  if (!hisTheta) continue;
990  if (his->GetEntries()<100) continue;
991  Double_t xref=fXIO;
992  if (isec0<36&&isec1<36) xref=fXIROC;
993  if (isec0>=36&&isec1>=36) xref=fXOROC;
994  Double_t meanTPC=his->GetMean();
995  Double_t meanTPCPhi=hisPhi->GetMean();
996  Double_t meanTPCTheta=hisTheta->GetMean();
997  Double_t meanTPCZ=hisZ->GetMean();
998  //
999  //
1000  Double_t kalman0= vecAlign(isec0,0);
1001  Double_t kalman1= vecAlign(isec1,0);
1002  Double_t kalmanY0= vecAlignY(isec0,0);
1003  Double_t kalmanY1= vecAlignY(isec1,0);
1004  Double_t kalmanTheta0= vecAlignTheta(isec0,0);
1005  Double_t kalmanTheta1= vecAlignTheta(isec1,0);
1006  Double_t kalmanZ0= vecAlignZ(isec0,0);
1007  Double_t kalmanZ1= vecAlignZ(isec1,0);
1008  //
1009  //
1010  AliTPCkalmanAlign::UpdateAlign1D((meanTPC)/xref, weightTPC, isec0,isec1, vecAlign,covAlign);
1011  AliTPCkalmanAlign::UpdateAlign1D(meanTPCPhi, weightTPC*10,isec0,isec1, vecAlign,covAlign);
1012  //
1013  AliTPCkalmanAlign::UpdateAlign1D(meanTPCTheta, weightTPCT, isec0,isec1, vecAlignTheta,covAlignTheta);
1014  if (side0==side1) AliTPCkalmanAlign::UpdateAlign1D(meanTPCZ, weightTPCT*100., isec0,isec1, vecAlignZ,covAlignZ);
1015  //printf("isec0\t%d\tisec1\t%d\t%f\t%f\t%f\n",isec0,isec1, meanTPC, meanITS0,meanITS1);
1016 
1017  if (iter>=0) (*pcstream)<<"align"<<
1018  "iter="<<iter<<
1019  "xref="<<xref<< // reference X
1020  "isec0="<<isec0<< // sector number
1021  "isec1="<<isec1<<
1022  "side0="<<side0<<
1023  "side1="<<side1<<
1024  //TPC
1025  "mTPC="<<meanTPC<< // delta Y / xref
1026  "mTPCPhi="<<meanTPCPhi<< // delta Phi
1027  "mTPCZ="<<meanTPCZ<< // delta Z
1028  "mTPCTheta="<<meanTPCTheta<< // delta Theta
1029  //ITS
1030  "mITS0="<<meanITS0<<
1031  "mITS1="<<meanITS1<<
1032  "mITSPhi0="<<meanITSPhi0<<
1033  "mITSPhi1="<<meanITSPhi1<<
1034  //
1035  "mITSPM0="<<meanITSPM0<<
1036  "mITSPM1="<<meanITSPM1<<
1037  "mITSPhiPM0="<<meanITSPhiPM0<<
1038  "mITSPhiPM1="<<meanITSPhiPM1<<
1039  //
1040  "mITSTheta0="<<meanITSTheta0<<
1041  "mITSTheta1="<<meanITSTheta1<<
1042  //Vertex
1043  "mVY0="<<meanVY0<<
1044  "mVY1="<<meanVY1<<
1045  "mVPhi0="<<meanVPhi0<<
1046  "mVPhi1="<<meanVPhi1<<
1047  "mVTheta0="<<meanVTheta0<<
1048  "mVTheta1="<<meanVTheta1<<
1049  // Vertex+ITS mean
1050  "mPhi0="<<meanPhi0<<
1051  "mPhi1="<<meanPhi1<<
1052  //Kalman
1053  "kY0="<<kalmanY0<<
1054  "kY1="<<kalmanY1<<
1055  "kPhi0="<<kalman0<<
1056  "kPhi1="<<kalman1<<
1057  "kTheta0="<<kalmanTheta0<<
1058  "kTheta1="<<kalmanTheta1<<
1059  "kZ0="<<kalmanZ0<<
1060  "kZ1="<<kalmanZ1<<
1061  "\n";
1062  }
1063  }
1064  }
1065  }
1066  pcstream->GetFile()->cd();
1067  vecAlign.Write("alignPhiMean");
1068  covAlign.Write("alingPhiCovar");
1069  vecAlignTheta.Write("alignThetaMean");
1070  covAlignTheta.Write("alingThetaCovar");
1071  vecAlignZ.Write("alignZMean");
1072  covAlignZ.Write("alingZCovar");
1073  delete pcstream;
1074  TFile f("combAlign.root");
1075  TTree * treeA = (TTree*)f.Get("align");
1076  treeA->SetMarkerStyle(25);
1077  treeA->SetMarkerSize(0.5);
1078 }
1079 
1080 
1081 
1082 
1083 
1084 void UpdateOCDBAlign(){
1091 
1092  AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
1093  TClonesArray * array = (TClonesArray*)entry->GetObject();
1094  Int_t entries = array->GetEntries();
1095  TClonesArray * arrayNew=(TClonesArray*)array->Clone();
1096  TFile f("combAlign.root");
1097  TMatrixD *matPhi=(TMatrixD*)f.Get("alignPhiMean");
1098  //
1099  //
1100  { for (Int_t i=0;i<entries; i++){
1101  //
1102  //
1103  AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
1104  AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i);
1105  Int_t imod;
1106  AliGeomManager::ELayerID ilayer;
1107  alignP->GetVolUID(ilayer, imod);
1108  if (ilayer==AliGeomManager::kInvalidLayer) continue;
1109  Int_t sector=imod;
1110  if (ilayer==AliGeomManager::kTPC2) sector+=36;
1111  Double_t transOld[3], rotOld[3];
1112  alignP->GetTranslation(transOld); // in cm
1113  alignP->GetAngles(rotOld); // in degrees
1114  printf("%d\t%d\t%d\t",ilayer, imod,sector);
1115  printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*matPhi)(sector,0)));
1116  alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]+(*matPhi)(sector,0)*TMath::RadToDeg());
1117  alignPNew->Print("kokot");
1118  alignP->Print("kokot");
1119  }
1120  }
1121 
1122 
1123 
1124  TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate";
1125  TString userName=gSystem->GetFromPipe("echo $USER");
1126  TString date=gSystem->GetFromPipe("date");
1127 
1128  AliCDBMetaData *metaData= new AliCDBMetaData();
1129  metaData->SetObjectClassName("TObjArray");
1130  metaData->SetResponsible("Marian Ivanov");
1131  metaData->SetBeamPeriod(1);
1132  metaData->SetAliRootVersion("05-26-02"); //root version
1133  metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
1134 
1135  AliCDBId* id1=NULL;
1136  id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity());
1137  AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1138  gStorage->Put(arrayNew, (*id1), metaData);
1139 
1140  TFile falign("falign.root","recreate");
1141  arrayNew->Write("new");
1142  array->Write("old");
1143  falign.Close();
1144 }
1145 
1146 
1147 
1148 void UpdateOCDBAlign0(){
1155 
1156  AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
1157  TClonesArray * array = (TClonesArray*)entry->GetObject();
1158  Int_t entries = array->GetEntries();
1159  TClonesArray * arrayNew=(TClonesArray*)array->Clone();
1160  TFile f("fitAlignCombined.root");
1161  TVectorD *vecDPhi = (TVectorD*) f.Get("vecDPhi");
1162  TVectorD *vecDGX = (TVectorD*) f.Get("vecDGX");
1163  TVectorD *vecDGY = (TVectorD*) f.Get("vecDGY");
1164 
1165  //
1166  //
1167  { for (Int_t i=0;i<entries; i++){
1168  //
1169  //
1170  AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
1171  AliAlignObjParams *alignPNew = (AliAlignObjParams*)arrayNew->UncheckedAt(i);
1172  Int_t imod;
1173  AliGeomManager::ELayerID ilayer;
1174  alignP->GetVolUID(ilayer, imod);
1175  if (ilayer==AliGeomManager::kInvalidLayer) continue;
1176  Int_t sector=imod;
1177  if (ilayer==AliGeomManager::kTPC2) sector+=36;
1178  Double_t transOld[3], rotOld[3];
1179  alignP->GetTranslation(transOld); // in cm
1180  alignP->GetAngles(rotOld); // in degrees
1181  printf("%d\t%d\t%d\t",ilayer, imod,sector);
1182  printf("%f mrad \t%f mrad\n",1000*rotOld[2]*TMath::DegToRad(), 1000*((*vecDPhi)[sector%36]));
1183  alignPNew->SetRotation(rotOld[0],rotOld[1], rotOld[2]-(*vecDPhi)[sector%36]*TMath::RadToDeg());
1184  alignPNew->SetTranslation(transOld[0]-(*vecDGX)[sector%36],transOld[1]-(*vecDGY)[sector%36], transOld[2]);
1185  alignPNew->Print("kokot");
1186  alignP->Print("kokot");
1187  }
1188  }
1189 
1190 
1191 
1192  TString ocdbStorage="local://"+gSystem->GetFromPipe("pwd")+"/OCDBUpdate";
1193  TString userName=gSystem->GetFromPipe("echo $USER");
1194  TString date=gSystem->GetFromPipe("date");
1195 
1196  AliCDBMetaData *metaData= new AliCDBMetaData();
1197  metaData->SetObjectClassName("TObjArray");
1198  metaData->SetResponsible("Marian Ivanov");
1199  metaData->SetBeamPeriod(1);
1200  metaData->SetAliRootVersion("05-26-02"); //root version
1201  metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
1202 
1203  AliCDBId* id1=NULL;
1204  id1=new AliCDBId("TPC/Align/Data", 0, AliCDBRunRange::Infinity());
1205  AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(ocdbStorage);
1206  gStorage->Put(arrayNew, (*id1), metaData);
1207 
1208  TFile falign("falign.root","recreate");
1209  arrayNew->Write("new");
1210  array->Write("old");
1211  falign.Close();
1212 }
AliTPCCalibGlobalMisalignment::AddAlign
void AddAlign(const AliTPCCalibGlobalMisalignment &add)
Definition: AliTPCCalibGlobalMisalignment.cxx:397
AliTPCCalibGlobalMisalignment::CreateMeanAlign
static AliTPCCalibGlobalMisalignment * CreateMeanAlign(const AliTPCCalibGlobalMisalignment *alignIn)
Definition: AliTPCCalibGlobalMisalignment.cxx:499
printf
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
vdyM
TTree * vdyM
Definition: FitAlignCombined.C:106
FitAlignCombined0
void FitAlignCombined0()
Definition: FitAlignCombined.C:639
TObjArray
#define TObjArray
Definition: AliMUONReconstructor.h:19
align
AliTPCcalibAlign align
Definition: CalibAlign.C:43
alignRot0
AliTPCCalibGlobalMisalignment * alignRot0
Definition: RegisterCorrection.C:146
alignTrans2
AliTPCCalibGlobalMisalignment * alignTrans2
Definition: RegisterCorrection.C:151
RegisterAlignFunction
void RegisterAlignFunction()
Definition: FitAlignCombined.C:114
cutS
TCut cutS
Definition: FitAlignCombined.C:111
alignRot2
AliTPCCalibGlobalMisalignment * alignRot2
Definition: RegisterCorrection.C:148
AliTPCROC::GetNRows
UInt_t GetNRows(UInt_t sector) const
Definition: AliTPCROC.h:28
f
TFile f("CalibObjects.root")
itsdyP
TTree * itsdyP
Definition: FitAlignCombined.C:93
vdz
TTree * vdz
Definition: FitAlignCombined.C:108
AliTPCROC.h
trddy
TTree * trddy
Definition: FitAlignCombined.C:102
UpdateOCDBAlign0
void UpdateOCDBAlign0()
Definition: FitAlignCombined.C:1148
vds
TTree * vds
Definition: FitAlignCombined.C:107
pcstream
TTreeSRedirector * pcstream
Definition: FitAlignCombined.C:87
LoadTrees
void LoadTrees()
Definition: FitAlignCombined.C:311
AliTPCcalibDB.h
AliTPCParam.h
AliTPCCalPad.h
itsdt
TTree * itsdt
Definition: FitAlignCombined.C:100
npoints
npoints
Definition: driftITSTPC.C:85
combAlignOCDBOld
AliTPCCalibGlobalMisalignment * combAlignOCDBOld
Definition: FitAlignCombined.C:81
MakeAlignFunctionSector
AliTPCCalibGlobalMisalignment * MakeAlignFunctionSector(TVectorD paramYLocal)
Definition: FitAlignCombined.C:220
AliTPCCalibGlobalMisalignment::Print
virtual void Print(Option_t *option="") const
Definition: AliTPCCalibGlobalMisalignment.cxx:372
array
TObjArray * array
Definition: AnalyzeLaser.C:12
AliTPCCalibGlobalMisalignment::SetAlignSectors
void SetAlignSectors(const TObjArray *arraySector)
Definition: AliTPCCalibGlobalMisalignment.cxx:234
UpdateOCDBAlign
void UpdateOCDBAlign()
Definition: FitAlignCombined.C:1084
alignTrans1
AliTPCCalibGlobalMisalignment * alignTrans1
Definition: RegisterCorrection.C:150
AliTPCCalROC.h
AliTPCCalibGlobalMisalignment.h
AliTPCPreprocessorOnline.h
itsdz
TTree * itsdz
Definition: FitAlignCombined.C:99
AliTPCExBEffectiveSector.h
chi2
Double_t chi2
Definition: AnalyzeLaser.C:7
vdy
TTree * vdy
Definition: FitAlignCombined.C:104
FitAlignCombinedCorr
void FitAlignCombinedCorr()
Definition: FitAlignCombined.C:382
combAlignLocal
AliTPCCalibGlobalMisalignment * combAlignLocal
Definition: FitAlignCombined.C:84
AliTPCComposedCorrection.h
proc
AliTPCROC * proc
Definition: FitAlignCombined.C:80
AliTPCROC
Geometry class for a single ROC.
Definition: AliTPCROC.h:14
alignRot1
AliTPCCalibGlobalMisalignment * alignRot1
Definition: RegisterCorrection.C:147
tofdy
TTree * tofdy
Definition: FitAlignCombined.C:101
vdyP
TTree * vdyP
Definition: FitAlignCombined.C:105
chainZ
TChain * chainZ
Definition: FitAlignCombined.C:90
MakeAlignFunctionGlobal
AliTPCCalibGlobalMisalignment * MakeAlignFunctionGlobal(TVectorD paramYGlobal)
Definition: FitAlignCombined.C:171
toolkit
TStatToolkit toolkit
Definition: gainCalib.C:18
tree
TTree * tree
Definition: FitAlignCombined.C:91
DrawFitQA
void DrawFitQA()
Definition: FitAlignCombined.C:526
AliTPCCorrection::AddVisualCorrection
static void AddVisualCorrection(AliTPCCorrection *corr, Int_t position)
Definition: AliTPCCorrection.cxx:3115
itsdy
TTree * itsdy
Definition: FitAlignCombined.C:92
vdt
TTree * vdt
Definition: FitAlignCombined.C:109
alignTrans0
AliTPCCalibGlobalMisalignment * alignTrans0
Definition: RegisterCorrection.C:149
AliTPCCalibGlobalMisalignment::SetAlignGlobal
void SetAlignGlobal(const TGeoMatrix *matrixGlobal)
Definition: AliTPCCalibGlobalMisalignment.cxx:216
itsdphiM
TTree * itsdphiM
Definition: FitAlignCombined.C:97
itsdp
TTree * itsdp
Definition: FitAlignCombined.C:95
FitAlignCombined
void FitAlignCombined()
Definition: FitAlignCombined.C:835
combAlignGlobal
AliTPCCalibGlobalMisalignment * combAlignGlobal
Definition: FitAlignCombined.C:83
combAlignOCDBNew
AliTPCCalibGlobalMisalignment * combAlignOCDBNew
Definition: FitAlignCombined.C:82
AliTPCROC::Instance
static AliTPCROC * Instance()
Definition: AliTPCROC.cxx:34
AliTPCROC::GetPadRowRadii
Float_t GetPadRowRadii(UInt_t isec, UInt_t irow) const
Definition: AliTPCROC.h:56
AliTPCCalibGlobalMisalignment::CreateOCDBAlign
static AliTPCCalibGlobalMisalignment * CreateOCDBAlign()
Definition: AliTPCCalibGlobalMisalignment.cxx:457
itsdyM
TTree * itsdyM
Definition: FitAlignCombined.C:94
chain
TChain * chain
Definition: FitAlignCombined.C:89
itsdphiP
TTree * itsdphiP
Definition: FitAlignCombined.C:96
AliMp::kY
direction in y
Definition: AliMpDirection.h:21