StRoot: StLaserAnalysisMaker/CheckLasers.C Source File

00001 #include "LaserBeams.h"
00002 #include "LASERINO.h"
00003 void CheckLasers() {
00004   cout << "NoLocals " << NoLocals << "\tNoBeams " << NoBeams << endl;
00005   const Int_t numSectors = 24;
00006   Double_t beta = 0;
00007   Double_t dBeta = 720./numSectors;
00008   Int_t sector = 0;
00009   Int_t half   = 0;
00010   TGeoHMatrix TpcHalf[2];
00011   Double_t rotHalfs[18] = {
00012     0, -1, 0, -1, 0, 0, 0, 0, -1,// sector  1-12
00013     0,  1, 0, -1, 0, 0, 0, 0,  1 // sector 13-24
00014   };
00015   for (half = 0; half <2; half++) TpcHalf[half].SetRotation(&rotHalfs[9*half]);
00016   TGeoHMatrix RotSec[24];
00017   for (sector = 1; sector <= numSectors; sector++) {
00018     if (sector > 12) beta = (numSectors-sector)*dBeta;
00019     else             beta = sector*dBeta;
00020     RotSec[sector-1].RotateZ(-beta);
00021   }
00022   
00023   for (Int_t raft = 1; raft <= 14; raft++) {
00024   //  for (Int_t raft = 9; raft <= 9; raft++) {
00025     TRVector mX;
00026     TRMatrix A(0,6);
00027     for (Int_t bundle = 1; bundle <= 6; bundle++) {
00028       for (Int_t mirror = 1; mirror <= 7; mirror++) {
00029         Int_t RBM = 100*raft + 10*bundle + mirror;
00030         if (RBM == 121 || // +1
00031             RBM == 122 || // +1
00032             RBM == 124 || // +1
00033             RBM == 126 || // +1
00034             RBM == 131 || // +1
00035             RBM == 132 )  // +1
00036         continue;
00037         if (RBM == 211 || // +1
00038             RBM == 221 || // +1
00039             RBM == 225 || // +1
00040             RBM == 231 || // +1
00041             RBM == 244 || // +1
00042             RBM == 246 || // +1
00043             RBM == 265)   // +1
00044         continue;
00045         if (RBM == 311 || // +1
00046             RBM == 331 || // +1
00047 //          RBM == 335 || // -1 
00048             RBM == 337 )  // +1
00049             continue;
00050         if (raft ==  5) continue;
00051         Int_t l = -1; 
00052         for (Int_t local = 0; local < NoLocals; local++) {
00053           if (Locals[local].Raft == raft && Locals[local].Bundle == bundle && Locals[local].Mirror == mirror) {l = local; break;}
00054         }
00055         Int_t k = -1;
00056         for (Int_t global = 0; global < NoBeams; global++) {
00057           if (LaserBeams[global].Raft == raft && LaserBeams[global].Bundle == bundle && LaserBeams[global].Mirror == mirror) {k = global; break;}
00058         } 
00059         if (l < 0 || k < 0) continue;
00060         //      cout << "Found match for sector " << LaserBeams[k].Sector << " raft " << raft << " bundle " << bundle << " mirror " << mirror << endl;
00061         Double_t *xg = &LaserBeams[k].X;// TRVector XG(3,xg); cout << " XG " << XG << endl;
00062         Double_t xG[3];
00063         sector = LaserBeams[k].Sector;
00064         half = 0;
00065         if (sector > 12) half = 1;
00066         TGeoHMatrix RotM = RotSec[sector-1] * TpcHalf[half];
00067         //      TGeoHMatrix RotM = TpcHalf[half] *  RotSec[sector-1];
00068         //      cout << " Raft " << raft << " Sector " << sector << " Bundle " << bundle << " Mirror " << mirror << endl; 
00069         //      RotM.Print();               //TRVector xgT(3,xg); cout << "xg " << xgT;
00070         RotM.MasterToLocal(xg,xG);  //TRVector xGT(3,xG); cout << "xG " << xGT;
00071         Double_t *xl = &Locals[l].X;
00072         Double_t xL[3] = {0.1*xl[0], 0.1*xl[1], 0.1*xl[2]};// TRVector XL(3,xL); cout << "XL " << XL;
00073         /*
00074                 (     1 -gamma  beta )
00075           DRT = ( gamma      1 -alpha); 
00076                 ( -beta  alpha     1 )
00077         */
00078         Double_t a[18] = {
00079               0.,  xL[2], -xL[1],  1., 0., 0.,
00080           -xL[2],     0.,  xL[0],  0., 1., 0.,
00081            xL[1], -xL[0],     0.,  0., 0., 1.
00082         };
00083         for (Int_t i = 0; i < 3; i++) {
00084           Double_t x = xG[i] - xL[i];
00085           mX.AddRow(&x);
00086           A.AddRow(&a[6*i]);
00087         }
00088       }
00089     }
00090     //    cout << "X "; mX.Print();
00091     //    cout << "A "; A.Print();
00092     if (A.GetNrows() < 1) {
00093       cout << "Raft " << raft << " is empty" << endl; 
00094       continue;
00095     } 
00096     cout << "Raft " << raft << " Sector = " << sector << "   =========================" << endl;  
00097     TRVector AmX(A,TRArray::kATxB,mX);      // cout << "AmX\t" << AmX << endl; 
00098     TRSymMatrix S(A,TRArray::kATxA);        // cout << "S\t" << S << endl;
00099     TRSymMatrix SInv(S,TRArray::kInverted); // cout << "SInv " << SInv << endl;
00100     TRVector  X(SInv,TRArray::kSxA,AmX);    // cout << "X " << X << endl;
00101     TGeoHMatrix T(Form("Raft%0i",raft));
00102     // normolize
00103     Double_t *xx = X.GetArray();
00104     T.RotateX(xx[0]*180/TMath::Pi());
00105     T.RotateY(xx[1]*180/TMath::Pi());
00106     T.RotateZ(xx[2]*180/TMath::Pi());
00107     T.SetTranslation(&xx[3]);
00108     cout << Form("{%i,%i",raft,sector);
00109     for (Int_t i = 0; i < 6; i++) cout << Form(",%15.5f",xx[i]);
00110     cout << "}," << endl;
00111     // cout << "Determinant " << T.Determinant() << endl;
00112     T.Print();
00113     TGeoHMatrix RotM = RotSec[sector-1] * TpcHalf[half] * T;
00114     //    TGeoHMatrix RotM = TpcHalf[half] * RotSec[sector-1] *  T;
00115     RotM.Print();
00116     // Check table
00117     for (Int_t bundle = 1; bundle <= 6; bundle++) {
00118       for (Int_t mirror = 1; mirror <= 7; mirror++) {
00119         Int_t l = -1; 
00120         for (Int_t local = 0; local < NoLocals; local++) {
00121           if (Locals[local].Raft == raft && Locals[local].Bundle == bundle && Locals[local].Mirror == mirror) {l = local; break;}
00122         }
00123         Int_t k = -1;
00124         for (Int_t global = 0; global < NoBeams; global++) {
00125           if (LaserBeams[global].Raft == raft && LaserBeams[global].Bundle == bundle && LaserBeams[global].Mirror == mirror) {k = global; break;}
00126         } 
00127         if (l < 0 || k < 0) continue;
00128         sector = LaserBeams[k].Sector;
00129         //      cout << "Found match for raft " << raft << " bundle " << bundle << " mirror " << mirror << endl;
00130         Double_t *xg = &LaserBeams[k].X; // TRVector XG(3,xg); cout << " XG " << XG;
00131         Double_t *xl = &Locals[l].X;     
00132         Double_t xL[3] = {0.1*xl[0],0.1*xl[1],0.1*xl[2]}; // TRVector XL(3,xL); cout << " xL " << XL;
00133         Double_t xG[3];
00134         
00135         RotM.LocalToMaster(xL,xG);     //  TRVector XGG(3,xG); cout << " xG " << XGG;
00136         Double_t dev = 0;
00137         for (Int_t i = 0; i < 3; i++) {
00138           Double_t dif = (xG[i] - xg[i]);
00139           dev += dif*dif;
00140         }
00141         dev = TMath::Sqrt(dev);
00142         if (dev > 0.0000) {
00143           cout << Form("S%02iR%02iB%iM%i dev = %f",sector,raft,bundle,mirror,dev);
00144           cout << Form(" dX/dY/dZ = %f/%f/%f",xG[0] - xg[0],xG[1] - xg[1],xG[2] - xg[2]);
00145           if (dev > 0.005) cout << " ===========";
00146           cout << endl;
00147         }
00148       }
00149     }
00150   }
00151 }