Changeset 7666 for trunk/MagicSoft/Mars
- Timestamp:
- 04/30/06 15:41:17 (19 years ago)
- Location:
- trunk/MagicSoft/Mars/mjtrain
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc
r7665 r7666 141 141 142 142 c.cd(1); 143 gPad->SetBorderMode(0); 144 gPad->SetFrameBorderMode(0); 143 145 gPad->SetLogx(); 144 146 h.DrawCopy(); … … 151 153 152 154 c.cd(2); 155 gPad->SetBorderMode(0); 156 gPad->SetFrameBorderMode(0); 153 157 gPad->SetLogx(); 154 158 MH::SetPalette("pretty"); … … 156 160 157 161 c.cd(4); 162 gPad->SetBorderMode(0); 163 gPad->SetFrameBorderMode(0); 158 164 gPad->SetLogx(); 159 165 MH::SetPalette("pretty"); … … 265 271 } 266 272 267 Double_t MJTrainSeparation::GetDataRate(MDataSet &set ) const273 Double_t MJTrainSeparation::GetDataRate(MDataSet &set, Double_t &num) const 268 274 { 269 275 TChain chain1("Events"); 270 276 set.AddFilesOff(chain1); 271 277 278 num = chain1.GetEntries(); 279 if (num<100) 280 { 281 *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl; 282 return -1; 283 } 284 285 TChain chain("EffectiveOnTime"); 286 set.AddFilesOff(chain); 287 288 chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff"); 289 290 TH1 *h = dynamic_cast<TH1*>(gROOT->FindObject("htemp")); 291 if (!h) 292 { 293 *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl; 294 return -1; 295 } 296 297 const Double_t ontime = h->Integral(); 298 delete h; 299 300 if (ontime<1) 301 { 302 *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl; 303 return -1; 304 } 305 306 *fLog << inf << "Found " << num << " background events in " << ontime << "s" << endl; 307 308 return num/ontime; 309 } 310 311 Double_t MJTrainSeparation::GetNumMC(MDataSet &set) const 312 { 313 TChain chain1("Events"); 314 set.AddFilesOn(chain1); 315 272 316 const Double_t num = chain1.GetEntries(); 273 317 if (num<100) … … 277 321 } 278 322 279 TChain chain("EffectiveOnTime"); 280 set.AddFilesOff(chain); 281 282 chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff"); 283 284 TH1 *h = dynamic_cast<TH1*>(gROOT->FindObject("htemp")); 285 if (!h) 286 { 287 *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl; 288 return -1; 289 } 290 291 const Double_t ontime = h->Integral(); 292 delete h; 293 294 if (ontime<1) 295 { 296 *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl; 297 return -1; 298 } 299 300 *fLog << inf << "Found " << num << " background events in " << ontime << "s" << endl; 301 302 return num/ontime; 323 return num; 303 324 } 304 325 … … 312 333 313 334 // Target spectrum 314 TF1 flx("Flux", "[0] *(x/1000)^(-2.6)", min, max);315 flx.SetParameter(0, 1e- 7);335 TF1 flx("Flux", "[0]/1000*(x/1000)^(-2.6)", min, max); 336 flx.SetParameter(0, 1e-6); 316 337 317 338 // Number n0 of events this spectrum would produce per s and m^2 … … 324 345 const Double_t R = n0*A; //[Hz] 325 346 347 *fLog << "Gamma rate from the source inside the MC production area: " << R << "Hz" << endl; 348 326 349 // Number N of events produced (in trainings sample) 327 350 const Double_t N = num; //[#] 328 351 352 *fLog << "Events produced by MC inside the production area: " << num << endl; 353 329 354 // This correponds to an observation time T [s] 330 355 const Double_t T = N/R; //[s] 331 356 357 *fLog << "Total time produced by the Monte Carlo: " << T << "s" << endl; 358 332 359 // With an average data rate after star of 333 const Double_t r = GetDataRate(fDataSetTrain); //[Hz] 360 Double_t data=0; 361 const Double_t r = GetDataRate(fDataSetTrain, data); //[Hz] 362 363 *fLog << "Events measured per second effective on time: " << r << "Hz" << endl; 334 364 335 365 // this yields a number of n events to be read for training 336 366 const Double_t n = r*T; //[#] 337 367 368 *fLog << "Events to be read from the data sample: " << n << endl; 369 *fLog << "Events available in data sample: " << data << endl; 370 338 371 if (r<0) 339 372 return kFALSE; 340 373 341 *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl; 342 *fLog << "For a data rate of " << r << "Hz this corresponds to " << TMath::Nint(n) << " data events." << endl; 343 344 fNumTrainOn = (UInt_t)-1; 345 fNumTrainOff = TMath::Nint(n); 374 Double_t nummc = GetNumMC(fDataSetTrain); 375 376 *fLog << "Events available in MC sample: " << nummc << endl; 377 378 *fLog << "MC read probability: " << data/n << endl; 379 380 // more data requested than available => Scale down num MC events 381 Double_t on, off; 382 if (data<n) 383 { 384 on = TMath::Nint(nummc*data/n); //(UInt_t)-1; 385 off = TMath::Nint(data); 386 *fLog << "Not enough data events available... scaling by " << data/n << endl; 387 } 388 else 389 { 390 on = TMath::Nint(nummc); 391 off = TMath::Nint(n); 392 } 393 394 if (fNumTrainOn>0 && fNumTrainOn<on) 395 { 396 fNumTrainOff = TMath::Nint(off*fNumTrainOn/on); 397 *fLog << "Less MC events requested... scaling by " << fNumTrainOn/on << endl; 398 } 399 else 400 { 401 fNumTrainOn = TMath::Nint(on); 402 fNumTrainOff = TMath::Nint(off); 403 } 404 405 *fLog << "Target number of MC events: " << fNumTrainOn << endl; 406 *fLog << "Target number of data events: " << fNumTrainOff << endl; 346 407 347 408 /* … … 359 420 MF f("f * MEventRate.fRate < rand"); 360 421 */ 361 362 422 363 423 return kTRUE; -
trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h
r7665 r7666 30 30 31 31 Bool_t GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const; 32 Double_t GetDataRate(MDataSet &set) const; 32 Double_t GetDataRate(MDataSet &set, Double_t &num) const; 33 Double_t GetNumMC(MDataSet &set) const; 33 34 Bool_t AutoTrain(); 34 35
Note:
See TracChangeset
for help on using the changeset viewer.