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 }