Script supersceeded by AliForwarddNdetaTask.C and
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibGlobalMisalignment.cxx
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
17 ////////////////////////////////////////////////////////////////////////////
18 //                                                                        //
19 // AliTPCCalibGlobalMisalignment class                                        //
20 // The class calculates the space point distortions due to simple         // 
21 // misalignments like shifts in caresian coordinates or a rotation        //
22 // of the TPC read out planes (A and C side)                              //
23 // Optionaly possible to use it for visualization of the alignemnt form the Alignment OCDB //
24 // fUseGeoManager has to be set to kTRUE to enable this option
25 //                                                                        //
26 // date: 06/05/2010                                                       //
27 // Authors: Stefan Rossegger, Jim Thomas, Magnus Mager                    //
28 ////////////////////////////////////////////////////////////////////////////
29
30 #include "AliTPCCalibGlobalMisalignment.h"
31 #include "TMath.h"
32 #include "TGeoMatrix.h"
33 #include "AliTPCROC.h"
34 #include "AliTPCcalibDB.h"
35 #include "AliTPCParam.h"
36 #include <TGeoPhysicalNode.h>
37 //
38 #include "AliAlignObjParams.h"
39 #include "AliGeomManager.h"
40 #include "AliCDBManager.h"
41 #include "AliCDBEntry.h"
42
43 AliTPCCalibGlobalMisalignment::AliTPCCalibGlobalMisalignment()
44   : AliTPCCorrection("mialign","Misalignment"),
45     fXShift(0.),fYShift(0.),fZShift(0.),
46     fRotPhiA(0.),fRotPhiC(0.),
47     fdRPhiOffsetA(0.), 
48     fdRPhiOffsetC(0.), 
49     fQuadrantQ0(0),   //OROC medium pads -delta ly+ - ly - shift
50     fQuadrantRQ0(0),  //OROC medium pads -delta ly+ - ly - rotation 
51     fQuadrantQ1(0),   //OROC long   pads -delta ly+ - ly - shift
52     fQuadrantQ2(0),   //OROC long   pads -shift
53     fQuadrantRQ1(0),  //OROC long   pads -delta ly+ - ly - rotation
54     fQuadrantRQ2(0),  //OROC long   pads -rotation
55     fMatrixGlobal(0), // global Alignment common
56     fArraySector(0)   // fArraySector
57 {
58   //
59   // default constructor
60   //
61 }
62
63 AliTPCCalibGlobalMisalignment::~AliTPCCalibGlobalMisalignment() {
64   //
65   // destructor
66   //
67   delete fQuadrantQ0;   //OROC medium pads -delta ly+ - ly - shift
68   delete fQuadrantRQ0;  //OROC medium pads -delta ly+ - ly - rotation 
69   delete fQuadrantQ1;   //OROC long   pads -delta ly+ - ly - shift
70   delete fQuadrantQ2;   //OROC long   pads -shift
71   delete fQuadrantRQ1;  //OROC long   pads -delta ly+ - ly - rotation
72   delete fQuadrantRQ2;  //OROC long   pads -rotation
73   delete fMatrixGlobal; // global matrix
74   delete fArraySector;  // sector matrices
75 }
76  
77 void AliTPCCalibGlobalMisalignment::SetQuadranAlign(const TVectorD *quadrantQ0, const TVectorD *quadrantRQ0, const TVectorD *quadrantQ1,const TVectorD *quadrantRQ1,  const TVectorD *quadrantQ2,  const TVectorD *quadrantRQ2){
78   //
79   // Set quadrant alignment
80   // 6 vectors for 36 (super) sectors
81   //
82   if (quadrantQ0) fQuadrantQ0   = new TVectorD(*quadrantQ0);
83   if (quadrantRQ0) fQuadrantRQ0 = new TVectorD(*quadrantRQ0);
84   //
85   if (quadrantQ1) fQuadrantQ1   = new TVectorD(*quadrantQ1);
86   if (quadrantQ1) fQuadrantRQ1  = new TVectorD(*quadrantRQ1);
87   if (quadrantQ2) fQuadrantQ2   = new TVectorD(*quadrantQ2);
88   if (quadrantQ2) fQuadrantRQ2  = new TVectorD(*quadrantRQ2);
89 }
90
91
92
93 void AliTPCCalibGlobalMisalignment::SetAlignGlobal(const TGeoMatrix * matrixGlobal){
94   //
95   // Set global misalignment
96   // Object is OWNER 
97   // 
98   if (fMatrixGlobal) delete fMatrixGlobal;
99   fMatrixGlobal=0;
100   if (matrixGlobal) fMatrixGlobal = new TGeoHMatrix(*matrixGlobal);
101 }
102
103 void AliTPCCalibGlobalMisalignment::SetAlignSectors(const TObjArray *arraySector){
104   //
105   // Set misalignment TObjArray of TGeoMatrices  - for each sector
106   // Object is OWNER
107   // 
108   if (fArraySector) delete fArraySector;
109   fArraySector=0;
110   if (arraySector) fArraySector = (TObjArray*)arraySector->Clone();
111 }
112
113
114 //void AliTPCCalibGlobalMisalignment::Init() {
115 //  //
116 // // Initialization funtion
117 //  //
118
119 //  // nothing to be initialized, results of this calibration class will go to the global aligment structure
120
121 //}
122
123 //void AliTPCCalibGlobalMisalignment::Update(const TTimeStamp &/*timeStamp*/) {
124 //  //
125 //  // Update function 
126 //  //
127 //
128 //  // nothing to be updated, results of this calibration class will go to the global aligment structure
129 //
130 //}
131
132
133
134 void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
135   //
136   // Calculates the simple correction due to a shift (in x,y,z) or an rotation of the TPC (around z)
137   //  
138   static AliTPCROC *tpcRoc =AliTPCROC::Instance();  
139   Double_t xref  = ( tpcRoc->GetPadRowRadii(0,0)+tpcRoc->GetPadRowRadii(36,tpcRoc->GetNRows(36)-1))*0.5;
140     
141   Double_t r=0, phi=0;
142   r   = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] );
143   phi = TMath::ATan2(x[1],x[0]);
144   // Getsector number
145   Double_t sec=TMath::Nint(-0.5+(phi*9./TMath::Pi()));
146   if (sec<0) sec+=18;
147   Int_t isec = TMath::Nint(sec);    
148   if  (roc%36>=18) isec+=18;
149   //
150   // Get the point on the local coordiante frame
151   //
152   Double_t alpha=(sec+0.5)*TMath::Pi()/9;
153   Double_t pos[3]={0,0,x[2]};
154   pos[0]=  TMath::Cos(alpha)*x[0]+TMath::Sin(alpha)*x[1];
155   pos[1]= -TMath::Sin(alpha)*x[0]+TMath::Cos(alpha)*x[1];
156   if (pos[0]>tpcRoc->GetPadRowRadiiUp(0))  isec+=36;
157
158   //
159   // apply quadrant alignment if available - in local coordinate frame
160   //
161   //
162   Double_t posQG[3]={x[0],x[1],x[2]};
163   {
164     Double_t dly=0;
165     Bool_t isQ0 = TMath::Abs(pos[0]-161)<28;
166     Bool_t isQ1 = (pos[0]>189);
167     Double_t  sign = (pos[1]>0)? 1.: -1.;
168     if (isQ0){
169       if (fQuadrantQ0)  dly+=sign*(*fQuadrantQ0)[isec%36];  // shift in cm
170       if (fQuadrantRQ0) dly+=sign*(*fQuadrantRQ0)[isec%36]*(pos[0]-xref);      
171     }
172     if (isQ1){
173       if (fQuadrantQ1)  dly+=sign*(*fQuadrantQ1)[isec%36];  // shift in cm
174       if (fQuadrantRQ1) dly+=sign*(*fQuadrantRQ1)[isec%36]*(pos[0]-xref);      
175       if (fQuadrantQ2)  dly+=(*fQuadrantQ2)[isec%36];  // shift in cm
176       if (fQuadrantRQ2) dly+=(*fQuadrantRQ2)[isec%36]*(pos[0]-xref);      
177     }
178     // Tranform the corrected point to the global frame
179     posQG[0]=  TMath::Cos(alpha)*pos[0]-TMath::Sin(alpha)*(pos[1]+dly);
180     posQG[1]=  TMath::Sin(alpha)*pos[0]+TMath::Cos(alpha)*(pos[1]+dly);
181   }
182   //
183   // rotation of the read-out planes
184   if  (roc%36<18) // A side
185     phi += fRotPhiA;
186   else         // C side
187     phi += fRotPhiC;
188   
189   // Simply adding a constant dRPHi residual. PURELY FOR CALIBRATION PURPOSES
190   if  (roc%36<18) // A side
191     phi += fdRPhiOffsetA/r;
192   else         // C side
193     phi += fdRPhiOffsetC/r;
194   
195   dx[0] = r * TMath::Cos(phi) - x[0];
196   dx[1] = r * TMath::Sin(phi) - x[1]; 
197   dx[2] = 0.; 
198
199   // Simple shifts
200   dx[0] -= fXShift;
201   dx[1] -= fYShift;
202   dx[2] -= fZShift;
203   // quadrant shifts
204   dx[0] += (posQG[0]-x[0]);
205   dx[1] += (posQG[1]-x[1]);
206   //
207   //
208   if (fMatrixGlobal){
209     // apply global alignment matrix
210     Double_t ppos[3]={x[0],x[1],x[2]};
211     Double_t pposC[3]={x[0],x[1],x[2]};
212     fMatrixGlobal->LocalToMaster(ppos,pposC);
213     dx[0]+=pposC[0]-ppos[0];
214     dx[1]+=pposC[1]-ppos[1];
215     dx[2]+=pposC[2]-ppos[2];
216   }
217
218   if (fArraySector){
219     // apply global alignment matrix
220     TGeoMatrix  *mat = (TGeoMatrix*)fArraySector->At(isec);
221     if (mat){
222       Double_t ppos[3]={x[0],x[1],x[2]};
223       Double_t pposC[3]={x[0],x[1],x[2]};
224       mat->LocalToMaster(ppos,pposC);
225       dx[0]+=pposC[0]-ppos[0];
226       dx[1]+=pposC[1]-ppos[1];
227       dx[2]+=pposC[2]-ppos[2];
228     }
229   }
230 }
231
232 void AliTPCCalibGlobalMisalignment::Print(Option_t* /*option*/ ) const {
233   //
234   // Print function to check the settings 
235   //
236   printf("%s",GetTitle());  
237   printf(" - Trivial Misalignments for calibration purposes: \n");
238   printf(" - X-Shift: %1.3f cm, Y-Shift: %1.3f cm, Z-Shift: %1.3f cm \n",fXShift,fYShift,fZShift);
239   printf(" - Phi-Rotations: A side: %1.5f rad, C side: %1.5f rad\n",fRotPhiA,fRotPhiC);
240   printf(" - dRPhi offsets: A side: %1.5f cm, C side: %1.5f cm\n",fdRPhiOffsetA,fdRPhiOffsetC);
241  
242  
243 }
244
245 void AliTPCCalibGlobalMisalignment::AddAlign(const  AliTPCCalibGlobalMisalignment & add){
246   //
247   // Add the alignmnet to current object
248   //
249   fXShift+=add.fXShift;               // Shift in global X [cm]
250   fYShift+=add.fYShift;               // Shift in global Y [cm]
251   fZShift+=add.fZShift;               // Shift in global Z [cm]
252
253   fRotPhiA+=add.fRotPhiA;      // simple rotation of A side read-out plane around the Z axis [rad]
254   fRotPhiC+=add.fRotPhiC;      // simple rotation of C side read-out plane around the Z axis [rad]
255   fdRPhiOffsetA+=add.fdRPhiOffsetA;  // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes!
256   fdRPhiOffsetC+=add.fdRPhiOffsetC;  // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes!
257   //
258   // Quadrant alignment
259   //
260   if (add.fQuadrantQ0) {
261     if (fQuadrantQ0) fQuadrantQ0->Add(*(add.fQuadrantQ0));
262     if (!fQuadrantQ0) fQuadrantQ0 = (TVectorD*)(add.fQuadrantQ0->Clone());
263   }
264   if (add.fQuadrantRQ0) {
265     if (fQuadrantRQ0) fQuadrantRQ0->Add(*(add.fQuadrantRQ0));
266     if (!fQuadrantRQ0) fQuadrantRQ0 = (TVectorD*)(add.fQuadrantRQ0->Clone());
267   }
268   //
269   if (add.fQuadrantQ1) {
270     if (fQuadrantQ1) fQuadrantQ1->Add(*(add.fQuadrantQ1));
271     if (!fQuadrantQ1) fQuadrantQ1 = (TVectorD*)(add.fQuadrantQ1->Clone());
272   }
273   if (add.fQuadrantRQ1) {
274     if (fQuadrantRQ1) fQuadrantRQ1->Add(*(add.fQuadrantRQ1));
275     if (!fQuadrantRQ1) fQuadrantRQ1 = (TVectorD*)(add.fQuadrantRQ1->Clone());
276   }
277   //
278   if (add.fQuadrantQ2) {
279     if (fQuadrantQ2) fQuadrantQ2->Add(*(add.fQuadrantQ2));
280     if (!fQuadrantQ2) fQuadrantQ2 = (TVectorD*)(add.fQuadrantQ2->Clone());
281   }
282   if (add.fQuadrantRQ2) {
283     if (fQuadrantRQ2) fQuadrantRQ2->Add(*(add.fQuadrantRQ2));
284     if (!fQuadrantRQ2) fQuadrantRQ2 = (TVectorD*)(add.fQuadrantRQ2->Clone());
285   }
286   //
287   // Global alignment - use native ROOT representation
288   //
289   if (add.fMatrixGlobal){
290     if (!fMatrixGlobal)  fMatrixGlobal = new TGeoHMatrix(*(add.fMatrixGlobal)); 
291     if (fMatrixGlobal)   ((TGeoHMatrix*)fMatrixGlobal)->Multiply(add.fMatrixGlobal); 
292   }
293   if (add.fArraySector){
294     if (!fArraySector) {SetAlignSectors(add.fArraySector);
295     }else{
296       for (Int_t isec=0; isec<72; isec++){
297         TGeoHMatrix *mat0= (TGeoHMatrix*)fArraySector->At(isec);
298         TGeoHMatrix *mat1= (TGeoHMatrix*)add.fArraySector->At(isec);
299         if (mat0&&mat1) mat0->Multiply(mat1);
300       }
301     }
302   }
303 }
304
305
306 AliTPCCalibGlobalMisalignment *  AliTPCCalibGlobalMisalignment::CreateOCDBAlign(){
307   //
308   // Create  AliTPCCalibGlobalMisalignment from OCDB Alignment entry
309   // OCDB has to be initialized before in user code
310   // All storages (defualt and specific)  and run number 
311   //
312   AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
313   if (!entry){
314     printf("Missing alignmnet entry. OCDB not initialized?\n");
315     return 0;
316   }
317   TClonesArray * array = (TClonesArray*)entry->GetObject();
318   Int_t entries = array->GetEntries();
319   TGeoHMatrix matrixGlobal;
320   TObjArray *alignArrayOCDB= new TObjArray(73);  // sector misalignment + global misalignment
321   //                                            // global is number 72
322   //
323   { for (Int_t i=0;i<entries; i++){
324       //
325       //
326       TGeoHMatrix matrix;
327       AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
328       alignP->GetMatrix(matrix);
329       Int_t imod;
330       AliGeomManager::ELayerID ilayer;
331       alignP->GetVolUID(ilayer, imod);
332       if (ilayer==AliGeomManager::kInvalidLayer) {
333         alignArrayOCDB->AddAt(matrix.Clone(),72);
334         alignP->GetMatrix(matrixGlobal);
335       }else{
336         Int_t sector=imod;
337         if (ilayer==AliGeomManager::kTPC2) sector+=36;
338         alignArrayOCDB->AddAt(matrix.Clone(),sector);
339       }
340     }
341   }
342   AliTPCCalibGlobalMisalignment *align = new  AliTPCCalibGlobalMisalignment;
343   align->SetAlignGlobal(&matrixGlobal);
344   align->SetAlignSectors(alignArrayOCDB);
345   return align;
346 }
347
348
349 AliTPCCalibGlobalMisalignment *  AliTPCCalibGlobalMisalignment::CreateMeanAlign(const AliTPCCalibGlobalMisalignment *alignIn){
350   //
351   // Create new object, disantangle common mean alignmnet and sector alignment
352   //
353   // 1. Try to get mean alignment
354   // 2. Remove mean alignment from sector alignment
355   // 3. Create new object
356   //
357   TObjArray * array = alignIn->GetAlignSectors();
358   TObjArray * arrayNew = new TObjArray(72);
359   //
360   //Get mean transformation
361   TGeoHMatrix matrix;  
362   {for (Int_t isec=0; isec<72; isec++){
363       const TGeoMatrix* cmatrix=(TGeoMatrix*)array->At(isec);
364       if (!cmatrix) continue;
365       matrix.Multiply(cmatrix);
366     }}
367   TGeoHMatrix matrixMean(matrix);
368   matrixMean.SetDx(matrix.GetTranslation()[0]/72.);
369   matrixMean.SetDy(matrix.GetTranslation()[1]/72.);
370   matrixMean.SetDz(matrix.GetTranslation()[2]/72.);
371   Double_t rotation[12];
372   {for (Int_t i=0; i<12; i++) {
373       rotation[i]=1.0;
374       if (TMath::Abs(matrix.GetRotationMatrix()[i]-1.)>0.1){
375       rotation[i]=matrix.GetRotationMatrix()[i]/72.;
376       }
377     }}
378   matrixMean.SetRotation(rotation);
379   TGeoHMatrix matrixInv = matrixMean.Inverse();
380   //
381   {for (Int_t isec=0; isec<72; isec++){
382       TGeoHMatrix* amatrix=(TGeoHMatrix*)(array->At(isec)->Clone());
383       if (!amatrix) continue;
384       amatrix->Multiply(&matrixInv);
385       arrayNew->AddAt(amatrix,isec);
386     }}
387   if (alignIn->GetAlignGlobal()) matrixMean.Multiply((alignIn->GetAlignGlobal()));
388   AliTPCCalibGlobalMisalignment *alignOut = new  AliTPCCalibGlobalMisalignment;
389   alignOut->SetAlignGlobal(&matrixMean);  
390   alignOut->SetAlignSectors(arrayNew);  
391   /*
392     Checks transformation:   
393     AliTPCCalibGlobalMisalignment *  alignIn =  AliTPCCalibGlobalMisalignment::CreateOCDBAlign()
394     AliTPCCalibGlobalMisalignment * alignOut =  AliTPCCalibGlobalMisalignment::CreateMeanAlign(alignIn)
395     alignOutM= (AliTPCCalibGlobalMisalignment*)alignOut->Clone();
396     alignOutS= (AliTPCCalibGlobalMisalignment*)alignOut->Clone();
397     alignOutS->SetAlignGlobal(0);  
398     alignOutM->SetAlignSectors(0);  
399     //
400     AliTPCCorrection::AddVisualCorrection(alignOut,0);
401     AliTPCCorrection::AddVisualCorrection(alignOutM,1);
402     AliTPCCorrection::AddVisualCorrection(alignOutS,2);
403     AliTPCCorrection::AddVisualCorrection(alignIn,3);
404     
405     TF1 f0("f0","AliTPCCorrection::GetCorrSector(x,85,0.9,1,0)",0,18);
406     TF1 f1("f1","AliTPCCorrection::GetCorrSector(x,85,0.9,1,1)",0,18);
407     TF1 f2("f2","AliTPCCorrection::GetCorrSector(x,85,0.9,1,2)",0,18);
408     TF1 f3("f3","AliTPCCorrection::GetCorrSector(x,85,0.9,1,3)",0,18);
409     f0->SetLineColor(1);
410     f1->SetLineColor(2);
411     f2->SetLineColor(3);
412     f3->SetLineColor(4);
413     f0->Draw();
414     f1->Draw("same");
415     f2->Draw("same");
416     f3->Draw("same");
417
418     TF2 f2D("f2D","AliTPCCorrection::GetCorrSector(x,y,0.9,1,0)-AliTPCCorrection::GetCorrSector(x,y,0.9,1,3)",0,18,85,245);
419   */
420   return alignOut;
421 }
422
423
424 void AliTPCCalibGlobalMisalignment::DumpAlignment( AliTPCCalibGlobalMisalignment* align, TTreeSRedirector *pcstream, const char *name){
425   //
426   // Dump alignment per sector into tree
427   //
428   TObjArray * array = align->GetAlignSectors();
429   if (!array) return;
430   //
431   //Get mean transformation
432   TGeoHMatrix matrix;  
433   {for (Int_t isec=0; isec<72; isec++){
434       TGeoHMatrix* cmatrix=(TGeoHMatrix*)array->At(isec);
435       TGeoHMatrix* cmatrixDown=(TGeoHMatrix*)array->At(isec%36);
436       TGeoHMatrix* cmatrixUp=(TGeoHMatrix*)array->At(isec%36+36);
437       TGeoHMatrix diff(*cmatrixDown);
438       diff.Multiply(&(cmatrixUp->Inverse()));
439       (*pcstream)<<name<<
440         "isec="<<isec<<
441         "m0.="<<cmatrix<<
442         "diff.="<<&diff<<
443         "\n";
444     }
445   }
446   
447 }