Global fit of alignment with E field distortion map (440 parameters)
[u/mrichter/AliRoot.git] / TPC / CalibMacros / FitRodShift.C
CommitLineData
32b5322b 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
81TTreeSRedirector *pcWorkspace= 0; // Workspace to store the fit parameters and QA histograms
82TTree * treeDY= 0; //distortion tree - Dy
83TTree * treeDZ= 0; //distortion tree - Dz
84// fit models ...
85AliTPCFCVoltError3D *rotOFC =0; // fit models
86AliTPCFCVoltError3D *rodOFC1 =0;
87AliTPCFCVoltError3D *rodOFC2 =0;
88AliTPCFCVoltError3D *rotIFC =0;
89AliTPCFCVoltError3D *rodIFC1 =0;
90AliTPCFCVoltError3D *rodIFC2 =0;
91AliTPCFCVoltError3D *rodAll =0; //all lookup initialized
92//
93AliTPCComposedCorrection *corFit0 = 0; // composed correction
94void FitFunctionQA();
95void MakeQA();
96void DrawAlignParam();
97
98void 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
115void 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
225void 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
284void 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
324Double_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
335void 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
393AliTPCCalibGlobalMisalignment *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
500AliTPCCalibGlobalMisalignment *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
581AliTPCFCVoltError3D* 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
651void 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
829void 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
846void 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
894void 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}