MUON/MUONSurveyCh5.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 5
  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 MUONSurveyCh5() {
  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 = 4;
  99   Int_t nChs = 2;
 100   Int_t nDetElems = 18;
 101   Int_t nDetElemsI = 9;
 102   //  Int_t nDetElemsO = 9;
 103   Int_t iDetElemToDetElemId[18] = {514,506,516,508,500,510,502,512,504,513,503,511,501,509,517,507,515,505};
 104   Int_t iDetElemPseudoIdToDetElem[18] = {4,12,6,10,8,17,1,15,3,13,5,11,7,9,0,16,2,14};
 105   Int_t iDetElemsOfChamber[2][9] = {{0,16,2,14,4,12,6,10,8},
 106                                     {9,7,11,5,13,3,15,1,17}};
 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("../Reports/AliceSt3_TC5_1381b.txt");
 119   myChamberO = (AliMUONSurveyChamber*)myChamberArray->At(1);
 120   myChamberO->GetSurveyObj()->FillFromLocalFile("../Reports/AliceSt3_TC5_1381b.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(Form("%dC_1000",chId+1),4);
 131   myChamberI->AddStickerTargets(Form("%dA_702",chId+1),9);
 132   myChamberI->AddStickerTargets(Form("%dA_703",chId+1),9);
 133   myChamberI->AddStickerTargets(Form("%dC_700",chId+1),9);
 134   myChamberI->AddStickerTargets(Form("%dC_701",chId+1),9);
 135   myChamberI->AddLButtonTargets(lSO->GetData(),Form("%d_1000",chId+1),4);
 136   myChamberO->AddGButtonTargets(Form("%dC_1010",chId+1),4);
 137   myChamberO->AddStickerTargets(Form("%dA_700",chId+1),9);
 138   myChamberO->AddStickerTargets(Form("%dA_701",chId+1),9);
 139   myChamberO->AddStickerTargets(Form("%dC_702",chId+1),9);
 140   myChamberO->AddStickerTargets(Form("%dC_703",chId+1),9);
 141   myChamberO->AddLButtonTargets(lSO->GetData(),Form("%d_1010",chId+1),4);
 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         if(iDetElem<9) baseName = Form("%dA_50%d",chId+1,iDetElem+1);
 153         else baseName = Form("%dC_50%d",chId+1,iDetElem+1-9);
 154         myDetElem->AddStickerTargets(baseName,9);
 155         if(iDetElem<9) baseName = Form("%dA_60%d",chId+1,iDetElem+1);
 156         else baseName = Form("%dC_60%d",chId+1,iDetElem+1-9);
 157         myDetElem->AddGButtonTargets(baseName,2);
 158         if(iDetElemToDetElemId[iDetElem]-500<10) baseName = Form("%d_50%d",chId+1,iDetElemToDetElemId[iDetElem]-500);
 159         else baseName = Form("%d_5%d",chId+1,iDetElemToDetElemId[iDetElem]-500);
 160         myDetElem->AddLButtonTargets(lSO->GetData(),baseName,2);
 161       }
 162     }
 163   }
 164   printf("All targets added! \n");
 165
 166   Double_t **lCenDetElem = new Double_t*[nDetElems+nChs];
 167   Double_t **lRotDetElem = new Double_t*[nDetElems+nChs];
 168   Double_t **lDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
 169   Double_t **lDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
 170   Double_t **lDiffThCenDetElem0 = new Double_t*[nDetElems+nChs];
 171   Double_t **lDiffThRotDetElem0 = new Double_t*[nDetElems+nChs];
 172   Double_t **lDeltaDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
 173   Double_t **lDeltaDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
 174
 175   for (int iDetElem=0; iDetElem<nDetElems+nChs; iDetElem++) {
 176     lCenDetElem[iDetElem] = new Double_t[3];
 177     lRotDetElem[iDetElem] = new Double_t[3];
 178     lDiffCenDetElem0[iDetElem] = new Double_t[3];
 179     lDiffRotDetElem0[iDetElem] = new Double_t[3];
 180     lDiffThCenDetElem0[iDetElem] = new Double_t[3];
 181     lDiffThRotDetElem0[iDetElem] = new Double_t[3];
 182     lDeltaDiffCenDetElem0[iDetElem] = new Double_t[3];
 183     lDeltaDiffRotDetElem0[iDetElem] = new Double_t[3];
 184   }
 185
 186   TGeoCombiTrans dtrfDetElem[nDetElems+nChs];
 187   TGeoCombiTrans localTrfDetElem[nDetElems+nChs];
 188   TGeoCombiTrans localTrfThDetElem[nDetElems+nChs];
 189
 190   // Import TGeo geometry
 191   char* geoFilename = "geometry.root";
 192   if ( ! AliGeomManager::GetGeometry() ) {
 193     AliGeomManager::LoadGeometry(geoFilename);
 194     if (! AliGeomManager::GetGeometry() ) {
 195       printf("MUONSurveyCh%d: getting geometry from file %s failed\n", chId+1,geoFilename);
 196       return;
 197     }
 198     cout << "geometry imported" << endl;
 199   }
 200
 201   AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
 202   transform->LoadGeometryData();
 203
 204   TGeoCombiTrans trfThChamber;
 205   TGeoCombiTrans trfThDetElem;
 206
 207   for (int iCh =0; iCh<=1; iCh++) {
 208     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 209
 210     trfThChamber = TGeoCombiTrans(*transform->GetModuleTransformerByDEId(iDetElemToDetElemId[iCh*nDetElemsI])->GetTransformation());
 211     trfThChamber.SetTranslation(trfThChamber.GetTranslation()[0]*10,
 212                                 trfThChamber.GetTranslation()[1]*10,
 213                                 trfThChamber.GetTranslation()[2]*10);
 214     trfThChamber.Print();
 215     myChamber->SetLocalTransformation(new TGeoCombiTrans(trfThChamber),kTRUE);
 216
 217     // Set Chamber plane function
 218     cout << "Setting plane for Chamber" << iCh+1 << " ..." << endl;
 219     myChamber->SetPlane(Form("fChamber%d",iCh+1));
 220     myChamber->SetPlaneParameters(0.,0.,0.);
 221
 222     // Fit a plane to sticker targets
 223     Double_t lCChi2 = myChamber->FitPlane();
 224     cout << "... chi2 = " << lCChi2 << " ..." << endl;
 225
 226     // Get best transformation from fitting method
 227     // (local to global transformation)
 228     cout << "Trying fitting method for chamber " << iCh << endl;
 229     myChamber->SurveyToAlign();
 230     //    myChamber->SurveyToAlign(myChamber->GetPlane()->GetParameter(0),myChamber->GetPlane()->GetParameter(1),myChamber->GetPlane()->GetParError(0),myChamber->GetPlane()->GetParError(1));
 231
 232     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 233       myDetElem =  myChamber->GetDetElem(iDetElemI);
 234       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 235
 236       trfThDetElem.Clear();
 237       trfThDetElem = TGeoCombiTrans(*transform->GetDetElement(iDetElemToDetElemId[iDetElem])->GetLocalTransformation());
 238       trfThDetElem.SetTranslation(trfThDetElem.GetTranslation()[0]*10,
 239                                   trfThDetElem.GetTranslation()[1]*10,
 240                                   trfThDetElem.GetTranslation()[2]*10);
 241       trfThDetElem.Print();
 242
 243       myDetElem->SetLocalTransformation(new TGeoCombiTrans(trfThDetElem),kTRUE);
 244
 245       for (int iCor=0; iCor<3; iCor++){
 246         lCenDetElem[iDetElem][iCor] = 0;
 247         lRotDetElem[iDetElem][iCor] = 0;
 248       }
 249
 250       if (bMonitor){
 251         // MONITOR: Draw graph with sticker targets for plane fitting
 252         myDetElem->DrawSTargets();
 253         gPad->Update();
 254       }
 255
 256       // Get DetElem transformation.
 257       // Psi and Theta are obtained by fitting a plane to the sticker targets.
 258       // Then Xc, Yc, Zc and Phi are obtained by solving the equations to the ref.
 259       // syst. transformation of the button targets
 260
 261       // Set DetElem plane function
 262       cout << "Setting plane for DetElem" << iDetElem+1 << " ..." << endl;
 263       myDetElem->SetPlane(Form("fDetElem%d",iDetElem+1));
 264       myDetElem->SetPlaneParameters(0.,0.,30.)
 265 ;
 266       // Fit a plane to sticker targets
 267       Double_t lChi2 = myDetElem->FitPlane();
 268       cout << "... chi2 = " << lChi2 << " ..." << endl;
 269
 270       lRotDetElem[iDetElem][0]=TMath::ATan(myDetElem->GetPlane()->GetParameter(0));
 271       lRotDetElem[iDetElem][1]=TMath::ATan(myDetElem->GetPlane()->GetParameter(1));
 272
 273       // Calculate Mean transformation using previous plane fit
 274       // and pairs of button targets
 275       myDetElem->CalculateMeanTransf(lCenDetElem[iDetElem],lRotDetElem[iDetElem]);
 276
 277       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;
 278
 279
 280       // Get best transformation from fitting method
 281       // (local to global transformation)
 282       cout << "Trying fitting method for DetElemId " << iDetElemToDetElemId[iDetElem] << endl;
 283       //      myDetElem->SurveyToAlign();
 284       myDetElem->SurveyToAlign(lRotDetElem[iDetElem][0],lRotDetElem[iDetElem][1],myDetElem->GetPlane()->GetParError(0),myDetElem->GetPlane()->GetParError(1));
 285     }
 286   }
 287
 288   // Print found transformation of Detection Element (plane fit + loop)
 289   for (int iCh =0; iCh<=1; iCh++) {
 290     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 291       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 292       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;
 293     }
 294   }
 295
 296   // Print Theoretical transformation of Detection Element
 297   for (int iCh =0; iCh<=1; iCh++) {
 298     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 299     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 300       myChamber->GetDetElem(iDetElemI)->PrintLocalTrf();
 301     }
 302   }
 303
 304   // Print found delta transformation of Detection Element
 305   for (int iCh =0; iCh<=1; iCh++) {
 306     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 307     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 308       myChamber->GetDetElem(iDetElemI)->PrintAlignTrf();
 309     }
 310   }
 311
 312   //
 313   // Compare transformations to expected ones
 314   //
 315   Int_t iDetElemToPos[18] = {0, 16, 2, 14, 4, 12, 6, 10, 8, 9, 7, 11, 5, 13, 3, 15, 1, 17};
 316
 317   TGraphErrors *gDeltaDiffCenXDetElem0 = new TGraphErrors(nDetElems);
 318   TGraphErrors *gDeltaDiffCenYDetElem0 = new TGraphErrors(nDetElems);
 319   TGraphErrors *gDeltaDiffCenZDetElem0 = new TGraphErrors(nDetElems);
 320   TGraphErrors *gDeltaDiffPsiDetElem0 = new TGraphErrors(nDetElems);
 321   TGraphErrors *gDeltaDiffThtDetElem0 = new TGraphErrors(nDetElems);
 322   TGraphErrors *gDeltaDiffPhiDetElem0 = new TGraphErrors(nDetElems);
 323
 324   for (int iCh =0; iCh<=1; iCh++) {
 325     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 326     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 327       myChamber->GetDetElem(iDetElemI)->GetAlignTrf()->Print();
 328     }
 329   }
 330
 331   for (int iCh =0; iCh<=1; iCh++) {
 332     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 333     // Store delta transformations to use for CDB entry creation
 334     dtrfDetElem[nDetElems+iCh].Clear();
 335     dtrfDetElem[nDetElems+iCh] = *(myChamber->GetAlignTrf());
 336
 337     // Those are for checks and visualizations
 338     localTrfDetElem[nDetElems+iCh].Clear();
 339     localTrfDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf())*(*(myChamber->GetAlignTrf())));
 340     localTrfDetElem[nDetElems+iCh].Print();
 341     lDiffCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfDetElem[nDetElems+iCh].GetTranslation();
 342     AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffRotDetElem0[nDetElems+iCh]);
 343
 344     localTrfThDetElem[nDetElems+iCh].Clear();
 345     localTrfThDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf()));
 346     localTrfThDetElem[nDetElems+iCh].Print();
 347     lDiffThCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfThDetElem[nDetElems+iCh].GetTranslation();
 348     AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffThRotDetElem0[nDetElems+iCh]);
 349
 350     for (int iCor=0; iCor<3; iCor++){
 351       lDeltaDiffCenDetElem0[nDetElems+iCh][iCor] = lDiffCenDetElem0[nDetElems+iCh][iCor]-lDiffThCenDetElem0[nDetElems+iCh][iCor];
 352       lDeltaDiffRotDetElem0[nDetElems+iCh][iCor] = lDiffRotDetElem0[nDetElems+iCh][iCor]-lDiffThRotDetElem0[nDetElems+iCh][iCor];
 353       if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]-=TMath::TwoPi();
 354       if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]+=TMath::TwoPi();
 355     }
 356
 357     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 358       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 359       myDetElem =  myChamber->GetDetElem(iDetElemI);
 360       // Store delta transformations to use for CDB entry creation
 361       dtrfDetElem[iDetElem].Clear();
 362       dtrfDetElem[iDetElem] = *(myDetElem->GetAlignTrf());
 363
 364       // Those are for checks and visualizations
 365       localTrfDetElem[iDetElem].Clear();
 366       localTrfDetElem[iDetElem] = (*(myDetElem->GetLocalTrf())*(*(myDetElem->GetAlignTrf())));
 367       //      localTrfDetElem[iDetElem] = (*(myDetElem->GetBaseTrf())*(*(myDetElem->GetAlignTrf())));
 368       localTrfDetElem[iDetElem].Print();
 369       lDiffCenDetElem0[iDetElem] = (Double_t*)localTrfDetElem[iDetElem].GetTranslation();
 370       AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[iDetElem].GetRotationMatrix(),lDiffRotDetElem0[iDetElem]);
 371 //       lDiffCenDetElem0[iDetElem] = lCenDetElem[iDetElem];
 372 //       lDiffRotDetElem0[iDetElem] = lRotDetElem[iDetElem];
 373
 374       localTrfThDetElem[iDetElem].Clear();
 375       localTrfThDetElem[iDetElem] = (*(myDetElem->GetLocalTrf()));
 376       //      localTrfThDetElem[iDetElem] = (*(myDetElem->GetBaseTrf()));
 377       localTrfThDetElem[iDetElem].Print();
 378       lDiffThCenDetElem0[iDetElem] = (Double_t*)localTrfThDetElem[iDetElem].GetTranslation();
 379       AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[iDetElem].GetRotationMatrix(),lDiffThRotDetElem0[iDetElem]);
 380
 381       for (int iCor=0; iCor<3; iCor++){
 382         lDeltaDiffCenDetElem0[iDetElem][iCor] = lDiffCenDetElem0[iDetElem][iCor]-lDiffThCenDetElem0[iDetElem][iCor];
 383         lDeltaDiffRotDetElem0[iDetElem][iCor] = lDiffRotDetElem0[iDetElem][iCor]-lDiffThRotDetElem0[iDetElem][iCor];
 384         if (lDeltaDiffRotDetElem0[iDetElem][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]-=TMath::TwoPi();
 385         if (lDeltaDiffRotDetElem0[iDetElem][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]+=TMath::TwoPi();
 386       }
 387
 388       gDeltaDiffCenXDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
 389       gDeltaDiffCenYDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
 390       gDeltaDiffCenZDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
 391       gDeltaDiffPsiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
 392       gDeltaDiffThtDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
 393       gDeltaDiffPhiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
 394       gDeltaDiffCenXDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(0),0.);
 395       gDeltaDiffCenYDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(1),0.);
 396       gDeltaDiffCenZDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(2),0.);
 397       gDeltaDiffPsiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(3),0.);
 398       gDeltaDiffThtDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(4),0.);
 399       gDeltaDiffPhiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(5),0.);
 400     }
 401   }
 402
 403   // Apply the found misalignments to the geometry
 404   AliMUONGeometryTransformer *newTransform = AliMUONSurveyUtil::ReAlign(transform,chId,nDetElems,iDetElemPseudoIdToDetElem,dtrfDetElem,true);
 405   newTransform->WriteTransformations(Form("transform2ReAlignSurveyCh%d.dat",chId+1));
 406
 407   if(bOCDB){
 408     // Generate realigned data in local cdb
 409     const TClonesArray* array = newTransform->GetMisAlignmentData();
 410
 411     // Set the alignment resolution in the align objects for this chamber
 412     Double_t chResX = 0.1*myChamberI->GetAlignResX();
 413     Double_t chResY = 0.1*myChamberI->GetAlignResY();
 414     Double_t deResX = 0.1*(myChamberI->GetMeanDetElemAlignResX()+myChamberO->GetMeanDetElemAlignResX())/2.;
 415     Double_t deResY = 0.1*(myChamberI->GetMeanDetElemAlignResY()+myChamberO->GetMeanDetElemAlignResY())/2.;
 416     printf("Chamber alignment resolutions: resX=%f , resY=%f\n",chResX,chResY);
 417     printf("Detection Elements alignment resolutions: resX=%f , resY=%f\n",deResX,deResY);
 418     chResX = TMath::Sqrt(0.1*0.1+chResX*chResX);
 419     chResY = TMath::Sqrt(0.1*0.1+chResY*chResY);
 420     AliMUONSurveyUtil::SetAlignmentResolution(array,chId,chResX,chResY,deResX,deResY);
 421
 422     // CDB manager
 423     AliCDBManager* cdbManager = AliCDBManager::Instance();
 424     cdbManager->SetDefaultStorage(Form("local://ReAlignSurveyCh%dCDB",chId+1));
 425
 426     AliCDBMetaData* cdbData = new AliCDBMetaData();
 427     cdbData->SetResponsible("Dimuon Offline project");
 428     cdbData->SetComment("MUON alignment objects with survey realignment");
 429     AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
 430     cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
 431   }
 432
 433   for(Int_t iCor=0; iCor<3; iCor++){
 434     for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
 435       cout << lDeltaDiffCenDetElem0[iDetElem][iCor] << " " << lDiffCenDetElem0[iDetElem][iCor] << " " << lDiffThCenDetElem0[iDetElem][iCor] << endl;
 436     }
 437     cout << endl;
 438   }
 439   cout << endl;
 440   for(Int_t iCor=0; iCor<3; iCor++){
 441     for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
 442       cout << lDeltaDiffRotDetElem0[iDetElem][iCor] << " " << lDiffRotDetElem0[iDetElem][iCor] << " " << lDiffThRotDetElem0[iDetElem][iCor] << endl;
 443     }
 444     cout << endl;
 445   }
 446
 447   TH1F *myDetElemDeltaDiffCenX = new TH1F("myDetElemDeltaDiffCenX","myDetElemDeltaDiffCenX",100,-10,10);
 448   myDetElemDeltaDiffCenX->SetMaximum(nDetElems+1);
 449   myDetElemDeltaDiffCenX->SetMinimum(0);
 450   TH1F *myDetElemDeltaDiffCenY = new TH1F("myDetElemDeltaDiffCenY","myDetElemDeltaDiffCenY",100,-10,10);
 451   myDetElemDeltaDiffCenY->SetMaximum(nDetElems+1);
 452   myDetElemDeltaDiffCenY->SetMinimum(0);
 453   TH1F *myDetElemDeltaDiffCenZ = new TH1F("myDetElemDeltaDiffCenZ","myDetElemDeltaDiffCenZ",100,-40,40);
 454   myDetElemDeltaDiffCenZ->SetMaximum(nDetElems+1);
 455   myDetElemDeltaDiffCenZ->SetMinimum(0);
 456
 457   TH1F *myDetElemDeltaDiffRotX = new TH1F("myDetElemDeltaDiffRotX","myDetElemDeltaDiffRotX",100,-40,40);
 458   myDetElemDeltaDiffRotX->SetMaximum(nDetElems+1);
 459   myDetElemDeltaDiffRotX->SetMinimum(0);
 460   TH1F *myDetElemDeltaDiffRotY = new TH1F("myDetElemDeltaDiffRotY","myDetElemDeltaDiffRotY",100,-40,40);
 461   myDetElemDeltaDiffRotY->SetMaximum(nDetElems+1);
 462   myDetElemDeltaDiffRotY->SetMinimum(0);
 463   TH1F *myDetElemDeltaDiffRotZ = new TH1F("myDetElemDeltaDiffRotZ","myDetElemDeltaDiffRotZ",100,-10,10);
 464   myDetElemDeltaDiffRotZ->SetMaximum(nDetElems+1);
 465   myDetElemDeltaDiffRotZ->SetMinimum(0);
 466
 467   //
 468   // ******** Starting plots
 469   //
 470   TCanvas *canvas;
 471   TPad *pad;
 472   TPaveLabel *theTitle;
 473   gStyle->SetPalette(1);
 474
 475   TPostScript *ps = 0;
 476
 477   if( saveps ){
 478     ps = new TPostScript(Form("surveyChamber%d",chId+1),filetype);
 479     ps->NewPage();
 480   }
 481
 482   // Observed misalignments
 483   str = Form("Chamber %d",chId+1);
 484   TCanvas *cvn2 = new TCanvas("cvn2",str,cWidth,cHeight);
 485   canvas = cvn2;
 486   canvas->Range(0,0,21,29);
 487
 488   title = Form(" MisAlignments Chamber %d - PL2G - In Frame ",chId+1);
 489   TPaveLabel *theTitle2 = new TPaveLabel(3,27.0,18,28.5,title,"br");
 490   theTitle = theTitle2;
 491   theTitle->SetFillColor(18);
 492   theTitle->SetTextFont(32);
 493   theTitle->SetTextSize(0.4);
 494   theTitle->SetTextColor(1);
 495   theTitle->Draw();
 496
 497   TPad *pad2 = new TPad("pad2","pad2",0.01,0.01,0.98,0.91,0);
 498   pad = pad2;
 499   pad->Draw();
 500   TLine *ch0Line = new TLine(0,1,0,2);
 501   TLine *ch1Line = new TLine(0,1,0,2);
 502   ch1Line->SetLineStyle(2);
 503   pad->Divide(3,2);
 504
 505   pad->cd(1);
 506   myDetElemDeltaDiffCenX->Draw();
 507   myDetElemDeltaDiffCenX->SetXTitle("#Delta[xc_{i}^{m}-xc_{i}^{th}] (mm)");
 508   myDetElemDeltaDiffCenX->SetYTitle("DetElem arbitrary ordering");
 509   gDeltaDiffCenXDetElem0->SetMarkerStyle(20);
 510   gDeltaDiffCenXDetElem0->Draw("P");
 511   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][0],0.5,lDeltaDiffCenDetElem0[nDetElems+0][0],9.5);
 512   ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][0],9.5,lDeltaDiffCenDetElem0[nDetElems+1][0],18.5);
 513
 514   pad->cd(2);
 515   myDetElemDeltaDiffCenY->Draw();
 516   myDetElemDeltaDiffCenY->SetXTitle("#Delta[yc_{i}^{m}-yc_{i}^{th}] (mm)");
 517   myDetElemDeltaDiffCenY->SetYTitle("DetElem arbitrary ordering");
 518   gDeltaDiffCenYDetElem0->SetMarkerStyle(20);
 519   gDeltaDiffCenYDetElem0->Draw("P");
 520   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][1],0.5,lDeltaDiffCenDetElem0[nDetElems+0][1],9.5);
 521   ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][1],9.5,lDeltaDiffCenDetElem0[nDetElems+1][1],18.5);
 522
 523   pad->cd(3);
 524   myDetElemDeltaDiffCenZ->Draw();
 525   myDetElemDeltaDiffCenZ->SetXTitle("#Delta[zc_{i}^{m}-zc_{i}^{th}] (mm)");
 526   myDetElemDeltaDiffCenZ->SetYTitle("DetElem arbitrary ordering");
 527   gDeltaDiffCenZDetElem0->SetMarkerStyle(20);
 528   gDeltaDiffCenZDetElem0->Draw("P");
 529   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][2],0.5,lDeltaDiffCenDetElem0[nDetElems+0][2],9.5);
 530   ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][2],9.5,lDeltaDiffCenDetElem0[nDetElems+1][2],18.5);
 531
 532   pad->cd(4);
 533   myDetElemDeltaDiffRotX->Draw();
 534   myDetElemDeltaDiffRotX->SetXTitle("#Delta[#psi_{i}^{m}-#psi_{i}^{th}] (mrad)");
 535   myDetElemDeltaDiffRotX->SetYTitle("DetElem arbitrary ordering");
 536   gDeltaDiffPsiDetElem0->SetMarkerStyle(20);
 537   gDeltaDiffPsiDetElem0->Draw("P");
 538   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],9.5);
 539   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],9.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],18.5);
 540
 541   pad->cd(5);
 542   myDetElemDeltaDiffRotY->Draw();
 543   myDetElemDeltaDiffRotY->SetXTitle("#Delta[#theta_{i}^{m}-#theta_{i}^{th}] (mrad)");
 544   myDetElemDeltaDiffRotY->SetYTitle("DetElem arbitrary ordering");
 545   gDeltaDiffThtDetElem0->SetMarkerStyle(20);
 546   gDeltaDiffThtDetElem0->Draw("P");
 547   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],9.5);
 548   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],9.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],18.5);
 549
 550   pad->cd(6);
 551   myDetElemDeltaDiffRotZ->Draw();
 552   myDetElemDeltaDiffRotZ->SetXTitle("#Delta[#phi_{i}^{m}-#phi_{i}^{th}] (mrad)");
 553   myDetElemDeltaDiffRotZ->SetYTitle("DetElem arbitrary ordering");
 554   gDeltaDiffPhiDetElem0->SetMarkerStyle(20);
 555   gDeltaDiffPhiDetElem0->Draw("P");
 556   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],9.5);
 557   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],9.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],18.5);
 558
 559   pad->Update();
 560
 561   if(bMonitor){
 562     // MONITOR: Histograms for monitoring
 563     TH2F *hCPSTa = new TH2F("hCPSTa","hCPSTa",100,-2000,2000,100,-2000,2000);
 564     TH2F *hCPSTc = new TH2F("hCPSTc","hCPSTc",100,-2000,2000,100,-2000,2000);
 565     TH2F *hSSTa = new TH2F("hSSTa","hSSTa",100,-2000,2000,100,-2000,2000);
 566     TH2F *hSSTc = new TH2F("hSSTc","hSSTc",100,-2000,2000,100,-2000,2000);
 567     TH2F *hSSTap = new TH2F("hSSTap","hSSTap",800,-2000,2000,800,-2000,2000);
 568     TH2F *hSSTcp = new TH2F("hSSTcp","hSSTcp",800,-2000,2000,800,-2000,2000);
 569
 570     // MONITOR: Fill histograms with chambers and slats sticker target positions
 571     for (int iCh =0; iCh<=1; iCh++) {
 572       myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
 573       cout << "Filling histograms of sticker target for chamber" << iCh+1 << " ..." << endl;
 574       myChamber->FillCPSTHistograms(TString("C"),hCPSTc,TString("A"),hCPSTa);
 575       myChamber->FillDESTHistograms(TString("C"),hSSTc,TString("A"),hSSTa);
 576
 577       for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
 578         //      Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
 579         myDetElem =  myChamber->GetDetElem(iDetElemI);
 580         // MONITOR: Fill slat plane for fit monitor.
 581         Double_t pl[3] = {0};
 582         Double_t pg[3] = {0};
 583         AliSurveyPoint *pointSBT0 = myDetElem->GetLButtonTarget(0);
 584         AliSurveyPoint *pointSBT1 = myDetElem->GetLButtonTarget(1);
 585         if(pointSBT0&&pointSBT1) {
 586           if (pointSBT0->GetX()>pointSBT1->GetX()){
 587             pointSBT0=myDetElem->GetLButtonTarget(1);
 588             pointSBT1=myDetElem->GetLButtonTarget(0);
 589           }
 590           Double_t lX = pointSBT0->GetX();
 591           while(lX<pointSBT1->GetX()){
 592             Double_t lY = pointSBT0->GetY()-200;
 593             while(lY<pointSBT0->GetY()+200){
 594               pl[0] = lX;  pl[1] = lY;  pl[2] = 0.;
 595               (TGeoCombiTrans((*(myDetElem->GetBaseTrf()))*(*(myDetElem->GetAlignTrf())))).LocalToMaster(pl,pg);
 596               if(myDetElem->GetGButtonTarget(0)->GetPointName().Contains("A")){
 597                 hSSTap->Fill(pg[0],pg[1],-pg[2]);
 598               }
 599               else {
 600                 hSSTcp->Fill(pg[0],pg[1],-pg[2]);
 601               }
 602               lY+=hSSTap->GetYaxis()->GetBinWidth(1);
 603             }
 604             lX+=hSSTap->GetXaxis()->GetBinWidth(1);
 605           }
 606         }
 607       }
 608     }
 609
 610     if( saveps ){
 611       ps->NewPage();
 612     }
 613
 614     // View from side A
 615     str = Form("Chamber %d - side A",chId+1);
 616     TCanvas *cvn0 = new TCanvas("cvn0",str,cWidth,cHeight);
 617     canvas = cvn0;
 618     canvas->Range(0,0,21,29);
 619
 620     title = Form(" Deformations of chamber %d - View from side A ",chId+1);
 621     TPaveLabel *theTitle0 = new TPaveLabel(3,27.0,18,28.5,title,"br");
 622     theTitle = theTitle0;
 623     theTitle->SetFillColor(18);
 624     theTitle->SetTextFont(32);
 625     theTitle->SetTextSize(0.4);
 626     theTitle->SetTextColor(1);
 627     theTitle->Draw();
 628
 629     TPad *pad0 = new TPad("pad0","pad0",0.01,0.01,0.98,0.91,0);
 630     pad = pad0;
 631     pad->Draw();
 632     pad->Divide(2,2);
 633
 634     pad->cd(1);
 635     hCPSTa->SetMinimum(9500);
 636     hCPSTa->SetMaximum(9750);
 637     hCPSTa->Draw("lego2z");
 638
 639     pad->cd(2);
 640     hSSTa->SetMinimum(9500);
 641     hSSTa->SetMaximum(9800);
 642     hSSTa->Draw("lego2z");
 643
 644     pad->cd(3);
 645
 646     pad->cd(4);
 647     hSSTap->SetMinimum(9500);
 648     hSSTap->SetMaximum(9800);
 649     hSSTap->Draw("surf2z");
 650
 651     pad->Update();
 652     if(saveps){
 653       ps->NewPage();
 654     }
 655
 656     // Inv Mass, Multiplicity
 657     str = Form("Chamber %d - side C",chId+1);
 658     TCanvas *cvn1 = new TCanvas("cvn1",str,cWidth,cHeight);
 659     canvas = cvn1;
 660     canvas->Range(0,0,21,29);
 661
 662     title = Form(" Deformations of chamber %d - View from side C ",chId+1);
 663     TPaveLabel *theTitle1 = new TPaveLabel(3,27.0,18,28.5,title,"br");
 664     theTitle = theTitle1;
 665     theTitle->SetFillColor(18);
 666     theTitle->SetTextFont(32);
 667     theTitle->SetTextSize(0.4);
 668     theTitle->SetTextColor(1);
 669     theTitle->Draw();
 670
 671     TPad *pad1 = new TPad("pad1","pad1",0.01,0.01,0.98,0.91,0);
 672     pad = pad1;
 673     pad->Draw();
 674     pad->Divide(2,2);
 675
 676     pad->cd(1);
 677     gStyle->SetPalette(1);
 678     hCPSTc->SetMinimum(9550);
 679     hCPSTc->SetMaximum(9850);
 680     hCPSTc->Draw("lego2z");
 681
 682     pad->cd(2);
 683     gStyle->SetPalette(1);
 684     hSSTc->SetMinimum(9600);
 685     hSSTc->SetMaximum(9850);
 686     hSSTc->Draw("lego2z");
 687
 688     pad->cd(3);
 689     gStyle->SetPalette(1);
 690
 691     pad->cd(4);
 692     gStyle->SetPalette(1);
 693     hSSTcp->SetMinimum(9600);
 694     hSSTcp->SetMaximum(9850);
 695     hSSTcp->Draw("surf2z");
 696   }
 697
 698   pad->Update();
 699   if( saveps ){
 700     ps->Close();
 701   }
 702
 703   TFile *hFile = new TFile(Form("spCH%d_PL2G_IF.root",chId+1),"RECREATE");
 704   gDeltaDiffCenXDetElem0->Write("gDeltaDiffCenXDetElem0");
 705   gDeltaDiffCenYDetElem0->Write("gDeltaDiffCenYDetElem0");
 706   gDeltaDiffCenZDetElem0->Write("gDeltaDiffCenZDetElem0");
 707   gDeltaDiffPsiDetElem0->Write("gDeltaDiffPsiDetElem0");
 708   gDeltaDiffThtDetElem0->Write("gDeltaDiffThtDetElem0");
 709   gDeltaDiffPhiDetElem0->Write("gDeltaDiffPhiDetElem0");
 710   myDetElemDeltaDiffCenX->Write();
 711   myDetElemDeltaDiffCenY->Write();
 712   myDetElemDeltaDiffCenZ->Write();
 713   myDetElemDeltaDiffRotX->Write();
 714   myDetElemDeltaDiffRotY->Write();
 715   myDetElemDeltaDiffRotZ->Write();
 716   hFile->Close();
 717   hFile->Delete();
 718 }