StRoot: macros/calib/alignTPC.C Source File

00001 void alignTPC(){
00002 
00003   gSystem->Load("libPhysics");
00004   gStyle->SetFrameBorderMode(0);
00005   gStyle->SetPadBorderMode(0);
00006 
00007   TVector3 EAO[25], EBO[25], EDO[25]; 
00008   TVector3 WAO[25], WBO[25], WDO[25]; 
00009 
00010   TVector3 tempVect(0,0,0);
00011   TVector3 crossVect(0,0,0);
00012   TVector3 ortVect(0,0,0);
00013   TVector3 uniVect(0,0,0);
00014   TVector3 nullVect(0,0,0);
00015 
00016   TVector3 axisSurveyX(1.,0.,0.);
00017   TVector3 axisSurveyY(0.,1.,0.);
00018   TVector3 axisSurveyZ(0.,0.,1.);
00019 
00020   double tAngle =0.;
00021 
00022   EAO[13].SetXYZ(-0.27465, -2.31703, 0.47066);
00023   EAO[14].SetXYZ(-0.47431, -2.31705, 0.27091);
00024   EAO[15].SetXYZ(-0.54726, -2.31713, -0.00182);
00025   EAO[16].SetXYZ(-0.47394, -2.31718, -0.27460);
00026   EAO[17].SetXYZ(-0.27426, -2.31711, -0.47420);
00027   EAO[18].SetXYZ(-0.00134, -2.31699, -0.54718);
00028   EAO[19].SetXYZ(0.27141, -2.31675, -0.47399);
00029   EAO[20].SetXYZ(0.47104, -2.31663, -0.27427);
00030   EAO[21].SetXYZ(0.54403, -2.31668, -0.00137);
00031   EAO[22].SetXYZ(0.47082, -2.31674, 0.27141);
00032   EAO[23].SetXYZ(0.27101, -2.31683, 0.47091);
00033   EAO[24].SetXYZ(-0.00182, -2.31695, 0.54387);
00034 
00035   EBO[13].SetXYZ( -0.59588, -2.31725, 1.02658);
00036   EBO[14].SetXYZ(-1.03053, -2.31719, 0.59175);
00037   EBO[15].SetXYZ(-1.18926, -2.31731, -0.00211);
00038   EBO[16].SetXYZ(-1.03017, -2.31692, -0.59588);
00039   EBO[17].SetXYZ(-0.59496, -2.31720, -1.03037);
00040   EBO[18].SetXYZ(-0.00110, -2.31688, -1.18911);
00041   EBO[19].SetXYZ(0., 0., 0.);
00042   EBO[20].SetXYZ(1.02724, -2.31631, -0.59508);
00043   EBO[21].SetXYZ(1.18602, -2.31632, -0.00115);
00044   EBO[22].SetXYZ(1.02673, -2.31650, 0.59253);
00045   EBO[23].SetXYZ(0.59184, -2.31677, 1.02697);
00046   EBO[24].SetXYZ(-0.00208, -2.31696, 1.18566);
00047 
00048   EDO[13].SetXYZ(-0.95525, -2.31701, 1.64841);
00049   EDO[14].SetXYZ(0., 0., 0.);
00050   EDO[15].SetXYZ(0., 0., 0.);
00051   EDO[16].SetXYZ(-1.65192, -2.31734, -0.95513);
00052   EDO[17].SetXYZ(-0.95400, -2.31674, -1.65246);
00053   EDO[18].SetXYZ(0., 0., 0.);
00054   EDO[19].SetXYZ(0.95194, -2.31608, -1.65172);
00055   EDO[20].SetXYZ(1.64927, -2.31585, -0.95392);
00056   EDO[21].SetXYZ(0., 0., 0.);
00057   EDO[22].SetXYZ(0., 0., 0.);
00058   EDO[23].SetXYZ(0.95070, -2.31622, 1.64905);
00059   EDO[24].SetXYZ(-0.00234, -2.31674, 1.90405);
00060 
00061   WAO[1].SetXYZ(0., 0., 0.);
00062   WAO[2].SetXYZ(0.46923, 2.31308, 0.27104);
00063   WAO[3].SetXYZ(0.54225, 2.31310, -0.00181);
00064   WAO[4].SetXYZ(0.46903, 2.31310, -0.27453);
00065   WAO[5].SetXYZ(0.26930, 2.31289, -0.47419);
00066   WAO[6].SetXYZ(-0.00342, 2.31277, -0.54732);
00067   WAO[7].SetXYZ(-0.27622, 2.31275, -0.47430);
00068   WAO[8].SetXYZ(-0.47585, 2.31267, -0.27457);
00069   WAO[9].SetXYZ(-0.54889, 2.31278, -0.00186);
00070   WAO[10].SetXYZ(-0.47583, 2.31282, 0.27090);
00071   WAO[11].SetXYZ(-0.27613, 2.31298, 0.47050);
00072   WAO[12].SetXYZ(-0.00335, 2.31298, 0.54363);
00073 
00074   WBO[1].SetXYZ(0.59054, 2.31298, 1.02659);
00075   WBO[2].SetXYZ(1.02529, 2.31314, 0.59191);
00076   WBO[3].SetXYZ(1.18429, 2.31314, -0.00186);
00077   WBO[4].SetXYZ(1.02510, 2.31299, -0.59558);
00078   WBO[5].SetXYZ(0.59031, 2.31274, -1.03022);
00079   WBO[6].SetXYZ(-0.00346, 2.31251, -1.18937);
00080   WBO[7].SetXYZ(-0.59727, 2.31251, -1.03034);
00081   WBO[8].SetXYZ(-1.03189, 2.31226, -0.59552);
00082   WBO[9].SetXYZ(-1.19092, 2.31227, -0.00178);
00083   WBO[10].SetXYZ(-1.03191, 2.31243, 0.59190);
00084   WBO[11].SetXYZ(-0.59707, 2.31270, 1.02648);
00085   WBO[12].SetXYZ(-0.00330, 2.31279, 1.18557);
00086 
00087   WDO[1].SetXYZ(0.94962, 2.31289, 1.64844);
00088   WDO[2].SetXYZ(1.64729, 2.31318, 0.95101);
00089   WDO[3].SetXYZ(1.90240, 2.31331, -0.00176);
00090   WDO[4].SetXYZ(0., 0., 0.);
00091   WDO[5].SetXYZ(0.94945, 2.31280, -1.65202);
00092   WDO[6].SetXYZ(-0.00344, 2.31223, -1.90736);
00093   WDO[7].SetXYZ(-0.95633, 2.31209, -1.65220);
00094   WDO[8].SetXYZ(0., 0., 0.);
00095   WDO[9].SetXYZ(0., 0., 0.);
00096   WDO[10].SetXYZ(-1.65380, 2.31169, 0.95098);
00097   WDO[11].SetXYZ(0., 0., 0.);
00098   WDO[12].SetXYZ(0., 0., 0.);
00099 
00100   TH3F *tpcCenterInSurvey = new TH3F("tpcCenterInSurvey","tpcCenterInSurvey",101,-10.1,10.1,101,-10.1,10.1,101,-10.1,10.1);
00101   TH1F *tpcCenterInSurveyX = new TH1F("tpcCenterInSurveyX","tpcCenterInSurveyX",201,-10.05,10.05);
00102   tpcCenterInSurveyX->SetXTitle("X_{survey} (mm)");
00103   TH1F *tpcCenterInSurveyY = new TH1F("tpcCenterInSurveyY","tpcCenterInSurveyY",201,-10.05,10.05);
00104   tpcCenterInSurveyY->SetXTitle("Y_{survey} (mm)");
00105   TH1F *tpcCenterInSurveyZ = new TH1F("tpcCenterInSurveyZ","tpcCenterInSurveyZ",201,-10.05,10.05);
00106   tpcCenterInSurveyZ->SetXTitle("Z_{survey} (mm)");
00107 
00108   TH3F *tpcWTiltInSurvey = new TH3F("tpcWTiltInSurvey","tpcWTiltInSurvey",101,-0.00101,0.00101,101,-0.00101,0.00101,101,-0.00101,0.00101);
00109   TH1F *tpcWTiltInSurveyX = new TH1F("tpcWTiltInSurveyX","tpcWTiltInSurveyX",201,-0.001005,0.001005);
00110   TH1F *tpcWTiltInSurveyY = new TH1F("tpcWTiltInSurveyY","tpcWTiltInSurveyY",201,-0.001005,0.001005);
00111   TH1F *tpcWTiltInSurveyZ = new TH1F("tpcWTiltInSurveyZ","tpcWTiltInSurveyZ",201,-0.001005,0.001005);
00112 
00113   TH3F *tpcETiltInSurvey = new TH3F("tpcETiltInSurvey","tpcETiltInSurvey",101,-1010,1010,101,-1010,1010,101,-1010,1010);
00114   TH1F *tpcETiltInSurveyX = new TH1F("tpcETiltInSurveyX","tpcETiltInSurveyX",201,-1005,1005);
00115   TH1F *tpcETiltInSurveyY = new TH1F("tpcETiltInSurveyY","tpcETiltInSurveyY",201,-1005,1005);
00116   TH1F *tpcETiltInSurveyZ = new TH1F("tpcETiltInSurveyZ","tpcETiltInSurveyZ",201,-1005,1005);
00117 
00118   TH1F *tpcWXInSurveyXY = new TH1F("tpcWXInSurveyXY","tpcWXInSurveyXY",201,-1.005,1.005);
00119   tpcWXInSurveyXY->SetXTitle("Angle (mrad)");
00120   TH1F *tpcWXInSurveyXZ = new TH1F("tpcWXInSurveyXZ","tpcWXInSurveyXZ",201,-1.005,1.005);
00121   tpcWXInSurveyXZ->SetXTitle("Angle (mrad)");
00122 
00123   TH1F *tpcWYInSurveyXY = new TH1F("tpcWYInSurveyXY","tpcWYInSurveyXY",201,-1.005,1.005);
00124   tpcWYInSurveyXY->SetXTitle("Angle (mrad)");
00125   TH1F *tpcWYInSurveyYZ = new TH1F("tpcWYInSurveyYZ","tpcWYInSurveyYZ",201,-1.005,1.005);
00126   tpcWYInSurveyYZ->SetXTitle("Angle (mrad)");
00127 
00128   TH1F *tpcWZInSurveyXZ = new TH1F("tpcWZInSurveyXZ","tpcWZInSurveyXZ",201,-1.005,1.005);
00129   TH1F *tpcWZInSurveyYZ = new TH1F("tpcWZInSurveyYZ","tpcWZInSurveyYZ",201,-1.005,1.005);
00130 
00131   TH1F *tpcEXInSurveyXY = new TH1F("tpcEXInSurveyXY","tpcEXInSurveyXY",201,-1.005,1.005);
00132   TH1F *tpcEXInSurveyXZ = new TH1F("tpcEXInSurveyXZ","tpcEXInSurveyXZ",201,-1.005,1.005);
00133 
00134   TH1F *tpcEYInSurveyXY = new TH1F("tpcEYInSurveyXY","tpcEYInSurveyXY",201,-1.005,1.005);
00135   tpcEYInSurveyXY->SetXTitle("Angle (mrad)");
00136   TH1F *tpcEYInSurveyYZ = new TH1F("tpcEYInSurveyYZ","tpcEYInSurveyYZ",201,-1.005,1.005);
00137   tpcEYInSurveyYZ->SetXTitle("Angle (mrad)");
00138 
00139   TH1F *tpcEZInSurveyXZ = new TH1F("tpcEZInSurveyXZ","tpcEZInSurveyXZ",201,-1.005,1.005);
00140   TH1F *tpcEZInSurveyYZ = new TH1F("tpcEZInSurveyYZ","tpcEZInSurveyYZ",201,-1.005,1.005);
00141 
00142 
00143 
00144   //
00145   // To calculate the center of the TPC (perfect cylinder, 12 perfectly symetric sectors at each end)
00146   // Calculate the midle of a line going from one end to the other through the center of the TPC
00147   // W1-E17, W2-E16, W3-E15, W4-E14, W5-E13
00148   // W6-E24, W7-E23, W8-E22, W9-E21, W10-E20, W11-E19, W12-E18
00149   for(int iSec=2; iSec<=5; iSec++){ // No WAO[1] data
00150     tempVect.SetXYZ(0.,0.,0.);  
00151     tempVect += WAO[iSec]; tempVect += EAO[18-iSec]; tempVect *= 0.5; tempVect *= 1000;
00152     tpcCenterInSurvey->Fill(tempVect.X(),tempVect.Y(),tempVect.Z());
00153     tpcCenterInSurveyX->Fill(tempVect.X());
00154     tpcCenterInSurveyY->Fill(tempVect.Y());
00155     tpcCenterInSurveyZ->Fill(tempVect.Z());
00156     printf("%f\n",tempVect.X());
00157   }
00158   for(int iSec=6; iSec<=12; iSec++){
00159     tempVect.SetXYZ(0.,0.,0.);  
00160     tempVect += WAO[iSec]; tempVect += EAO[30-iSec]; tempVect *= 0.5; tempVect *= 1000;
00161     tpcCenterInSurvey->Fill(tempVect.X(),tempVect.Y(),tempVect.Z());
00162     tpcCenterInSurveyX->Fill(tempVect.X());
00163     tpcCenterInSurveyY->Fill(tempVect.Y());
00164     tpcCenterInSurveyZ->Fill(tempVect.Z());
00165     printf("%f\n",tempVect.X());
00166   }
00167 
00168   for(int iSec=1; iSec<=5; iSec++){ 
00169     tempVect.SetXYZ(0.,0.,0.);  
00170     tempVect += WBO[iSec]; tempVect += EBO[18-iSec]; tempVect *= 0.5; tempVect *= 1000;
00171     tpcCenterInSurvey->Fill(tempVect.X(),tempVect.Y(),tempVect.Z());
00172     tpcCenterInSurveyX->Fill(tempVect.X());
00173     tpcCenterInSurveyY->Fill(tempVect.Y());
00174     tpcCenterInSurveyZ->Fill(tempVect.Z());
00175     printf("%f\n",tempVect.X());
00176   }
00177   for(int iSec=6; iSec<=12; iSec++){
00178     if(iSec!=11){ // No EBO[19] data
00179       tempVect.SetXYZ(0.,0.,0.);  
00180       tempVect += WBO[iSec]; tempVect += EBO[30-iSec]; tempVect *= 0.5; tempVect *= 1000;
00181       tpcCenterInSurvey->Fill(tempVect.X(),tempVect.Y(),tempVect.Z());
00182       tpcCenterInSurveyX->Fill(tempVect.X());
00183       tpcCenterInSurveyY->Fill(tempVect.Y());
00184       tpcCenterInSurveyZ->Fill(tempVect.Z());
00185       printf("%f\n",tempVect.X());
00186     }
00187   }
00188 
00189   for(int iSec=1; iSec<=5; iSec++){ 
00190     if(iSec!=3&&iSec!=4){ // No EDO[19] data
00191       tempVect.SetXYZ(0.,0.,0.);  
00192       tempVect += WDO[iSec]; tempVect += EDO[18-iSec]; tempVect *= 0.5; tempVect *= 1000;
00193       tpcCenterInSurvey->Fill(tempVect.X(),tempVect.Y(),tempVect.Z());
00194       tpcCenterInSurveyX->Fill(tempVect.X());
00195       tpcCenterInSurveyY->Fill(tempVect.Y());
00196       tpcCenterInSurveyZ->Fill(tempVect.Z());
00197       printf("%f\n",tempVect.X());
00198     }
00199   }
00200   for(int iSec=6; iSec<=12; iSec++){
00201     if(iSec!=8&&iSec!=9&&iSec!=11&&iSec!=12){ // No EDO[19] data
00202       tempVect.SetXYZ(0.,0.,0.);  
00203       tempVect += WDO[iSec]; tempVect += EDO[30-iSec]; tempVect *= 0.5; tempVect *= 1000;
00204       tpcCenterInSurvey->Fill(tempVect.X(),tempVect.Y(),tempVect.Z());
00205       tpcCenterInSurveyX->Fill(tempVect.X());
00206       tpcCenterInSurveyY->Fill(tempVect.Y());
00207       tpcCenterInSurveyZ->Fill(tempVect.Z());
00208       printf("%f\n",tempVect.X());
00209     }
00210   }
00211 
00212   //
00213   // To calculate the angle TPC with respect to the Surveynet,
00214   // Calculate the normal vector of a line in the endcap plane
00215   // 
00216   // A targets
00217   for(int iSec=1; iSec<=10; iSec++){
00218     if(WAO[iSec]!=nullVect){
00219       for(int jSec=iSec+1; jSec<=11; jSec++){
00220         if(WAO[jSec]!=nullVect){
00221           tempVect.SetXYZ(0.,0.,0.);  
00222           tempVect += WAO[iSec];  tempVect -= WAO[jSec];
00223           for(int kSec=jSec+1; kSec<=12; kSec++){
00224             if(WAO[kSec]!=nullVect){    
00225               crossVect.SetXYZ(0.,0.,0.);  
00226               crossVect += WAO[jSec];  crossVect -= WAO[kSec];
00227               crossVect = crossVect.Cross(tempVect);
00228               uniVect = crossVect.Unit();
00229               //printf("(%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00230 
00231               tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00232               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00233               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00234               tAngle *= 1e3;
00235               if (tAngle>10000) printf("WAO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00236               //printf("In XY plane, Angle to Y %f\n",tAngle);
00237               tpcWYInSurveyXY->Fill(tAngle);
00238               tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00239               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00240               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00241               tAngle *= 1e3;
00242               if (tAngle>10000) printf("WAO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00243               //printf("In YZ plane, Angle to Y %f\n",tAngle);
00244               tpcWYInSurveyYZ->Fill(tAngle);
00245             }
00246           }
00247         }
00248       }
00249     }
00250   }
00251 
00252   for(int iSec=1; iSec<=9; iSec++){
00253     if(WAO[iSec]!=nullVect){
00254       for(int jSec=iSec+1; jSec<=10; jSec++){
00255         if(WAO[jSec]!=nullVect){
00256           tempVect.SetXYZ(0.,0.,0.);  
00257           tempVect += WAO[iSec];  tempVect -= WAO[jSec];
00258           for(int kSec=jSec+1; kSec<=11; kSec++){
00259             if(WAO[kSec]!=nullVect){    
00260               for(int lSec=kSec+1; lSec<=12; lSec++){
00261                 if(WAO[lSec]!=nullVect){        
00262                   crossVect.SetXYZ(0.,0.,0.);  
00263                   crossVect += WAO[kSec];  crossVect -= WAO[lSec];
00264                   crossVect = crossVect.Cross(tempVect);
00265                   uniVect = crossVect.Unit();
00266                   
00267                   tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00268                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00269                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00270                   tAngle *= 1e3;
00271                   if (tAngle>10000) printf("WAO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00272                   //printf("In XY plane, Angle to Y %f\n",tAngle);
00273                   tpcWYInSurveyXY->Fill(tAngle);
00274                   tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00275                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00276                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00277                   tAngle *= 1e3;
00278                   if (tAngle>10000) printf("WAO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00279                   //printf("In YZ plane, Angle to Y %f\n",tAngle);
00280                   tpcWYInSurveyYZ->Fill(tAngle);
00281                 }
00282               }
00283             }
00284           }
00285         }
00286       }
00287     }
00288   }
00289   
00290   // B targets
00291   for(int iSec=1; iSec<=10; iSec++){
00292     if(WBO[iSec]!=nullVect){
00293       for(int jSec=iSec+1; jSec<=11; jSec++){
00294         if(WBO[jSec]!=nullVect){
00295           tempVect.SetXYZ(0.,0.,0.);  
00296           tempVect += WBO[iSec];  tempVect -= WBO[jSec];
00297           for(int kSec=jSec+1; kSec<=12; kSec++){
00298             if(WBO[kSec]!=nullVect){    
00299               crossVect.SetXYZ(0.,0.,0.); 
00300               crossVect += WBO[jSec];  crossVect -= WBO[kSec];
00301               crossVect = crossVect.Cross(tempVect);
00302               uniVect = crossVect.Unit();
00303 
00304               tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00305               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00306               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00307               tAngle *= 1e3;
00308               if (tAngle>10000) printf("WBO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00309               //printf("In XY plane, Angle to Y %f\n",tAngle);
00310               tpcWYInSurveyXY->Fill(tAngle);
00311               tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00312               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00313               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00314               tAngle *= 1e3;
00315               if (tAngle>10000) printf("WBO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00316               //printf("In YZ plane, Angle to Y %f\n",tAngle);
00317               tpcWYInSurveyYZ->Fill(tAngle);
00318             }
00319           }
00320         }
00321       }
00322     }
00323   }
00324 
00325   for(int iSec=1; iSec<=9; iSec++){
00326     if(WBO[iSec]!=nullVect){
00327       for(int jSec=iSec+1; jSec<=10; jSec++){
00328         if(WBO[jSec]!=nullVect){
00329           tempVect.SetXYZ(0.,0.,0.);  
00330           tempVect += WBO[iSec];  tempVect -= WBO[jSec];
00331           for(int kSec=jSec+1; kSec<=11; kSec++){
00332             if(WBO[kSec]!=nullVect){    
00333               for(int lSec=kSec+1; lSec<=12; lSec++){
00334                 if(WBO[lSec]!=nullVect){        
00335                   crossVect.SetXYZ(0.,0.,0.);  
00336                   crossVect += WBO[kSec];  crossVect -= WBO[lSec];
00337                   crossVect = crossVect.Cross(tempVect);
00338                   uniVect = crossVect.Unit();
00339 
00340                   tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00341                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00342                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00343                   tAngle *= 1e3;
00344                   if (tAngle>10000) printf("WBO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00345                   //printf("In XY plane, Angle to Y %f\n",tAngle);
00346                   tpcWYInSurveyXY->Fill(tAngle);
00347                   tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00348                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00349                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00350                   tAngle *= 1e3;
00351                   if (tAngle>10000) printf("WBO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00352                   //printf("In YZ plane, Angle to Y %f\n",tAngle);
00353                   tpcWYInSurveyYZ->Fill(tAngle);
00354                 }
00355               }
00356             }
00357           }
00358         }
00359       }
00360     }
00361   }
00362 
00363   // D targets
00364   for(int iSec=1; iSec<=10; iSec++){
00365     if(WDO[iSec]!=nullVect){
00366       for(int jSec=iSec+1; jSec<=11; jSec++){
00367         if(WDO[jSec]!=nullVect){
00368           tempVect.SetXYZ(0.,0.,0.);  
00369           tempVect += WDO[iSec];  tempVect -= WDO[jSec];
00370           for(int kSec=jSec+1; kSec<=12; kSec++){
00371             if(WDO[kSec]!=nullVect){    
00372               crossVect.SetXYZ(0.,0.,0.);  
00373               crossVect += WDO[jSec];  crossVect -= WDO[kSec];
00374               crossVect = crossVect.Cross(tempVect);
00375               uniVect = crossVect.Unit();
00376 
00377               tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00378               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00379               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00380               tAngle *= 1e3;
00381               if (tAngle>10000) printf("WDO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00382               //printf("In XY plane, Angle to Y %f\n",tAngle);
00383               tpcWYInSurveyXY->Fill(tAngle);
00384               tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00385               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00386               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00387               tAngle *= 1e3;
00388               if (tAngle>10000) printf("WDO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00389               //printf("In YZ plane, Angle to Y %f\n",tAngle);
00390               tpcWYInSurveyYZ->Fill(tAngle);
00391             }
00392           }
00393         }
00394       }
00395     }
00396   }
00397 
00398   for(int iSec=1; iSec<=9; iSec++){
00399     if(WDO[iSec]!=nullVect){
00400       for(int jSec=iSec+1; jSec<=10; jSec++){
00401         if(WDO[jSec]!=nullVect){
00402           tempVect.SetXYZ(0.,0.,0.);  
00403           tempVect += WDO[iSec];  tempVect -= WDO[jSec];
00404           for(int kSec=jSec+1; kSec<=11; kSec++){
00405             if(WDO[kSec]!=nullVect){    
00406               for(int lSec=kSec+1; lSec<=12; lSec++){
00407                 if(WDO[lSec]!=nullVect){        
00408                   crossVect.SetXYZ(0.,0.,0.);  
00409                   crossVect += WDO[kSec];  crossVect -= WDO[lSec];
00410                   crossVect = crossVect.Cross(tempVect);
00411                   uniVect = crossVect.Unit();
00412                   
00413                   tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00414                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00415                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00416                   tAngle *= 1e3;
00417                   if (tAngle>10000) printf("WDO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00418                   //printf("In XY plane, Angle to Y %f\n",tAngle);
00419                   tpcWYInSurveyXY->Fill(tAngle);
00420                   tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00421                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00422                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00423                   tAngle *= 1e3;
00424                   if (tAngle>10000) printf("WDO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00425                   //printf("In YZ plane, Angle to Y %f\n",tAngle);
00426                   tpcWYInSurveyYZ->Fill(tAngle);
00427                 }
00428               }
00429             }
00430           }
00431         }
00432       }
00433     }
00434   }
00435 
00436 
00437   // East
00438   // A targets
00439   for(int iSec=13; iSec<=22; iSec++){
00440     if(EAO[iSec]!=nullVect){
00441       for(int jSec=iSec+1; jSec<=23; jSec++){
00442         if(EAO[jSec]!=nullVect){
00443           tempVect.SetXYZ(0.,0.,0.);  
00444           tempVect += EAO[iSec];  tempVect -= EAO[jSec];
00445           for(int kSec=jSec+1; kSec<=24; kSec++){
00446             if(EAO[kSec]!=nullVect){    
00447               crossVect.SetXYZ(0.,0.,0.);  
00448               crossVect += EAO[jSec];  crossVect -= EAO[kSec];
00449               crossVect = crossVect.Cross(tempVect);
00450               uniVect = crossVect.Unit();
00451 
00452               tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00453               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00454               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00455               tAngle *= 1e3;
00456               if (tAngle>10000) printf("EAO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00457               //printf("In XY plane, Angle to Y %f\n",tAngle);
00458               tpcEYInSurveyXY->Fill(tAngle);
00459               tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00460               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00461               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00462               tAngle *= 1e3;
00463               if (tAngle>10000) printf("EAO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00464               //printf("In YZ plane, Angle to Y %f\n",tAngle);
00465               tpcEYInSurveyYZ->Fill(tAngle);
00466             }
00467           }
00468         }
00469       }
00470     }
00471   }
00472 
00473   for(int iSec=13; iSec<=21; iSec++){
00474     if(EAO[iSec]!=nullVect){
00475       for(int jSec=iSec+1; jSec<=22; jSec++){
00476         if(EAO[jSec]!=nullVect){
00477           tempVect.SetXYZ(0.,0.,0.);  
00478           tempVect += EAO[iSec];  tempVect -= EAO[jSec];
00479           for(int kSec=jSec+1; kSec<=23; kSec++){
00480             if(EAO[kSec]!=nullVect){    
00481               for(int lSec=kSec+1; lSec<=24; lSec++){
00482                 if(EAO[lSec]!=nullVect){        
00483                   crossVect.SetXYZ(0.,0.,0.);  
00484                   crossVect += EAO[kSec];  crossVect -= EAO[lSec];
00485                   crossVect = crossVect.Cross(tempVect);
00486                   uniVect = crossVect.Unit();
00487                   
00488                   tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00489                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00490                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00491                   tAngle *= 1e3;
00492                   if (tAngle>10000) printf("EAO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00493                   //printf("In XY plane, Angle to Y %f\n",tAngle);
00494                   tpcEYInSurveyXY->Fill(tAngle);
00495                   tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00496                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00497                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00498                   tAngle *= 1e3;
00499                   if (tAngle>10000) printf("EAO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00500                   //printf("In YZ plane, Angle to Y %f\n",tAngle);
00501                   tpcEYInSurveyYZ->Fill(tAngle);
00502                 }
00503               }
00504             }
00505           }
00506         }
00507       }
00508     }
00509   }
00510 
00511   // B targets
00512   for(int iSec=13; iSec<=22; iSec++){
00513     if(EBO[iSec]!=nullVect){
00514       for(int jSec=iSec+1; jSec<=23; jSec++){
00515         if(EBO[jSec]!=nullVect){
00516           tempVect.SetXYZ(0.,0.,0.);  
00517           tempVect += EBO[iSec];  tempVect -= EBO[jSec];
00518           for(int kSec=jSec+1; kSec<=24; kSec++){
00519             if(EBO[kSec]!=nullVect){    
00520               crossVect.SetXYZ(0.,0.,0.);  
00521               crossVect += EBO[jSec];  crossVect -= EBO[kSec];
00522               crossVect = crossVect.Cross(tempVect);
00523               uniVect = crossVect.Unit();
00524 
00525               tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00526               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00527               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00528               tAngle *= 1e3;
00529               if (tAngle>10000) printf("EBO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00530               //printf("In XY plane, Angle to Y %f\n",tAngle);
00531               tpcEYInSurveyXY->Fill(tAngle);
00532               tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00533               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00534               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00535               tAngle *= 1e3;
00536               if (tAngle>10000) printf("EBO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00537               //printf("In YZ plane, Angle to Y %f\n",tAngle);
00538               tpcEYInSurveyYZ->Fill(tAngle);
00539             }
00540           }
00541         }
00542       }
00543     }
00544   }
00545 
00546   for(int iSec=13; iSec<=21; iSec++){
00547     if(EBO[iSec]!=nullVect){
00548       for(int jSec=iSec+1; jSec<=22; jSec++){
00549         if(EBO[jSec]!=nullVect){
00550           tempVect.SetXYZ(0.,0.,0.);  
00551           tempVect += EBO[iSec];  tempVect -= EBO[jSec];
00552           for(int kSec=jSec+1; kSec<=23; kSec++){
00553             if(EBO[kSec]!=nullVect){    
00554               for(int lSec=kSec+1; lSec<=24; lSec++){
00555                 if(EBO[lSec]!=nullVect){        
00556                   crossVect.SetXYZ(0.,0.,0.);  
00557                   crossVect += EBO[kSec];  crossVect -= EBO[lSec];
00558                   crossVect = crossVect.Cross(tempVect);
00559                   uniVect = crossVect.Unit();
00560                   
00561                   tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00562                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00563                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00564                   tAngle *= 1e3;
00565                   if (tAngle>10000) printf("EBO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00566                   //printf("In XY plane, Angle to Y %f\n",tAngle);
00567                   tpcEYInSurveyXY->Fill(tAngle);
00568                   tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00569                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00570                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00571                   tAngle *= 1e3;
00572                   if (tAngle>10000) printf("EBO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00573                   //printf("In YZ plane, Angle to Y %f\n",tAngle);
00574                   tpcEYInSurveyYZ->Fill(tAngle);
00575                 }
00576               }
00577             }
00578           }
00579         }
00580       }
00581     }
00582   }
00583 
00584   // D targets
00585   for(int iSec=13; iSec<=22; iSec++){
00586     if(EDO[iSec]!=nullVect){
00587       for(int jSec=iSec+1; jSec<=23; jSec++){
00588         if(EDO[jSec]!=nullVect){
00589           tempVect.SetXYZ(0.,0.,0.);  
00590           tempVect += EDO[iSec];  tempVect -= EDO[jSec];
00591           for(int kSec=jSec+1; kSec<=24; kSec++){
00592             if(EDO[kSec]!=nullVect){    
00593               crossVect.SetXYZ(0.,0.,0.);  
00594               crossVect += EDO[jSec];  crossVect -= EDO[kSec];
00595               crossVect = crossVect.Cross(tempVect);
00596               uniVect = crossVect.Unit();
00597 
00598               tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00599               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00600               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00601               tAngle *= 1e3;
00602               if (tAngle>10000) printf("EDO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00603               //printf("In XY plane, Angle to Y %f\n",tAngle);
00604               tpcEYInSurveyXY->Fill(tAngle);
00605               tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00606               if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00607               if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00608               tAngle *= 1e3;
00609               if (tAngle>10000) printf("EDO (%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00610               //printf("In YZ plane, Angle to Y %f\n",tAngle);
00611               tpcEYInSurveyYZ->Fill(tAngle);
00612             }
00613           }
00614         }
00615       }
00616     }
00617   }
00618 
00619   for(int iSec=13; iSec<=21; iSec++){
00620     if(EDO[iSec]!=nullVect){
00621       for(int jSec=iSec+1; jSec<=22; jSec++){
00622         if(EDO[jSec]!=nullVect){
00623           tempVect.SetXYZ(0.,0.,0.);  
00624           tempVect += EDO[iSec];  tempVect -= EDO[jSec];
00625           for(int kSec=jSec+1; kSec<=23; kSec++){
00626             if(EDO[kSec]!=nullVect){    
00627               for(int lSec=kSec+1; lSec<=24; lSec++){
00628                 if(EDO[lSec]!=nullVect){        
00629                   crossVect.SetXYZ(0.,0.,0.);  
00630                   crossVect += EDO[kSec];  crossVect -= EDO[lSec];
00631                   crossVect = crossVect.Cross(tempVect);
00632                   uniVect = crossVect.Unit();
00633                   
00634                   tAngle = TMath::ATan2(uniVect.X(),uniVect.Y());
00635                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00636                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00637                   tAngle *= 1e3;
00638                   if (tAngle>10000) printf("EDO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00639                   //printf("In XY plane, Angle to Y %f\n",tAngle);
00640                   tpcEYInSurveyXY->Fill(tAngle);
00641                   tAngle = TMath::ATan2(uniVect.Z(),uniVect.Y());
00642                   if (tAngle>TMath::Pi()/2.) tAngle -= TMath::Pi();
00643                   if (tAngle<-TMath::Pi()/2.) tAngle += TMath::Pi();
00644                   tAngle *= 1e3;
00645                   if (tAngle>10000) printf("EDO (%d,%d,%d,%d) (%f,%f,%f)\n",iSec,jSec,kSec,lSec,uniVect.X(),uniVect.Y(),uniVect.Z());
00646                   //printf("In YZ plane, Angle to Y %f\n",tAngle);
00647                   tpcEYInSurveyYZ->Fill(tAngle);
00648                 }
00649               }
00650             }
00651           }
00652         }
00653       }
00654     }
00655   }
00656 
00657 
00658 
00659 /*   for(int iSec=13; iSec<=15; iSec++){ */
00660 /*     tempVect.SetXYZ(0.,0.,0.);   */
00661 /*     tempVect += EAO[iSec]; tempVect -= EAO[iSec+6]; */
00662 /*     crossVect.SetXYZ(0.,0.,0.);   */
00663 /*     crossVect += EAO[iSec+3]; crossVect -= EAO[iSec+3+6]; */
00664 /*     crossVect = crossVect.Cross(tempVect); */
00665     
00666 /*     uniVect = crossVect.Unit(); */
00667 /*     uniVect.SetX(0.); */
00668 /*     tAngle = uniVect.Angle(axisSurveyY); */
00669 /*     tAngle *= 1e3; */
00670 /*     printf("In YZ plane, Angle to Y %f\n",tAngle); */
00671 /*     uniVect = crossVect.Unit(); */
00672 /*     uniVect.SetZ(0.); */
00673 /*     tAngle = uniVect.Angle(axisSurveyY); */
00674 /*     tAngle *= 1e3; */
00675 /*     printf("In XY plane, Angle to Y %f\n",tAngle); */
00676 
00677 /*   } */
00678 
00679 }