- Timestamp:
- 10/26/12 09:23:07 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/rootmacros/PulseTemplates/pulse.C
r14479 r14532 124 124 Pulse::InitMembers() 125 125 { 126 mName = "pulse"; 127 mHisto = NULL; 128 mModel = NULL; 129 mOptions = "S"; 130 mBsl = 0; 131 mHeight = 0; 132 mT0 = 0; 133 mT1 = 0; 134 mTau1 = 0; 135 mTau2 = 0; 136 mIntegral = 0; 137 mAmplitude = 0; 138 mPhE = 0; 139 mType = 0; 140 mFitProb = 0; 141 mFitNCalls = 0; 142 mFitNdf = 0; 143 mChi2 = 0; 126 mName = "pulse"; 127 mHisto = NULL; 128 mModel = NULL; 129 mOptions = "S"; 130 mBsl.first = 0; 131 mHeight.first = 0; 132 mT0.first = 0; 133 mT1.first = 0; 134 mTau1.first = 0; 135 mTau2.first = 0; 136 mIntegral.first = 0; 137 mAmplitude.first = 0; 138 mPhE.first = 0; 139 mBsl.second = 0; 140 mHeight.second = 0; 141 mT0.second = 0; 142 mT1.second = 0; 143 mTau1.second = 0; 144 mTau2.second = 0; 145 mIntegral.second = 0; 146 mAmplitude.second = 0; 147 mPhE.second = 0; 148 mOrder = 0; 149 mType = 0; 150 mFitProb = 0; 151 mFitNCalls = 0; 152 mFitNdf = 0; 153 mChi2 = 0; 144 154 } 145 155 … … 154 164 mFitMin = fitMin; 155 165 mFitMax = fitMax; 156 mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, 6);166 mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, 7 ); 157 167 158 168 // ====================================================================== … … 168 178 bsl = bsl/30; 169 179 170 double stop = mHisto->GetMaximumBin(); //pos of max 171 double height = mHisto->GetBinContent(stop); 172 double start = stop-10; //pos 10 slices before maximum 180 double stop = mHisto->GetMaximumBin(); //pos of max 181 double amplitude = mHisto->GetBinContent(stop); 182 double start = stop-10; //pos 10 slices before maximum 183 int phe = amplitude/10; 173 184 // ====================================================================== 174 185 175 double para[] = {bsl, height, start, stop, tau, tau};186 double para[] = {bsl, amplitude, start, stop, tau, tau, phe}; 176 187 177 188 178 189 179 190 mModel->SetParameters(para); 180 mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2" );191 mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2", "PhE"); 181 192 mModel->SetLineColor(kRed); 182 193 183 mFitResultPtr = mHisto->Fit(mModel, fitOptions); 184 185 mBsl = mModel->GetParameter(0); 186 mHeight = mModel->GetParameter(1); 187 mT0 = mModel->GetParameter(2); 188 mT1 = mModel->GetParameter(3); 189 mTau1 = mModel->GetParameter(4); 190 mTau2 = mModel->GetParameter(5); 194 mFitResultPtr = mHisto->Fit(mModel, fitOptions); 195 mBsl.first = mModel->GetParameter(0); 196 mBsl.second = mModel->GetParError(0); 197 mHeight.first = mModel->GetParameter(1); 198 mHeight.second = mModel->GetParError(1); 199 mT0.first = mModel->GetParameter(2); 200 mT0.second = mModel->GetParError(2); 201 mT1.first = mModel->GetParameter(3); 202 mT1.second = mModel->GetParameter(3); 203 mTau1.first = mModel->GetParameter(4); 204 mTau1.second = mModel->GetParError(4); 205 mTau2.first = mModel->GetParameter(5); 206 mTau2.second = mModel->GetParError(5); 207 mPhE.first = mModel->GetParameter(6); 208 mPhE.second = mModel->GetParError(6); 191 209 } 192 210 … … 201 219 mFitMin = fitMin; 202 220 mFitMax = fitMax; 203 mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, 5);221 mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, 6 ); 204 222 205 223 // ====================================================================== … … 218 236 double amplitude= mHisto->GetBinContent(stop); 219 237 double start = stop-10; //pos 10 slices before maximum 238 int phe = amplitude/10; 220 239 // ====================================================================== 221 240 222 241 223 double para[] = {bsl, amplitude, start, tau, tau };242 double para[] = {bsl, amplitude, start, tau, tau, phe}; 224 243 mModel->SetParameters(para); 225 mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2" );244 mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2", "PhE"); 226 245 mModel->SetLineColor(kBlue); 227 246 228 mFitResultPtr = mHisto->Fit(mModel, fitOptions); 229 230 mBsl = mModel->GetParameter(0); 231 mHeight = mModel->GetParameter(1); 232 mT0 = mModel->GetParameter(2); 233 mTau1 = mModel->GetParameter(3); 234 mTau2 = mModel->GetParameter(4); 247 mFitResultPtr = mHisto->Fit(mModel, fitOptions); 248 mBsl.first = mModel->GetParameter(0); 249 mBsl.second = mModel->GetParError(0); 250 mHeight.first = mModel->GetParameter(1); 251 mHeight.second = mModel->GetParError(1); 252 mT0.first = mModel->GetParameter(2); 253 mT0.second = mModel->GetParError(2); 254 mTau1.first = mModel->GetParameter(3); 255 mTau1.second = mModel->GetParError(3); 256 mTau2.first = mModel->GetParameter(4); 257 mTau2.second = mModel->GetParError(4); 258 mPhE.first = mModel->GetParameter(5); 259 mPhE.second = mModel->GetParError(5); 260 261 cout << "test" << mBsl.first << "\t +/-" << mBsl.second << endl; 235 262 } 236 263 … … 238 265 Pulse::CalculateParameters() 239 266 { 240 mIntegral = mModel->Integral(mFitMin, mFitMax);241 mAmplitude = mModel->GetMaximum() - mBsl;242 mFitProb = mFitResultPtr->Prob();243 mFitNCalls = mFitResultPtr->NCalls();244 mFitNdf = mFitResultPtr->Ndf();245 mChi2 = mFitResultPtr->Chi2();267 mIntegral.first = mModel->Integral(mFitMin, mFitMax); 268 mAmplitude.first = mModel->GetMaximum() - mBsl.first; 269 mFitProb = mFitResultPtr->Prob(); 270 mFitNCalls = mFitResultPtr->NCalls(); 271 mFitNdf = mFitResultPtr->Ndf(); 272 mChi2 = mFitResultPtr->Chi2(); 246 273 } 247 274 … … 251 278 252 279 TString Pulse::GetName(){ return mName;} 253 double Pulse::GetBsl(){ return mBsl;} 254 double Pulse::GetHeight(){ return mHeight;} 255 double Pulse::GetT0(){ return mT0;} 256 double Pulse::GetT1(){ return mT1;} 257 double Pulse::GetTau1(){ return mTau1;} 258 double Pulse::GetTau2(){ return mTau2;} 259 double Pulse::GetIntegral(){ return mIntegral;} 260 double Pulse::GetAmplitude(){ return mAmplitude;} 261 int Pulse::GetPhe(){ return mPhE;} 280 double Pulse::GetBsl(){ return mBsl.first;} 281 double Pulse::GetHeight(){ return mHeight.first;} 282 double Pulse::GetT0(){ return mT0.first;} 283 double Pulse::GetT1(){ return mT1.first;} 284 double Pulse::GetTau1(){ return mTau1.first;} 285 double Pulse::GetTau2(){ return mTau2.first;} 286 double Pulse::GetIntegral(){ return mIntegral.first;} 287 double Pulse::GetAmplitude(){ return mAmplitude.first;} 288 int Pulse::GetPhE(){ return mPhE.first;} 289 double Pulse::GetBslErr(){ return mBsl.second;} 290 double Pulse::GetHeightErr(){ return mHeight.second;} 291 double Pulse::GetT0Err(){ return mT0.second;} 292 double Pulse::GetT1Err(){ return mT1.second;} 293 double Pulse::GetTau1Err(){ return mTau1.second;} 294 double Pulse::GetTau2Err(){ return mTau2.second;} 295 double Pulse::GetIntegralErr(){ return mIntegral.second;} 296 double Pulse::GetAmplitudeErr(){ return mAmplitude.second;} 297 int Pulse::GetPhEErr(){ return mPhE.second;} 298 int Pulse::GetOrder(){ return mOrder;} 262 299 int Pulse::GetType(){return mType;} 263 300 int Pulse::GetFitMin(){return mFitMin;} … … 293 330 double start = par[2]; 294 331 // double rising = par[3]; 332 double stop = par[3]; 295 333 double tau1 = par[4]; 296 334 double tau2 = par[5]; 297 double stop = par[3];335 int phe = par[6]; 298 336 299 337 // helper variables … … 305 343 returnval += Heaviside(t[0]-start)*e1; 306 344 returnval += Heaviside(t[0]-stop)*e2; 307 returnval *= height;345 returnval *= phe*height; 308 346 returnval += bsl; 309 347 // if (t[0] > start) … … 362 400 double tau1 = par[3]; 363 401 double tau2 = par[4]; 402 int phe = par[5]; 364 403 365 404 // helper variables … … 369 408 // calculate return value 370 409 returnval += bsl; 371 returnval += gain*e1*e2*Heaviside(t[0]-t_0);410 returnval += phe*gain*e1*e2*Heaviside(t[0]-t_0); 372 411 373 412 return returnval;
Note:
See TracChangeset
for help on using the changeset viewer.