MUON/MUONSurveyCh1.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id:$ */
  17
  18 /// \ingroup macros
  19 /// \file MUONSurveyCh8L.C
  20 /// \brief Macro to process survey and photogrammetry data of chamber 1
  21 ///
  22 /// Macro loads the survey data from .txt file using AliSurveyObj.
  23 /// Macro uses AliMUONSurvey... classes.
  24 /// The transformations of the detection elements or chambers can be obtained
  25 /// in three ways:
  26 /// - A: Fit of local to global transformation using the fixed button targets.
  27 /// - B: Fit a plane to the sticker targets -> psi, theta
  28 ///      Use above psi and theta and fit remaining 4 parameters using the fixed
  29 ///      button targets
  30 /// - C: Fit a plane to the sticker targets -> psi, theta
  31 ///      Use above psi and theta to calculate xc, yc, zc and phi by solving
  32 ///      the equations from a local to global transformation of the fixed button
  33 ///      targets
  34 ///
  35 /// Methods A and B are prefered to C, and B is better if sticker targets are
  36 /// available and lie on a plane!
  37 /// For slats only methods B and C can be used.
  38 /// Various histograms are filled and printed for monitoring.
  39 /// MisAlignment object is then created.
  40 ///
  41 /// \author Javier Castillo
  42
  43 #if !defined(__CINT__) || defined(__MAKECINT__)
  44
  45 #include "AliMUONGeometryTransformer.h"
  46 #include "AliMUONGeometryModuleTransformer.h"
  47 #include "AliMUONGeometryDetElement.h"
  48 #include "AliMUONGeometryMisAligner.h"
  49 #include "AliMUONSurveyChamber.h"
  50 #include "AliMUONSurveyDetElem.h"
  51 #include "AliMUONSurveyUtil.h"
  52
  53 #include "AliSurveyObj.h"
  54 #include "AliSurveyPoint.h"
  55 #include "AliGeomManager.h"
  56 #include "AliCDBManager.h"
  57 #include "AliCDBMetaData.h"
  58 #include "AliCDBId.h"
  59 #include "AliAlignObjMatrix.h"
  60 #include "AliAlignObj.h"
  61
  62 #include <TROOT.h>
  63 #include <TGeoManager.h>
  64 #include <TClonesArray.h>
  65 #include <TObjArray.h>
  66 #include <TArrayD.h>
  67 #include <TMath.h>
  68 #include <TString.h>
  69 #include <TFitter.h>
  70 #include <TH2.h>
  71 #include <TF2.h>
  72 #include <TGraphErrors.h>
  73 #include <TCanvas.h>
  74 #include <TPad.h>
  75 #include <TLine.h>
  76 #include <TPostScript.h>
  77 #include <TPaveLabel.h>
  78 #include <TStyle.h>
  79 #include <TFile.h>
  80 #include <TMatrixDSym.h>
  81
  82 #include <fstream>
  83
  84 #endif
  85
  86 void MUONSurveyCh1() {
  87
  88   TString str;
  89   TString title;
  90
  91   Bool_t bMonitor = kTRUE;
  92   Bool_t bOCDB = kTRUE;
  93   Bool_t saveps = kFALSE;
  94   const int cWidth = (int)(700*(29./21.));
  95   const int cHeight = 700;
  96   const int filetype  = 112; // landscape
  97
  98   Int_t chId = 0;
  99   Int_t nChs = 2;
 100   Int_t nDetElems = 4;
 101   Int_t nDetElemsI = 2;
 102   //  Int_t nDetElemsO = 2;
 103   Int_t iDetElemToDetElemId[4] = {100,101,102,103};
 104   Int_t iDetElemPseudoIdToDetElem[4] = {0,1,2,3};
 105   Int_t iDetElemsOfChamber[2][2] = {{0,3},
 106                                     {1,2}};
 107
 108   TObjArray *myChamberArray = new TObjArray();
 109   myChamberArray->Add(new AliMUONSurveyChamber(chId));
 110   myChamberArray->Add(new AliMUONSurveyChamber(chId));
 111
 112   AliMUONSurveyChamber *myChamberI = 0x0;
 113   AliMUONSurveyChamber *myChamberO = 0x0;
 114   AliMUONSurveyChamber *myChamber = 0x0;
 115   AliMUONSurveyDetElem *myDetElem = 0x0;
 116
 117   myChamberI = (AliMUONSurveyChamber*)myChamberArray->At(0);
 118   myChamberI->GetSurveyObj()->FillFromLocalFile("~/phgwrk/Reports/AliceSt1_TC1_1171d.txt");
 119   myChamberO = (AliMUONSurveyChamber*)myChamberArray->At(1);
 120   myChamberO->GetSurveyObj()->FillFromLocalFile("~/phgwrk/Reports/AliceSt1_TC1_1171d.txt");
 121
 122   myChamber = myChamberI;
 123   myChamber->PrintSurveyReport();
 124
 125   // Chamber & DetElems button targets local coordinates
 126   AliSurveyObj *lSO = new AliSurveyObj();
 127   lSO->FillFromLocalFile("$ALICE_ROOT/MUON/data/MUONTrackingTargetsLocal.txt");
 128
 129   // Set survey targets for chambers
 130   myChamberI->AddGButtonTargets(TString("T4_520"),3);
 131   myChamberI->AddStickerTargets(TString("T4_620"),6);
 132   myChamberI->AddLButtonTargets(lSO->GetData(),Form("%d_520",chId+1),3);
 133   myChamberI->AddGButtonTargets(TString("T4_510"),3);
 134   myChamberI->AddStickerTargets(TString("T4_610"),6);
 135   myChamberI->AddLButtonTargets(lSO->GetData(),Form("%d_510",chId+1),3);
 136   myChamberO->AddGButtonTargets(TString("T4_520"),3);
 137   myChamberO->AddStickerTargets(TString("T4_620"),6);
 138   myChamberO->AddLButtonTargets(lSO->GetData(),Form("%d_520",chId+1),3);
 139   myChamberO->AddGButtonTargets(TString("T4_510"),3);
 140   myChamberO->AddStickerTargets(TString("T4_610"),6);
 141   myChamberO->AddLButtonTargets(lSO->GetData(),Form("%d_510",chId+1),3);
 142
 143   // Set survey targets for detection elements
 144   for (int iCh =0; iCh<=1; iCh++) {
 145     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 146     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 147       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 148       myChamber->AddSurveyDetElem(iDetElemToDetElemId[iDetElem]);
 149       TString baseName;
 150       myDetElem =  myChamber->GetDetElem(iDetElemI);
 151       if (myDetElem) {
 152         baseName = Form("T4_53%d",iDetElem+1);
 153         myDetElem->AddGButtonTargets(baseName,3);
 154         baseName = Form("%d_53%d",chId+1,iDetElem+1);
 155         myDetElem->AddLButtonTargets(lSO->GetData(),baseName,3);
 156       }
 157     }
 158   }
 159   printf("All targets added! \n");
 160
 161   Double_t **lCenDetElem = new Double_t*[nDetElems+nChs];
 162   Double_t **lRotDetElem = new Double_t*[nDetElems+nChs];
 163   Double_t **lDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
 164   Double_t **lDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
 165   Double_t **lDiffThCenDetElem0 = new Double_t*[nDetElems+nChs];
 166   Double_t **lDiffThRotDetElem0 = new Double_t*[nDetElems+nChs];
 167   Double_t **lDeltaDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
 168   Double_t **lDeltaDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
 169
 170   for (int iDetElem=0; iDetElem<nDetElems+nChs; iDetElem++) {
 171     lCenDetElem[iDetElem] = new Double_t[3];
 172     lRotDetElem[iDetElem] = new Double_t[3];
 173     lDiffCenDetElem0[iDetElem] = new Double_t[3];
 174     lDiffRotDetElem0[iDetElem] = new Double_t[3];
 175     lDiffThCenDetElem0[iDetElem] = new Double_t[3];
 176     lDiffThRotDetElem0[iDetElem] = new Double_t[3];
 177     lDeltaDiffCenDetElem0[iDetElem] = new Double_t[3];
 178     lDeltaDiffRotDetElem0[iDetElem] = new Double_t[3];
 179   }
 180
 181   TGeoCombiTrans dtrfDetElem[nDetElems+nChs];
 182   TGeoCombiTrans localTrfDetElem[nDetElems+nChs];
 183   TGeoCombiTrans localTrfThDetElem[nDetElems+nChs];
 184
 185   // Import TGeo geometry
 186   char* geoFilename = "geometry.root";
 187   if ( ! AliGeomManager::GetGeometry() ) {
 188     AliGeomManager::LoadGeometry(geoFilename);
 189     if (! AliGeomManager::GetGeometry() ) {
 190       printf("MUONSurveyCh%d: getting geometry from file %s failed\n", chId+1,geoFilename);
 191       return;
 192     }
 193     cout << "geometry imported" << endl;
 194   }
 195
 196   AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
 197   transform->LoadGeometryData();
 198
 199   TGeoCombiTrans trfThChamber;
 200   TGeoCombiTrans trfThDetElem;
 201
 202   for (int iCh =0; iCh<=1; iCh++) {
 203     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 204
 205     trfThChamber = TGeoCombiTrans(*transform->GetModuleTransformerByDEId(iDetElemToDetElemId[iCh*nDetElemsI])->GetTransformation());
 206     trfThChamber.Print();
 207     myChamber->SetLocalTransformation(new TGeoCombiTrans(trfThChamber),kTRUE);
 208
 209     // Set Chamber plane function
 210     cout << "Setting plane for Chamber" << iCh+1 << " ..." << endl;
 211     myChamber->SetPlane(Form("fChamber%d",iCh+1));
 212     myChamber->SetPlaneParameters(0.,0.,0.);
 213
 214     // Fit a plane to sticker targets
 215     Double_t lCChi2 = myChamber->FitPlane();
 216     cout << "... chi2 = " << lCChi2 << " ..." << endl;
 217
 218     // Get best transformation from fitting method
 219     // (local to global transformation)
 220     cout << "Trying fitting method for chamber " << iCh << endl;
 221     myChamber->SurveyToAlign();
 222     //    myChamber->SurveyToAlign(myChamber->GetPlane()->GetParameter(0),myChamber->GetPlane()->GetParameter(1),myChamber->GetPlane()->GetParError(0),myChamber->GetPlane()->GetParError(1));
 223
 224     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 225       myDetElem =  myChamber->GetDetElem(iDetElemI);
 226       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 227
 228       trfThDetElem.Clear();
 229       trfThDetElem = TGeoCombiTrans(*transform->GetDetElement(iDetElemToDetElemId[iDetElem])->GetLocalTransformation());
 230       trfThDetElem.Print();
 231
 232       myDetElem->SetLocalTransformation(new TGeoCombiTrans(trfThDetElem),kTRUE);
 233
 234       for (int iCor=0; iCor<3; iCor++){
 235         lCenDetElem[iDetElem][iCor] = 0;
 236         lRotDetElem[iDetElem][iCor] = 0;
 237       }
 238
 239       if (bMonitor){
 240         // MONITOR: Draw graph with sticker targets for plane fitting
 241         myDetElem->DrawSTargets();
 242         gPad->Update();
 243       }
 244
 245       // Get DetElem transformation.
 246       // Psi and Theta are obtained by fitting a plane to the sticker targets.
 247       // Then Xc, Yc, Zc and Phi are obtained by solving the equations to the ref.
 248       // syst. transformation of the button targets
 249
 250       // Set DetElem plane function
 251       cout << "Setting plane for DetElem" << iDetElem+1 << " ..." << endl;
 252       myDetElem->SetPlane(Form("fDetElem%d",iDetElem+1));
 253       myDetElem->SetPlaneParameters(0.,0.,3.);
 254
 255       // Fit a plane to sticker targets
 256       Double_t lChi2 = myDetElem->FitPlane();
 257       cout << "... chi2 = " << lChi2 << " ..." << endl;
 258
 259       lRotDetElem[iDetElem][0]=TMath::ATan(myDetElem->GetPlane()->GetParameter(0));
 260       lRotDetElem[iDetElem][1]=TMath::ATan(myDetElem->GetPlane()->GetParameter(1));
 261
 262       // Calculate Mean transformation using previous plane fit
 263       // and pairs of button targets
 264       myDetElem->CalculateMeanTransf(lCenDetElem[iDetElem],lRotDetElem[iDetElem]);
 265
 266       cout << "DetElem" << iDetElem+1 << " : mycen(" << lCenDetElem[iDetElem][0] << "," << lCenDetElem[iDetElem][1] << "," << lCenDetElem[iDetElem][2] << "); rot(" << lRotDetElem[iDetElem][0] << "," << lRotDetElem[iDetElem][1] << "," << lRotDetElem[iDetElem][2] << ")"  << endl;
 267
 268
 269       // Get best transformation from fitting method
 270       // (local to global transformation)
 271       cout << "Trying fitting method for DetElemId " << iDetElemToDetElemId[iDetElem] << endl;
 272       myDetElem->SurveyToAlign();
 273       //      myDetElem->SurveyToAlign(lRotDetElem[iDetElem][0],lRotDetElem[iDetElem][1],myDetElem->GetPlane()->GetParError(0),myDetElem->GetPlane()->GetParError(1));
 274     }
 275   }
 276
 277   // Print found transformation of Detection Element (plane fit + loop)
 278   for (int iCh =0; iCh<=1; iCh++) {
 279     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 280       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 281       cout << "DetElem" << iDetElem+1 << " : mycen(" << lCenDetElem[iDetElem][0] << "," << lCenDetElem[iDetElem][1] << "," << lCenDetElem[iDetElem][2] << "); rot(" << lRotDetElem[iDetElem][0] << "," << lRotDetElem[iDetElem][1] << "," << lRotDetElem[iDetElem][2] << ")"  << endl;
 282     }
 283   }
 284
 285   // Print Theoretical transformation of Detection Element
 286   for (int iCh =0; iCh<=1; iCh++) {
 287     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 288     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 289       myChamber->GetDetElem(iDetElemI)->PrintLocalTrf();
 290     }
 291   }
 292
 293   // Print found delta transformation of Detection Element
 294   for (int iCh =0; iCh<=1; iCh++) {
 295     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 296     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 297       myChamber->GetDetElem(iDetElemI)->PrintAlignTrf();
 298     }
 299   }
 300
 301   //
 302   // Compare transformations to expected ones
 303   //
 304   Int_t iDetElemToPos[18] = {0, 1, 2, 3};
 305
 306   TGraphErrors *gDeltaDiffCenXDetElem0 = new TGraphErrors(nDetElems);
 307   TGraphErrors *gDeltaDiffCenYDetElem0 = new TGraphErrors(nDetElems);
 308   TGraphErrors *gDeltaDiffCenZDetElem0 = new TGraphErrors(nDetElems);
 309   TGraphErrors *gDeltaDiffPsiDetElem0 = new TGraphErrors(nDetElems);
 310   TGraphErrors *gDeltaDiffThtDetElem0 = new TGraphErrors(nDetElems);
 311   TGraphErrors *gDeltaDiffPhiDetElem0 = new TGraphErrors(nDetElems);
 312
 313   for (int iCh =0; iCh<=1; iCh++) {
 314     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 315     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 316       myChamber->GetDetElem(iDetElemI)->GetAlignTrf()->Print();
 317     }
 318   }
 319
 320   for (int iCh =0; iCh<=1; iCh++) {
 321     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 322     // Store delta transformations to use for CDB entry creation
 323     dtrfDetElem[nDetElems+iCh].Clear();
 324     dtrfDetElem[nDetElems+iCh] = *(myChamber->GetAlignTrf());
 325
 326     // Those are for checks and visualizations
 327     localTrfDetElem[nDetElems+iCh].Clear();
 328     localTrfDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf())*(*(myChamber->GetAlignTrf())));
 329     localTrfDetElem[nDetElems+iCh].Print();
 330     lDiffCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfDetElem[nDetElems+iCh].GetTranslation();
 331     AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffRotDetElem0[nDetElems+iCh]);
 332
 333     localTrfThDetElem[nDetElems+iCh].Clear();
 334     localTrfThDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf()));
 335     localTrfThDetElem[nDetElems+iCh].Print();
 336     lDiffThCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfThDetElem[nDetElems+iCh].GetTranslation();
 337     AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffThRotDetElem0[nDetElems+iCh]);
 338
 339     for (int iCor=0; iCor<3; iCor++){
 340       lDeltaDiffCenDetElem0[nDetElems+iCh][iCor] = lDiffCenDetElem0[nDetElems+iCh][iCor]-lDiffThCenDetElem0[nDetElems+iCh][iCor];
 341       lDeltaDiffRotDetElem0[nDetElems+iCh][iCor] = lDiffRotDetElem0[nDetElems+iCh][iCor]-lDiffThRotDetElem0[nDetElems+iCh][iCor];
 342       if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]-=TMath::TwoPi();
 343       if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]+=TMath::TwoPi();
 344     }
 345
 346     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 347       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 348       myDetElem =  myChamber->GetDetElem(iDetElemI);
 349       // Store delta transformations to use for CDB entry creation
 350       dtrfDetElem[iDetElem].Clear();
 351       dtrfDetElem[iDetElem] = *(myDetElem->GetAlignTrf());
 352
 353       // Those are for checks and visualizations
 354       localTrfDetElem[iDetElem].Clear();
 355       localTrfDetElem[iDetElem] = (*(myDetElem->GetLocalTrf())*(*(myDetElem->GetAlignTrf())));
 356       //      localTrfDetElem[iDetElem] = (*(myDetElem->GetBaseTrf())*(*(myDetElem->GetAlignTrf())));
 357       localTrfDetElem[iDetElem].Print();
 358       lDiffCenDetElem0[iDetElem] = (Double_t*)localTrfDetElem[iDetElem].GetTranslation();
 359       AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[iDetElem].GetRotationMatrix(),lDiffRotDetElem0[iDetElem]);
 360 //       lDiffCenDetElem0[iDetElem] = lCenDetElem[iDetElem];
 361 //       lDiffRotDetElem0[iDetElem] = lRotDetElem[iDetElem];
 362
 363       localTrfThDetElem[iDetElem].Clear();
 364       localTrfThDetElem[iDetElem] = (*(myDetElem->GetLocalTrf()));
 365       //      localTrfThDetElem[iDetElem] = (*(myDetElem->GetBaseTrf()));
 366       localTrfThDetElem[iDetElem].Print();
 367       lDiffThCenDetElem0[iDetElem] = (Double_t*)localTrfThDetElem[iDetElem].GetTranslation();
 368       AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[iDetElem].GetRotationMatrix(),lDiffThRotDetElem0[iDetElem]);
 369
 370       for (int iCor=0; iCor<3; iCor++){
 371         lDeltaDiffCenDetElem0[iDetElem][iCor] = lDiffCenDetElem0[iDetElem][iCor]-lDiffThCenDetElem0[iDetElem][iCor];
 372         lDeltaDiffRotDetElem0[iDetElem][iCor] = lDiffRotDetElem0[iDetElem][iCor]-lDiffThRotDetElem0[iDetElem][iCor];
 373         if (lDeltaDiffRotDetElem0[iDetElem][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]-=TMath::TwoPi();
 374         if (lDeltaDiffRotDetElem0[iDetElem][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]+=TMath::TwoPi();
 375       }
 376
 377       gDeltaDiffCenXDetElem0->SetPoint(iDetElem,1e1*lDeltaDiffCenDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
 378       gDeltaDiffCenYDetElem0->SetPoint(iDetElem,1e1*lDeltaDiffCenDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
 379       gDeltaDiffCenZDetElem0->SetPoint(iDetElem,1e1*lDeltaDiffCenDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
 380       gDeltaDiffPsiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
 381       gDeltaDiffThtDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
 382       gDeltaDiffPhiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
 383       gDeltaDiffCenXDetElem0->SetPointError(iDetElem,1e1*myDetElem->GetFitter()->GetParError(0),0.);
 384       gDeltaDiffCenYDetElem0->SetPointError(iDetElem,1e1*myDetElem->GetFitter()->GetParError(1),0.);
 385       gDeltaDiffCenZDetElem0->SetPointError(iDetElem,1e1*myDetElem->GetFitter()->GetParError(2),0.);
 386       gDeltaDiffPsiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(3),0.);
 387       gDeltaDiffThtDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(4),0.);
 388       gDeltaDiffPhiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(5),0.);
 389     }
 390   }
 391
 392   // Apply the found misalignments to the geometry
 393   AliMUONGeometryTransformer *newTransform = AliMUONSurveyUtil::ReAlign(transform,chId,nDetElems,iDetElemPseudoIdToDetElem,dtrfDetElem,true);
 394   newTransform->WriteTransformations(Form("transform2ReAlignSurveyCh%d.dat",chId+1));
 395
 396   if(bOCDB){
 397     // Generate realigned data in local cdb
 398     const TClonesArray* array = newTransform->GetMisAlignmentData();
 399
 400     // Set the alignment resolution in the align objects for this chamber
 401     Double_t chResX = myChamberI->GetAlignResX();
 402     Double_t chResY = myChamberI->GetAlignResY();
 403     Double_t deResX = (myChamberI->GetMeanDetElemAlignResX()+myChamberO->GetMeanDetElemAlignResX())/2.;
 404     Double_t deResY = (myChamberI->GetMeanDetElemAlignResY()+myChamberO->GetMeanDetElemAlignResY())/2.;
 405
 406     printf("Chamber alignment resolutions: resX=%f , resY=%f\n",chResX,chResY);
 407     printf("Detection Elements alignment resolutions: resX=%f , resY=%f\n",deResX,deResY);
 408     chResX = TMath::Sqrt(0.1*0.1+chResX*chResX);
 409     chResY = TMath::Sqrt(0.1*0.1+chResY*chResY);
 410     AliMUONSurveyUtil::SetAlignmentResolution(array,chId,chResX,chResY,deResX,deResY);
 411
 412     // CDB manager
 413     AliCDBManager* cdbManager = AliCDBManager::Instance();
 414     cdbManager->SetDefaultStorage(Form("local://ReAlignSurveyCh%dCDB",chId+1));
 415
 416     AliCDBMetaData* cdbData = new AliCDBMetaData();
 417     cdbData->SetResponsible("Dimuon Offline project");
 418     cdbData->SetComment("MUON alignment objects with survey realignment");
 419     AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
 420     cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
 421   }
 422
 423   for(Int_t iCor=0; iCor<3; iCor++){
 424     for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
 425       cout << lDeltaDiffCenDetElem0[iDetElem][iCor] << " " << lDiffCenDetElem0[iDetElem][iCor] << " " << lDiffThCenDetElem0[iDetElem][iCor] << endl;
 426     }
 427     cout << endl;
 428   }
 429   cout << endl;
 430   for(Int_t iCor=0; iCor<3; iCor++){
 431     for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
 432       cout << lDeltaDiffRotDetElem0[iDetElem][iCor] << " " << lDiffRotDetElem0[iDetElem][iCor] << " " << lDiffThRotDetElem0[iDetElem][iCor] << endl;
 433     }
 434     cout << endl;
 435   }
 436
 437   // Plot the found misalignements (Units here are mm and mrad!)
 438   TH1F *myDetElemDeltaDiffCenX = new TH1F("myDetElemDeltaDiffCenX","myDetElemDeltaDiffCenX",100,-10,10);
 439   myDetElemDeltaDiffCenX->SetMaximum(nDetElems+1);
 440   myDetElemDeltaDiffCenX->SetMinimum(0);
 441   TH1F *myDetElemDeltaDiffCenY = new TH1F("myDetElemDeltaDiffCenY","myDetElemDeltaDiffCenY",100,-10,10);
 442   myDetElemDeltaDiffCenY->SetMaximum(nDetElems+1);
 443   myDetElemDeltaDiffCenY->SetMinimum(0);
 444   TH1F *myDetElemDeltaDiffCenZ = new TH1F("myDetElemDeltaDiffCenZ","myDetElemDeltaDiffCenZ",100,-20,20);
 445   myDetElemDeltaDiffCenZ->SetMaximum(nDetElems+1);
 446   myDetElemDeltaDiffCenZ->SetMinimum(0);
 447
 448   TH1F *myDetElemDeltaDiffRotX = new TH1F("myDetElemDeltaDiffRotX","myDetElemDeltaDiffRotX",100,-10,10);
 449   myDetElemDeltaDiffRotX->SetMaximum(nDetElems+1);
 450   myDetElemDeltaDiffRotX->SetMinimum(0);
 451   TH1F *myDetElemDeltaDiffRotY = new TH1F("myDetElemDeltaDiffRotY","myDetElemDeltaDiffRotY",100,-10,10);
 452   myDetElemDeltaDiffRotY->SetMaximum(nDetElems+1);
 453   myDetElemDeltaDiffRotY->SetMinimum(0);
 454   TH1F *myDetElemDeltaDiffRotZ = new TH1F("myDetElemDeltaDiffRotZ","myDetElemDeltaDiffRotZ",100,-6,6);
 455   myDetElemDeltaDiffRotZ->SetMaximum(nDetElems+1);
 456   myDetElemDeltaDiffRotZ->SetMinimum(0);
 457
 458   //
 459   // ******** Starting plots
 460   //
 461   TCanvas *canvas;
 462   TPad *pad;
 463   TPaveLabel *theTitle;
 464   gStyle->SetPalette(1);
 465
 466   TPostScript *ps = 0;
 467
 468   if( saveps ){
 469     ps = new TPostScript(Form("surveyChamber%d",chId+1),filetype);
 470     ps->NewPage();
 471   }
 472
 473   // Observed misalignments
 474   str = Form("Chamber %d",chId+1);
 475   TCanvas *cvn2 = new TCanvas("cvn2",str,cWidth,cHeight);
 476   canvas = cvn2;
 477   canvas->Range(0,0,21,29);
 478
 479   title = Form(" MisAlignments Chamber %d - L2G - In Single Frame ",chId+1);
 480   TPaveLabel *theTitle2 = new TPaveLabel(3,27.0,18,28.5,title,"br");
 481   theTitle = theTitle2;
 482   theTitle->SetFillColor(18);
 483   theTitle->SetTextFont(32);
 484   theTitle->SetTextSize(0.4);
 485   theTitle->SetTextColor(1);
 486   theTitle->Draw();
 487
 488   TPad *pad2 = new TPad("pad2","pad2",0.01,0.01,0.98,0.91,0);
 489   pad = pad2;
 490   pad->Draw();
 491   TLine *ch0Line = new TLine(0,1,0,2);
 492   TLine *ch1Line = new TLine(0,1,0,2);
 493   ch1Line->SetLineStyle(2);
 494   pad->Divide(3,2);
 495
 496   pad->cd(1);
 497   myDetElemDeltaDiffCenX->Draw();
 498   myDetElemDeltaDiffCenX->SetXTitle("#Delta[xc_{i}^{m}-xc_{i}^{th}] (mm)");
 499   myDetElemDeltaDiffCenX->SetYTitle("DetElem arbitrary ordering");
 500   gDeltaDiffCenXDetElem0->SetMarkerStyle(20);
 501   gDeltaDiffCenXDetElem0->Draw("P");
 502   ch0Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+0][0],0.5,1e1*lDeltaDiffCenDetElem0[nDetElems+0][0],1.5);
 503   ch0Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+0][0],3.5,1e1*lDeltaDiffCenDetElem0[nDetElems+0][0],4.5);
 504   ch1Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+1][0],1.5,1e1*lDeltaDiffCenDetElem0[nDetElems+1][0],3.5);
 505
 506   pad->cd(2);
 507   myDetElemDeltaDiffCenY->Draw();
 508   myDetElemDeltaDiffCenY->SetXTitle("#Delta[yc_{i}^{m}-yc_{i}^{th}] (mm)");
 509   myDetElemDeltaDiffCenY->SetYTitle("DetElem arbitrary ordering");
 510   gDeltaDiffCenYDetElem0->SetMarkerStyle(20);
 511   gDeltaDiffCenYDetElem0->Draw("P");
 512   ch0Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+0][1],0.5,1e1*lDeltaDiffCenDetElem0[nDetElems+0][1],1.5);
 513   ch0Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+0][1],3.5,1e1*lDeltaDiffCenDetElem0[nDetElems+0][1],4.5);
 514   ch1Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+1][1],1.5,1e1*lDeltaDiffCenDetElem0[nDetElems+1][1],3.5);
 515
 516   pad->cd(3);
 517   myDetElemDeltaDiffCenZ->Draw();
 518   myDetElemDeltaDiffCenZ->SetXTitle("#Delta[zc_{i}^{m}-zc_{i}^{th}] (mm)");
 519   myDetElemDeltaDiffCenZ->SetYTitle("DetElem arbitrary ordering");
 520   gDeltaDiffCenZDetElem0->SetMarkerStyle(20);
 521   gDeltaDiffCenZDetElem0->Draw("P");
 522   ch0Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+0][2],0.5,1e1*lDeltaDiffCenDetElem0[nDetElems+0][2],1.5);
 523   ch0Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+0][2],3.5,1e1*lDeltaDiffCenDetElem0[nDetElems+0][2],4.5);
 524   ch1Line->DrawLine(1e1*lDeltaDiffCenDetElem0[nDetElems+1][2],1.5,1e1*lDeltaDiffCenDetElem0[nDetElems+1][2],3.5);
 525
 526   pad->cd(4);
 527   myDetElemDeltaDiffRotX->Draw();
 528   myDetElemDeltaDiffRotX->SetXTitle("#Delta[#psi_{i}^{m}-#psi_{i}^{th}] (mrad)");
 529   myDetElemDeltaDiffRotX->SetYTitle("DetElem arbitrary ordering");
 530   gDeltaDiffPsiDetElem0->SetMarkerStyle(20);
 531   gDeltaDiffPsiDetElem0->Draw("P");
 532   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],1.5);
 533   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],3.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],4.5);
 534   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],1.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],3.5);
 535
 536   pad->cd(5);
 537   myDetElemDeltaDiffRotY->Draw();
 538   myDetElemDeltaDiffRotY->SetXTitle("#Delta[#theta_{i}^{m}-#theta_{i}^{th}] (mrad)");
 539   myDetElemDeltaDiffRotY->SetYTitle("DetElem arbitrary ordering");
 540   gDeltaDiffThtDetElem0->SetMarkerStyle(20);
 541   gDeltaDiffThtDetElem0->Draw("P");
 542   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],1.5);
 543   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],3.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],4.5);
 544   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],1.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],3.5);
 545
 546   pad->cd(6);
 547   myDetElemDeltaDiffRotZ->Draw();
 548   myDetElemDeltaDiffRotZ->SetXTitle("#Delta[#phi_{i}^{m}-#phi_{i}^{th}] (mrad)");
 549   myDetElemDeltaDiffRotZ->SetYTitle("DetElem arbitrary ordering");
 550   gDeltaDiffPhiDetElem0->SetMarkerStyle(20);
 551   gDeltaDiffPhiDetElem0->Draw("P");
 552   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],1.5);
 553   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],3.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],4.5);
 554   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],1.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],3.5);
 555
 556   pad->Update();
 557
 558   if(bMonitor){
 559     // MONITOR: Histograms for monitoring
 560     TH2F *hCPSTa = new TH2F("hCPSTa","hCPSTa",100,-200,200,100,-200,200);
 561     TH2F *hCPSTc = new TH2F("hCPSTc","hCPSTc",100,-200,200,100,-200,200);
 562     TH2F *hSSTa = new TH2F("hSSTa","hSSTa",100,-200,200,100,-200,200);
 563     TH2F *hSSTc = new TH2F("hSSTc","hSSTc",100,-200,200,100,-200,200);
 564     TH2F *hSSTap = new TH2F("hSSTap","hSSTap",800,-200,200,800,-200,200);
 565     TH2F *hSSTcp = new TH2F("hSSTcp","hSSTcp",800,-200,200,800,-200,200);
 566
 567     // MONITOR: Fill histograms with chambers and slats sticker target positions
 568     for (int iCh =0; iCh<=1; iCh++) {
 569       myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 570       cout << "Filling histograms of sticker target for chamber" << iCh+1 << " ..." << endl;
 571       myChamber->FillCPSTHistograms(Form("T4_6%d",2-iCh),hCPSTc);
 572       myChamber->FillDESTHistograms(TString("T4_46"),hSSTc);
 573
 574       for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 575         //      Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 576         myDetElem =  myChamber->GetDetElem(iDetElemI);
 577         // MONITOR: Fill slat plane for fit monitor.
 578         Double_t pl[3] = {0};
 579         Double_t pg[3] = {0};
 580         AliSurveyPoint *pointSBT0 = myDetElem->GetLButtonTarget(0);
 581         AliSurveyPoint *pointSBT1 = myDetElem->GetLButtonTarget(1);
 582         if(pointSBT0&&pointSBT1) {
 583           if (pointSBT0->GetX()>pointSBT1->GetX()){
 584             pointSBT0=myDetElem->GetLButtonTarget(1);
 585             pointSBT1=myDetElem->GetLButtonTarget(0);
 586           }
 587           Double_t lX = pointSBT0->GetX();
 588           while(lX<pointSBT1->GetX()){
 589             Double_t lY = pointSBT0->GetY()-20;
 590             while(lY<pointSBT0->GetY()+20){
 591               pl[0] = lX;  pl[1] = lY;  pl[2] = 0.;
 592               (TGeoCombiTrans((*(myDetElem->GetBaseTrf()))*(*(myDetElem->GetAlignTrf())))).LocalToMaster(pl,pg);
 593               if(myDetElem->GetGButtonTarget(0)->GetPointName().Contains("A")){
 594                 if (hSSTap->GetBinContent(hSSTap->FindBin(pg[0],pg[1]))==0)
 595                   hSSTap->Fill(pg[0],pg[1],-pg[2]);
 596               }
 597               else {
 598                 if (hSSTcp->GetBinContent(hSSTcp->FindBin(pg[0],pg[1]))==0)
 599                   hSSTcp->Fill(pg[0],pg[1],-pg[2]);
 600               }
 601               lY+=hSSTap->GetYaxis()->GetBinWidth(1);
 602             }
 603             lX+=hSSTap->GetXaxis()->GetBinWidth(1);
 604           }
 605         }
 606       }
 607     }
 608
 609     if( saveps ){
 610       ps->NewPage();
 611     }
 612
 613     // View from side A
 614     str = Form("Chamber %d - side A",chId+1);
 615     TCanvas *cvn0 = new TCanvas("cvn0",str,cWidth,cHeight);
 616     canvas = cvn0;
 617     canvas->Range(0,0,21,29);
 618
 619     title = Form(" Deformations of chamber %d - View from side A ",chId+1);
 620     TPaveLabel *theTitle0 = new TPaveLabel(3,27.0,18,28.5,title,"br");
 621     theTitle = theTitle0;
 622     theTitle->SetFillColor(18);
 623     theTitle->SetTextFont(32);
 624     theTitle->SetTextSize(0.4);
 625     theTitle->SetTextColor(1);
 626     theTitle->Draw();
 627
 628     TPad *pad0 = new TPad("pad0","pad0",0.01,0.01,0.98,0.91,0);
 629     pad = pad0;
 630     pad->Draw();
 631     pad->Divide(2,2);
 632
 633     pad->cd(1);
 634     hCPSTa->SetMinimum(950);
 635     hCPSTa->SetMaximum(975);
 636     hCPSTa->Draw("lego2z");
 637
 638     pad->cd(2);
 639     hSSTa->SetMinimum(950);
 640     hSSTa->SetMaximum(980);
 641     hSSTa->Draw("lego2z");
 642
 643     pad->cd(3);
 644
 645     pad->cd(4);
 646     hSSTap->SetMinimum(950);
 647     hSSTap->SetMaximum(980);
 648     hSSTap->Draw("surf2z");
 649
 650     pad->Update();
 651     if(saveps){
 652       ps->NewPage();
 653     }
 654
 655     // Inv Mass, Multiplicity
 656     str = Form("Chamber %d - side C",chId+1);
 657     TCanvas *cvn1 = new TCanvas("cvn1",str,cWidth,cHeight);
 658     canvas = cvn1;
 659     canvas->Range(0,0,21,29);
 660
 661     title = Form(" Deformations of chamber %d - View from side C ",chId+1);
 662     TPaveLabel *theTitle1 = new TPaveLabel(3,27.0,18,28.5,title,"br");
 663     theTitle = theTitle1;
 664     theTitle->SetFillColor(18);
 665     theTitle->SetTextFont(32);
 666     theTitle->SetTextSize(0.4);
 667     theTitle->SetTextColor(1);
 668     theTitle->Draw();
 669
 670     TPad *pad1 = new TPad("pad1","pad1",0.01,0.01,0.98,0.91,0);
 671     pad = pad1;
 672     pad->Draw();
 673     pad->Divide(2,2);
 674
 675     Double_t lMin, lMax;
 676
 677     pad->cd(1);
 678     lMin = hCPSTc->GetMinimum(0);
 679     hCPSTc->SetMinimum(TMath::Floor(lMin));
 680     lMax = hCPSTc->GetMaximum();
 681     hCPSTc->SetMaximum(TMath::Ceil(lMax));
 682     hCPSTc->Draw("lego2z");
 683
 684     pad->cd(2);
 685 //     lMin = hSSTc->GetMinimum(0);
 686 //     hSSTc->SetMinimum(TMath::Floor(lMin));
 687 //     lMax = hSSTc->GetMaximum();
 688 //     hSSTc->SetMaximum(TMath::Ceil(lMax));
 689 //    hSSTc->Draw("lego2z");
 690
 691
 692     pad->cd(3);
 693
 694     pad->cd(4);
 695     lMin = hSSTcp->GetMinimum(0);
 696     hSSTcp->SetMinimum(TMath::Floor(lMin));
 697     lMax = hSSTcp->GetMaximum();
 698     hSSTcp->SetMaximum(TMath::Ceil(lMax));
 699     hSSTcp->Draw("surf2z");
 700   }
 701
 702   pad->Update();
 703   if( saveps ){
 704     ps->Close();
 705   }
 706
 707   TFile *hFile = new TFile(Form("spCH%d_L2G_ISF.root",chId+1),"RECREATE");
 708   gDeltaDiffCenXDetElem0->Write("gDeltaDiffCenXDetElem0");
 709   gDeltaDiffCenYDetElem0->Write("gDeltaDiffCenYDetElem0");
 710   gDeltaDiffCenZDetElem0->Write("gDeltaDiffCenZDetElem0");
 711   gDeltaDiffPsiDetElem0->Write("gDeltaDiffPsiDetElem0");
 712   gDeltaDiffThtDetElem0->Write("gDeltaDiffThtDetElem0");
 713   gDeltaDiffPhiDetElem0->Write("gDeltaDiffPhiDetElem0");
 714   myDetElemDeltaDiffCenX->Write();
 715   myDetElemDeltaDiffCenY->Write();
 716   myDetElemDeltaDiffCenZ->Write();
 717   myDetElemDeltaDiffRotX->Write();
 718   myDetElemDeltaDiffRotY->Write();
 719   myDetElemDeltaDiffRotZ->Write();
 720   hFile->Close();
 721   hFile->Delete();
 722 }