00001 #include "tofcalgsec/calib_barrel_common.h"
00002 #include "TF1.h"
00003
00004 calib_barrel_common::calib_barrel_common( const unsigned int nzbin ):TofCalibFit( true, nBarrelCommon ) {
00005
00006 nKind = 4;
00007 nBinPerCounter = nzbin;
00008
00009 nHistPerCounter = 0;
00010 nCanvasPerCounter = 0;
00011 nHistogram = nKind*nBinPerCounter+1;
00012 nCanvas = 4;
00013 CanvasName.push_back( static_cast<string>("Offset") );
00014 CanvasName.push_back( static_cast<string>("Offset-PlusMinus") );
00015 CanvasName.push_back( static_cast<string>("Sigma") );
00016 CanvasName.push_back( static_cast<string>("Sigma-TCorrelation") );
00017 nGraphPerCanvas.push_back(2);
00018 nGraphPerCanvas.push_back(2);
00019 nGraphPerCanvas.push_back(2);
00020 nGraphPerCanvas.push_back(3);
00021
00022 int numGraphs = 0;
00023 std::vector<unsigned int>::iterator iter = nGraphPerCanvas.begin();
00024 for( ; iter!=nGraphPerCanvas.end(); iter++ ) {
00025 numGraphs = numGraphs + (*iter);
00026 }
00027 if( numGraphs != nGraphTotalCommon ) {
00028 cout << "tofcalgsec::calib_barrel_common: the number of Graphs is NOT reasonable!!!" << endl;
00029 exit(0);
00030 }
00031
00032 m_name = string("calib_barrel_common");
00033
00034 const int tbin = 150;
00035 const double tbegin = -1.5;
00036 const double tend = 1.5;
00037
00038 char hname[256];
00039
00040 for( unsigned int j=0; j<nKind; j++ ) {
00041 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00042 if( j==0 ) { sprintf( hname, "tleft-z%i", k); }
00043 else if( j==1 ) { sprintf( hname, "tright-z%i", k); }
00044 else if( j==2 ) { sprintf( hname, "tplus-z%i", k); }
00045 else if( j==3 ) { sprintf( hname, "tminus-z%i", k); }
00046 m_fitresult.push_back( HepVector( nParCommon,0 ) );
00047 m_histograms.push_back( new TH1F( hname, hname, tbin, tbegin, tend ) );
00048 }
00049 }
00050 sprintf( hname, "deltaT" );
00051 m_fitresult.push_back( HepVector( nParCommon,0 ) );
00052 m_histograms.push_back( new TH1F( hname, hname, tbin, tbegin, tend ) );
00053
00054 zpos.resize( nBinPerCounter );
00055 zposerr.resize( nBinPerCounter );
00056 zstep = ( zend - zbegin )/nBinPerCounter;
00057 for( unsigned int i=0; i<nBinPerCounter; i++ ) {
00058 zpos[i] = zbegin + ( i+0.5 )*zstep;
00059 zposerr[i] = 0.5*zstep;
00060 }
00061
00062 }
00063
00064
00065 calib_barrel_common::~calib_barrel_common() {
00066 m_fitresult.clear();
00067 zpos.clear();
00068 zposerr.clear();
00069 }
00070
00071
00072 void calib_barrel_common::calculate( RecordSet*& data, unsigned int icounter ) {
00073
00074 std::cout << setiosflags(ios::left) << setw(10) << icounter << setw(8) << data->size() << setw(30) << name() << std::endl;
00075
00076 if( data->size() > 0 ) {
00077 std::vector<Record*>::iterator iter = data->begin();
00078 for( ; iter!=data->end(); iter++ ) {
00079 fillRecord( (*iter) );
00080 }
00081 }
00082
00083 if( icounter==(NBarrel-1) ) {
00084 fitHistogram();
00085 fillGraph();
00086 fitGraph();
00087 }
00088
00089 return;
00090 }
00091
00092
00093
00094 void calib_barrel_common::fillRecord( const Record* r ) {
00095 double zhit = r->zrhit();
00096 if( zhit<zbegin || zhit>zend ) return;
00097 int zbin = static_cast<int>((zhit-zbegin)/zstep);
00098 if( ( zbin<0 ) || ( zbin>static_cast<int>(nBinPerCounter-1) ) ) {
00099 cout << "tofcalgsec::calib_barrel_common: zhit is out of range, zhit=" << zhit << " zbin=" << zbin << endl;
00100 return;
00101 }
00102
00103 std::vector<TH1F*>::iterator iter = m_histograms.begin();
00104 (*(iter+zbin))->Fill( r->tleft() );
00105 (*(iter+nBinPerCounter+zbin))->Fill( r->tright() );
00106 (*(iter+2*nBinPerCounter+zbin))->Fill( (r->tleft()+r->tright())/2.0 );
00107 (*(iter+3*nBinPerCounter+zbin))->Fill( (r->tleft()-r->tright())/2.0 );
00108 (*(iter+nKind*nBinPerCounter))->Fill( r->tleft() );
00109 (*(iter+nKind*nBinPerCounter))->Fill( r->tright() );
00110
00111 return;
00112 }
00113
00114
00115 void calib_barrel_common::fitHistogram() {
00116 TF1* g = new TF1("g", "gaus");
00117 g->SetLineColor(2);
00118 g->SetLineWidth(1);
00119
00120 std::vector<TH1F*>::iterator iter1 = m_histograms.begin();
00121 std::vector<HepVector>::iterator iter2 = m_fitresult.begin();
00122 for( ; iter1!=m_histograms.end(); iter1++, iter2++ ) {
00123 (*iter1)->Fit( g, "Q");
00124 (*iter2)[0] = g->GetParameter(1);
00125 (*iter2)[1] = g->GetParError(1);
00126 (*iter2)[2] = g->GetParameter(2);
00127 (*iter2)[3] = g->GetParError(2);
00128 }
00129
00130 iter2 = m_fitresult.end() - 1;
00131 X[2] = (*iter2)[0];
00132 X[3] = (*iter2)[1];
00133
00134 return;
00135
00136 }
00137
00138
00139 void calib_barrel_common::fillGraph() {
00140
00141 char gname1[256], gname2[256];
00142 TH1F *gra[nGraphTotalCommon];
00143
00144
00145
00146
00147
00148
00149
00150
00151 std::vector<double> toffset, toffseterr;
00152 std::vector<double> tsigma, tsigmaerr;
00153 toffset.resize( nBinPerCounter );
00154 toffseterr.resize( nBinPerCounter );
00155 tsigma.resize( nBinPerCounter );
00156 tsigmaerr.resize( nBinPerCounter );
00157
00158 unsigned int number = 0;
00159 std::vector<HepVector>::iterator iter = m_fitresult.begin();
00160 for( unsigned int j=0; j<nKind; j++ ) {
00161 if( j==0 ) {
00162 sprintf( gname1, "tlefttright-offset" );
00163 sprintf( gname2, "tlefttright-sigma" );
00164 }
00165 else if( j==1 ) {
00166 sprintf( gname1, "tcommon-offset" );
00167 sprintf( gname2, "tcommon-sigma" );
00168 }
00169 else if( j==2 ) {
00170 sprintf( gname1, "tplusminus-offset" );
00171 sprintf( gname2, "tplusminus-sigma" );
00172 }
00173 else if( j==3 ) {
00174 sprintf( gname1, "tcorrelation-offset" );
00175 sprintf( gname2, "tcorrelation-sigma" );
00176 }
00177
00178 gra[j] = new TH1F( gname1, gname1, nBinPerCounter, zbegin, zend );
00179 gra[j+4] = new TH1F( gname2, gname2, nBinPerCounter, zbegin, zend );
00180
00181 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00182 number = j*nBinPerCounter + k;
00183 toffset[k] = (*(iter+number))[0];
00184 toffseterr[k] = (*(iter+number))[1];
00185 tsigma[k] = (*(iter+number))[2];
00186 tsigmaerr[k] = (*(iter+number))[3];
00187 gra[j]->SetBinContent( k+1, toffset[k] );
00188 gra[j]->SetBinError( k+1, toffseterr[k] );
00189 gra[j+4]->SetBinContent( k+1, tsigma[k] );
00190 gra[j+4]->SetBinError( k+1, tsigmaerr[k] );
00191 }
00192 }
00193
00194 gra[nGraphTotalCommon-1] = new TH1F( "Sigma", "Sigma", nBinPerCounter, zbegin, zend );
00195 for( unsigned int k=0; k<nBinPerCounter; k++ ) {
00196 number = (nKind-1)*nBinPerCounter + k;
00197 double sigPlus = (*(iter+number-nBinPerCounter))[2];
00198 double sigMinus = (*(iter+number))[2];
00199 double sigErrPlus = (*(iter+number-nBinPerCounter))[3];
00200 double sigErrMinus = (*(iter+number))[3];
00201 if( sigPlus>sigMinus ) {
00202 tsigma[k] = sqrt( sigPlus*sigPlus - sigMinus*sigMinus );
00203 }
00204 else {
00205 tsigma[k] = 0.0 - sqrt( sigMinus*sigMinus - sigPlus*sigPlus );
00206 }
00207 tsigmaerr[k] = sqrt( sigErrPlus*sigErrPlus + sigErrMinus*sigErrMinus );
00208
00209 gra[nGraphTotalCommon-1]->SetBinContent( k+1, tsigma[k] );
00210 gra[nGraphTotalCommon-1]->SetBinError( k+1, tsigmaerr[k] );
00211 }
00212
00213 for( int j=0; j<nGraphTotalCommon; j++, j++ ) {
00214 gra[j]->SetMarkerSize(1.5);
00215 gra[j]->SetMarkerStyle(20);
00216 gra[j]->SetMarkerColor(2);
00217 if( j==4 ) {
00218 gra[j]->SetMaximum( 0.22 );
00219 gra[j]->SetMinimum( 0.07 );
00220 }
00221 else if( j==6 ) {
00222 gra[j]->SetMaximum( 0.20 );
00223 gra[j]->SetMinimum( -0.02 );
00224 }
00225 }
00226 for( int j=1; j<nGraphTotalCommon; j++, j++ ) {
00227 gra[j]->SetMarkerSize(1.5);
00228 gra[j]->SetMarkerStyle(21);
00229 gra[j]->SetMarkerColor(4);
00230 }
00231 gra[nGraphTotalCommon-1]->SetMarkerStyle(4);
00232 gra[nGraphTotalCommon-1]->SetMarkerColor(6);
00233
00234 for( int j=0; j<nGraphTotalCommon; j++ ) {
00235 m_graphs.push_back( gra[j] );
00236 }
00237
00238 return;
00239 }
00240
00241
00242 void calib_barrel_common::fitGraph() {
00243 TF1* p0 = new TF1("p0", "pol0");
00244 p0->SetLineColor(1);
00245 p0->SetLineWidth(1);
00246
00247 std::vector<TH1F*>::iterator iter=m_graphs.end()-1;
00248 (*iter)->Fit( "p0", "Q" );
00249 X[0] = p0->GetParameter(0);
00250 X[1] = p0->GetParError(0);
00251
00252 m_result.push_back( X );
00253 return;
00254 }