Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / TPC / CalibMacros / FitRodShift.C
1 /// \file FitRodShift.C
2 ///
3 /// \author marian.ivanov@cern.ch, Stefan.Rossegger@cern.ch
4 ///
5 /// Macro to fit  alignment/E field distortion maps
6 /// As a input output cluster distortion maps (produced by  AliTPCcalibAling class)
7 /// are used. Cluster distrotion maps granularity (180 *phi x44 *theta x 53 R)
8 ///
9 /// In total 440 parameters to fit using global fit:
10 ///
11 /// 1. Rotation and translation for each sector (72x2)
12 /// 2. Rotation and translation for each quadrant of OROC (36x4x2)
13 /// 3. Rod/strip shifts in IFC and OFC (18 sectors x 2 sides x 2 ) 
14 /// 4. Rotated clips in IFC and OFC 
15 ///
16 ///
17 /// Input file mean.root with distortion tree expected to be in directory:
18 /// ../mergeField0/clusterDY.root
19 /// The ouput file fitRodShift.root contains:
20 /// 1. Resulting (residual) AliTPCCalibMisalignment  and AliTPCFCVoltError3D classes
21 /// 2. All important temporary results are stored in the workspace (TTreeSRedirector)  associated
22 ///    with the file - fitrodShift
23 ///    2.a  Fit parameters with errors
24 ///    2.b  QA default graphs
25 ///    2.c  QA defualt histograms
26 ///
27 ///
28 ///
29 /// Functions:
30 /// 1. LoadModels()                   - load models to fit - Rod shift (10,20)
31 /// 2. RegisterAlignFunction()        - register align functions (10-16)
32 /// 3. LoadTrees                      - load trees and make aliases
33 /// 4. PrintFit                       - helper function to print substring of the fit 
34 /// 5. DeltaLookup                    - function to calulate the distortion for given distortion cunction
35 ///                                   -
36 /// 6. MakeAliases                    - Make tree aliases shortcuts for the fitting function
37 ///
38 /// 7. MakeAlignCorrection            - Crete alignment entry to be stored in the OCDB
39 ///                                   - Dump fit parameters and errors
40 /// 8. MakeQuadrantCorrection         - 
41 ///                                   - Dump fit parameters and errors
42 ///
43 /// 9. FitRodShifts                   - main minimization function
44 ///
45 /// ~~~
46 /// .x ~/rootlogon.C
47 /// .x ~/NimStyle.C
48 /// .L $ALICE_ROOT/TPC/CalibMacros/FitRodShift.C+
49 /// FitRodShift(kFALSE);
50 /// ~~~
51
52 #if !defined(__CINT__) || defined(__MAKECINT__)
53 #include "TFile.h"
54 #include "TTreeStream.h"
55 #include "TMath.h" 
56 #include "TGraph.h" 
57 #include "TRandom.h"
58 #include "TTree.h"
59 #include "TF1.h"
60 #include "TH2F.h"
61 #include "TH3F.h"
62 #include "TAxis.h"
63 #include "TPad.h"
64 #include "TCanvas.h"
65 #include "TStyle.h"
66 #include "AliTPCParamSR.h"
67 #include "TDatabasePDG.h"
68 #include "AliTPCFCVoltError3D.h"
69 #include "AliTPCROCVoltError3D.h"
70 #include "AliTPCComposedCorrection.h"
71 #include "AliTPCBoundaryVoltError.h"
72 #include "AliCDBEntry.h"
73 #include "AliTPCROC.h"
74 #include <TStatToolkit.h>
75 #include "TCut.h"
76 #include "TGraphErrors.h"
77 #include "AliTPCCalibGlobalMisalignment.h"
78 #endif
79
80 TTreeSRedirector *pcWorkspace= 0;    // Workspace to store the fit parameters and QA histograms
81 TTree * treeDY= 0;                       //distortion tree - Dy
82 TTree * treeDZ= 0;                       //distortion tree - Dz
83 // fit models ...
84 AliTPCFCVoltError3D *rotOFC  =0;          // fit models
85 AliTPCFCVoltError3D *rodOFC1 =0;
86 AliTPCFCVoltError3D *rodOFC2 =0;
87 AliTPCFCVoltError3D *rotIFC  =0;
88 AliTPCFCVoltError3D *rodIFC1 =0;
89 AliTPCFCVoltError3D *rodIFC2 =0;
90 AliTPCFCVoltError3D *rodAll  =0;        //all lookup initialized
91 //
92 AliTPCComposedCorrection *corFit0 = 0;  // composed correction
93 void FitFunctionQA();
94 void MakeQA();
95 void DrawAlignParam();
96
97 void PrintFit(TString fitString, TString filter){
98   /// helper function to print substring of the fit
99   /// TString fitString =*strFitGA;
100   /// TString filter="rot";
101
102   TObjArray *arr = fitString.Tokenize("++");
103   Int_t entries=arr->GetEntries();
104   for (Int_t i=0; i<entries; i++){
105     TString aaa = arr->At(i)->GetName();
106     if (aaa.Contains(filter)==1){
107       printf("%d\t%s\n",i,aaa.Data());
108     }    
109   }
110 }
111
112
113 void LoadModels(){
114   /// Load the models from the file
115   /// Or create it
116
117   Int_t volt = 1;
118   TFile f("OCDB-RodShifts.root");
119   AliCDBEntry *entry = (AliCDBEntry*) f.Get("AliCDBEntry");
120   if (entry) { // load file    
121     AliTPCComposedCorrection *cor = (AliTPCComposedCorrection*) entry->GetObject();    
122     cor->Print();
123     TCollection *iter = cor->GetCorrections();    
124     rotOFC = (AliTPCFCVoltError3D*)iter->FindObject("rotOFC");
125     rodOFC1 = (AliTPCFCVoltError3D*)iter->FindObject("rodOFC1");
126     rodOFC2 = (AliTPCFCVoltError3D*)iter->FindObject("rodOFC2");
127     rotIFC = (AliTPCFCVoltError3D*)iter->FindObject("rotIFC");
128     rodIFC1 = (AliTPCFCVoltError3D*)iter->FindObject("rodIFC1");
129     rodIFC2 = (AliTPCFCVoltError3D*)iter->FindObject("rodIFC2");
130     rodAll = (AliTPCFCVoltError3D*)iter->FindObject("rodAll");
131     // rotOFC->CreateHistoDZinXY(1,500,500)->Draw("surf2");
132   } else {    
133     // OFC 
134     rotOFC = new AliTPCFCVoltError3D();
135     rotOFC->SetOmegaTauT1T2(0,1,1);
136     rotOFC->SetRotatedClipVoltA(1,volt);
137     rotOFC->SetRotatedClipVoltC(1,volt);
138     //
139     rodOFC1 = new AliTPCFCVoltError3D();
140     rodOFC1->SetOmegaTauT1T2(0,1,1);
141     rodOFC1->SetRodVoltShiftA(18,volt);
142     rodOFC1->SetRodVoltShiftC(18,volt);
143     //
144     rodOFC2 = new AliTPCFCVoltError3D();
145     rodOFC2->SetOmegaTauT1T2(0,1,1);
146     rodOFC2->SetCopperRodShiftA(18,volt);
147     rodOFC2->SetCopperRodShiftC(18,volt);    
148     // IFC     
149     rotIFC = new AliTPCFCVoltError3D();
150     rotIFC->SetOmegaTauT1T2(0,1,1);
151     rotIFC->SetRotatedClipVoltA(0,volt);
152     rotIFC->SetRotatedClipVoltC(0,volt);
153     //
154     rodIFC1 = new AliTPCFCVoltError3D();
155     rodIFC1->SetOmegaTauT1T2(0,1,1);
156     rodIFC1->SetRodVoltShiftA(0,volt);
157     rodIFC1->SetRodVoltShiftC(0,volt);
158     //
159     rodIFC2 = new AliTPCFCVoltError3D();
160     rodIFC2->SetOmegaTauT1T2(0,1,1);
161     rodIFC2->SetCopperRodShiftA(0,volt);
162     rodIFC2->SetCopperRodShiftC(0,volt);
163     // dummy object with all inits
164     rodAll = new AliTPCFCVoltError3D();
165     Double_t volt0=0.0000000001;
166     rodAll->SetRotatedClipVoltA(1,volt0);
167     rodAll->SetRotatedClipVoltC(1,volt0);
168     rodAll->SetRodVoltShiftA(18,volt0);
169     rodAll->SetRodVoltShiftC(18,volt0);
170     rodAll->SetCopperRodShiftA(18,volt0);
171     rodAll->SetCopperRodShiftC(18,volt0);    
172     rodAll->SetRotatedClipVoltA(0,volt0);
173     rodAll->SetRotatedClipVoltC(0,volt0);
174     rodAll->SetRodVoltShiftA(0,volt0);
175     rodAll->SetRodVoltShiftC(0,volt0);
176     rodAll->SetCopperRodShiftA(0,volt0);
177     rodAll->SetCopperRodShiftC(0,volt0);
178     rodAll->SetOmegaTauT1T2(0,1,1);
179     //
180     //
181     // Initialization of the lookup tables
182     //
183     printf(" ------- OFC rotated clip:\n"); rotOFC->InitFCVoltError3D();
184     printf(" ------- OFC rod & strip:\n");  rodOFC1->InitFCVoltError3D();
185     printf(" ------- OFC copper rod:\n");   rodOFC2->InitFCVoltError3D();
186     printf(" ------- IFC rotated clip:\n"); rotIFC->InitFCVoltError3D();
187     printf(" ------- IFC rod & strip:\n");  rodIFC1->InitFCVoltError3D();
188     printf(" ------- IFC copper rod:\n");   rodIFC2->InitFCVoltError3D();
189     printf(" ------- Dummy all:\n");        rodAll->InitFCVoltError3D();
190     // give names
191     rotOFC->SetName("rotOFC");rotOFC->SetTitle("rotOFC");
192     rodOFC1->SetName("rodOFC1");rodOFC1->SetTitle("rodOFC1");
193     rodOFC2->SetName("rodOFC2");rodOFC2->SetTitle("rodOFC2");
194     rotIFC->SetName("rotIFC");rotIFC->SetTitle("rotIFC");
195     rodIFC1->SetName("rodIFC1");rodIFC1->SetTitle("rodIFC1");
196     rodIFC2->SetName("rodIFC2");rodIFC2->SetTitle("rodIFC2");
197     rodAll->SetName("rodAll");rodAll->SetTitle("rodAll");
198     //
199     // save in file
200     AliTPCComposedCorrection *corFit0 = new AliTPCComposedCorrection();
201     TObjArray *cs = new TObjArray();
202     cs->Add(rotIFC); cs->Add(rotOFC);
203     cs->Add(rodIFC1); cs->Add(rodOFC1);
204     cs->Add(rodIFC2); cs->Add(rodOFC2);
205     cs->Add(rodAll);
206     corFit0->SetCorrections(cs);
207     corFit0->SetOmegaTauT1T2(0,1.,1.);
208     corFit0->Print();    
209     corFit0->StoreInOCDB(0,1000000,"testForFit");    
210   }
211   //
212   AliTPCCorrection::AddVisualCorrection(rotOFC,0); 
213   AliTPCCorrection::AddVisualCorrection(rodOFC1,1); 
214   AliTPCCorrection::AddVisualCorrection(rodOFC2,2); 
215   AliTPCCorrection::AddVisualCorrection(rotIFC,3); 
216   AliTPCCorrection::AddVisualCorrection(rodIFC1,4); 
217   AliTPCCorrection::AddVisualCorrection(rodIFC2,5); 
218   AliTPCCorrection::AddVisualCorrection(rodAll,6); 
219 }
220
221
222 void RegisterAlignFunction(){
223   /// Register primitive alignment components.
224   /// Linear conbination of primitev forulas used for fit
225   /// The nominal delta 1 mm in shift and 1 mrad in rotation
226   /// Primitive formulas registeren in AliTPCCoreection::AddvisualCorrection
227   /// 10 - deltaX
228   /// 11 - deltaY
229   /// 12 - deltaZ
230   /// 13 - rot0 (phi)
231   /// 14 - rot1 (theta)
232   /// 15 - rot2
233
234   TGeoHMatrix matrixX;
235   TGeoHMatrix matrixY;
236   TGeoHMatrix matrixZ;
237   TGeoRotation rot0;
238   TGeoRotation rot1;
239   TGeoRotation rot2;  //transformation matrices
240   TGeoRotation rot90;  //transformation matrices
241   matrixX.SetDx(0.1); matrixY.SetDy(0.1); matrixZ.SetDz(0.1); //1 mm translation
242   rot0.SetAngles(0.001*TMath::RadToDeg(),0,0);
243   rot1.SetAngles(0,0.001*TMath::RadToDeg(),0);
244   rot2.SetAngles(0,0,0.001*TMath::RadToDeg());
245   //how to get rot02 ?
246   rot90.SetAngles(0,90,0);
247   rot2.MultiplyBy(&rot90,kTRUE);
248   rot90.SetAngles(0,-90,0);
249   rot2.MultiplyBy(&rot90,kFALSE);
250   AliTPCCalibGlobalMisalignment *alignRot0  =new  AliTPCCalibGlobalMisalignment;
251   alignRot0->SetAlignGlobal(&rot0);
252   AliTPCCalibGlobalMisalignment *alignRot1=new  AliTPCCalibGlobalMisalignment;
253   alignRot1->SetAlignGlobal(&rot1);
254   AliTPCCalibGlobalMisalignment *alignRot2=new  AliTPCCalibGlobalMisalignment;
255   alignRot2->SetAlignGlobal(&rot2);
256   //
257   AliTPCCalibGlobalMisalignment *alignTrans0  =new  AliTPCCalibGlobalMisalignment;
258   alignTrans0->SetAlignGlobal(&matrixX);
259   AliTPCCalibGlobalMisalignment *alignTrans1=new  AliTPCCalibGlobalMisalignment;
260   alignTrans1->SetAlignGlobal(&matrixY);
261   AliTPCCalibGlobalMisalignment *alignTrans2=new  AliTPCCalibGlobalMisalignment;
262   alignTrans2->SetAlignGlobal(&matrixZ);
263   AliTPCCorrection::AddVisualCorrection(alignTrans0  ,10);
264   AliTPCCorrection::AddVisualCorrection(alignTrans1  ,11);
265   AliTPCCorrection::AddVisualCorrection(alignTrans2  ,12);
266   AliTPCCorrection::AddVisualCorrection(alignRot0    ,13);
267   AliTPCCorrection::AddVisualCorrection(alignRot1    ,14);
268   AliTPCCorrection::AddVisualCorrection(alignRot2    ,15);
269   TObjArray arrAlign(6);
270   arrAlign.AddAt(alignTrans0->Clone(),0);
271   arrAlign.AddAt(alignTrans1->Clone(),1);
272   arrAlign.AddAt(alignTrans2->Clone(),2);
273   arrAlign.AddAt(alignRot0->Clone(),3);
274   arrAlign.AddAt(alignRot1->Clone(),4);
275   arrAlign.AddAt(alignRot2->Clone(),5);
276   //combAlign.SetCorrections((TObjArray*)arrAlign.Clone());
277 }
278
279
280 void LoadTrees(){
281   /// 1. Load trees
282   /// 2. Make standard cuts - filter the tree
283   /// 3. makeStandard aliases
284
285   TFile *fy = new TFile("clusterDY.root");
286   TFile *fz = new TFile("clusterDZ.root");
287   treeDY= (TTree*)fy->Get("delta");
288   treeDZ= (TTree*)fz->Get("delta");
289   TCut cutAll = "entries>3000&&abs(kZ)<1&&localX>80&&localX<246&&abs(sector-int(sector)-0.5)<0.41&&abs(localX-134)>2&&rmsG<0.5&&abs(localX-189)>3";
290
291   treeDY->Draw(">>outListY",cutAll,"entryList");
292   TEntryList *elistOutY = (TEntryList*)gDirectory->Get("outListY");
293   treeDY->SetEntryList(elistOutY);
294   treeDZ->Draw(">>outListZ",cutAll,"entryList");
295   TEntryList *elistOutZ = (TEntryList*)gDirectory->Get("outListZ");
296   treeDZ->SetEntryList(elistOutZ);
297   //
298   // Make aliases
299   //
300   treeDY->SetAlias("dsec","(sector-int(sector)-0.5)");
301   treeDY->SetAlias("dsec0","(sector-int(sector))");
302   treeDY->SetAlias("signy","(-1.+2*(sector-int(sector)>0.5))");
303   treeDY->SetAlias("dx","(localX-165.5)");   // distance X to ref. plane
304   treeDY->SetAlias("dxm","0.001*(localX-165.5)");
305   treeDY->SetAlias("rx","(localX/166.5)");   //ratio distrnace to reference plane
306   treeDY->SetAlias("dq1","(((q1==0)*(-rx)+q1*(1-rx))*signy)");
307   //
308   {for (Int_t isec=0; isec<18; isec++){ // sectors
309       treeDY->SetAlias(Form("sec%d",isec), Form("(abs(sector-%3.1lf)<0.5)",isec+0.5));
310       treeDZ->SetAlias(Form("sec%d",isec), Form("(abs(sector-%3.1lf)<0.5)",isec+0.5));
311     }}
312   treeDY->SetAlias("gy","localX*sin(pi*sector/9)");
313   treeDY->SetAlias("gx","localX*cos(pi*sector/9)");
314   treeDY->SetAlias("side","(-1.+2.*(kZ>0))");  
315   treeDY->SetAlias("drphi","mean");
316   treeDY->SetMarkerStyle(25);
317 }
318
319 Double_t DeltaLookup(Double_t sector, Double_t localX, Double_t kZ, Double_t xref, Int_t value, Int_t corr){
320   /// Distortion maps are calculated relative to the reference X plane
321   /// The same procedure applied for given correction
322   /// fd(sector, localX, kZ) ==>  fd(sector, localX, kZ)-fd(sector, xref, kZ)*localX/xref
323
324   Double_t distortion    = AliTPCCorrection::GetCorrSector(sector,localX,kZ,value,corr);
325   Double_t distortionRef = AliTPCCorrection::GetCorrSector(sector,xref,kZ,value,corr)*localX/xref;
326   return distortion-distortionRef;
327 }
328
329 void MakeAliases(){
330   /// make alias names
331
332   AliTPCROC * roc = AliTPCROC::Instance();
333   Double_t xref  = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
334   treeDY->SetAlias("iroc","(localX<134)");  // IROC
335   treeDY->SetAlias("oroc","(localX>134)");  // OROC
336   treeDY->SetAlias("q1","abs(localX-161)<28"); // OROC first haf
337   treeDY->SetAlias("q2","(localX>189)");       // OROC second half
338   //
339   treeDY->SetAlias("dIFC","abs(localX-83)");    // distance to IFC
340   treeDY->SetAlias("dOFC","abs(localX-250)");   // distance to OFC
341   treeDY->SetAlias("errY","(1.+1/(1+(dIFC/2.))+1/(1+dOFC/2.))");
342
343   //
344   // rotated clip - IFC --------------------
345   treeDY->SetAlias("rotClipI","DeltaLookup(sector,localX,kZ,165.6,1,3+0)");
346   // rotated clip - OFC --------------------
347   treeDY->SetAlias("rotClipO","DeltaLookup(sector,localX,kZ,165.6,1,0+0)");
348   for (Int_t isec=0;isec<18; isec++){    
349     //
350     // alignment IROC - shift cm - rotation mrad
351     treeDY->SetAlias(Form("dyI_%d",isec),Form("(iroc*sec%d)",isec));
352     treeDY->SetAlias(Form("drotI_%d",isec),Form("(iroc*sec%d*(localX-%3.1lf)/1000.)",isec,xref));
353     // alignment OROC
354     treeDY->SetAlias(Form("dyO_%d",isec),Form("(sec%d*oroc-sec%d*localX/165.6)",isec,isec)); // fix local Y shift - do not use it
355     treeDY->SetAlias(Form("drotO_%d",isec),Form("(sec%d*oroc*(localX-%3.1lf)/1000.)",isec,xref));
356     //
357     // quadrant shift OROC 
358     treeDY->SetAlias(Form("dqO0_%d",isec),Form("(sec%d*sign(dsec)*(q1-localX/165.6))",isec));
359     // quadrant delta rotation OFC inner
360     treeDY->SetAlias(Form("dqOR0_%d",isec),Form("(sec%d*sign(dsec)*q1*(localX-165.6)/1000.)",isec));
361     //
362     treeDY->SetAlias(Form("dqO1_%d",isec),Form("(q2*sec%d*sign(dsec))",isec));  // delta
363     treeDY->SetAlias(Form("dqOR1_%d",isec),Form("(q2*sec%d*sign(dsec)*(localX-165.6)/1000.)",isec));
364     treeDY->SetAlias(Form("dqO2_%d",isec),Form("(q2*sec%d)",isec)); //common
365     treeDY->SetAlias(Form("dqOR2_%d",isec),Form("(q2*sec%d*(localX-165.6)/1000.)",isec));
366     //
367     //
368     // shifted rod -  OFC
369     treeDY->SetAlias(Form("rodStripO_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,1+0)",isec));
370     // Shifted cooper rod OFC 
371     treeDY->SetAlias(Form("coppRodO_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,2+0)",isec)); 
372     // shifted rod - IFC 
373     treeDY->SetAlias(Form("rodStripI_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,4+0)",isec)); 
374     // Shifted cooper rod IFC 
375     treeDY->SetAlias(Form("coppRodI_%d",isec),Form("DeltaLookup(sector-%d,localX,kZ,165.6,1,5+0)",isec)); 
376   }
377   //
378   // fitted correction - Using the saved coposed corrections
379   //
380   treeDY->SetAlias("dAlign","DeltaLookup(sector,localX,kZ,165.6,1,1000+0)");
381   treeDY->SetAlias("dQuadrant","DeltaLookup(sector,localX,kZ,165.6,1,1100+0)");
382   treeDY->SetAlias("dField","DeltaLookup(sector,localX,kZ,165.6,1,100+0)");
383   treeDY->SetAlias("dAll","(dAlign+dQuadrant+dField)");
384 }
385
386 AliTPCCalibGlobalMisalignment *MakeAlignCorrection(TVectorD paramA, TVectorD paramC,  TMatrixD covar, Double_t chi2){
387   /// Make a global alignmnet
388   /// Take a fit parameters and make a combined correction:
389   /// Only delta local Y and delta phi are fitted - not sensitivity for other parameters
390   /// GX and GY shift extracted per side.
391   ///
392   /// Algorithm:
393   ///   1. Loop over sectors
394   ///   2. Make combined AliTPCCalibGlobalMisalignment
395
396   AliTPCCalibGlobalMisalignment *alignLocal  =new  AliTPCCalibGlobalMisalignment;  
397   Int_t offset=3;
398   AliTPCROC * roc = AliTPCROC::Instance();
399   Double_t xref  = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
400   // 
401   pcWorkspace->GetFile()->cd();
402   TObjArray * array = new TObjArray(72);
403   
404   for (Int_t side=-1; side<2; side+=2){
405     TVectorD &param = (side==1) ? paramA : paramC;
406     TGeoHMatrix matrixGX;
407     TGeoHMatrix matrixGY;
408     matrixGX.SetDx(0.1*param[1]); 
409     matrixGY.SetDy(0.1*param[2]); 
410     //  
411     for (Int_t isec=0; isec<18; isec++){ 
412       TGeoHMatrix matrixOROC;
413       TGeoHMatrix matrixIROC; 
414       TGeoRotation rotIROC;
415       TGeoRotation rotOROC;
416       Double_t phi= (Double_t(isec)+0.5)*TMath::Pi()/9.;
417       //
418       Double_t drotIROC = -param[offset+isec+18]*0.001;
419       Double_t drotOROC = -param[offset+isec+36]*0.001;
420       Double_t dlyIROC  = param[offset+isec];
421       Double_t dlyIROC0 = param[offset+isec]+drotIROC*xref;     // shift at x=0
422       Double_t dlyOROC  = 0;
423       Double_t dlyOROC0 = 0+drotOROC*xref;  // shift at x=0
424       //
425       Double_t dgxIROC = TMath::Cos(phi)*0+TMath::Sin(phi)*dlyIROC0;
426       Double_t dgyIROC = TMath::Sin(phi)*0-TMath::Cos(phi)*dlyIROC0;
427       Double_t dgxOROC = TMath::Cos(phi)*0+TMath::Sin(phi)*dlyOROC0;
428       Double_t dgyOROC = TMath::Sin(phi)*0-TMath::Cos(phi)*dlyOROC0;    
429       Double_t errYIROC=TMath::Sqrt(covar(offset+isec, offset+isec)*chi2);
430       Double_t errPhiIROC=TMath::Sqrt(covar(offset+isec+18, offset+isec+18)*chi2);
431       Double_t errPhiOROC=TMath::Sqrt(covar(offset+isec+36, offset+isec+36)*chi2);
432       matrixIROC.SetDx(dgxIROC);
433       matrixIROC.SetDy(dgyIROC);
434       matrixOROC.SetDx(dgxOROC);
435       matrixOROC.SetDy(dgyOROC);
436       rotIROC.SetAngles(drotIROC*TMath::RadToDeg(),0,0);
437       matrixIROC.Multiply(&matrixGX);
438       matrixIROC.Multiply(&matrixGY);
439       matrixIROC.Multiply(&rotIROC);
440       rotOROC.SetAngles(drotOROC*TMath::RadToDeg(),0,0);
441       matrixOROC.Multiply(&matrixGX);
442       matrixOROC.Multiply(&matrixGY);
443       matrixOROC.Multiply(&rotOROC);
444       if (side>0){
445         array->AddAt(matrixIROC.Clone(),isec);
446         array->AddAt(matrixOROC.Clone(),isec+36);
447       }
448       if (side<0){
449         array->AddAt(matrixIROC.Clone(),isec+18);
450         array->AddAt(matrixOROC.Clone(),isec+36+18);
451       }
452       Double_t rms=TMath::Sqrt(chi2); //chi2 normalized 1/npoints before
453       (*pcWorkspace)<<"align"<<
454         "isec="<<isec<<           // sector
455         "side="<<side<<           // side
456         "phi="<<phi<<             // phi
457         "rms="<<rms<<             // rms in cm
458         "cov.="<<&covar<<
459         // errors
460         "errYIROC="<<errYIROC<<      //error in local y
461         "errPhiIROC="<<errPhiIROC<<  //error in phi IROC
462         "errPhiOROC="<<errPhiOROC<<  //error in phi OROC
463         //
464         "dlyIROC="<<dlyIROC<<        //delta Local Y at refX=165 -IROC
465         "dlyIROC0="<<dlyIROC0<<      //delta Local Y at refX=0   -IROC
466         "dgxIROC="<<dgxIROC<<        //
467         "dgyIROC="<<dgyIROC<<
468         "drotIROC="<<drotIROC<<
469         //
470         "dlyOROC="<<dlyOROC<<       // assuming 0 at 
471         "dlyOROC0="<<dlyOROC0<<
472         "dgxOROC="<<dgxOROC<<
473         "dgyOROC="<<dgyOROC<<
474         "drotOROC="<<drotOROC<<
475         "\n";
476     } 
477   }
478   alignLocal->SetAlignSectors(array);
479   AliTPCCorrection::AddVisualCorrection(alignLocal,1000);
480   /*
481     Check:
482     alignLocal->CreateHistoDRPhiinXY(10,100,100)->Draw("colz");  // OK
483     treeDY->Draw("DeltaLookup(sector,localX,kZ,165.6,1,1001):localX:int(sector)","kZ>0","colz");
484     //
485     (( (*pcWorkspace)<<"align"))->GetTree()->Draw("dlyOROC-dlyOROC:isec","","*")
486     treeDY->Draw("tfitA:DeltaLookup(sector,localX,kZ,165.6,1,1001):sector","kZ>0","colz");
487
488   */
489   return alignLocal;
490 }
491
492 AliTPCCalibGlobalMisalignment *MakeQuadrantCorrection(TVectorD paramA, TVectorD paramC, TMatrixD covar, Double_t chi2){
493   /// Make a global alignmnet
494   /// side= 1 - A side
495   /// side=-1 - C side
496   /// Take a fit parameters and make a combined correction:
497   /// Only delta local Y and delta phi are fitted - not sensitivity for other parameters
498   /// GX and GY shift extracted per side.
499   ///
500   /// Algorithm:
501   ///   1. Loop over sectors
502   ///   2. Make combined AliTPCCalibGlobalMisalignment
503
504   AliTPCCalibGlobalMisalignment *alignLocalQuadrant  =new  AliTPCCalibGlobalMisalignment;  
505   //
506   Int_t offset=3+3*18;
507   TVectorD quadrantQ0(36);   //dly+=sign*(*fQuadrantQ0)[isec];  // shift in cm
508   TVectorD quadrantRQ0(36);  //dly+=sign*(*fQuadrantRQ0)[isec]*(pos[0]-xref);      
509   TVectorD quadrantQ1(36);   //dly+=sign*(*fQuadrantQ1)[isec];  // shift in cm
510   TVectorD quadrantRQ1(36);  //dly+=sign*(*fQuadrantRQ1)[isec]*(pos[0]-xref);      
511   TVectorD quadrantQ2(36);   //dly+=(*fQuadrantQ2)[isec];  // shift in cm
512   TVectorD quadrantRQ2(36);  //dly+=(*fQuadrantRQ2)[isec]*(pos[0]-xref);      
513
514   for (Int_t side=-1; side<2; side+=2){
515     Int_t offsetSec= (side==1) ? 0:18;
516     TVectorD &param = (side==1) ? paramA : paramC;
517     for (Int_t isec=0; isec<18; isec++){ 
518       Double_t  q0=-param[offset+isec+0*18];
519       Double_t rq0=-param[offset+isec+1*18]*0.001;
520       Double_t  q1=-param[offset+isec+2*18];
521       Double_t rq1=-param[offset+isec+3*18]*0.001;
522       Double_t  q2=-param[offset+isec+4*18];
523       Double_t rq2=-param[offset+isec+5*18]*0.001;
524       //
525       Double_t  sq0=TMath::Sqrt(covar(offset+isec+0*18,offset+isec+0*18)*chi2);
526       Double_t srq0=TMath::Sqrt(covar(offset+isec+1*18,offset+isec+1*18)*chi2)*0.001;
527       Double_t  sq1=TMath::Sqrt(covar(offset+isec+2*18,offset+isec+2*18)*chi2);
528       Double_t srq1=TMath::Sqrt(covar(offset+isec+3*18,offset+isec+3*18)*chi2)*0.001;
529       Double_t  sq2=TMath::Sqrt(covar(offset+isec+4*18,offset+isec+4*18)*chi2);
530       Double_t srq2=TMath::Sqrt(covar(offset+isec+5*18,offset+isec+5*18)*chi2)*0.001;
531       //
532       quadrantQ0[offsetSec+isec]=q0;
533       quadrantRQ0[offsetSec+isec]=rq0;
534       quadrantQ1[offsetSec+isec]=q1;
535       quadrantRQ1[offsetSec+isec]=rq1;
536       quadrantQ2[offsetSec+isec]=q2;
537       quadrantRQ2[offsetSec+isec]=rq2;
538       Double_t rms=TMath::Sqrt(chi2); //chi2 normalized 1/npoints before
539       (*pcWorkspace)<<"quadrant"<<
540         "isec="<<isec<<           // sector
541         "side="<<side<<           // side
542         "rms="<<rms<<             // rms in cm
543         "cov.="<<&covar<<
544         //
545         "q0="<<q0<<               // quadrant alignment
546         "rq0="<<rq0<<
547         "q1="<<q1<<
548         "rq1="<<rq1<<
549         "q2="<<q2<<
550         "rq2="<<rq2<<
551         "sq0="<<sq0<<             //rms of quadrant parameters 
552         "srq0="<<srq0<<
553         "sq1="<<sq1<<
554         "srq1="<<srq1<<
555         "sq2="<<sq2<<
556         "srq2="<<srq2<<
557         "\n";    
558     }
559   }
560   alignLocalQuadrant->SetQuadranAlign(&quadrantQ0, &quadrantRQ0, &quadrantQ1, &quadrantRQ1, &quadrantQ2, &quadrantRQ2);
561   AliTPCCorrection::AddVisualCorrection(alignLocalQuadrant,1100);  
562   /*
563     (( (*pcWorkspace)<<"quadrant"))->GetTree()->Draw("q0:isec","","*")
564
565     alignLocalQuadrant->CreateHistoDRPhiinXY(10,100,100)->Draw("colz");  //OK
566     
567     treeDY->Draw("tfitA-DeltaLookup(sector,localX,kZ,165.6,1,1000):DeltaLookup(sector,localX,kZ,165.6,1,1100):int(sector)","kZ>0&&oroc","colz");
568   */
569   return alignLocalQuadrant;
570 }
571
572 AliTPCFCVoltError3D* MakeEfieldCorrection(TVectorD paramA, TVectorD paramC, TMatrixD covar, Double_t chi2){
573   /// Make a global  AliTPCFCVoltError3D object
574   ///
575   /// Take a fit parameters and make a combined correction:
576   /// Only delta local Y and delta phi are fitted - not sensitivity for other parameters
577   /// GX and GY shift extracted per side.
578   ///
579   /// Algorithm:
580   ///   1. Loop over sectors
581   ///   2. Make combined AliTPCCalibGlobalMisalignment
582
583   Int_t offset=3+9*18;
584   //
585   AliTPCFCVoltError3D* corrField = new AliTPCFCVoltError3D;
586   //
587   Double_t rotIROCA=paramA[offset+0];
588   Double_t rotIROCC=paramC[offset+0];
589   Double_t rotOROCA=paramA[offset+1];
590   Double_t rotOROCC=paramC[offset+1];
591   //
592   corrField->SetRotatedClipVoltA(0,paramA[offset+0]);
593   corrField->SetRotatedClipVoltC(0,paramC[offset+0]);
594   corrField->SetRotatedClipVoltA(1,paramA[offset+1]);
595   corrField->SetRotatedClipVoltC(1,paramC[offset+1]);
596   {for (Int_t isec=0; isec<18; isec++){
597       corrField->SetRodVoltShiftA(isec,paramA[offset+2+36+isec]);     // rod shift IFC
598       corrField->SetRodVoltShiftA(18+isec,paramA[offset+2+isec]);     // rod shift OFC
599       corrField->SetRodVoltShiftC(isec,paramC[offset+2+36+isec]);    // rod shift IFC
600       corrField->SetRodVoltShiftC(18+isec,paramC[offset+2+isec]);    // rod shift OFC
601       Double_t rodIROCA=paramA[offset+2+36+isec];
602       Double_t rodIROCC=paramC[offset+2+36+isec];
603       Double_t rodOROCA=paramA[offset+2+isec];
604       Double_t rodOROCC=paramC[offset+2+isec];
605       //
606       Double_t srodIROCA=TMath::Sqrt(covar(offset+2+36+isec,offset+2+36+isec)*chi2);
607       Double_t srodOROCA=TMath::Sqrt(covar(offset+2+isec,offset+2+isec)*chi2);
608       Double_t phi= (Double_t(isec)+0.5)*TMath::Pi()/9.;
609       Double_t rms=TMath::Sqrt(chi2); //chi2 normalized 1/npoints before
610       (*pcWorkspace)<<"field"<<
611         "isec="<<isec<<           // sector
612         //"side="<<side<<           // side
613         "phi="<<phi<<             // phi
614         "rms="<<rms<<             // rms in cm
615         "cov.="<<&covar<<
616         //
617         "rotIROCA="<<rotIROCA<<  
618         "rotIROCC="<<rotIROCC<<
619         "rotOROCA="<<rotOROCA<<
620         "rotOROCC="<<rotOROCC<<
621         //
622         "rodIROCA="<<rodIROCA<<
623         "rodIROCC="<<rodIROCC<<
624         "rodOROCA="<<rodOROCA<<
625         "rodOROCC="<<rodOROCC<<
626         //
627         "srodIROCA="<<srodIROCA<<
628         "srodOROCA="<<srodOROCA<<
629         "srodIROCC="<<srodIROCA<<
630         "srodOROCC="<<srodOROCA<<
631         "\n";
632     }    
633   }
634   corrField->SetOmegaTauT1T2(0,1.,1.);
635   corrField->InitFCVoltError3D();
636   AliTPCCorrection::AddVisualCorrection(corrField,100);  
637   return corrField;
638 }
639
640
641 void FitRodShift(Bool_t flagIFCcopper = kTRUE) {
642   /// Main fit function
643   /// In total 440 parameters to fit using global fit:
644   ///   1. Rotation and translation for each sector (72x2)
645   ///   2. Rotation and translation for each quadrant of OROC (36x4x2)
646   ///   3. Rod/strip shifts in IFC and OFC (18 sectors x 2 sides x 2 )
647   ///   4. Rotated clips in IFC and OFC
648
649   LoadTrees();
650   LoadModels();
651   RegisterAlignFunction();
652   MakeAliases();
653
654   pcWorkspace= new TTreeSRedirector("fitAlignLookup.root");
655   TFormula::SetMaxima(10000);
656   //
657   // 1. fit Global
658   //
659   Double_t chi2G=0;   Int_t npointsG=0;  TVectorD paramG;  TMatrixD covarG;
660   TString fstringGlobal="";  
661   fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,10-0)++"; // deltaGX
662   fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,11-0)++"; // deltaGY
663   fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,10-0)*side++"; // deltaGX - sides
664   fstringGlobal+="DeltaLookup(sector,localX,kZ,165.6,1,11-0)*side++"; // deltaGY -sides
665   //
666   TString *strFitGlobal = TStatToolkit::FitPlane(treeDY,"drphi", fstringGlobal.Data(),"1", chi2G,npointsG,paramG,covarG,-1,0, 10000000, kTRUE);
667   treeDY->SetAlias("fitYGlobal",strFitGlobal->Data());
668   strFitGlobal->Tokenize("++")->Print();
669   printf("chi2=%f\n",TMath::Sqrt(chi2G/npointsG));
670   // testplots - if necessary - e.g.
671   // treeDY->Draw("AliTPCCorrection::GetCorrSector(sector,localX,kZ,1,2):sector:localX","Cut","colz")
672
673   // FIT PREPARATION +++++++++++++++++++++++++++++++++
674   // define the fit function
675
676   TString  fstring="";         
677
678   //
679   // Alignment
680   //
681   {
682     fstring+="DeltaLookup(sector,localX,kZ,165.6,1,10-0)++"; // deltaGX
683     fstring+="DeltaLookup(sector,localX,kZ,165.6,1,11-0)++"; // deltaGY
684     // alignment IROC
685     for (Int_t i=0;i<18;i++){
686       fstring+=Form("dyI_%d++",i);  // alignment - linear shift (in cm)
687     }
688     for (Int_t i=0;i<18;i++){
689       fstring+=Form("drotI_%d++",i);  // alignment - rotation (in mrad)
690     }
691     // alignment OROC
692     for (Int_t i=0;i<18;i++){ // shift of OROC is not allowed
693       //fstring+=Form("dyO_%d++",i);  // alignment - linear shift (in cm)
694     }
695     for (Int_t i=0;i<18;i++){
696       fstring+=Form("drotO_%d++",i);  // alignment - rotation (in mrad)
697     }
698     //
699     for (Int_t i=0;i<18;i++){
700       fstring+=Form("dqO0_%d++",i);  // alignment - quadrant shift OROC in
701     }
702     for (Int_t i=0;i<18;i++){
703       fstring+=Form("dqOR0_%d++",i);  // alignment - quadrant rotation OROC in
704     }
705     for (Int_t i=0;i<18;i++){
706       fstring+=Form("dqO1_%d++",i);  // alignment - quadrant shift OROC out
707     }
708     for (Int_t i=0;i<18;i++){
709       fstring+=Form("dqOR1_%d++",i);  // alignment - quadrant rotation OROC out
710     }
711     for (Int_t i=0;i<18;i++){
712       fstring+=Form("dqO2_%d++",i);  // alignment - quadrant shift OROC out
713     }
714     for (Int_t i=0;i<18;i++){
715       fstring+=Form("dqOR2_%d++",i);  // alignment - quadrant rotation OROC out
716     }
717   }
718   //
719   // Field distortion
720   //
721   {
722     // rotated clip - IFC --------------------
723     treeDY->SetAlias("rotClipI","DeltaLookup(sector,localX,kZ,165.6,1,3+0)");
724     fstring+="rotClipI++";
725     // rotated clip - OFC --------------------
726     treeDY->SetAlias("rotClipO","DeltaLookup(sector,localX,kZ,165.6,1,0+0)");
727     fstring+="rotClipO++";
728     
729     // shifted rod -  OFC
730     for (Int_t i=0;i<18;i++){
731       fstring+=Form("rodStripO_%d++",i);
732     }    
733     // Shifted cooper rod OFC     
734     for (Int_t i=0;i<18;i++){
735       fstring+=Form("coppRodO_%d++",i);
736     }    
737     // shifted rod - IFC 
738     for (Int_t i=0;i<18;i++){
739       fstring+=Form("rodStripI_%d++",i);
740     }    
741     // Shifted cooper rod IFC 
742     if (flagIFCcopper) for (Int_t i=0;i<18;i++){
743         fstring+=Form("coppRodI_%d++",i);
744     }  
745   }
746   //fstring+=fstringGlobal;
747  
748   // FIT ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
749
750   Double_t chi2A=0;   Int_t npointsA=0;  TVectorD paramA;  TMatrixD covarA;
751   Double_t chi2C=0;   Int_t npointsC=0;  TVectorD paramC;  TMatrixD covarC;
752   
753   
754   printf("Fitting A side\n");
755   TCut cutAllA = "Cut&&kZ>0";
756   TString *strFitGA = TStatToolkit::FitPlane(treeDY,"drphi:errY", fstring.Data(),"kZ>0", chi2A,npointsA,paramA,covarA,-1,0, 10000000, kTRUE);
757   treeDY->SetAlias("tfitA",strFitGA->Data());
758   // strFitGA->Tokenize("++")->Print();
759   printf("chi2=%f\n",TMath::Sqrt(chi2A/npointsA));
760   printf("Sigma Y:\t%f (mm)\n",10.*TMath::Sqrt(covarA(3,3)*chi2A/npointsA));
761   printf("IROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarA(3+18,3+18)*chi2A/npointsA));
762   printf("OROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarA(3+36,3+36)*chi2A/npointsA));
763
764   printf("Fitting C side\n");
765   TCut cutAllC = "Cut&&kZ<0";
766   TString *strFitGC = TStatToolkit::FitPlane(treeDY,"drphi:errY", fstring.Data(),"kZ<0", chi2C,npointsC,paramC,covarC,-1,0, 10000000, kTRUE);
767   treeDY->SetAlias("tfitC",strFitGC->Data());
768   //  strFitGC->Tokenize("++")->Print();
769   printf("chi2=%f\n",TMath::Sqrt(chi2C/npointsC));
770   printf("Sigma Y:\t%f (mm)\n",10.*TMath::Sqrt(covarC(3,3)*chi2C/npointsC));
771   printf("IROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarC(3+18,3+18)*chi2C/npointsC));
772   printf("OROC Sigma Angle:\t%f (mrad)\n",TMath::Sqrt(covarC(3+36,3+36)*chi2C/npointsC));
773
774   //
775   // make combined correction
776   //
777   TH2 *hdist =0;
778   //
779   AliTPCCalibGlobalMisalignment * alignLocal = MakeAlignCorrection(paramA, paramC, covarA, chi2A/npointsA);
780   alignLocal->SetName("alignLocal");
781   alignLocal->Write("alignLocal");
782   hdist = alignLocal->CreateHistoDRPhiinXY(10,250,250);
783   hdist->SetName("AlignLocalAside"); hdist->SetTitle("Alignment map (A side)");
784   hdist->Write("AlignLocalAside");
785   hdist = alignLocal->CreateHistoDRPhiinXY(-10,250,250);
786   hdist->SetName("AlignLocalCside"); hdist->SetTitle("Alignment map (C side)");
787   hdist->Write("AlignLocalCside");
788   //
789   AliTPCCalibGlobalMisalignment *alignQuadrant = MakeQuadrantCorrection(paramA, paramC,  covarA, chi2A/npointsA);
790   alignQuadrant->SetName("alignQuadrant");
791   alignQuadrant->Write("alignQuadrant");
792   hdist = alignQuadrant->CreateHistoDRPhiinXY(10,250,250);
793   hdist->SetName("AlignQuadrantAside"); hdist->SetTitle("Quadrant Alignment map (A side)");
794   hdist->Write("AlignQuadrantAside");
795   hdist = alignQuadrant->CreateHistoDRPhiinXY(-10,250,250);
796   hdist->SetName("AlignQuadrantCside"); hdist->SetTitle("Quadrant Alignment map (C side)");
797   hdist->Write("AlignQuadrantCside");
798   //
799   AliTPCFCVoltError3D* corrField = MakeEfieldCorrection(paramA, paramC, covarA, chi2A/npointsA);
800   corrField->SetName("corrField");
801   corrField->Write("corrField");
802
803   hdist = corrField->CreateHistoDRPhiinXY(10,250,250);
804   hdist->SetName("AlignEfieldAside"); hdist->SetTitle("Efield Alignment map (A side)");
805   hdist->Write("AlignEfieldAside");
806   hdist = corrField->CreateHistoDRPhiinXY(-10,250,250);
807   hdist->SetName("AlignEfieldCside"); hdist->SetTitle("Efield Alignment map (C side)");
808   hdist->Write("AlignEfieldCside");
809
810
811   //
812   delete pcWorkspace;
813   MakeQA();
814   return;  
815 }
816
817
818 void MakeQA(){
819   LoadTrees();
820   LoadModels();
821   RegisterAlignFunction();
822   MakeAliases();
823   TFile f("fitAlignLookup.root","update");  
824   AliTPCCalibGlobalMisalignment*     alignLocal =  (AliTPCCalibGlobalMisalignment*)f.Get("alignLocal");
825   AliTPCCalibGlobalMisalignment*  alignQuadrant =  (AliTPCCalibGlobalMisalignment*)f.Get("alignQuadrant");
826   AliTPCFCVoltError3D      *corrField= (AliTPCFCVoltError3D*)f.Get("corrField");
827   AliTPCCorrection::AddVisualCorrection(alignLocal,1000);
828   AliTPCCorrection::AddVisualCorrection(alignQuadrant,1100);
829   AliTPCCorrection::AddVisualCorrection(corrField,100);
830   FitFunctionQA();
831   DrawAlignParam();
832   f.Close();
833 }
834
835 void FitFunctionQA(){
836   ///
837
838   TH1 *his=0;
839   TCanvas *canvasDist= new TCanvas("FitQA","fitQA",1200,800);
840   canvasDist->Divide(2,2);
841   //
842   canvasDist->cd(1)->SetLogy(kFALSE); canvasDist->cd(1)->SetRightMargin(0.15);
843   treeDY->Draw("10*mean:sector:localX","kZ<0&&localX<134","colz"); 
844   his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi1");
845   his->GetXaxis()->SetTitle("Sector"); 
846   his->GetZaxis()->SetTitle("R (cm)"); 
847   his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
848   his->Draw("cont2"); 
849   his->Draw("colz");
850   //
851   canvasDist->cd(2)->SetLogy(kFALSE); canvasDist->cd(2)->SetRightMargin(0.15);
852   treeDY->Draw("10*mean:sector:localX","kZ<0&&localX>160","colz"); 
853   his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi2");
854   his->GetXaxis()->SetTitle("Sector"); 
855   his->GetZaxis()->SetTitle("R (cm)"); 
856   his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)");
857   his->Draw("cont2"); 
858   his->Draw("colz");
859   //
860   canvasDist->cd(3)->SetLogy(kFALSE); canvasDist->cd(3)->SetRightMargin(0.15);
861   treeDY->Draw("10*(mean-dAll):sector:localX","kZ<0&&localX<134","colz"); 
862   his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi3"); his->SetTitle("Delta #RPhi");
863   his->GetXaxis()->SetTitle("Sector"); 
864   his->GetZaxis()->SetTitle("R (cm)"); 
865   his->GetYaxis()->SetTitle("#Delta_{r#phifit}-#Delta_{r#phi} (mm)"); 
866   his->Draw("cont2"); 
867   his->Draw("colz");
868   //
869   canvasDist->cd(4)->SetLogy(kFALSE); canvasDist->cd(4)->SetRightMargin(0.15);
870   treeDY->SetMarkerColor(1);
871   treeDY->Draw("10*mean:10*dAll:localX","","colz"); 
872   his= (TH1*)treeDY->GetHistogram()->Clone();   his->SetName("DeltaRPhi4");
873   his->GetXaxis()->SetTitle("Fit value #Delta_{r#phi} (mm)"); 
874   his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); 
875   his->GetZaxis()->SetTitle("R (cm)"); 
876   his->Draw("cont2");
877   his->Draw("colz");
878   //
879   canvasDist->Write("FitQA");
880 }
881
882 void DrawAlignParam(){
883   ///
884
885   TFile f("fitAlignLookup.root");
886   TTree * treeAlign=(TTree*)f.Get("align");
887   TTree * treeQuadrant=(TTree*)f.Get("quadrant");
888   TCanvas *canvasAlign=new TCanvas("align","align",1000,800); 
889   TH1 *his=0;
890   canvasAlign->Divide(2,2);
891   treeAlign->SetMarkerStyle(25);
892   treeQuadrant->SetMarkerStyle(25);
893   gStyle->SetOptStat(kTRUE);
894   //
895   canvasAlign->cd(1); canvasAlign->cd(1)->SetRightMargin(0.15);
896   treeAlign->Draw("10*dlyIROC*side:isec:1+side","","colz");
897   his= (TH1*)treeAlign->GetHistogram()->Clone();   his->SetName("dlyIROC");his->SetTitle("IROC Alignment #Delta_{r#phi}");
898   his->GetXaxis()->SetTitle("sector"); 
899   his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); 
900   his->GetZaxis()->SetTitle("side");
901   his->Draw("cont2"); 
902   his->Draw("colz");
903   //
904   canvasAlign->cd(2); canvasAlign->cd(2)->SetRightMargin(0.15);
905   treeAlign->Draw("1000*drotIROC*side:isec:1+side","","colz");
906   his= (TH1*)treeAlign->GetHistogram()->Clone();   his->SetName("drotIROC");his->SetTitle("IROC Angular Alignment #Delta_{r#phi}");
907   his->GetXaxis()->SetTitle("sector"); 
908   his->GetYaxis()->SetTitle("#Delta_{r#phi} (mrad)"); 
909   his->GetZaxis()->SetTitle("side");
910   his->Draw("cont2");
911   his->Draw("colz");
912   //
913   canvasAlign->cd(4);canvasAlign->cd(4)->SetRightMargin(0.15);
914   treeAlign->Draw("1000*drotOROC*side:isec:1+side","","colz");
915   his= (TH1*)treeAlign->GetHistogram()->Clone();   his->SetName("drotOROC");his->SetTitle("OROC Angular Alignment #Delta_{r#phi}");
916   his->GetXaxis()->SetTitle("sector"); 
917   his->GetYaxis()->SetTitle("#Delta_{r#phi} (mrad)"); 
918   his->GetZaxis()->SetTitle("side");
919   his->Draw("cont2");
920   his->Draw("colz");
921   //
922   canvasAlign->cd(3);canvasAlign->cd(3)->SetRightMargin(0.15);
923   treeQuadrant->Draw("10*q2*side:isec:1+side","","colz");
924   his= (TH1*)treeQuadrant->GetHistogram()->Clone();   his->SetName("drphiOROC");his->SetTitle("OROC Alignment  Outer Quadrant #Delta_{r#phi}");
925   his->GetXaxis()->SetTitle("sector"); 
926   his->GetYaxis()->SetTitle("#Delta_{r#phi} (mm)"); 
927   his->GetZaxis()->SetTitle("side");
928   his->Draw("cont2");
929   his->Draw("colz");
930
931
932 }