00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #define INTEGRATE_OVER_HOURS
00027
00028 #if !defined(__CINT__) || defined(__MAKECINT__)
00029
00030 #include "Riostream.h"
00031 #include <stdio.h>
00032 #include "TROOT.h"
00033 #include "TSystem.h"
00034 #include "TMath.h"
00035 #include "TH1.h"
00036 #include "TH2.h"
00037 #include "TH3.h"
00038 #include "TStyle.h"
00039 #include "TF1.h"
00040 #include "TProfile.h"
00041 #include "TTree.h"
00042 #include "TChain.h"
00043 #include "TFile.h"
00044 #include "TNtuple.h"
00045 #include "TCanvas.h"
00046 #include "TString.h"
00047 #include "TList.h"
00048 #include "TLegend.h"
00049 #include "TFileSet.h"
00050 #include "TDataSetIter.h"
00051 #include "TDirIter.h"
00052 #include "TTreeIter.h"
00053 #endif
00054 static Int_t run = 0;
00055 static Int_t date = 0;
00056 static Int_t Time = 0;
00057
00058 static Double_t DVAll[2][3];
00059 static Double_t dDVAll[2][3];
00060 static Int_t _debug = 0;
00061 TH2D *dv = 0;
00062 TH2D *slope = 0;
00063 TNtuple *runNT = 0;
00064 struct Run_t {
00065 Float_t run, date, time, events, day, dvAll, ddvAll, dvWest, ddvWest, dvEast, ddvEast, slAll, dslAll, slWest, dslWest, slEast, dslEast;
00066 };
00067 const Char_t *vRun = "run:date:time:events:day:dvAll:ddvAll:dvWest:ddvWest:dvEast:ddvEast:slAll:dslAll:slWest:dslWest:slEast:dslEast";
00068 Run_t Run;
00069
00070 Double_t ScaleE2W(Double_t day) {
00071
00072 static Double_t par[2] = {1.57361e-01,-4.59752e-03};
00073 return par[0] + par[1]*day;
00074 }
00075
00076 void MakeTable() {
00077 #ifndef SeparateWestandEast
00078 if (! (DVAll[0][0] > 5.3 && DVAll[0][0] < 5.9 && dDVAll[0][0] > 0 && dDVAll[0][0]< 5e-4)) {
00079 cout << "Run " << run << " fails =============================" << endl;
00080 return;
00081 }
00082 #else
00083 if (! (DVAll[0][1] > 5.3 && DVAll[0][1] < 5.9 && dDVAll[0][1] > 0 && dDVAll[0][1]< 4e-5 ||
00084 DVAll[0][2] > 5.3 && DVAll[0][2] < 5.9 && dDVAll[0][2] > 0 && dDVAll[0][2]< 4e-5)) {
00085 cout << "Run " << run << " fails =============================" << endl;
00086 return;
00087 }
00088 #endif
00089 TString fOut = Form("tpcDriftVelocity.%8i.%06i.C",date,Time);
00090 ofstream out;
00091 cout << "Create " << fOut << endl;
00092 out.open(fOut.Data());
00093 out << "TDataSet *CreateTable() {" << endl;
00094 out << " if (!gROOT->GetClass(\"St_tpcDriftVelocity\")) return 0;" << endl;
00095 out << " St_tpcDriftVelocity *tableSet = new St_tpcDriftVelocity(\"tpcDriftVelocity\",1);" << endl;
00096 out << " tpcDriftVelocity_st row;// Laser Run " << run << endl;
00097 Double_t dvEast = DVAll[0][2];
00098 Double_t ddvEast = dDVAll[0][2];
00099 Double_t dvWest = DVAll[0][1];
00100 Double_t ddvWest = dDVAll[0][1];
00101 #ifdef SeparateWestandEast
00102 if (! (dvWest > 5.5 && dvWest < 5.9 && ddvWest > 0 && ddvWest< 4e-5) ) {
00103 ddvWest = -1;
00104 dvWest = dvEast/(1 + 1e-3*ScaleE2W(Run.day));
00105 }
00106 if (! (dvEast > 5.5 && dvEast < 5.9 && ddvEast > 0 && ddvEast< 4e-5) ) {
00107 ddvEast = -1;
00108 dvEast = dvWest*(1 + 1e-3*ScaleE2W(Run.day));
00109 }
00110 out << " row.laserDriftVelocityEast = " << dvEast << "; // +/- " << ddvEast
00111 << " cm/us East: Slope = " << DVAll[1][2] << " +/- " << dDVAll[1][2] << " DV = " << DVAll[0][2] << " +/- " << dDVAll[0][2]<< endl;
00112 out << " row.laserDriftVelocityWest = " << dvWest << "; // +/- " << ddvWest
00113 << " cm/us West: Slope = " << DVAll[1][1] << " +/- " << dDVAll[1][1] << " DV = " << DVAll[0][1] << " +/- " << dDVAll[0][1]<< endl;
00114 out << " row.cathodeDriftVelocityEast = 0; // cm/us : from cathode emission ;" << endl;
00115 out << " row.cathodeDriftVelocityWest = 0; // cm/us : from cathode emission ;" << endl;
00116 out << " tableSet->AddAt(&row); " << endl;
00117 out << " return (TDataSet *)tableSet;" << endl;
00118 #else
00119 out << " row.laserDriftVelocityEast = " << DVAll[0][0] << "; // +/- " << dDVAll[0][0]
00120 << " cm/us All: East = " << DVAll[1][2] << " +/- " << dDVAll[1][2] << endl;
00121 out << " row.laserDriftVelocityWest = " << DVAll[0][0] << "; // +/- " << dDVAll[0][0]
00122 << " cm/us All: West = " << DVAll[1][1] << " +/- " << dDVAll[1][1] << endl;
00123 out << " row.cathodeDriftVelocityEast = 0; // cm/us : from cathode emission ;" << endl;
00124 out << " row.cathodeDriftVelocityWest = 0; // cm/us : from cathode emission ;" << endl;
00125 out << " tableSet->AddAt(&row);// 1e3*Delta: All = " << DVAll[0][0] << " +/- " << dDVAll[0][0] << endl;
00126 out << " return (TDataSet *)tableSet;//"
00127 << " West = " << DVAll[0][1] << " +/- " << dDVAll[0][1]
00128 << " East = " << DVAll[0][2] << " +/- " << dDVAll[0][2] << endl;
00129 #endif
00130 out << "};" << endl;
00131 }
00132
00133 void Fit() {
00134 static TDatime t0(2007,1,1,0,0,0);
00135 static UInt_t u0 = t0.Convert();
00136 date = 20000000 + (Int_t) Run.date;
00137 Time = (Int_t) Run.time;
00138 TDatime t(date, Time);
00139 UInt_t u = t.Convert();
00140 Run.day = 1. + (u - u0)/(24.*60.*60.);
00141 run = (Int_t) (1000000*((Int_t) (Run.date/1000000)) + Run.run);
00142 memset(&DVAll[0][0], 0, 6*sizeof(Double_t));
00143 memset(&dDVAll[0][0], 0, 6*sizeof(Double_t));
00144 Run.events = slope->GetEntries();
00145 cout << "Run " << run << " has " << Run.events << " entries" << endl;
00146 TH2D *hist[2] = {dv, slope};
00147 Float_t *par = &Run.dvAll;
00148 for (Int_t l = 0; l < 2; l++) {
00149 #ifdef ADJUSTABLE_BINNING
00150 hist[l]->BufferEmpty();
00151 #endif
00152 hist[l]->Write();
00153 hist[l]->FitSlicesY(0,1,0,10,"QNRI");
00154 TString fitN(hist[l]->GetName());
00155 fitN += "_1";
00156 TH1D *fit = (TH1D *) gDirectory->Get(fitN);
00157 if (! fit) continue;
00158 fit->SetMarkerStyle(20);
00159 Double_t xmin = fit->GetXaxis()->GetXmin();
00160 Double_t xmax = fit->GetXaxis()->GetXmax();
00161 TF1 *pol0 = (TF1*) gROOT->GetFunction("pol0");
00162 if (pol0) {
00163 for (Int_t k = 0; k < 3; k++) {
00164 pol0->SetLineColor(k+1);
00165 if (k == 0) fit->Fit(pol0,"er","",xmin,xmax);
00166 if (k == 1) {
00167 if (xmin >= 13) continue;
00168 fit->Fit(pol0,"er+","",xmin,12.5);
00169 }
00170 if (k == 2) {
00171 if (xmax <= 15) continue;
00172 fit->Fit(pol0,"er+","",13.5,xmax);
00173 }
00174 DVAll[l][k] = pol0->GetParameter(0);
00175 dDVAll[l][k] = pol0->GetParError(0);
00176 par[2*(k+3*l)] = DVAll[l][k];
00177 par[2*(k+3*l)+1] = dDVAll[l][k];
00178 }
00179 }
00180 fit->Write();
00181 }
00182 runNT->Fill(&Run.run);
00183 MakeTable();
00184 }
00185
00186 void LoopOverLaserTrees(const Char_t *files="./st_laser_*.laser.root") {
00187 TDirIter Dir(files);
00188 TTreeIter iter("laser");
00189 iter.AddFile(files);
00190
00191 const Int_t& fEvtHdr_fRun = iter("fEvtHdr.fRun");
00192 const Int_t& fEvtHdr_fDate = iter("fEvtHdr.fDate");
00193 const Int_t& fEvtHdr_fTime = iter("fEvtHdr.fTime");
00194
00195 const Float_t& fEvtHdr_fDriVel = iter("fEvtHdr.fDriVel");
00196 const Float_t& fEvtHdr_fClock = iter("fEvtHdr.fClock");
00197
00198
00199
00200
00201
00202 const Float_t& fEvtHdr_fOnlClock = iter("fEvtHdr.fOnlClock");
00203 const Int_t*& fFit_N = iter("fFit.N");
00204 const Int_t*& fFit_Sector = iter("fFit.Sector");
00205 #ifdef __REFIT__
00206 const Int_t*& fFit_Bundle = (Int_t **) iter("fFit.Bundle[42]");
00207 const Int_t*& fFit_Mirror = (Int_t **) iter("fFit.Mirror[42]");
00208 #endif
00209 const Double32_t*& fFit_offset = iter("fFit.offset");
00210 const Double32_t*& fFit_slope = iter("fFit.slope");
00211 const Double32_t*& fFit_doffset = iter("fFit.doffset");
00212 const Double32_t*& fFit_dslope = iter("fFit.dslope");
00213 const Double32_t*& fFit_chisq = iter("fFit.chisq");
00214 #ifdef __REFIT__
00215 const Double32_t**& fFit_X = (Double_t **) iter("fFit.X[42]");
00216 const Double32_t**& fFit_Y = (Double_t **) iter("fFit.Y[42]");
00217 #endif
00218 const Double32_t*& fFit_Prob = iter("fFit.Prob");
00219 const Int_t*& fFit_ndf = iter("fFit.ndf");
00220
00221 TFile *fOut = new TFile("LaserPlots.root","recreate");
00222 runNT = new TNtuple("RunNT","Run date time",vRun);
00223 static const Double_t EastWRTWestDiff = 0;
00224
00225 Double_t OnlFreq = 0;
00226 Int_t oldRun = -1;
00227 Double_t oldDate = -1;
00228 TDatime t0(2007,1,1,0,0,0);
00229 UInt_t ut0 = t0.Convert();
00230 while (iter.Next()) {
00231 #if 0
00232 Int_t iok = -1;
00233 for (Int_t i = 0; i < NZFR; i++) {if (ZFieldRuns[i] == fEvtHdr_fRun) {iok = i; break;}}
00234 if (iok >= 0) continue;
00235 #endif
00236 if (fEvtHdr_fRun != oldRun) {
00237 OnlFreq = fEvtHdr_fOnlClock;
00238 fOut->cd();
00239 TDatime t(fEvtHdr_fDate,fEvtHdr_fTime);
00240 UInt_t ut = t.Convert();
00241 Double_t Date = 1. + (ut - ut0)/(24.*60.*60.);
00242 #ifdef INTEGRATE_OVER_HOURS
00243 if (Date - oldDate > 0.125) {
00244 #endif
00245 if (oldRun != -1) Fit();
00246 cout << "New run " << fEvtHdr_fRun << " Date Old/New " << oldDate << "/" << Date << endl;
00247 Run.run = fEvtHdr_fRun%1000000;
00248 Run.date = fEvtHdr_fDate%1000000;
00249 Run.time = fEvtHdr_fTime;
00250 Run.events = 0;
00251 #ifdef ADJUSTABLE_BINNING
00252 dv = new TH2D(Form("DV%i",fEvtHdr_fRun%1000000),Form("Drift Velocity for run %i",fEvtHdr_fRun%1000000),12,1,25,400,1,-1);
00253 dv->SetBuffer(100000);
00254 #else
00255 dv = new TH2D(Form("DV%i",fEvtHdr_fRun%1000000),Form("Drift Velocity for run %i",fEvtHdr_fRun%1000000),12,1,25,2000,5.3,5.9);
00256 #endif
00257 dv->SetXTitle("Sector");
00258 dv->SetYTitle("Drift Velocity ");
00259 #ifdef ADJUSTABLE_BINNING
00260 slope = new TH2D(Form("SL%i",fEvtHdr_fRun%1000000),Form("Slope for run %i",fEvtHdr_fRun%1000000),12,1,25,400,1,-1);
00261 slope->SetBuffer(100000);
00262 #else
00263 slope = new TH2D(Form("SL%i",fEvtHdr_fRun%1000000),Form("Slope for run %i",fEvtHdr_fRun%1000000),12,1,25,2000,-10.,10.);
00264 #endif
00265 slope->SetXTitle("Sector");
00266 slope->SetYTitle("Difference wrt reference Drift Velocity in pemill");
00267 oldRun = (Int_t) fEvtHdr_fRun;
00268 oldDate = Date;
00269 #ifdef INTEGRATE_OVER_HOURS
00270 }
00271 #endif
00272 }
00273 Double_t dt = fEvtHdr_fDate%100000 + ((Double_t) fEvtHdr_fTime)*1e-6;
00274 Double_t DT = Run.date + Run.time*1e-6;
00275 if (dt < DT) {
00276 Run.date = fEvtHdr_fDate%1000000;
00277 Run.time = fEvtHdr_fTime;
00278 }
00279 for (Int_t k = 0; k < 12; k++) {
00280 #ifndef __REFIT__
00281 if (fFit_ndf[k] > 2 && fFit_Prob[k] > 1.e-2) {
00282 Double_t Slope = 1e3*fFit_slope[k];
00283 if (fFit_Sector[k] > 12) Slope -= EastWRTWestDiff;
00284 slope->Fill(fFit_Sector[k], Slope);
00285 Double_t Vm = OnlFreq/fEvtHdr_fClock*fEvtHdr_fDriVel;
00286 Double_t V = 1e-6*Vm/(1.+1e-3*Slope-Vm/TMath::Ccgs());
00287 dv->Fill(fFit_Sector[k],V);
00288 }
00289 #else
00290 Int_t N = fFit_N[k];
00291 if (N > 2) {
00292 Double_t Flag[42];
00293 memset (Flag, 0, 42*sizeof(Double_t));
00294 Double_t YA = 0;
00295 for (Int_t i = 0; i < N; i++) {
00296 YA += fFit_Y[42*k+i];
00297 }
00298 YA /= N;
00299 for (Int_t i = 0; i < N; i++) {if (TMath::Abs(fFit_Y[42*k+i] - YA) > 10) Flag[i] = 1;}
00300 for (;;) {
00301 TRVector Y;
00302 TRMatrix A(0,2);
00303 for (Int_t i = 0; i < N; i++) {
00304 if (! fFit_Flag[42*k+i]) {
00305 Double_t dev = fFit_Y[42*k+i]/sigma;
00306 Y.AddRow(&dev);
00307 Double_t a[2] = {1./sigma, fFit_X[42*k+i]/sigma};
00308 A.AddRow(a);
00309 }
00310 }
00311 Int_t ndf = A.GetNrows() - 2;
00312 if (ndf < 1) break;
00313 TRSymMatrix S(A,TRArray::kATxA); if (_debug) cout << "S: " << S << endl;
00314 TRVector B(A,TRArray::kATxB,Y); if (_debug) cout << "B: " << B << endl;
00315 TRSymMatrix SInv(S,TRArray::kInverted); if (_debug) cout << "SInv: " << SInv << endl;
00316 TRVector X(SInv,TRArray::kSxA,B); if (_debug) cout << "X: " << X << endl;
00317 TRVector R(Y);
00318 R -= TRVector(A,TRArray::kAxB,X); if (_debug) cout << "R: " << R << endl;
00319 Double_t chisq = R*R;
00320 Double_t prob = TMath::Prob(chisq,ndf);
00321 if (prob > 0.01) {
00322 Double_t offset = X[0];
00323 Double_t slope = X[1];
00324 Double_t chisq = chisq;
00325 Double_t dslope = SInv[2];
00326 Double_t doffset = SInv[0];
00327 if (ndf > 4 && prob > 1.e-2) {
00328 Double_t Slope = 1e3*slope;
00329 if (fFit_Sector[k] > 12) Slope -= EastWRTWestDiff;
00330 slope->Fill(fFit_Sector[k], Slope);
00331 Double_t Vm = OnlFreq/fEvtHdr_fClock*fEvtHdr_fDriVel;
00332 Double_t V = 1e-6*Vm/(1.+1e-3*Slope-Vm/TMath::Ccgs());
00333 dv->Fill(fFit_Sector[k],V);
00334 }
00335 }
00336 break;
00337 } else {
00338 Int_t j = -1;
00339 Int_t imax = -1;
00340 Double_t Rmax = 0;
00341 for (Int_t i = 0; i < N; i++) {
00342 if (! fit->Flag[i]) {
00343 j++;
00344 Double_t RR = TMath::Abs(R[j]);
00345 if (RR > Rmax) {
00346 imax = i;
00347 Rmax = RR;
00348 }
00349 }
00350 }
00351 if (imax < 0) break;
00352 Flag[imax] = 1;
00353 }
00354 }
00355 #endif
00356 }
00357 }
00358 Fit();
00359 runNT->Write();
00360 delete fOut;
00361 }
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376 Double_t OffSets(Double_t *x, Double_t *par = 0) {
00377
00378
00379 static Double_t Toffsets[12] = {10.33, 3.34, 6.14, 13.11, -1., 17.31,
00380 19.88, 12.95, 5.97, 3.18, 10.17, 17.14};
00381 static Double_t DV = 55.42840e-4;
00382 Int_t sector2 = (Int_t) ((x[0]-1)/2.);
00383 Double_t off = -1;
00384 if (sector2 < 0 || sector2 > 12) return off;
00385 return DV*Toffsets[sector2];
00386 }
