00001 #include "tofcalgsec/calib_endcap_atten.h"
00002 #include "TF1.h"
00003
00004
00005 static double endcapQFunc(double *x, double *par) {
00006 double xx = x[0];
00007 double f = par[0] + par[1]*(xx-44.5) + par[2]*(xx-44.5)*(xx-44.5);
00008
00009 return f;
00010 }
00011
00012
00013 calib_endcap_atten::calib_endcap_atten( const unsigned int nrbin ):TofCalibFit( false, nEndcapAtten ) {
00014
00015 nKind = 1;
00016 nBinPerCounter = nrbin + 1;
00017
00018 nHistPerCounter = nKind*nBinPerCounter;
00019 nCanvasPerCounter = 2;
00020 CanvasPerCounterName.push_back( static_cast<string>("Pulse Height Most Probable Value") );
00021 CanvasPerCounterName.push_back( static_cast<string>("Pulse Height Sigma") );
00022 nGraphPerCanvasPerCounter.push_back(1);
00023 nGraphPerCanvasPerCounter.push_back(1);
00024
00025 nHistogram = 0;
00026 nCanvas = 2;
00027 CanvasName.push_back( static_cast<string>("Pulse Height Most Probable Value vs TOF Counter Number") );
00028 CanvasName.push_back( static_cast<string>("Pulse Height Sigma vs TOF Counter Number") );
00029 nGraphPerCanvas.push_back(1);
00030 nGraphPerCanvas.push_back(1);
00031
00032 int numGraphs1 = 0;
00033 std::vector<unsigned int>::iterator iter = nGraphPerCanvasPerCounter.begin();
00034 for( ; iter!=nGraphPerCanvasPerCounter.end(); iter++ ) {
00035 numGraphs1 = numGraphs1 + (*iter);
00036 }
00037 if( numGraphs1 != nGraphEcAtten ) {
00038 cout << "tofcalgsec::calib_endcap_atten: the number of Graphs is NOT reasonable!!!" << endl;
00039 exit(0);
00040 }
00041 int numGraphs2 = 0;
00042 iter = nGraphPerCanvas.begin();
00043 for( ; iter!=nGraphPerCanvas.end(); iter++ ) {
00044 numGraphs2 = numGraphs2 + (*iter);
00045 }
00046 if( numGraphs2 != nGraphEcAtten ) {
00047 cout << "tofcalgsec::calib_endcap_atten: the number of Graphs is NOT reasonable!!!" << endl;
00048 exit(0);
00049 }
00050
00051 m_name = string("calib_endcap_atten");
00052
00053 const int qbin = 100;
00054 const double qbegin = 0.0;
00055 const double qend = 5000.0;
00056
00057
00058 char hname[256];
00059 for( unsigned int i=0; i<NEndcap; i++ ) {
00060 m_result.push_back( HepVector(nEndcapAtten,0) );
00061 for( unsigned int k=0; k<nrbin; k++ ) {
00062 sprintf( hname, "Q-id%i-r%i", i, k);
00063 m_histograms.push_back( new TH1F( hname, hname, qbin, qbegin, qend ) );
00064 m_fitresult.push_back( HepVector(nParEcAtten,0) );
00065 }
00066 sprintf( hname, "Q0-id%i", i );
00067 m_histograms.push_back( new TH1F( hname, hname, qbin, qbegin, qend ) );
00068 m_fitresult.push_back( HepVector(nParEcAtten,0) );
00069 }
00070
00071 rpos.resize( nrbin );
00072 rposerr.resize( nrbin );
00073 rstep = ( rend - rbegin )/nrbin;
00074 for( unsigned int i=0; i<nrbin; i++ ) {
00075 rpos[i] = rbegin + ( i+0.5 )*rstep;
00076 rposerr[i] = 0.5*rstep;
00077 }
00078 itofid.resize( NEndcap );
00079 itofiderr.resize( NEndcap );
00080 itofidstep = 1.0;
00081 for( unsigned int i=0; i<NEndcap; i++ ) {
00082 itofid[i] = i*1.0;
00083 itofiderr[i] = 0.5;
00084 }
00085
00086 }
00087
00088
00089 calib_endcap_atten::~calib_endcap_atten() {
00090 m_fitresult.clear();
00091 rpos.clear();
00092 rposerr.clear();
00093 itofid.clear();
00094 itofiderr.clear();
00095 }
00096
00097
00098 void calib_endcap_atten::calculate( RecordSet*& data, unsigned int icounter ) {
00099
00100 std::cout << setiosflags(ios::left) << setw(10) << icounter << setw(8) << data->size() << setw(30) << name() << std::endl;
00101
00102 if( data->size() > 0 ) {
00103 std::vector<Record*>::iterator iter = data->begin();
00104 for( ; iter!=data->end(); iter++ ) {
00105 fillRecord( (*iter), icounter );
00106 }
00107 fitHistogram( icounter );
00108 fillGraph( icounter );
00109 fitGraph( icounter );
00110 }
00111 else {
00112 fillGraph( icounter );
00113 }
00114
00115 if( data->size() > 0 ) {
00116 std::vector<Record*>::iterator iter = data->begin();
00117 for( ; iter!=data->end(); iter++ ) {
00118 updateData( (*iter), icounter );
00119 fillRecordQ0( (*iter), icounter );
00120 }
00121 fitHistogramQ0( icounter );
00122 }
00123
00124 if( icounter==(NEndcap-1) ) {
00125 fillGraphQ0();
00126 }
00127
00128 return;
00129 }
00130
00131
00132 void calib_endcap_atten::fillRecord( const Record* r, unsigned int icounter ) {
00133
00134 double rhit = r->zrhit();
00135 if( rhit<rbegin || rhit>rend ) return;
00136 int rbin = static_cast<int>((rhit-rbegin)/rstep);
00137 if( rbin<0 ) { rbin = 0; }
00138 else if( rbin>static_cast<int>(nBinPerCounter-1-1) ) {
00139 cout << "tofcalgsec::calib_endcap_atten:fillRecord: rhit is out of range, rhit=" << rhit << " rbin=" << rbin << endl;
00140 return;
00141 }
00142
00143 std::vector<TH1F*>::iterator iter = m_histograms.begin() + icounter*nKind*nBinPerCounter + rbin;
00144 (*iter)->Fill( r->qleft()*abs(r->theta()) );
00145
00146 return;
00147 }
00148
00149
00150 void calib_endcap_atten::fitHistogram( unsigned int icounter ) {
00151 TF1* ld = new TF1("ld", "landau");
00152 ld->SetLineColor(2);
00153 ld->SetLineWidth(1);
00154
00155 std::vector<TH1F*>::iterator iter1 = m_histograms.begin() + icounter*nKind*nBinPerCounter;
00156 std::vector<HepVector>::iterator iter2 = m_fitresult.begin() + icounter*nKind*nBinPerCounter;
00157 for( unsigned int j=0; j<nBinPerCounter-1; j++, iter1++, iter2++ ) {
00158 (*iter1)->Fit( ld, "Q");
00159 (*iter2)[0] = ld->GetParameter(1);
00160 (*iter2)[1] = ld->GetParError(1);
00161 (*iter2)[2] = ld->GetParameter(2);
00162 (*iter2)[3] = ld->GetParError(2);
00163 }
00164
00165 return;
00166
00167 }
00168
00169
00170 void calib_endcap_atten::fillGraph( unsigned int icounter ) {
00171
00172 char gname1[256], gname2[256];
00173
00174
00175
00176
00177
00178 std::vector<double> toffset, toffseterr;
00179 std::vector<double> tsigma, tsigmaerr;
00180 toffset.resize( nBinPerCounter-1 );
00181 toffseterr.resize( nBinPerCounter-1 );
00182 tsigma.resize( nBinPerCounter-1 );
00183 tsigmaerr.resize( nBinPerCounter-1 );
00184
00185 std::vector<HepVector>::iterator iter = m_fitresult.begin() + icounter*nKind*nBinPerCounter;
00186 for( unsigned int k=0; k<nBinPerCounter-1; k++ ) {
00187 toffset[k] = log((*(iter+k))[0]);
00188 toffseterr[k] = log((*(iter+k))[0])*((*(iter+k))[1])/((*(iter+k))[0]);
00189 tsigma[k] = (*(iter+k))[2];
00190 tsigmaerr[k] = (*(iter+k))[3];
00191 }
00192
00193 sprintf( gname1, "Q MPV-tofid-%i", icounter );
00194 m_graphs.push_back( new TH1F( gname1, gname1, nBinPerCounter-1, rbegin, rend ) );
00195 std::vector<TH1F*>::iterator itgraph = m_graphs.end() - 1;
00196 (*itgraph)->SetMarkerSize(1.5);
00197 (*itgraph)->SetMarkerStyle(20);
00198 (*itgraph)->SetMarkerColor(2);
00199 for( unsigned int k=0; k<nBinPerCounter-1; k++ ) {
00200 (*itgraph)->SetBinContent( k+1, toffset[k] );
00201 (*itgraph)->SetBinError( k+1, toffseterr[k] );
00202 }
00203
00204 sprintf( gname2, "Q sigma-tofid-%i", icounter );
00205 m_graphs.push_back( new TH1F( gname2, gname2, nBinPerCounter-1, rbegin, rend ) );
00206 itgraph = m_graphs.end() - 1;
00207 (*itgraph)->SetMarkerSize(1.5);
00208 (*itgraph)->SetMarkerStyle(21);
00209 (*itgraph)->SetMarkerColor(4);
00210 for( unsigned int k=0; k<nBinPerCounter-1; k++ ) {
00211 (*itgraph)->SetBinContent( k+1, tsigma[k] );
00212 (*itgraph)->SetBinError( k+1, tsigmaerr[k] );
00213 }
00214
00215 return;
00216 }
00217
00218
00219 void calib_endcap_atten::fitGraph( unsigned int icounter ) {
00220
00221 TF1 *fsingleq = new TF1( "fsingleq", endcapQFunc, rbegin, rend, 3 );
00222 fsingleq->SetLineColor(1);
00223 fsingleq->SetLineWidth(1);
00224 fsingleq->SetParameters( 6.5, 0.0, 0.0 );
00225
00226 std::vector<TH1F*>::iterator itgraph = m_graphs.begin() + icounter*nGraphEcAtten;
00227
00228 (*itgraph)->Fit( "fsingleq", "Q", "", rbegin, rend );
00229 X = HepVector( m_npar, 0 );
00230 X[0] = fsingleq->GetParameter(0);
00231 X[1] = fsingleq->GetParameter(1);
00232 X[2] = fsingleq->GetParameter(2);
00233 X[3] = 0.;
00234 X[4] = 0.;
00235 X[5] = 0.;
00236
00237 std::vector<HepVector>::iterator iter = m_result.begin() + icounter;
00238 (*iter) = X;
00239
00240 return;
00241 }
00242
00243
00244 void calib_endcap_atten::updateData( Record* r, unsigned int icounter ) {
00245 double rhit = r->zrhit();
00246 double q = r->qleft();
00247 double costheta = abs(r->theta());
00248
00249 double par[3];
00250 std::vector<HepVector>::iterator iter = m_result.begin() + icounter;
00251 for( unsigned int i=0; i<3; i++ ) {
00252 par[i] = (*iter)[i];
00253 }
00254
00255 par[0] = 0.;
00256 double q0 = q*costheta/exp(endcapQFunc(&rhit,par));
00257 r->setQ0( q0 );
00258
00259 return;
00260 }
00261
00262
00263 void calib_endcap_atten::fillRecordQ0( const Record* r, unsigned int icounter ) {
00264 std::vector<TH1F*>::iterator iter = m_histograms.begin() + icounter*nKind*nBinPerCounter + nBinPerCounter - 1;
00265 (*iter)->Fill( r->q0() );
00266
00267 return;
00268 }
00269
00270
00271 void calib_endcap_atten::fitHistogramQ0( unsigned int icounter ) {
00272 TF1* ld = new TF1("ld", "landau");
00273 ld->SetLineColor(2);
00274 ld->SetLineWidth(1);
00275
00276 std::vector<TH1F*>::iterator iter1 = m_histograms.begin() + icounter*nKind*nBinPerCounter + nBinPerCounter - 1;
00277 std::vector<HepVector>::iterator iter2 = m_fitresult.begin() + icounter*nKind*nBinPerCounter + nBinPerCounter - 1;
00278 (*iter1)->Fit( ld, "Q");
00279 (*iter2)[0] = ld->GetParameter(1);
00280 (*iter2)[1] = ld->GetParError(1);
00281 (*iter2)[2] = ld->GetParameter(2);
00282 (*iter2)[3] = ld->GetParError(2);
00283
00284 return;
00285 }
00286
00287
00288 void calib_endcap_atten::fillGraphQ0() {
00289 char gname1[256], gname2[256];
00290
00291
00292
00293
00294
00295 std::vector<double> toffset, toffseterr;
00296 std::vector<double> tsigma, tsigmaerr;
00297 toffset.resize( NEndcap );
00298 toffseterr.resize( NEndcap );
00299 tsigma.resize( NEndcap );
00300 tsigmaerr.resize( NEndcap );
00301
00302 unsigned int number = 0;
00303 std::vector<HepVector>::iterator iter1 = m_fitresult.begin() + nBinPerCounter - 1;
00304 std::vector<HepVector>::iterator iter2 = m_result.begin();
00305 for( unsigned int i=0; i<NEndcap; i++ ) {
00306 number = i*nBinPerCounter;
00307 toffset[i] = (*(iter1+number))[0];
00308 toffseterr[i] = (*(iter1+number))[1];
00309 tsigma[i] = (*(iter1+number))[2];
00310 tsigmaerr[i] = (*(iter1+number))[3];
00311
00312
00313 (*(iter2+i))[3] = toffset[i]/toffset[0];
00314 (*(iter2+i))[4] = toffset[i];
00315 }
00316
00317 sprintf( gname1, "Q0 MPV vs TOF Counter Number" );
00318 m_graphs.push_back( new TH1F( gname1, gname1, NEndcap, -0.5, NEndcap-0.5 ) );
00319 std::vector<TH1F*>::iterator itgraph = m_graphs.end() - 1;
00320 (*itgraph)->SetMarkerSize(1.5);
00321 (*itgraph)->SetMarkerStyle(20);
00322 (*itgraph)->SetMarkerColor(2);
00323 for( unsigned int k=0; k<nBinPerCounter-1; k++ ) {
00324 (*itgraph)->SetBinContent( k+1, tsigma[k] );
00325 (*itgraph)->SetBinError( k+1, tsigmaerr[k] );
00326 }
00327 for( unsigned int i=0; i<NEndcap; i++ ) {
00328 (*itgraph)->SetBinContent( i+1, toffset[i] );
00329 (*itgraph)->SetBinError( i+1, toffseterr[i] );
00330 }
00331
00332 sprintf( gname2, "Q0 Sigma vs TOF Counter Number" );
00333 m_graphs.push_back( new TH1F( gname2, gname2, NEndcap, -0.5, NEndcap-0.5 ) );
00334 itgraph = m_graphs.end() - 1;
00335 (*itgraph)->SetTitle(gname2);
00336 (*itgraph)->SetMarkerSize(1.5);
00337 (*itgraph)->SetMarkerStyle(21);
00338 (*itgraph)->SetMarkerColor(4);
00339 for( unsigned int i=0; i<NEndcap; i++ ) {
00340 (*itgraph)->SetBinContent( i+1, tsigma[i] );
00341 (*itgraph)->SetBinError( i+1, tsigmaerr[i] );
00342 }
00343
00344 return;
00345 }