00001
00002 #if !defined(__CINT__) || defined(__MAKECINT__)
00003 #include "Riostream.h"
00004 #include <stdio.h>
00005 #include "TSystem.h"
00006 #include "TMath.h"
00007 #include "TH1.h"
00008 #include "TH2.h"
00009 #include "TH3.h"
00010 #include "TProfile.h"
00011 #include "TStyle.h"
00012 #include "TF1.h"
00013 #include "TTree.h"
00014 #include "TChain.h"
00015 #include "TFile.h"
00016 #include "TNtuple.h"
00017 #include "TCanvas.h"
00018 #include "TMinuit.h"
00019 #include "TSpectrum.h"
00020 #include "TString.h"
00021 #include "TLine.h"
00022 #include "TText.h"
00023 #include "TROOT.h"
00024 #include "TList.h"
00025 #include "TPolyMarker.h"
00026 #include "StBichsel/Bichsel.h"
00027 #include "BetheBloch.h"
00028 #include "TDirIter.h"
00029 #include "TTreeIter.h"
00030 #else
00031 class TSystem;
00032 class TMath;
00033 class TH1;
00034 class TH2;
00035 class TH3;
00036 class TProfile;
00037 class TStyle;
00038 class TF1;
00039 class TTree;
00040 class TChain;
00041 class TFile;
00042 class TNtuple;
00043 class TCanvas;
00044 class TMinuit;
00045 class TSpectrum;
00046 class TString;
00047 class TLine;
00048 class TText;
00049 class TROOT;
00050 class TList;
00051 class TPolyMarker;
00052 class Bichsel;
00053 class BetheBloch;
00054 class TDirIter;
00055 class TTreeIter;
00056 #endif
00057 #include "StLaserAnalysisMaker/laserino.h"
00058 void CheckMirrors(Int_t sector, Int_t var = 1) {
00059 static const Char_t *Names[3] = {"X","Y","Z"};
00060 TTree *laser = (TTree *) gDirectory->Get("laser");
00061 if (! laser ) return;
00062 TCanvas *c = new TCanvas(Form("Sector%02i%s",sector,Names[var-1]),Form("Sector%02i%s",sector,Names[var-1]),600,800);
00063 c->Divide(2,3);
00064 for (Int_t Bundle = 1; Bundle <= 6; Bundle++) {
00065 c->cd(Bundle);
00066 laser->Draw(Form("fTracks.Laser.XyzB.mX%i:fTracks.Laser.Mirror",var),
00067 Form("fTracks.Laser.Sector==%i&&fTracks.Laser.Bundle==%i",sector,Bundle),"colz",10000);
00068 }
00069 }
00070
00071 void CheckMirrors(const Char_t *files = "./laser_8102110.root") {
00072 TDirIter Dir(files);
00073 TTreeIter iter("laser");
00074 iter.AddFile(files);
00075 const Int_t& fNtrack = iter("fNtrack");
00076 const UShort_t*& fTracks_mNumberOfFitPointsTpc = iter("fTracks.mNumberOfFitPointsTpc");
00077 const Double32_t*& fTracks_fpTInv = iter("fTracks.fpTInv");
00078 const Int_t*& fTracks_Flag = iter("fTracks.Flag");
00079 const Double_t *&fTracks_dirPU_mX1 = iter("fTracks.dirPU.mX1");
00080 const Double_t*& fTracks_dirPM_mX1 = iter("fTracks.dirPM.mX1");
00081 const Double_t*& fTracks_dirPM_mX2 = iter("fTracks.dirPM.mX2");
00082 const Int_t *&fTracks_Laser_Sector = iter("fTracks.Laser.Sector");
00083 const Int_t *&fTracks_Laser_Bundle = iter("fTracks.Laser.Bundle");
00084 const Int_t *&fTracks_Laser_Mirror = iter("fTracks.Laser.Mirror");
00085 const Double_t*& fTracks_XyzPB_mX1 = iter("fTracks.XyzPB.mX1");
00086 const Double_t*& fTracks_XyzPB_mX2 = iter("fTracks.XyzPB.mX2");
00087 const Double_t*& fTracks_XyzPB_mX3 = iter("fTracks.XyzPB.mX3");
00088 const Double_t*& fTracks_XyzPL_mX3 = iter("fTracks.XyzPL.mX3");
00089 const Double_t*& fTracks_Laser_XyzB_mX1 = iter("fTracks.Laser.XyzB.mX1");
00090 const Double_t*& fTracks_Laser_XyzB_mX2 = iter("fTracks.Laser.XyzB.mX2");
00091 const Double_t*& fTracks_Laser_XyzB_mX3 = iter("fTracks.Laser.XyzB.mX3");
00092 const Int_t NS = 12;
00093 const Int_t NB = 6;
00094 const Int_t NM = 7;
00095 TH2D *X = new TH2D("X","dX versus mirror, bundle, sector",NS*NB*NM,0,NS*NB*NM,500,-.5,0.5);
00096 TH2D *Y = new TH2D("Y","dY versus mirror, bundle, sector",NS*NB*NM,0,NS*NB*NM,500,-.5,0.5);
00097 TH2D *Z = new TH2D("Z","dZ versus mirror, bundle, sector",NS*NB*NM,0,NS*NB*NM,500,-.5,0.5);
00098 while (iter.Next()) {
00099 for (Int_t k = 0; k < fNtrack; k++) {
00100 if (fTracks_Flag[k] < 2) continue;
00101 if (fTracks_dirPU_mX1[k] > -0.5) continue;
00102 if (fTracks_mNumberOfFitPointsTpc[k] < 35) continue;
00103 static const Double_t pTInv0 = 4.78815e-03;
00104 static const Double_t DpTInv0 = 9.75313e-03;
00105 if (TMath::Abs(fTracks_fpTInv[k] - pTInv0) > 3*DpTInv0) continue;
00106 Int_t s = fTracks_Laser_Sector[k]/2 - 1;
00107 if (s < 0 || s >= NS) continue;
00108 Int_t b = fTracks_Laser_Bundle[k] - 1;
00109 if (b < 0 || b >= NB) continue;
00110 Int_t m = fTracks_Laser_Mirror[k] - 1;
00111 if (m < 0 || m >= NM) continue;
00112 Double_t x = fTracks_XyzPB_mX1[k] - fTracks_Laser_XyzB_mX1[k];
00113 Double_t y = fTracks_XyzPB_mX2[k] - fTracks_Laser_XyzB_mX2[k];
00114 Double_t z = fTracks_XyzPB_mX3[k] - fTracks_Laser_XyzB_mX3[k];
00115
00116 #if 0
00117 if (s < 6) z -= 6.39057;
00118 else z -= 6.59919;
00119 z -= 3.29319000000000009e-04*fTracks_XyzPL_mX3[k];
00120 #else
00121 if (s < 6) z -= 6.45043e+00+1.26835e-04*fTracks_XyzPL_mX3[k];
00122 else z -= 6.60302e+00+3.93078e-04*fTracks_XyzPL_mX3[k];
00123 #endif
00124 Double_t indx = m + NM*(b + NB*s) + 0.5;
00125 X->Fill(indx,x);
00126 Y->Fill(indx,y);
00127 Z->Fill(indx,z);
00128 }
00129 }
00130 X->FitSlicesY();
00131 Y->FitSlicesY();
00132 Z->FitSlicesY();
00133 TH1D *Fit[3][3];
00134 const Char_t *xyz[3] = {"X","Y","Z"};
00135 const Char_t *res[3] = {"1","2","chi2"};
00136 for (Int_t i = 0; i < 3; i++)
00137 for (Int_t j = 0; j < 3; j++) Fit[i][j] = (TH1D*) gDirectory->Get(Form("%s_%s",xyz[i],res[j]));
00138 ofstream out;
00139 TString fOut(files);
00140 fOut,ReplaceAll("*","");
00141 fOut.ReplaceAll(".root","");
00142 fOut,ReplaceAll(".","");
00143 fOut.ReplaceAll("/","");
00144 fOut += ".data";
00145 cout << "Create " << fOut << endl;
00146 out.open(fOut.Data());
00147 for (Int_t s = 0; s < 12; s++) {
00148 for (Int_t b = 0; b < 6; b++) {
00149 for (Int_t m = 0; m < 7; m++) {
00150 Int_t bin = m + NM*(b + NB*s) + 1;
00151 Double_t dxyz[3] = {0,0,0};
00152 Double_t ddxyz[3] = {0,0,0};
00153 Int_t iflag[3] = {0,0,0};
00154 for (Int_t i = 0; i < 3; i++) {
00155 if (Fit[i][0]->GetBinError(bin) > 0 && Fit[i][0]->GetBinError(bin) < 0.1) {
00156 if (TMath::Abs(Fit[i][0]->GetBinContent(bin)) > Fit[i][0]->GetBinError(bin)) iflag[i] = 1;
00157 dxyz[i] = Fit[i][0]->GetBinContent(bin);
00158 ddxyz[i] = Fit[i][0]->GetBinError(bin);
00159 }
00160 }
00161 if (iflag[0] + iflag[1] + iflag[2] || ddxyz[0] > 0 || dxyz[1] > 0 || dxyz[2] > 0)
00162 cout << Form("Sector: %2i",2*s + 2) << " Bundle: " << b+1 << " Mirror: " << m+1
00163 << " <x> = " << Form("%7.4f",dxyz[0]) << " +/- " << Form("%7.4f",ddxyz[0])
00164 << " <y> = " << Form("%7.4f",dxyz[1]) << " +/- " << Form("%7.4f",ddxyz[1])
00165 << " <z> = " << Form("%7.4f",dxyz[2]) << " +/- " << Form("%7.4f",ddxyz[2])
00166 << endl;
00167 out << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dX+dxyz[0],ddxyz[0])
00168 << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dY+dxyz[1],ddxyz[1])
00169 << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dZ+dxyz[2],ddxyz[2])
00170 << "}";
00171 if (m == 6) out << "}";
00172 if (b == 5 && m == 6) out << "}";
00173 out << ",";
00174 out << "// S/B/M " << 2*(s+1) << "/" << b+1 << "/" << m+1 << endl;
00175 }
00176 }
00177 }
00178 out.close();
00179 }
00180 #if 0
00181
00182 void CheckMirrors(const Char_t *files = "./laser_8102110.root") {
00183 TDirIter Dir(files);
00184 TTreeIter iter("laser");
00185 iter.AddFile(files);
00186 const Int_t& fNtrack = iter("fNtrack");
00187 const UShort_t*& fTracks_mNumberOfFitPointsTpc = iter("fTracks.mNumberOfFitPointsTpc");
00188 const Double32_t*& fTracks_fpTInv = iter("fTracks.fpTInv");
00189 const Int_t*& fTracks_Flag = iter("fTracks.Flag");
00190 const Double_t *&fTracks_dirPU_mX1 = iter("fTracks.dirPU.mX1");
00191 const Double_t*& fTracks_dirPM_mX1 = iter("fTracks.dirPM.mX1");
00192 const Double_t*& fTracks_dirPM_mX2 = iter("fTracks.dirPM.mX2");
00193 const Int_t *&fTracks_Laser_Sector = iter("fTracks.Laser.Sector");
00194 const Int_t *&fTracks_Laser_Bundle = iter("fTracks.Laser.Bundle");
00195 const Int_t *&fTracks_Laser_Mirror = iter("fTracks.Laser.Mirror");
00196 const Double_t*& fTracks_XyzPB_mX1 = iter("fTracks.XyzPB.mX1");
00197 const Double_t*& fTracks_XyzPB_mX2 = iter("fTracks.XyzPB.mX2");
00198 const Double_t*& fTracks_XyzPB_mX3 = iter("fTracks.XyzPB.mX3");
00199 const Double_t*& fTracks_Laser_XyzB_mX1 = iter("fTracks.Laser.XyzB.mX1");
00200 const Double_t*& fTracks_Laser_XyzB_mX2 = iter("fTracks.Laser.XyzB.mX2");
00201 const Double_t*& fTracks_Laser_XyzB_mX3 = iter("fTracks.Laser.XyzB.mX3");
00202 Int_t N = 0;
00203 struct Kxy_t {
00204 Int_t N;
00205 Double_t K, Kx, K2, K2x, Ky, x, y, x2, xy, y2, del, Kdel, z, z2;
00206 };
00207 Kxy_t Offset[12][6][8];
00208 const Int_t NW = (sizeof(Kxy_t) - sizeof(Int_t))/sizeof(Double_t);
00209 memset(Offset, 0, 12*6*8*sizeof(Kxy_t));
00210 while (iter.Next()) {
00211
00212 N++;
00213
00214 for (Int_t k = 0; k < fNtrack; k++) {
00215 if (fTracks_Flag[k] < 2) continue;
00216 if (fTracks_dirPU_mX1[k] > -0.5) continue;
00217 if (fTracks_mNumberOfFitPointsTpc[k] < 35) continue;
00218 static const Double_t pTInv0 = 4.78815e-03;
00219 static const Double_t DpTInv0 = 9.75313e-03;
00220 if (TMath::Abs(fTracks_fpTInv[k] - pTInv0) > 3*DpTInv0) continue;
00221 Int_t s = fTracks_Laser_Sector[k]/2 - 1;
00222 if (s < 0 || s > 11) continue;
00223 Int_t b = fTracks_Laser_Bundle[k] - 1;
00224 if (b < 0 || b > 5) continue;
00225 Int_t m = fTracks_Laser_Mirror[k] - 1;
00226 if (m < 0 || m > 6) continue;
00227 LOOP:
00228 Kxy_t *mirror = &Offset[s][b][m];
00229 Double_t K = fTracks_dirPM_mX2[k]/fTracks_dirPM_mX1[k];
00230 Double_t x = fTracks_XyzPB_mX1[k] - fTracks_Laser_XyzB_mX1[k];
00231 Double_t y = fTracks_XyzPB_mX2[k] - fTracks_Laser_XyzB_mX2[k];
00232 Double_t z = fTracks_XyzPB_mX3[k] - fTracks_Laser_XyzB_mX3[k];
00233 if (s < 6) z -= 6.46746;
00234 else z -= 6.56439;
00235 if (TMath::Abs(z) > 0.5) continue;
00236 if (TMath::Abs(x) > 1 || TMath::Abs(y) > 1) continue;
00237 Double_t del = y - K*x;
00238 mirror->N++;
00239 mirror->K += K;
00240 mirror->Kx += x*K;
00241 mirror->K2 += K*K;
00242 mirror->K2x += K*K*x;
00243 mirror->Ky += K*y;
00244 mirror->x += x;
00245 mirror->y += y;
00246 mirror->x2 += x*x;
00247 mirror->xy += x*y;
00248 mirror->y2 += y*y;
00249 mirror->del += del;
00250 mirror->Kdel+= K*del;
00251 mirror->z += z;
00252 mirror->z2 += z*z;
00253 if (m != 7) {
00254 m = 7; goto LOOP;
00255 }
00256 }
00257 }
00258 ofstream out;
00259 TString fOut = "LaserCorrection.data";
00260 cout << "Create " << fOut << endl;
00261 out.open(fOut.Data());
00262 for (Int_t s = 0; s < 12; s++) {
00263 for (Int_t b = 0; b < 6; b++) {
00264 for (Int_t m = 0; m < 8; m++) {
00265 Kxy_t *mirror = &Offset[s][b][m];
00266 Double_t *xx = &mirror->K;
00267 if (mirror->N > 2) for (Int_t i = 0; i < NW; i++) xx[i] /= mirror->N;
00268 else for (Int_t i = 0; i < NW; i++) xx[i] = 0;
00269 Double_t dx, dy, dz;
00270 dx = dy = dz = 0;
00271 if ( mirror->N > 2 ) {
00272 dx = TMath::Sqrt(mirror->x2 - mirror->x* mirror->x);
00273 dy = TMath::Sqrt(mirror->y2 - mirror->y* mirror->y);
00274 dz = TMath::Sqrt(mirror->z2 - mirror->z* mirror->z);
00275 cout << "Sector: " << 2*s + 2 << " Bundle: " << b+1 << " Mirror: " << m+1 << " N: " << Form("%6i",mirror->N)
00276 << " <x> = " << Form("%7.4f",mirror->x) << " +/- " << Form("%7.4f",dx)
00277 << " <y> = " << Form("%7.4f",mirror->y) << " +/- " << Form("%7.4f",dy)
00278 << " <z> = " << Form("%7.4f",mirror->z) << " +/- " << Form("%7.4f",dz)
00279 << endl;
00280 }
00281 if (m < 7) {
00282 out << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dX+mirror->x,dx)
00283 << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dY+mirror->y,dy)
00284 << Form(",%7.4f,%7.4f",Mirrors[s][b][m].dZ+mirror->z,dz)
00285 << "}";
00286 if (m == 6) out << "}";
00287 if (b == 5 && m == 6) out << "}";
00288 out << ",";
00289 out << "// S/B/M " << 2*(s+1) << "/" << b+1 << "/" << m+1 << " => " << mirror->N << endl;
00290 }
00291 #if 0
00292 else {
00293 Double_t K = mirror->K;
00294 Double_t K2 = mirror->K2;
00295 Double_t del = mirror->del;
00296 Double_t Kdel = mirror->Kdel;
00297 Double_t det = K2 - K*K;
00298 cout << "======================================================" << endl;
00299 cout << "Sector: " << 2*s + 2 << " Bundle: " << b+1 << " Mirror: " << m+1 << " N: " << mirror->N << " det: " << det << endl;
00300 if (TMath::Abs(det) < 1e-7) continue;
00301 Double_t x0 = - (Kdel - K*del)/det;
00302 Double_t y0 = K*x0 + del;
00303 cout << "Sector: " << 2*s + 2 << " Bundle: " << b+1 << " Mirror: " << m+1
00304 << " x0/y0 = " << x0 << " / " << y0 << endl;
00305 cout << "======================================================" << endl;
00306 }
00307 #endif
00308 }
00309 }
00310 }
00311 out.close();
00312 }
00313
00314 void dZ(const Char_t *files = "./laser_8102110.root") {
00315 TDirIter Dir(files);
00316 TTreeIter iter("laser");
00317 iter.AddFile(files);
00318 const Int_t &fNtrack = iter("fNtrack");
00319 const Int_t *&fTracks_Flag = iter("fTracks.Flag");
00320 const Double_t *&fTracks_dU_mX1 = iter("fTracks.dU.mX1");
00321 const Double_t *&fTracks_dU_mX2 = iter("fTracks.dU.mX2");
00322 const Double_t *&fTracks_dU_mX3 = iter("fTracks.dU.mX3");
00323 const Int_t *&fTracks_Laser_Sector = iter("fTracks.Laser.Sector");
00324 const Int_t *&fTracks_Laser_Bundle = iter("fTracks.Laser.Bundle");
00325 const Int_t *&fTracks_Laser_Mirror = iter("fTracks.Laser.Mirror");
00326 const Double_t *&fTracks_Laser_XyzU_mX3 = iter("fTracks.Laser.XyzU.mX3");
00327
00328 while (iter.Next()) {
00329 for (Int_t k = 0; k < fNtrack; k++) {
00330 if (fTracks_Flag[k] < 2) continue;
00331 if (TMath::Abs(fTracks_dU_mX1[k]) > 0.1 || TMath::Abs(fTracks_dU_mX2[k]) > 0.1) continue;
00332 }
00333 }
00334 #endif
