00001 #include "tofcalgsec/calib_barrel_sigma.h"
00002 #include "TF1.h"
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 calib_barrel_sigma::calib_barrel_sigma( const unsigned int nzbin ):TofCalibFit( true, nBarrelSigma ) {
00023
00024 nKind = 5;
00025 nBinPerCounter = nzbin;
00026
00027 nHistPerCounter = nKind*nBinPerCounter;
00028 nCanvasPerCounter = 4;
00029 CanvasPerCounterName.push_back( static_cast<string>("Barrel-offset") );
00030 CanvasPerCounterName.push_back( static_cast<string>("Offset-TimeCorrelation") );
00031 CanvasPerCounterName.push_back( static_cast<string>("Barrel-sigma") );
00032 CanvasPerCounterName.push_back( static_cast<string>("Sigma-TimeCorrelation") );
00033 nGraphPerCanvasPerCounter.push_back(3);
00034 nGraphPerCanvasPerCounter.push_back(2);
00035 nGraphPerCanvasPerCounter.push_back(3);
00036 nGraphPerCanvasPerCounter.push_back(3);
00037
00038 nHistogram = 0;
00039 nCanvas = 0;
00040
00041 int numGraphs = 0;
00042 std::vector<unsigned int>::iterator iter = nGraphPerCanvasPerCounter.begin();
00043 for( ; iter!=nGraphPerCanvasPerCounter.end(); iter++ ) {
00044 numGraphs = numGraphs + (*iter);
00045 }
00046 if( numGraphs != nGraphTotalSigma ) {
00047 cout << "tofcalgsec::calib_barrel_sigma: the number of Graphs is NOT reasonable!!!" << endl;
00048 exit(0);
00049 }
00050
00051 m_name = string("calib_barrel_sigma");
00052
00053 const int tbin = 150;
00054 const double tbegin = -1.5;
00055 const double tend = 1.5;
00056
00057
00058 char hname[256];
00059 for( unsigned int i=0; i<NBarrel; i++ ) {
00060 m_result.push_back( HepVector(nBarrelSigma,0) );
00061 for( unsigned int j=0; j<nKind; j++ ) {
00062 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00063 if( j==0 ) { sprintf( hname, "tleft-id%i-z%i", i, k); }
00064 else if( j==1 ) { sprintf( hname, "tright-id%i-z%i", i, k); }
00065 else if( j==2 ) { sprintf( hname, "t0-id%i-z%i", i, k); }
00066 else if( j==3 ) { sprintf( hname, "tplus-id%i-z%i", i, k); }
00067 else if( j==4 ) { sprintf( hname, "tminus-id%i-z%i", i, k); }
00068 m_histograms.push_back( new TH1F( hname, hname, tbin, tbegin, tend ) );
00069
00070 m_fitresult.push_back( HepVector(nParSigma,0) );
00071 }
00072 }
00073 }
00074
00075 zpos.resize( nBinPerCounter );
00076 zposerr.resize( nBinPerCounter );
00077 zstep = ( zend - zbegin )/nBinPerCounter;
00078 for( unsigned int i=0; i<nBinPerCounter; i++ ) {
00079 zpos[i] = zbegin + ( i+0.5 )*zstep;
00080 zposerr[i] = 0.5*zstep;
00081 }
00082
00083 }
00084
00085
00086 calib_barrel_sigma::~calib_barrel_sigma() {
00087 m_fitresult.clear();
00088 zpos.clear();
00089 zposerr.clear();
00090 }
00091
00092
00093 void calib_barrel_sigma::calculate( RecordSet*& data, unsigned int icounter ) {
00094
00095 std::cout << setiosflags(ios::left) << setw(10) << icounter << setw(8) << data->size() << setw(30) << name() << std::endl;
00096
00097 if( data->size() > 0 ) {
00098 std::vector<Record*>::iterator iter = data->begin();
00099 for( ; iter!=data->end(); iter++ ) {
00100 fillRecord( (*iter), icounter );
00101 }
00102 }
00103 fitHistogram( icounter );
00104 fillGraph( icounter );
00105 fitGraph( icounter );
00106
00107 if( data->size() > 0 ) {
00108 std::vector<Record*>::iterator iter = data->begin();
00109 for( ; iter!=data->end(); iter++ ) {
00110 updateData( (*iter), icounter );
00111 fillRecordT0( (*iter), icounter );
00112 }
00113 }
00114 fitHistogramT0( icounter );
00115 fillGraphT0( icounter );
00116 fitGraphT0( icounter );
00117
00118 return;
00119 }
00120
00121
00122 void calib_barrel_sigma::fillRecord( const Record* r, unsigned int icounter ) {
00123
00124 double zhit = r->zrhit();
00125 if( zhit<zbegin || zhit>zend ) return;
00126 int zbin = static_cast<int>((zhit-zbegin)/zstep);
00127 if( zbin<0 ) { zbin = 0; }
00128 else if( zbin>static_cast<int>(nBinPerCounter-1) ) {
00129 cout << "tofcalgsec::calib_barrel_sigma:fillRecord: zhit is out of range, zhit=" << zhit << " zbin=" << zbin << endl;
00130 return;
00131 }
00132
00133 std::vector<TH1F*>::iterator iter = m_histograms.begin();
00134 unsigned int number = icounter*nKind*nBinPerCounter + zbin;
00135 (*(iter+number))->Fill( r->tleft() );
00136 (*(iter+nBinPerCounter+number))->Fill( r->tright() );
00137 (*(iter+3*nBinPerCounter+number))->Fill( (r->tleft()+r->tright())/2.0 );
00138 (*(iter+4*nBinPerCounter+number))->Fill( (r->tleft()-r->tright())/2.0 );
00139
00140 return;
00141 }
00142
00143
00144 void calib_barrel_sigma::fitHistogram( unsigned int icounter ) {
00145 TF1* g = new TF1("g", "gaus");
00146 g->SetLineColor(2);
00147 g->SetLineWidth(1);
00148
00149 std::vector<TH1F*>::iterator iter1 = m_histograms.begin() + icounter*nKind*nBinPerCounter;
00150 std::vector<HepVector>::iterator iter2 = m_fitresult.begin() + icounter*nKind*nBinPerCounter;
00151 for( unsigned int i=0; i<nKind; i++ ) {
00152 for( unsigned int j=0; j<nBinPerCounter; j++ ) {
00153 if( i!=2 ) {
00154 (*iter1)->Fit( g, "Q");
00155 (*iter2)[0] = g->GetParameter(1);
00156 (*iter2)[1] = g->GetParError(1);
00157 (*iter2)[2] = g->GetParameter(2);
00158 (*iter2)[3] = g->GetParError(2);
00159 }
00160 iter1++;
00161 iter2++;
00162 }
00163 }
00164
00165 return;
00166
00167 }
00168
00169
00170 void calib_barrel_sigma::fillGraph( unsigned int icounter ) {
00171
00172 char gname1[256], gname2[256];
00173
00174
00175
00176
00177
00178
00179
00180 std::vector<double> toffset, toffseterr;
00181 std::vector<double> tsigma, tsigmaerr;
00182 toffset.resize( nBinPerCounter );
00183 toffseterr.resize( nBinPerCounter );
00184 tsigma.resize( nBinPerCounter );
00185 tsigmaerr.resize( nBinPerCounter );
00186
00187 unsigned int number = 0;
00188 std::vector<HepVector>::iterator iter = m_fitresult.begin() + icounter*nKind*nBinPerCounter;
00189 for( unsigned int j=0; j<nKind; j++ ) {
00190 if( j==0 ) { sprintf( gname1, "tleft-offset-tofid-%i", icounter ); }
00191 else if( j==1 ) { sprintf( gname1, "tright-offset-tofid-%i", icounter ); }
00192 else if( j==2 ) { sprintf( gname1, "t0-offset-tofid-%i", icounter ); }
00193 else if( j==3 ) { sprintf( gname1, "tplus-offset-tofid-%i", icounter ); }
00194 else if( j==4 ) { sprintf( gname1, "tminus-offset-tofid-%i", icounter ); }
00195
00196 m_graphs.push_back( new TH1F( gname1, gname1, nBinPerCounter, zbegin, zend ) );
00197 std::vector<TH1F*>::iterator itgraph = m_graphs.end() - 1;
00198
00199 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00200 number = j*nBinPerCounter + k;
00201 toffset[k] = (*(iter+number))[0];
00202 toffseterr[k] = (*(iter+number))[1];
00203 (*itgraph)->SetBinContent( k+1, toffset[k] );
00204 (*itgraph)->SetBinError( k+1, toffseterr[k] );
00205 }
00206 (*itgraph)->SetMarkerSize(1.5);
00207 if( j==0 || j==3) {
00208 (*itgraph)->SetMarkerStyle(20);
00209 (*itgraph)->SetMarkerColor(2);
00210 (*itgraph)->SetMaximum( 0.15 );
00211 (*itgraph)->SetMinimum(-0.15 );
00212 }
00213 else if( j==1 || j==4 ) {
00214 (*itgraph)->SetMarkerStyle(21);
00215 (*itgraph)->SetMarkerColor(4);
00216 }
00217 else if( j==2 ) {
00218 (*itgraph)->SetMarkerStyle(4);
00219 (*itgraph)->SetMarkerColor(2);
00220 }
00221 }
00222
00223 for( unsigned int j=0; j<nKind; j++ ) {
00224 if( j==0 ) { sprintf( gname2, "tleft-sigma-tofid-%i", icounter ); }
00225 else if( j==1 ) { sprintf( gname2, "tright-sigma-tofid-%i", icounter ); }
00226 else if( j==2 ) { sprintf( gname2, "t0-sigma-tofid-%i", icounter ); }
00227 else if( j==3 ) { sprintf( gname2, "tplus-sigma-tofid-%i", icounter ); }
00228 else if( j==4 ) { sprintf( gname2, "tminus-sigma-tofid-%i", icounter ); }
00229 m_graphs.push_back( new TH1F( gname2, gname2, nBinPerCounter, zbegin, zend ) );
00230 std::vector<TH1F*>::iterator itgraph = m_graphs.end() - 1;
00231
00232 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00233 number = j*nBinPerCounter + k;
00234 tsigma[k] = (*(iter+number))[2];
00235 tsigmaerr[k] = (*(iter+number))[3];
00236 (*itgraph)->SetBinContent( k+1, tsigma[k] );
00237 (*itgraph)->SetBinError( k+1, tsigmaerr[k] );
00238 }
00239 (*itgraph)->SetMarkerSize(1.5);
00240 if( j==0 || j==3 ) {
00241 (*itgraph)->SetMarkerStyle(20);
00242 (*itgraph)->SetMarkerColor(2);
00243 (*itgraph)->SetMaximum( 0.3 );
00244 (*itgraph)->SetMinimum( 0.0 );
00245 }
00246 else if( j==1 || j==4 ) {
00247 (*itgraph)->SetMarkerStyle(21);
00248 (*itgraph)->SetMarkerColor(4);
00249 }
00250 else if( j==2 ) {
00251 (*itgraph)->SetMarkerStyle(4);
00252 (*itgraph)->SetMarkerColor(2);
00253 }
00254 }
00255
00256 sprintf( gname2, "sigma-tofid-%i", icounter );
00257 m_graphs.push_back( new TH1F( gname2, gname2, nBinPerCounter, zbegin, zend ) );
00258 std::vector<TH1F*>::iterator itgraph = m_graphs.end() - 1;
00259 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00260 number = (nKind-1)*nBinPerCounter + k;
00261 double sigPlus = (*(iter+number-nBinPerCounter))[2];
00262 double sigMinus = (*(iter+number))[2];
00263 double sigErrPlus = (*(iter+number-nBinPerCounter))[3];
00264 double sigErrMinus = (*(iter+number))[3];
00265 tsigma[k] = sqrt( sigPlus*sigPlus - sigMinus*sigMinus );
00266 tsigmaerr[k] = sqrt( sigErrPlus*sigErrPlus + sigErrMinus*sigErrMinus );
00267 (*itgraph)->SetBinContent( k+1, tsigma[k] );
00268 (*itgraph)->SetBinError( k+1, tsigmaerr[k] );
00269 }
00270 (*itgraph)->SetMarkerSize(1.5);
00271 (*itgraph)->SetMarkerStyle(4);
00272 (*itgraph)->SetMarkerColor(2);
00273
00274 return;
00275 }
00276
00277
00278 void calib_barrel_sigma::fitGraph( unsigned int icounter ) {
00279
00280 TF1* fsingle = new TF1("fsingle", "pol4");
00281 fsingle->SetLineColor(1);
00282 fsingle->SetLineWidth(1);
00283
00284 std::vector<unsigned int>::iterator itnumber = nGraphPerCanvasPerCounter.begin();
00285 std::vector<TH1F*>::iterator itgraph = m_graphs.begin() + icounter*nGraphTotalSigma + (*itnumber) + (*(itnumber+1));
00286
00287 fsingle->SetParameter( 0, 0.14 );
00288 fsingle->SetParameter( 1, -4.0e-4 );
00289 (*itgraph)->Fit( "fsingle", "QR", "", zbegin, zend );
00290 X = HepVector( m_npar, 0 );
00291 for( unsigned int i=0; i<5; i++ ) {
00292 X[i] = fsingle->GetParameter(i);
00293 }
00294
00295 fsingle->SetParameter( 0, 0.14 );
00296 fsingle->SetParameter( 1, 4.0e-4 );
00297 (*(itgraph+1))->Fit( "fsingle", "QR", "", zbegin, zend );
00298 for( unsigned int i=0; i<5; i++ ) {
00299 X[i+5] = fsingle->GetParameter(i);
00300 }
00301
00302 std::vector<HepVector>::iterator iter = m_result.begin() + icounter;
00303 (*iter) = X;
00304
00305 return;
00306 }
00307
00308
00309 void calib_barrel_sigma::updateData( Record* r, unsigned int icounter ) {
00310 double zhit = r->zrhit();
00311 double t1 = r->tleft();
00312 double t2 = r->tright();
00313
00314 double par1[5], par2[5];
00315 std::vector<HepVector>::iterator iter = m_result.begin() + icounter;
00316 for( unsigned int i=0; i<5; i++ ) {
00317 par1[i] = (*iter)[i];
00318 par2[i] = (*iter)[i+5];
00319 }
00320
00321 double tsigma1 = par1[0]+par1[1]*zhit+par1[2]*pow(zhit,2)+par1[3]*pow(zhit,3) + par1[4]*pow(zhit,4);
00322 double tsigma2 = par2[0]+par2[1]*zhit+par2[2]*pow(zhit,2)+par2[3]*pow(zhit,3) + par2[4]*pow(zhit,4);
00323 double tc = m_tcorrelation[0];
00324
00325 double weight1 = (tsigma2*tsigma2-tc*tc)/(tsigma1*tsigma1+tsigma2*tsigma2-2.0*tc*tc);
00326 double weight2 = (tsigma1*tsigma1-tc*tc)/(tsigma1*tsigma1+tsigma2*tsigma2-2.0*tc*tc);
00327 double t0 = weight1*t1 + weight2*t2;
00328
00329 r->setT0( t0 );
00330
00331 return;
00332 }
00333
00334
00335 void calib_barrel_sigma::fillRecordT0( const Record* r, unsigned int icounter ) {
00336 double zhit = r->zrhit();
00337 if( zhit<zbegin || zhit>zend ) return;
00338 int zbin = static_cast<int>((zhit-zbegin)/zstep);
00339 if( zbin<0 ) { zbin = 0; }
00340 else if( zbin>static_cast<int>(nBinPerCounter-1) ) {
00341 cout << "tofcalgsec::calib_barrel_sigma:fillRecordT0: zhit is out of range, zhit=" << zhit << " zbin=" << zbin << endl;
00342 return;
00343 }
00344
00345 std::vector<TH1F*>::iterator iter = m_histograms.begin();
00346 unsigned int number = icounter*nKind*nBinPerCounter + 2*nBinPerCounter + zbin;
00347 (*(iter+number))->Fill( r->t0() );
00348
00349 return;
00350 }
00351
00352
00353 void calib_barrel_sigma::fitHistogramT0( unsigned int icounter ) {
00354 TF1* g = new TF1("g", "gaus");
00355 g->SetLineColor(2);
00356 g->SetLineWidth(1);
00357
00358 std::vector<TH1F*>::iterator iter1 = m_histograms.begin() + icounter*nKind*nBinPerCounter + 2*nBinPerCounter;
00359 std::vector<HepVector>::iterator iter2 = m_fitresult.begin() + icounter*nKind*nBinPerCounter + 2*nBinPerCounter;
00360 for( unsigned int j=0; j<nBinPerCounter; j++, iter1++, iter2++ ) {
00361 (*iter1)->Fit( g, "Q");
00362 (*iter2)[0] = g->GetParameter(1);
00363 (*iter2)[1] = g->GetParError(1);
00364 (*iter2)[2] = g->GetParameter(2);
00365 (*iter2)[3] = g->GetParError(2);
00366 }
00367
00368 return;
00369 }
00370
00371
00372 void calib_barrel_sigma::fillGraphT0( unsigned int icounter ) {
00373 char gname1[256], gname2[256];
00374
00375
00376
00377
00378
00379
00380
00381 std::vector<double> toffset, toffseterr;
00382 std::vector<double> tsigma, tsigmaerr;
00383 toffset.resize( nBinPerCounter );
00384 toffseterr.resize( nBinPerCounter );
00385 tsigma.resize( nBinPerCounter );
00386 tsigmaerr.resize( nBinPerCounter );
00387
00388 std::vector<HepVector>::iterator iter = m_fitresult.begin() + icounter*nKind*nBinPerCounter + 2*nBinPerCounter;
00389 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00390 toffset[k] = (*(iter+k))[0];
00391 toffseterr[k] = (*(iter+k))[1];
00392 tsigma[k] = (*(iter+k))[2];
00393 tsigmaerr[k] = (*(iter+k))[3];
00394 }
00395
00396 sprintf( gname1, "offset-tofid-%i", icounter );
00397 std::vector<TH1F*>::iterator itgraph = m_graphs.begin() + icounter*nGraphTotalSigma + 2;
00398 for( unsigned int i=0; i<nBinPerCounter; i++ ) {
00399
00400
00401 (*itgraph)->SetBinContent( i+1, toffset[i] );
00402 (*itgraph)->SetBinError( i+1, toffseterr[i] );
00403 }
00404
00405 sprintf( gname2, "sigma-tofid-%i", icounter );
00406 itgraph = m_graphs.begin() + icounter*nGraphTotalSigma + 7;
00407 for( unsigned int i=0; i<nBinPerCounter; i++ ) {
00408
00409
00410 (*itgraph)->SetBinContent( i+1, tsigma[i] );
00411 (*itgraph)->SetBinError( i+1, tsigmaerr[i] );
00412 }
00413
00414 return;
00415 }
00416
00417
00418 void calib_barrel_sigma::fitGraphT0( unsigned int icounter ) {
00419
00420
00421 TF1 *fdouble = new TF1( "fdouble", "pol4", zbegin, zend );
00422 fdouble->SetLineColor(1);
00423 fdouble->SetLineWidth(1);
00424
00425 std::vector<unsigned int>::iterator itnumber = nGraphPerCanvasPerCounter.begin();
00426 std::vector<TH1F*>::iterator itgraph = m_graphs.begin() + icounter*nGraphTotalSigma + (*itnumber) + (*(itnumber+1)) + 2;
00427 (*itgraph)->Fit( "fdouble", "Q", "", zbegin, zend );
00428
00429 std::vector<HepVector>::iterator iter = m_result.begin() + icounter;
00430 (*iter)[10] = fdouble->GetParameter(0);
00431 (*iter)[11] = fdouble->GetParameter(1);
00432 (*iter)[12] = fdouble->GetParameter(2);
00433 (*iter)[13] = fdouble->GetParameter(3);
00434 (*iter)[14] = fdouble->GetParameter(4);
00435
00436 return;
00437 }