Changeset 7539 for trunk/MagicSoft/Mars
- Timestamp:
- 02/27/06 14:06:48 (19 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r7538 r7539 63 63 - added GetForest 64 64 65 * mjtrain/MJTrainSeparation.[h,cc]: 66 - addedsome code for upcomming automatic event selection 67 65 68 66 69 -
trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc
r7420 r7539 18 18 ! Author(s): Thomas Bretz 11/2005 <mailto:tbretz@astro.uni-wuerzburg.de> 19 19 ! 20 ! Copyright: MAGIC Software Development, 200 520 ! Copyright: MAGIC Software Development, 2006 21 21 ! 22 22 ! … … 30 30 #include "MJTrainSeparation.h" 31 31 32 #include <TF1.h> 32 33 #include <TH2.h> 34 #include <TChain.h> 33 35 #include <TGraph.h> 34 36 #include <TVirtualPad.h> … … 136 138 } 137 139 140 /* 141 Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const 142 { 143 fLog->Separator("Initialize energy weighting"); 144 145 if (!CheckEnv(w)) 146 { 147 *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl; 148 return kFALSE; 149 } 150 151 TChain chain("RunHeaders"); 152 set.AddFilesOn(chain); 153 154 MMcCorsikaRunHeader *h=0; 155 chain.SetBranchAddress("MMcCorsikaRunHeader.", &h); 156 chain.GetEntry(1); 157 158 if (!h) 159 { 160 *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl; 161 return kFALSE; 162 } 163 164 if (!w.Set(*h)) 165 { 166 *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl; 167 return kFALSE; 168 } 169 170 w.Print(); 171 return kTRUE; 172 } 173 174 Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const 175 { 176 // Some debug output 177 fLog->Separator("Compiling original MC distribution"); 178 179 weight.SetNameMcEvt("MMcEvtBasic"); 180 const TString w(weight.GetFormulaWeights()); 181 weight.SetNameMcEvt(); 182 183 *fLog << inf << "Using weights: " << w << endl; 184 *fLog << "Please stand by, this may take a while..." << flush; 185 186 if (fDisplay) 187 fDisplay->SetStatusLine1("Compiling MC distribution..."); 188 189 // Create chain 190 TChain chain("OriginalMC"); 191 set.AddFilesOn(chain); 192 193 // Prepare histogram 194 h.Reset(); 195 196 // Fill histogram from chain 197 h.SetDirectory(gROOT); 198 if (h.InheritsFrom(TH2::Class())) 199 { 200 h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)"); 201 h.SetXTitle("\\Theta [\\circ]"); 202 h.SetYTitle("E [GeV]"); 203 h.SetZTitle("Counts"); 204 chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff"); 205 } 206 else 207 { 208 h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)"); 209 h.SetXTitle("\\Theta [\\circ]"); 210 h.SetYTitle("Counts"); 211 chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff"); 212 } 213 h.SetDirectory(0); 214 215 *fLog << "done." << endl; 216 if (fDisplay) 217 fDisplay->SetStatusLine2("done."); 218 219 if (h.GetEntries()>0) 220 return kTRUE; 221 222 *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl; 223 224 return h.GetEntries()>0; 225 } 226 */ 227 228 Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const 229 { 230 TChain chain("OriginalMC"); 231 set.AddFilesOn(chain); 232 233 min = chain.GetMinimum("MMcEvtBasic.fEnergy"); 234 max = chain.GetMaximum("MMcEvtBasic.fEnergy"); 235 236 num = chain.GetEntries(); 237 238 if (num<100) 239 *fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl; 240 241 return num>=100; 242 } 243 244 Double_t MJTrainSeparation::GetDataRate(MDataSet &set) const 245 { 246 TChain chain1("Events"); 247 set.AddFilesOff(chain1); 248 249 const Double_t num = chain1.GetEntries(); 250 if (num<100) 251 { 252 *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl; 253 return -1; 254 } 255 256 TChain chain("EffectiveOnTime"); 257 set.AddFilesOff(chain); 258 259 TH1F h; 260 h.SetDirectory(gROOT); 261 h.SetNameTitle("OnTime", "Effective on-time"); 262 chain.Draw("MEffectiveOnTime.fVal>>OnTime", "", "goff"); 263 h.SetDirectory(0); 264 265 if (h.Integral()<1) 266 { 267 *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl; 268 return -1; 269 } 270 271 *fLog << inf << "Found " << num << " events in " << h.Integral(); 272 *fLog << "s (" << num/h.Integral() << "Hz)" << endl; 273 274 return num/h.Integral(); 275 } 276 277 /* 278 Scale: 279 280 281 TF1 fold("old", "x^(-2.6)", emin, emax); 282 TF1 fnew("new", "x^(-4.0)", emin, emax); 283 284 TF1 q("q", "new/old", emin, emax); 285 286 Double_t scale = 1./q.GetMaximum(emin, emax); 287 288 // Anzahl produzierter Events vor MFEnergySlope: 289 Double_t nold = fold.Integral(emin, emax); 290 291 // Anzahl produzierter Events nach MFEnergySlope: 292 Double_t nnew = fnew.Integral(emin, emax)*scale; 293 294 class MFSpectrum : MMcSpectrumWeight 295 { 296 Double_t fScale; 297 Bool_t fResult; 298 299 MFSpectrum::MFSpectrum(const char *name, const char *title) 300 { 301 fName = name ? name : "MMcSpectrumWeight"; 302 fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; 303 304 Init(fName, fTitle); 305 306 } 307 308 Int_t PreProcess(MParList *pList) 309 { 310 Int_t rc = MFSpectrumWeight::PreProcess(pList); 311 if (rc!=kTRUE) 312 return rc; 313 314 fScale = fEval->GetMaximum(fEnergyMin, fEnergyMax); 315 316 return kTRUE; 317 } 318 319 Int_t Process() 320 { 321 const Double_t e = fMcEvt->GetEnergy(); 322 323 Double_t prob = fFunc->Eval(e)/fScale; 324 325 const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope); 326 const Float_t Nrnd = ; 327 328 fResult = Nexp >= gRandom->Uniform(); 329 } 330 331 } 332 333 334 */ 335 336 Bool_t MJTrainSeparation::AutoTrain() 337 { 338 Double_t num, min, max; 339 if (!GetEventsProduced(fDataSetTrain, num, min, max)) 340 return kFALSE; 341 342 *fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl; 343 344 // Target spectrum 345 TF1 flx("Flux", "[0]*(x/1000)^(-2.6)", min, max); 346 flx.SetParameter(0, 1e-5); 347 348 // Number n0 of events this spectrum would produce per s and m^2 349 const Double_t n0 = flx.Integral(min, max); //[#] 350 351 // Area produced in MC 352 const Double_t A = TMath::Pi()*300*300; //[m²] 353 354 // Rate R of events this spectrum would produce per s 355 const Double_t R = n0*A; //[Hz] 356 357 // Number N of events produced (in trainings sample) 358 const Double_t N = num; //[#] 359 360 // This correponds to an observation time T [s] 361 const Double_t T = N/R; //[s] 362 363 // With an average data rate after star of 364 const Double_t r = GetDataRate(fDataSetTrain); //[Hz] 365 366 // this yields a number of n events to be read for training 367 const Double_t n = r*T; //[#] 368 369 if (r<0) 370 return kFALSE; 371 372 *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl; 373 *fLog << "For a data rate of " << r << "Hz this corresponds to " << n << " data events." << endl; 374 375 fNumTrainOn = (UInt_t)-1; 376 fNumTrainOff = TMath::Nint(n); 377 378 /* 379 An event rate dependent selection? 380 ---------------------------------- 381 Total average data rate: R 382 Goal number of events: N 383 Number of data events: N0 384 Rate assigned to single evt: r 385 386 Selection probability: N/N0 * r/R 387 388 f := N/N0 * r 389 390 MF f("f * MEventRate.fRate < rand"); 391 */ 392 393 394 return kTRUE; 395 } 396 138 397 Bool_t MJTrainSeparation::Train(const char *out) 139 398 { … … 148 407 return kFALSE; 149 408 } 409 410 // ----------------------- Auto Train? ---------------------- 411 412 if (fAutoTrain) 413 if (!AutoTrain()) 414 return kFALSE; 150 415 151 416 // --------------------- Setup files -------------------- -
trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h
r7420 r7539 24 24 UInt_t fNumTestOff; 25 25 26 Bool_t fAutoTrain; 27 26 28 void DisplayResult(MH3 &h31, MH3 &h32); 29 30 Bool_t GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const; 31 Double_t GetDataRate(MDataSet &set) const; 32 Bool_t AutoTrain(); 27 33 28 34 public: … … 44 50 } 45 51 52 void EnableAutoTrain(Bool_t b=kTRUE) { fAutoTrain = kTRUE; } 53 46 54 Bool_t Train(const char *out); 47 55
Note:
See TracChangeset
for help on using the changeset viewer.