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