Changeset 1809 for trunk/MagicSoft/Mars
- Timestamp:
- 03/08/03 14:00:30 (22 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 18 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r1800 r1809 1 1 -*-*- END OF LINE -*-*- 2 3 4 2003/03/08: Wolfgang Wittek 5 6 * macros/AnalyseCT1.C 7 - for the CT1 data analysis 8 9 * mhist/MHMatrix.[h,cc] 10 - let refcolumn start at 1 (not at 0) 11 12 * mhist/MHSigmaTheta.[h,cc] 13 - Draw replaced by DrawCopy 14 - add SetDirectory(NULL) 15 16 * manalysis/ MSelBasic.[h,cc] 17 MSelStandard.[h,cc] 18 MSelFinal.[h,cc] 19 - more detailed output for errors 20 - bugs removed 21 22 * manalysis/ MPadSchweizer.[h,cc] 23 - add SetDirectory(NULL) 24 - add fErrors 25 26 * mfilter/ MFEventSelector.[h,cc] 27 - add fErrors 28 29 * manalysis/MMultiDimDistCalc.[h,cc] 30 - check division by zero 31 32 * mhist/MHHadronness.[h,cc] 33 - check division by zero 34 - normalize distributions of hadronness 35 36 * mfileio/MCT1ReadPreProc.[h,cc] 37 - add event number (event.isecs_since_midday) 38 - change definition of "fIsMcFile", 39 because outpars.bmontecarlo is set wrongly sometimes 40 - copy pedestalRMS for each event from the header information 41 - check for the presence of a footer record even after reading 42 a run header 43 44 * mmc/ MMcEvt.[hxx,cxx]] 45 - add GetEvtNumber() 46 47 2 48 2003/02/27: Abelardo Moralejo 3 49 … … 63 109 * mmc/MMcTrigHeader.hxx 64 110 - Added GetMultiplicity() and GetMeanThreshold() 65 66 111 67 112 -
trunk/MagicSoft/Mars/macros/AnalyseCT1.C
r1767 r1809 1 1 void AnalyseCT1() 2 2 { 3 const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; 4 const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; 3 //----------------------------------------------- 4 //const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; 5 //const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; 5 6 //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6"; 6 7 //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc"; 7 const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc"; 8 //const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc"; 9 10 //----------------------------------------------- 11 const char *offfile = "~/Mars200203/CT1data/offdata.preproc"; 12 13 const char *onfile = "~/Mars200203/CT1data/mkn421_on.preproc"; 14 //const char *onfile = "~/Mars200203/CT1data/Crab_preproc_daniel_0.6"; 15 16 //const char *mcfile = "~/Mars200203/CT1data/mc_gammas.preproc"; //Version 0.5 17 //const char *mcfile = "~/Mars200203/CT1data/mc_Gammas.preproc.0.6"; 18 const char *mcfile = "~/Mars200203/CT1data/mkn421_mc_pedrms_0.001.preproc"; 19 20 //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6"; 21 //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6_040303"; 22 23 //----------------------------------------------- 24 //const char *offfile = "~magican/home/ct1test/PreprocData/offdata.preproc"; 25 //const char *onfile = "~magican/home/ct1test/PreprocData/mkn421_on.preproc"; 26 //const char *mcfile = "~magican/home/ct1test/PreprocData/mc_gammas.preproc"; 27 28 8 29 const char *filename; 9 30 10 Bool_t Data = kFALSE; // loop over data ? 11 Bool_t WLook = kFALSE; // write out look-alike events for hadrons ? 12 Bool_t WHist = kFALSE; // write out histograms for padding ? 13 14 Bool_t MC = kTRUE; // loop over MC gamma data ? 15 Bool_t RHist = kTRUE; // read in histograms for padding ? 16 Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ? 17 18 19 cout << "" << endl; 20 cout << "Macro AnalyseCT1 : Data, WHist = " << Data << ", " 21 << WHist << endl; 22 23 cout << " RHist, MC = " << RHist << ", " 24 << MC << endl; 25 cout << "" << endl; 31 //************************************************************************ 32 // Job 1 : read OFF data 33 // (generate sigmabar vs. Theta plot; write root file for OFF data (OFF1); 34 // generate hadron matrix for g/h separation) 35 36 Bool_t Job1 = kFALSE; 37 Bool_t WHistOFF = kTRUE; // write out histogram sigmabar vs. Theta ? 38 Bool_t WLookOFF = kTRUE; // write out look-alike events for hadrons ? 39 Bool_t WOFF1 = kTRUE; // write out root file OFF1 ? 40 41 42 // Job 2 : read ON data 43 // (generate sigmabar vs. Theta plot; write root file for ON data (ON1)) 44 45 Bool_t Job2 = kFALSE; 46 Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ? 47 Bool_t WLookON = kFALSE; // write out look-alike events for hadrons ? 48 Bool_t WON1 = kFALSE; // write out root file ON1 ? 49 50 // Job 3 : read MC gamma data, read sigmabar vs. Theta plot from OFF data 51 // (do padding; write root file for MC gammas (MC1); 52 // generate gamma matrix for g/h separation) 53 54 Bool_t Job3 = kFALSE; 55 Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ? 56 Bool_t WMC1 = kTRUE; // write out root file MC1 ? 57 58 59 60 // Job 4 : read OFF1 and MC1 files, read hadron and gamma matrix 61 // (produce the Neyman-Pearson plot) 62 Bool_t Job4 = kTRUE; 63 64 65 // Job 5 : read MC1 file, read hadron and gamma matrix 66 // (do g/h separation; write root file for MC gammas after final cuts (MC2)) 67 Bool_t Job5 = kFALSE; 68 Bool_t WMC2 = kFALSE; // write out root file MC2 ? 69 70 71 // Job 6 : read ON data, read hadron and gamma matrix 72 // (do g/h separation; write root file for ON data after final cuts (ON2)) 73 Bool_t Job6 = kFALSE; 74 Bool_t WON2 = kTRUE; // write out root file ON2 ? 75 //************************************************************************ 76 77 78 26 79 27 80 //--------------------------------------------------------------------- 28 // start of section for ON data 29 // (in the CT1 analysis the ON data are also used as a sample representing 30 // the hadrons for the g/h separation) 31 32 // read ON data file 33 34 // - to produce the 2D-histogram "sigmabar versus Theta" 81 // Job 1 82 //====== 83 84 // read OFF data file 85 86 // - produce the 2D-histogram "sigmabar versus Theta" for OFF data 35 87 // (to be used for the padding of the MC gamma data) 36 88 37 // - to write a file of look-alike events (for hadrons)89 // - write a file of look-alike events (hadron matrix) 38 90 // (to be used later together with the corresponding MC gamma file 39 91 // for the g/h separation in the analysis of the ON data) 40 92 41 // - to write a file of ON events (after the standard cuts,42 // 93 // - write a file of OFF events (OFF1) 94 // (after the standard cuts, before the g/h separation) 43 95 // (to be used together with the corresponding MC gamma file 44 96 // for the optimization of the g/h separation) 45 97 46 if ( Data)98 if (Job1) 47 99 { 48 100 cout << "=====================================================" << endl; 49 cout << "Macro AnalyseCT1 : Start of section for ON data" << endl; 101 cout << "Macro AnalyseCT1 : Start of Job 1" << endl; 102 103 cout << "" << endl; 104 cout << "Macro AnalyseCT1 : Job1, WHistOFF, WLookOFF, WOFF1 = " 105 << Job1 << ", " << WHistOFF << ", " << WLookOFF << ", " 106 << WOFF1 << endl; 107 108 109 TString typeData = "OFF"; 110 cout << "typeData = " << typeData << endl; 50 111 51 112 MTaskList tliston; … … 112 173 // 113 174 114 filename = o nfile;115 readname = "ReadCT1data";116 MCT1ReadPreProc read(filename, readname);175 filename = offfile; 176 RName = "ReadCT1data"; 177 MCT1ReadPreProc read(filename, RName); 117 178 118 179 MSelBasic selbasic; … … 150 211 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); 151 212 152 MSelStandard selstand(fHil las, fHillasSrc);213 MSelStandard selstand(fHilName, fHilSrcName); 153 214 154 215 MFillH hfill1("MHHillas", fHilName); … … 168 229 169 230 const char* mtxName = "MatrixHadrons"; 170 Int_t howMany = 30000; 171 char* outName = "matrix_hadrons.root"; 231 Int_t howMany = 50000; 232 TString outName = "matrix_hadrons_"; 233 outName += typeData; 234 outName += ".root"; 172 235 173 236 cout << "" << endl; … … 182 245 183 246 // matrix limitation for look-alike events (approximate number) 184 MFEventSelector selector("", "", readname);247 MFEventSelector selector("", "", RName); 185 248 selector.SetNumSelectEvts(howMany); 186 249 MFilterList flist; … … 192 255 MHMatrix *pmatrix = new MHMatrix(mtxName); 193 256 MHMatrix& matrix = *pmatrix; 257 258 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); 259 matrix.AddColumn("MSigmabar.fSigmabar"); 260 matrix.AddColumn("log10(Hillas.fSize)"); 261 matrix.AddColumn("HillasSrc.fDist"); 194 262 matrix.AddColumn("Hillas.fWidth"); 195 263 matrix.AddColumn("Hillas.fLength"); 196 matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize"); 197 matrix.AddColumn("abs(Hillas.fAsym)"); 264 matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); 198 265 matrix.AddColumn("abs(Hillas.fM3Long)"); 199 matrix.AddColumn("abs(Hillas.fM3Trans)");200 matrix.AddColumn("abs(HillasSrc.fHeadTail)");201 266 matrix.AddColumn("Hillas.fConc"); 202 267 matrix.AddColumn("Hillas.fConc1"); 203 matrix.AddColumn("HillasSrc.fDist"); 204 matrix.AddColumn("log10(Hillas.fSize)"); 205 matrix.AddColumn("HillasSrc.fAlpha"); 206 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); 207 //matrix.AddColumn("MMcEvt.fTheta"); 268 208 269 pliston.AddToList(pmatrix); 209 270 … … 211 272 fillmatg.SetFilter(&flist); 212 273 274 275 if (WOFF1) 276 { 277 TString outNameImage = typeData; 278 outNameImage += "1.root"; 279 MWriteRootFile write(outNameImage); 280 write.AddContainer("MSigmabar", "Events"); 281 write.AddContainer("Hillas", "Events"); 282 write.AddContainer("MMcEvt", "Events"); 283 write.AddContainer("HillasSrc", "Events"); 284 write.AddContainer("MRawRunHeader", "RunHeaders"); 285 //write.AddContainer("MMcRunHeader","RunHeaders"); 286 write.AddContainer("MSrcPosCam", "RunHeaders"); 287 } 288 213 289 //+++++++++++++++++++++++++++++++++++++++++++++++++++ 214 290 215 MSelFinal selfinal(fHillas, fHillasSrc);216 291 217 292 //***************************** … … 229 304 230 305 tliston.AddToList(&selstand); 306 if (WOFF1) 307 tliston.AddToList(&write); 308 309 tliston.AddToList(&flist); 310 tliston.AddToList(&fillmatg); 311 231 312 tliston.AddToList(&hfill1); 232 313 tliston.AddToList(&hfill2); … … 234 315 tliston.AddToList(&hfill4); 235 316 236 tliston.AddToList(&flist);237 tliston.AddToList(&fillmatg);238 239 //tliston.AddToList(&selfinal);240 317 //***************************** 241 318 … … 243 320 // Execute event loop 244 321 // 322 MProgressBar bar; 245 323 MEvtLoop evtloop; 246 324 evtloop.SetParList(&pliston); 325 evtloop.SetProgressBar(&bar); 247 326 248 327 Int_t maxevents = 1000000000; 328 //Int_t maxevents = 1000; 249 329 if ( !evtloop.Eventloop(maxevents) ) 250 330 return; … … 255 335 // Display the histograms 256 336 // 257 pliston.FindObject("MHHillas")->DrawClone(); 258 pliston.FindObject("MHHillasSrc")->DrawClone(); 337 TObject *c; 338 c = pliston.FindObject("MHHillas")->DrawClone(); 339 340 c = pliston.FindObject("MHHillasSrc")->DrawClone(); 259 341 260 342 //pliston.FindObject("MHHillasExt")->DrawClone(); 261 pliston.FindObject(nxt)->DrawClone();262 263 pliston.FindObject("MHStarMap")->DrawClone();343 c = pliston.FindObject(nxt)->DrawClone(); 344 345 c = pliston.FindObject("MHStarMap")->DrawClone(); 264 346 265 347 … … 277 359 Int_t nmaxevts = 2000; 278 360 279 // target distribution for the variable in column refcolumn ;361 // target distribution for the variable in column refcolumn (start at 1); 280 362 // set refcolumn negative if distribution is not to be changed 281 Int_t refcolumn = 1 2;363 Int_t refcolumn = 1; 282 364 Int_t nbins = 10; 283 365 Float_t frombin = 0.5; … … 311 393 // write out look-alike events (for hadrons) 312 394 // 313 if (WLook )395 if (WLookOFF) 314 396 { 315 397 TFile *writejens = new TFile(outName, "RECREATE", ""); … … 324 406 //------------------------------------------- 325 407 // Write histograms onto a file 326 if (WHist )408 if (WHistOFF) 327 409 { 328 410 MHSigmaTheta *sigtheta = … … 330 412 if (!sigtheta) 331 413 { 332 *fLog<< "Object 'SigmaTheta' not found" << endl;414 cout << "Object 'SigmaTheta' not found" << endl; 333 415 return; 334 416 } … … 337 419 TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); 338 420 339 TFile outfile("SigmaTheta.root", "RECREATE"); 421 TString outNameSigTh = "SigmaTheta_"; 422 outNameSigTh += typeData; 423 outNameSigTh += ".root"; 424 425 TFile outfile(outNameSigTh, "RECREATE"); 340 426 fHSigTh->Write(); 341 427 fHSigPixTh->Write(); … … 343 429 outfile.Close(); 344 430 345 cout << "File 'SigmaTheta.root'was written out" << endl;431 cout << "File " << outNameSigTh << " was written out" << endl; 346 432 } 347 433 348 434 349 cout << "Macro AnalyseCT1 : End of section for ON data" << endl;435 cout << "Macro AnalyseCT1 : End of Job 1" << endl; 350 436 cout << "===================================================" << endl; 351 437 } … … 353 439 354 440 //--------------------------------------------------------------------- 355 // start of section for MC gamma data 441 // Job 2 442 //====== 443 // read ON data file 444 445 // - produce the 2D-histogram "sigmabar versus Theta" for ON data 446 // (to be used for the padding of the MC gamma data) 447 448 // - write a file of look-alike events (hadron matrix) 449 // (to be used later together with the corresponding MC gamma file 450 // for the g/h separation in the analysis of the ON data) 451 452 // - write a file of ON events (ON1) 453 // (after the standard cuts, before the g/h separation) 454 // (to be used together with the corresponding MC gamma file 455 // for the optimization of the g/h separation) 456 457 if (Job2) 458 { 459 cout << "=====================================================" << endl; 460 cout << "Macro AnalyseCT1 : Start of Job 2" << endl; 461 462 cout << "" << endl; 463 cout << "Macro AnalyseCT1 : Job2, WHistON, WLookON, WON1 = " 464 << Job2 << ", " << WHistON << ", " << WLookON << ", " 465 << WON1 << endl; 466 467 TString typeData = "ON"; 468 cout << "typeData = " << typeData << endl; 469 470 MTaskList tliston; 471 MParList pliston; 472 pliston.AddToList(&tliston); 473 474 475 //:::::::::::::::::::::::::::::::::::::::::: 476 477 MBinning binssize("BinningSize"); 478 binssize.SetEdgesLog(50, 10, 1.0e5); 479 pliston.AddToList(&binssize); 480 481 MBinning binsdistc("BinningDist"); 482 binsdistc.SetEdges(50, 0, 1.4); 483 pliston.AddToList(&binsdistc); 484 485 MBinning binswidth("BinningWidth"); 486 binswidth.SetEdges(50, 0, 1.0); 487 pliston.AddToList(&binswidth); 488 489 MBinning binslength("BinningLength"); 490 binslength.SetEdges(50, 0, 1.0); 491 pliston.AddToList(&binslength); 492 493 MBinning binsalpha("BinningAlpha"); 494 binsalpha.SetEdges(100, -100, 100); 495 pliston.AddToList(&binsalpha); 496 497 MBinning binsasym("BinningAsym"); 498 binsasym.SetEdges(50, -1.5, 1.5); 499 pliston.AddToList(&binsasym); 500 501 MBinning binsht("BinningHeadTail"); 502 binsht.SetEdges(50, -1.5, 1.5); 503 pliston.AddToList(&binsht); 504 505 MBinning binsm3l("BinningM3Long"); 506 binsm3l.SetEdges(50, -1.5, 1.5); 507 pliston.AddToList(&binsm3l); 508 509 MBinning binsm3t("BinningM3Trans"); 510 binsm3t.SetEdges(50, -1.5, 1.5); 511 pliston.AddToList(&binsm3t); 512 513 514 //..... 515 MBinning binsb("BinningSigmabar"); 516 binsb.SetEdges( 100, 0.0, 5.0); 517 pliston.AddToList(&binsb); 518 519 MBinning binth("BinningTheta"); 520 binth.SetEdges( 70, -0.5, 69.5); 521 pliston.AddToList(&binth); 522 523 MBinning binsdiff("BinningDiffsigma2"); 524 binsdiff.SetEdges(100, -5.0, 20.0); 525 pliston.AddToList(&binsdiff); 526 //:::::::::::::::::::::::::::::::::::::::::: 527 528 529 //------------------------------------------- 530 // create the tasks which should be executed 531 // 532 533 filename = onfile; 534 RName = "ReadCT1data"; 535 MCT1ReadPreProc read(filename, RName); 536 537 MSelBasic selbasic; 538 539 MBlindPixelCalc blind; 540 blind.SetUseInterpolation(); 541 blind.SetUseBlindPixels(); 542 // blind.SetUseCetralPixel(); 543 544 // create an object of class MSigmabar, 545 // because class MHSigmaTheta will look for it 546 MSigmabar sigbar; 547 pliston.AddToList(&sigbar); 548 MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); 549 fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); 550 551 MImgCleanStd clean; 552 553 // calculate also extended image parameters 554 TString fHilName = "Hillas"; 555 MHillasExt hext(fHilName); 556 pliston.AddToList(&hext); 557 MHillasExt *fHillas = &hext; 558 559 // name of output container for MHillasCalc 560 // = name of MHillas object 561 MHillasCalc hcalc(fHilName); 562 563 // name of output container for MHillasSrcCalc 564 // = name of MHillasSrc object 565 TString fHilSrcName = "HillasSrc"; 566 MHillasSrc hilsrc(fHilSrcName); 567 MHillasSrc *fHillasSrc = &hilsrc; 568 pliston.AddToList(&hilsrc); 569 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); 570 571 MSelStandard selstand(fHilName, fHilSrcName); 572 573 MFillH hfill1("MHHillas", fHilName); 574 MFillH hfill2("MHStarMap", fHilName); 575 576 TString nxt = "HillasExt"; 577 MHHillasExt hhilext(nxt, "", fHilName); 578 pliston.AddToList(&hhilext); 579 TString namext = nxt; 580 namext += "[MHHillasExt]"; 581 MFillH hfill3(namext, fHilSrcName); 582 583 MFillH hfill4("MHHillasSrc", fHilSrcName); 584 585 //+++++++++++++++++++++++++++++++++++++++++++++++++++ 586 // prepare writing of look-alike events (for hadrons) 587 588 const char* mtxName = "MatrixHadrons"; 589 Int_t howMany = 50000; 590 TString outName = "matrix_hadrons_"; 591 outName += typeData; 592 outName += ".root"; 593 594 595 cout << "" << endl; 596 cout << "========================================================" << endl; 597 cout << "Matrix for (hadrons)" << endl; 598 cout << "input file = " << filename<< endl; 599 cout << "matrix name = " << mtxName << endl; 600 cout << "max. no.of look-alike events = " << howMany << endl; 601 cout << "name of output root file = " << outName << endl; 602 cout << "" << endl; 603 604 605 // matrix limitation for look-alike events (approximate number) 606 MFEventSelector selector("", "", RName); 607 selector.SetNumSelectEvts(howMany); 608 MFilterList flist; 609 flist.AddToList(&selector); 610 611 // 612 // --- setup the matrix and add it to the parlist 613 // 614 MHMatrix *pmatrix = new MHMatrix(mtxName); 615 MHMatrix& matrix = *pmatrix; 616 617 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); 618 matrix.AddColumn("MSigmabar.fSigmabar"); 619 matrix.AddColumn("log10(Hillas.fSize)"); 620 matrix.AddColumn("HillasSrc.fDist"); 621 matrix.AddColumn("Hillas.fWidth"); 622 matrix.AddColumn("Hillas.fLength"); 623 matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); 624 matrix.AddColumn("abs(Hillas.fM3Long)"); 625 matrix.AddColumn("Hillas.fConc"); 626 matrix.AddColumn("Hillas.fConc1"); 627 628 pliston.AddToList(pmatrix); 629 630 MFillH fillmatg(mtxName); 631 fillmatg.SetFilter(&flist); 632 633 634 if (WON1) 635 { 636 TString outNameImage = typeData; 637 outNameImage += "1.root"; 638 MWriteRootFile write(outNameImage); 639 write.AddContainer("MSigmabar", "Events"); 640 write.AddContainer("Hillas", "Events"); 641 write.AddContainer("MMcEvt", "Events"); 642 write.AddContainer("HillasSrc", "Events"); 643 write.AddContainer("MRawRunHeader", "RunHeaders"); 644 //write.AddContainer("MMcRunHeader","RunHeaders"); 645 write.AddContainer("MSrcPosCam", "RunHeaders"); 646 } 647 648 //+++++++++++++++++++++++++++++++++++++++++++++++++++ 649 650 651 //***************************** 652 // tasks to be executed 653 654 tliston.AddToList(&read); 655 656 tliston.AddToList(&selbasic); 657 tliston.AddToList(&blind); 658 tliston.AddToList(&fillsigtheta); 659 tliston.AddToList(&clean); 660 661 tliston.AddToList(&hcalc); 662 tliston.AddToList(&csrc1); 663 664 tliston.AddToList(&selstand); 665 if (WON1) 666 tliston.AddToList(&write); 667 668 tliston.AddToList(&flist); 669 tliston.AddToList(&fillmatg); 670 671 tliston.AddToList(&hfill1); 672 tliston.AddToList(&hfill2); 673 tliston.AddToList(&hfill3); 674 tliston.AddToList(&hfill4); 675 676 677 //***************************** 678 679 //------------------------------------------- 680 // Execute event loop 681 // 682 MProgressBar bar; 683 MEvtLoop evtloop; 684 evtloop.SetParList(&pliston); 685 evtloop.SetProgressBar(&bar); 686 687 Int_t maxevents = 1000000000; 688 //Int_t maxevents = 1000; 689 if ( !evtloop.Eventloop(maxevents) ) 690 return; 691 692 tliston.PrintStatistics(0, kTRUE); 693 694 //------------------------------------------- 695 // Display the histograms 696 // 697 TObject *c; 698 c = pliston.FindObject("MHHillas")->DrawClone(); 699 700 c = pliston.FindObject("MHHillasSrc")->DrawClone(); 701 702 //pliston.FindObject("MHHillasExt")->DrawClone(); 703 c = pliston.FindObject(nxt)->DrawClone(); 704 705 c = pliston.FindObject("MHStarMap")->DrawClone(); 706 707 708 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 709 // 710 // ---------------------------------------------------------- 711 // Definition of the reference sample of look-alike events 712 // (this is a subsample of the original sample) 713 // ---------------------------------------------------------- 714 // 715 cout << "" << endl; 716 cout << "========================================================" << endl; 717 // Select a maximum of nmaxevts events from the sample of look-alike 718 // events. They will form the reference sample. 719 Int_t nmaxevts = 2000; 720 721 // target distribution for the variable in column refcolumn (start at 1); 722 // set refcolumn negative if distribution is not to be changed 723 Int_t refcolumn = 1; 724 Int_t nbins = 10; 725 Float_t frombin = 0.5; 726 Float_t tobin = 1.0; 727 TH1F *thsh = new TH1F("thsh","target distribution", 728 nbins, frombin, tobin); 729 thsh->SetXTitle("cos( \\Theta)"); 730 thsh->SetYTitle("Counts"); 731 Float_t dbin = (tobin-frombin)/nbins; 732 Float_t lbin = frombin +dbin*0.5; 733 for (Int_t j=1; j<=nbins; j++) 734 { 735 thsh->Fill(lbin,1.0); 736 lbin += dbin; 737 } 738 739 cout << "" << endl; 740 cout << "========================================================" << endl; 741 cout << "Macro AnalyseCT1 : define reference sample for hadrons" << endl; 742 cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = " 743 << refcolumn << ", " << nmaxevts << endl; 744 745 if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) 746 { 747 cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl; 748 return; 749 } 750 //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 751 752 //------------------------------------------- 753 // write out look-alike events (for hadrons) 754 // 755 if (WLookON) 756 { 757 TFile *writejens = new TFile(outName, "RECREATE", ""); 758 matrix.Write(); 759 cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file " 760 << outName << endl; 761 762 writejens->Close(); 763 delete writejens; 764 } 765 766 //------------------------------------------- 767 // Write histograms onto a file 768 if (WHistON) 769 { 770 MHSigmaTheta *sigtheta = 771 (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); 772 if (!sigtheta) 773 { 774 cout << "Object 'SigmaTheta' not found" << endl; 775 return; 776 } 777 TH2D *fHSigTh = sigtheta->GetSigmaTheta(); 778 TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta(); 779 TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); 780 781 TString outNameSigTh = "SigmaTheta_"; 782 outNameSigTh += typeData; 783 outNameSigTh += ".root"; 784 785 TFile outfile(outNameSigTh, "RECREATE"); 786 fHSigTh->Write(); 787 fHSigPixTh->Write(); 788 fHDifPixTh->Write(); 789 outfile.Close(); 790 791 cout << "File " << outNameSigTh << " was written out" << endl; 792 } 793 794 795 cout << "Macro AnalyseCT1 : End of Job 2" << endl; 796 cout << "===================================================" << endl; 797 } 798 799 800 //--------------------------------------------------------------------- 801 // Job 3 802 //====== 356 803 357 804 // read MC gamma data … … 360 807 // (using the 2D-histogram "sigmabar versus Theta" of the OFF data) 361 808 362 // - to write a file of look-alike events ( for gammas)809 // - to write a file of look-alike events (gammas matrix) 363 810 // (to be used later together with the corresponding hadron file 364 811 // for the g/h separation in the analysis of the ON data) 365 812 366 // - to write a file of padded MC gamma events 813 // - to write a file of padded MC gamma events (MC1) 367 814 // (after the standard cuts, before the g/h separation) 368 815 // (to be used together with the corresponding hadron file 369 816 // for the optimization of the g/h separation) 370 817 371 // - to write a file of padded MC gamma events (after the final cuts) 372 // (to be used for the optimization of the energy estimator 373 // and for calculating the effective collection areas) 374 375 if (MC) 818 819 if (Job3) 376 820 { 377 cout << "=============================================================" 378 << endl; 379 cout << "Macro : AnalyseCT1 : Start of section for MC gamma data" 380 << endl; 381 382 MTaskList tlist; 383 MParList plist; 384 385 plist.AddToList(&tlist); 821 cout << "=====================================================" << endl; 822 cout << "Macro AnalyseCT1 : Start of Job 3" << endl; 823 824 cout << "" << endl; 825 cout << "Macro AnalyseCT1 : Job3, WLookMC, WMC1 = " 826 << Job3 << ", " << WLookMC << ", " << WMC1 << endl; 827 828 829 TString typeMC = "MC"; 830 cout << "typeMC = " << typeMC << endl; 831 832 // use for padding sigmabar vs. Theta from ON or OFF data 833 TString typeData = "ON"; 834 //TString typeData = "OFF"; 835 cout << "typeData = " << typeData << endl; 386 836 387 837 //------------------------------------ … … 389 839 // and "3D-ThetaPixSigma" 390 840 // and "3D-ThetaPixDiff" 391 if (RHist) 392 { 393 cout << "Reading in file 'SigmaTheta.root' " << endl; 394 395 TFile *infile = new TFile("SigmaTheta.root"); 841 842 843 TString outNameSigTh = "SigmaTheta_"; 844 outNameSigTh += typeData; 845 outNameSigTh += ".root"; 846 847 848 cout << "Reading in file " << outNameSigTh << endl; 849 850 TFile *infile = new TFile(outNameSigTh); 396 851 infile->ls(); 397 852 … … 400 855 if (!fHSigmaTheta) 401 856 { 402 *fLog<< "Object '2D-ThetaSigmabar' not found on root file" << endl;857 cout << "Object '2D-ThetaSigmabar' not found on root file" << endl; 403 858 return; 404 859 } … … 409 864 if (!fHSigmaPixTheta) 410 865 { 411 *fLog<< "Object '3D-ThetaPixSigma' not found on root file" << endl;866 cout << "Object '3D-ThetaPixSigma' not found on root file" << endl; 412 867 return; 413 868 } … … 418 873 if (!fHSigmaPixTheta) 419 874 { 420 *fLog<< "Object '3D-ThetaPixDiff' not found on root file" << endl;875 cout << "Object '3D-ThetaPixDiff' not found on root file" << endl; 421 876 return; 422 877 } 423 878 cout << "Object '3D-ThetaPixDiff' was read in" << endl; 424 } 425 else 426 { 427 MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); 428 if (!sigtheta) 429 { 430 cout << "Object with name 'SigmaTheta' not found" << endl; 431 return; 432 } 433 434 TH2D *fHSigmaTheta = sigtheta->GetSigmaTheta(); 435 TH3D *fHSigmaPixTheta = sigtheta->GetSigmaPixTheta(); 436 TH3D *fHDiffPixTheta = sigtheta->GetDiffPixTheta(); 437 } 879 438 880 //------------------------------------ 439 881 882 MTaskList tlist; 883 MParList plist; 884 885 plist.AddToList(&tlist); 440 886 441 887 … … 539 985 MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); 540 986 541 MSelStandard selstand(fHil las, fHillasSrc);987 MSelStandard selstand(fHilName, fHilSrcName); 542 988 543 989 MFillH hfill1("MHHillas", fHilName); … … 557 1003 558 1004 const char* mtxName = "MatrixGammas"; 559 Int_t howMany = 30000; 560 char* outName = "matrix_gammas.root"; 1005 Int_t howMany = 50000; 1006 1007 TString outName = "matrix_gammas_"; 1008 outName += typeMC; 1009 outName += "_"; 1010 outName += typeData; 1011 outName += ".root"; 1012 561 1013 562 1014 cout << "" << endl; … … 581 1033 MHMatrix *pmatrix = new MHMatrix(mtxName); 582 1034 MHMatrix& matrix = *pmatrix; 1035 1036 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); 1037 matrix.AddColumn("MSigmabar.fSigmabar"); 1038 matrix.AddColumn("log10(Hillas.fSize)"); 1039 matrix.AddColumn("HillasSrc.fDist"); 583 1040 matrix.AddColumn("Hillas.fWidth"); 584 1041 matrix.AddColumn("Hillas.fLength"); 585 matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize"); 586 matrix.AddColumn("abs(Hillas.fAsym)"); 1042 matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); 587 1043 matrix.AddColumn("abs(Hillas.fM3Long)"); 588 matrix.AddColumn("abs(Hillas.fM3Trans)");589 matrix.AddColumn("abs(HillasSrc.fHeadTail)");590 1044 matrix.AddColumn("Hillas.fConc"); 591 1045 matrix.AddColumn("Hillas.fConc1"); 592 matrix.AddColumn("HillasSrc.fDist"); 593 matrix.AddColumn("log10(Hillas.fSize)"); 594 matrix.AddColumn("HillasSrc.fAlpha"); 595 matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); 596 //matrix.AddColumn("MMcEvt.fTheta"); 1046 597 1047 plist.AddToList(pmatrix); 598 1048 … … 600 1050 fillmatg.SetFilter(&flist); 601 1051 1052 1053 if (WMC1) 1054 { 1055 TString outNameImage = typeMC; 1056 outNameImage += "_"; 1057 outNameImage += typeData; 1058 outNameImage += "1.root"; 1059 MWriteRootFile write(outNameImage); 1060 write.AddContainer("MSigmabar", "Events"); 1061 write.AddContainer("Hillas", "Events"); 1062 write.AddContainer("MMcEvt", "Events"); 1063 write.AddContainer("HillasSrc", "Events"); 1064 write.AddContainer("MRawRunHeader", "RunHeaders"); 1065 //write.AddContainer("MMcRunHeader","RunHeaders"); 1066 write.AddContainer("MSrcPosCam", "RunHeaders"); 1067 } 1068 1069 602 1070 //+++++++++++++++++++++++++++++++++++++++++++++++++++ 603 1071 604 605 MSelFinal selfinal(fHillas, fHillasSrc);606 1072 607 1073 //***************************** … … 619 1085 620 1086 tlist.AddToList(&selstand); 1087 if (WMC1) 1088 tlist.AddToList(&write); 1089 1090 tlist.AddToList(&flist); 1091 tlist.AddToList(&fillmatg); 1092 621 1093 tlist.AddToList(&hfill1); 622 1094 tlist.AddToList(&hfill2); … … 624 1096 tlist.AddToList(&hfill4); 625 1097 626 tlist.AddToList(&flist); 627 tlist.AddToList(&fillmatg); 628 629 //tlist.AddToList(&selfinal); 1098 630 1099 //***************************** 631 1100 … … 634 1103 // Execute event loop 635 1104 // 1105 MProgressBar bar; 636 1106 MEvtLoop evtloop; 637 1107 evtloop.SetParList(&plist); 1108 evtloop.SetProgressBar(&bar); 638 1109 639 1110 Int_t maxevents = 1000000000; 1111 //Int_t maxevents = 100; 640 1112 if ( !evtloop.Eventloop(maxevents) ) 641 1113 return; … … 670 1142 // target distribution for the variable in column refcolumn; 671 1143 // set refcolumn negative if distribution is not to be changed 672 Int_t refcolumn = 1 2;1144 Int_t refcolumn = 1; 673 1145 Int_t nbins = 10; 674 1146 Float_t frombin = 0.5; … … 711 1183 } 712 1184 713 cout << "Macro AnalyseCT1 : End of section for MC gamma data"1185 cout << "Macro AnalyseCT1 : End of Job 3" 714 1186 << endl; 715 1187 cout << "=========================================================" 716 1188 << endl; 717 1189 } 1190 1191 1192 1193 //--------------------------------------------------------------------- 1194 // Job 4 1195 //====== 1196 1197 // read OFF1 and MC1 data files 1198 // 1199 // - produce Neyman-Pearson plot 1200 1201 if (Job4) 1202 { 1203 cout << "=====================================================" << endl; 1204 cout << "Macro AnalyseCT1 : Start of Job 4" << endl; 1205 1206 cout << "" << endl; 1207 cout << "Macro AnalyseCT1 : Job4 = " << Job4 << endl; 1208 1209 1210 TString typeMC = "MC"; 1211 cout << "typeMC = " << typeMC << endl; 1212 1213 // use hadron matrix from ON or OFF data 1214 TString typeData = "ON"; 1215 //TString typeData = "OFF"; 1216 cout << "typeData = " << typeData << endl; 1217 1218 //************************************************************************* 1219 // read in matrices of look-alike events 1220 1221 const char* mtxName = "MatrixGammas"; 1222 1223 1224 TString outName = "matrix_gammas_"; 1225 outName += typeMC; 1226 outName += "_"; 1227 outName += typeData; 1228 outName += ".root"; 1229 1230 cout << "" << endl; 1231 cout << "========================================================" << endl; 1232 cout << "Get matrix for (gammas)" << endl; 1233 cout << "matrix name = " << mtxName << endl; 1234 cout << "name of root file = " << outName << endl; 1235 cout << "" << endl; 1236 1237 1238 // read in the object with the name 'mtxName' from file 'outName' 1239 // 1240 TFile *fileg = new TFile(outName); 1241 1242 MHMatrix matrixg; 1243 matrixg.Read(mtxName); 1244 pmatrixg = &matrixg; 1245 1246 delete fileg; 1247 1248 1249 //***************************************************************** 1250 1251 const char* mtxName = "MatrixHadrons"; 1252 1253 TString outName = "matrix_hadrons_"; 1254 outName += typeData; 1255 outName += ".root"; 1256 1257 cout << "" << endl; 1258 cout << "========================================================" << endl; 1259 cout << " Get matrix for (hadrons)" << endl; 1260 cout << "matrix name = " << mtxName << endl; 1261 cout << "name of root file = " << outName << endl; 1262 cout << "" << endl; 1263 1264 1265 // read in the object with the name mtxName 1266 // 1267 TFile *fileh = new TFile(outName); 1268 1269 MHMatrix matrixh; 1270 matrixh.Read(mtxName); 1271 pmatrixh = &matrixh; 1272 1273 delete fileh; 1274 1275 //***************************************************************** 1276 1277 MTaskList tliston; 1278 MParList pliston; 1279 pliston.AddToList(&tliston); 1280 1281 pliston.AddToList(pmatrixg); 1282 pliston.AddToList(pmatrixh); 1283 1284 //:::::::::::::::::::::::::::::::::::::::::: 1285 1286 MBinning binssize("BinningSize"); 1287 binssize.SetEdgesLog(50, 10, 1.0e5); 1288 pliston.AddToList(&binssize); 1289 1290 MBinning binsdistc("BinningDist"); 1291 binsdistc.SetEdges(50, 0, 1.4); 1292 pliston.AddToList(&binsdistc); 1293 1294 MBinning binswidth("BinningWidth"); 1295 binswidth.SetEdges(50, 0, 1.0); 1296 pliston.AddToList(&binswidth); 1297 1298 MBinning binslength("BinningLength"); 1299 binslength.SetEdges(50, 0, 1.0); 1300 pliston.AddToList(&binslength); 1301 1302 MBinning binsalpha("BinningAlpha"); 1303 binsalpha.SetEdges(100, -100, 100); 1304 pliston.AddToList(&binsalpha); 1305 1306 MBinning binsasym("BinningAsym"); 1307 binsasym.SetEdges(50, -1.5, 1.5); 1308 pliston.AddToList(&binsasym); 1309 1310 MBinning binsht("BinningHeadTail"); 1311 binsht.SetEdges(50, -1.5, 1.5); 1312 pliston.AddToList(&binsht); 1313 1314 MBinning binsm3l("BinningM3Long"); 1315 binsm3l.SetEdges(50, -1.5, 1.5); 1316 pliston.AddToList(&binsm3l); 1317 1318 MBinning binsm3t("BinningM3Trans"); 1319 binsm3t.SetEdges(50, -1.5, 1.5); 1320 pliston.AddToList(&binsm3t); 1321 1322 1323 //..... 1324 MBinning binsb("BinningSigmabar"); 1325 binsb.SetEdges( 100, 0.0, 5.0); 1326 pliston.AddToList(&binsb); 1327 1328 MBinning binth("BinningTheta"); 1329 binth.SetEdges( 70, -0.5, 69.5); 1330 pliston.AddToList(&binth); 1331 1332 MBinning binsdiff("BinningDiffsigma2"); 1333 binsdiff.SetEdges(100, -5.0, 20.0); 1334 pliston.AddToList(&binsdiff); 1335 //:::::::::::::::::::::::::::::::::::::::::: 1336 1337 1338 //------------------------------------------- 1339 // create the tasks which should be executed 1340 // 1341 1342 TString filenameData = typeData; 1343 filenameData += "1.root"; 1344 1345 TString filenameMC = typeMC; 1346 filenameMC += "_"; 1347 filenameMC += typeData; 1348 filenameMC += "1.root"; 1349 1350 1351 cout << "filenameData = " << filenameData << endl; 1352 cout << "filenameMC = " << filenameMC << endl; 1353 1354 MReadMarsFile read("Events", filenameMC); 1355 read.AddFile(filenameData); 1356 read.DisableAutoScheme(); 1357 1358 //....................................................................... 1359 // g/h separation 1360 1361 MMultiDimDistCalc ghsep; 1362 ghsep.SetUseNumRows(25); 1363 ghsep.SetUseKernelMethod(kFALSE); 1364 1365 //....................................................................... 1366 1367 // geometry is needed in MHHillas... classes 1368 MGeomCam *fGeom = 1369 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1370 1371 TString fHilName = "Hillas"; 1372 TString fHilSrcName = "HillasSrc"; 1373 1374 Float_t hadcut = 0.2; 1375 MSelFinal selfinalgh(fHilName, fHilSrcName); 1376 selfinalgh.SetHadronnessCut(hadcut); 1377 selfinalgh.SetAlphaCut(100.0); 1378 1379 MFillH fillh("MHHadronness"); 1380 1381 MSelFinal selfinal(fHilName, fHilSrcName); 1382 selfinal.SetHadronnessCut(hadcut); 1383 selfinal.SetAlphaCut(20.0); 1384 1385 MFillH hfill1("MHHillas", fHilName); 1386 MFillH hfill2("MHStarMap", fHilName); 1387 1388 TString nxt = "HillasExt"; 1389 MHHillasExt hhilext(nxt, "", fHilName); 1390 pliston.AddToList(&hhilext); 1391 TString namext = nxt; 1392 namext += "[MHHillasExt]"; 1393 MFillH hfill3(namext, fHilSrcName); 1394 1395 MFillH hfill4("MHHillasSrc", fHilSrcName); 1396 1397 1398 //***************************** 1399 // tasks to be executed 1400 1401 tliston.AddToList(&read); 1402 1403 tliston.AddToList(&ghsep); 1404 tliston.AddToList(&fillh); 1405 1406 tliston.AddToList(&selfinalgh); 1407 tliston.AddToList(&hfill1); 1408 tliston.AddToList(&hfill2); 1409 tliston.AddToList(&hfill3); 1410 tliston.AddToList(&hfill4); 1411 1412 tliston.AddToList(&selfinal); 1413 1414 //***************************** 1415 1416 //------------------------------------------- 1417 // Execute event loop 1418 // 1419 MProgressBar bar; 1420 MEvtLoop evtloop; 1421 evtloop.SetParList(&pliston); 1422 evtloop.SetProgressBar(&bar); 1423 1424 Int_t maxevents = 1000000000; 1425 //Int_t maxevents = 1000; 1426 if ( !evtloop.Eventloop(maxevents) ) 1427 return; 1428 1429 tliston.PrintStatistics(0, kTRUE); 1430 1431 1432 //------------------------------------------- 1433 // Display the histograms 1434 // 1435 TObject *c; 1436 1437 c = pliston.FindObject("MHHadronness")->DrawClone(); 1438 c->Print(); 1439 1440 //c = pliston.FindObject("MHHillas")->DrawClone(); 1441 1442 c = pliston.FindObject("MHHillasSrc")->DrawClone(); 1443 1444 //pliston.FindObject("MHHillasExt")->DrawClone(); 1445 //c = pliston.FindObject(nxt)->DrawClone(); 1446 1447 c = pliston.FindObject("MHStarMap")->DrawClone(); 1448 1449 1450 1451 cout << "Macro AnalyseCT1 : End of Job 4" << endl; 1452 cout << "===================================================" << endl; 1453 } 1454 1455 1456 //--------------------------------------------------------------------- 1457 // Job 5 1458 //====== 1459 1460 // - read in the matrices of look-alike events for gammas and hadrons 1461 1462 // then read MC1 data file 1463 // 1464 // - perform the g/h separation 1465 // - apply the final cuts 1466 // 1467 // - write a file of MC gamma events (MC2) 1468 // (after the g/h separation and after the final cuts) 1469 1470 if (Job5) 1471 { 1472 cout << "=====================================================" << endl; 1473 cout << "Macro AnalyseCT1 : Start of Job 5" << endl; 1474 1475 cout << "" << endl; 1476 cout << "Macro AnalyseCT1 : Job5, WMC2 = " 1477 << Job5 << ", " << WMC2 << endl; 1478 1479 TString typeInput = "MC"; 1480 cout << "typeInput = " << typeInput << endl; 1481 1482 TString typeMC = "MC"; 1483 cout << "typeMC = " << typeMC << endl; 1484 1485 // use hadron matrix from ON or OFF data 1486 TString typeData = "ON"; 1487 //TString typeData = "OFF"; 1488 cout << "typeData = " << typeData << endl; 1489 1490 //************************************************************************* 1491 // read in matrices of look-alike events 1492 1493 const char* mtxName = "MatrixGammas"; 1494 1495 TString outName = "matrix_gammas_"; 1496 outName += typeMC; 1497 outName += "_"; 1498 outName += typeData; 1499 outName += ".root"; 1500 1501 1502 cout << "" << endl; 1503 cout << "========================================================" << endl; 1504 cout << "Get matrix for (gammas)" << endl; 1505 cout << "matrix name = " << mtxName << endl; 1506 cout << "name of root file = " << outName << endl; 1507 cout << "" << endl; 1508 1509 1510 // read in the object with the name 'mtxName' from file 'outName' 1511 // 1512 TFile *fileg = new TFile(outName); 1513 1514 MHMatrix matrixg; 1515 matrixg.Read(mtxName); 1516 pmatrixg = &matrixg; 1517 1518 delete fileg; 1519 1520 //***************************************************************** 1521 1522 const char* mtxName = "MatrixHadrons"; 1523 1524 TString outName = "matrix_hadrons_"; 1525 outName += typeData; 1526 outName += ".root"; 1527 1528 1529 cout << "" << endl; 1530 cout << "========================================================" << endl; 1531 cout << " Get matrix for (hadrons)" << endl; 1532 cout << "matrix name = " << mtxName << endl; 1533 cout << "name of root file = " << outName << endl; 1534 cout << "" << endl; 1535 1536 1537 // read in the object with the name mtxName 1538 // 1539 TFile *fileh = new TFile(outName); 1540 1541 MHMatrix matrixh; 1542 matrixh.Read(mtxName); 1543 pmatrixh = &matrixh; 1544 1545 delete fileh; 1546 1547 //************************************************************************* 1548 1549 1550 MTaskList tliston; 1551 MParList pliston; 1552 pliston.AddToList(&tliston); 1553 1554 pliston.AddToList(pmatrixg); 1555 pliston.AddToList(pmatrixh); 1556 1557 //:::::::::::::::::::::::::::::::::::::::::: 1558 1559 MBinning binssize("BinningSize"); 1560 binssize.SetEdgesLog(50, 10, 1.0e5); 1561 pliston.AddToList(&binssize); 1562 1563 MBinning binsdistc("BinningDist"); 1564 binsdistc.SetEdges(50, 0, 1.4); 1565 pliston.AddToList(&binsdistc); 1566 1567 MBinning binswidth("BinningWidth"); 1568 binswidth.SetEdges(50, 0, 1.0); 1569 pliston.AddToList(&binswidth); 1570 1571 MBinning binslength("BinningLength"); 1572 binslength.SetEdges(50, 0, 1.0); 1573 pliston.AddToList(&binslength); 1574 1575 MBinning binsalpha("BinningAlpha"); 1576 binsalpha.SetEdges(100, -100, 100); 1577 pliston.AddToList(&binsalpha); 1578 1579 MBinning binsasym("BinningAsym"); 1580 binsasym.SetEdges(50, -1.5, 1.5); 1581 pliston.AddToList(&binsasym); 1582 1583 MBinning binsht("BinningHeadTail"); 1584 binsht.SetEdges(50, -1.5, 1.5); 1585 pliston.AddToList(&binsht); 1586 1587 MBinning binsm3l("BinningM3Long"); 1588 binsm3l.SetEdges(50, -1.5, 1.5); 1589 pliston.AddToList(&binsm3l); 1590 1591 MBinning binsm3t("BinningM3Trans"); 1592 binsm3t.SetEdges(50, -1.5, 1.5); 1593 pliston.AddToList(&binsm3t); 1594 1595 1596 //..... 1597 MBinning binsb("BinningSigmabar"); 1598 binsb.SetEdges( 100, 0.0, 5.0); 1599 pliston.AddToList(&binsb); 1600 1601 MBinning binth("BinningTheta"); 1602 binth.SetEdges( 70, -0.5, 69.5); 1603 pliston.AddToList(&binth); 1604 1605 MBinning binsdiff("BinningDiffsigma2"); 1606 binsdiff.SetEdges(100, -5.0, 20.0); 1607 pliston.AddToList(&binsdiff); 1608 //:::::::::::::::::::::::::::::::::::::::::: 1609 1610 1611 //------------------------------------------- 1612 // create the tasks which should be executed 1613 // 1614 1615 TString filenameMC = typeInput; 1616 filenameMC += "_"; 1617 filenameMC += typeData; 1618 filenameMC += "1.root"; 1619 1620 readname = "ReadCT1MCdata"; 1621 MReadMarsFile read("Events", filenameMC); 1622 read.DisableAutoScheme(); 1623 1624 1625 //....................................................................... 1626 // g/h separation 1627 1628 MMultiDimDistCalc ghsep; 1629 ghsep.SetUseNumRows(25); 1630 ghsep.SetUseKernelMethod(kFALSE); 1631 1632 1633 1634 //....................................................................... 1635 1636 // geometry is needed in MHHillas... classes 1637 MGeomCam *fGeom = 1638 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1639 1640 TString fHilName = "Hillas"; 1641 TString fHilSrcName = "HillasSrc"; 1642 1643 Float_t hadcut = 0.2; 1644 MSelFinal selfinalgh(fHilName, fHilSrcName); 1645 selfinalgh.SetHadronnessCut(hadcut); 1646 selfinalgh.SetAlphaCut(100.0); 1647 1648 MFillH fillh("MHHadronness"); 1649 1650 MSelFinal selfinal(fHilName, fHilSrcName); 1651 selfinal.SetHadronnessCut(hadcut); 1652 selfinal.SetAlphaCut(20.0); 1653 1654 if (WMC2) 1655 { 1656 TString outNameImage = typeInput; 1657 outNameImage += "_"; 1658 outNameImage += typeData; 1659 outNameImage += "2.root"; 1660 MWriteRootFile write(outNameImage); 1661 write.AddContainer("MHadronness", "Events"); 1662 write.AddContainer("MSigmabar", "Events"); 1663 write.AddContainer("Hillas", "Events"); 1664 write.AddContainer("MMcEvt", "Events"); 1665 write.AddContainer("HillasSrc", "Events"); 1666 write.AddContainer("MRawRunHeader", "RunHeaders"); 1667 //write.AddContainer("MMcRunHeader","RunHeaders"); 1668 write.AddContainer("MSrcPosCam", "RunHeaders"); 1669 } 1670 1671 1672 MFillH hfill1("MHHillas", fHilName); 1673 MFillH hfill2("MHStarMap", fHilName); 1674 1675 TString nxt = "HillasExt"; 1676 MHHillasExt hhilext(nxt, "", fHilName); 1677 pliston.AddToList(&hhilext); 1678 TString namext = nxt; 1679 namext += "[MHHillasExt]"; 1680 MFillH hfill3(namext, fHilSrcName); 1681 1682 MFillH hfill4("MHHillasSrc", fHilSrcName); 1683 1684 1685 1686 //***************************** 1687 // tasks to be executed 1688 1689 tliston.AddToList(&read); 1690 1691 tliston.AddToList(&ghsep); 1692 tliston.AddToList(&fillh); 1693 1694 tliston.AddToList(&selfinalgh); 1695 tliston.AddToList(&hfill1); 1696 tliston.AddToList(&hfill2); 1697 tliston.AddToList(&hfill3); 1698 tliston.AddToList(&hfill4); 1699 1700 tliston.AddToList(&selfinal); 1701 if (WMC2) 1702 tliston.AddToList(&write); 1703 1704 //***************************** 1705 1706 //------------------------------------------- 1707 // Execute event loop 1708 // 1709 MProgressBar bar; 1710 MEvtLoop evtloop; 1711 evtloop.SetParList(&pliston); 1712 evtloop.SetProgressBar(&bar); 1713 1714 Int_t maxevents = 1000000000; 1715 //Int_t maxevents = 1000; 1716 if ( !evtloop.Eventloop(maxevents) ) 1717 return; 1718 1719 tliston.PrintStatistics(0, kTRUE); 1720 1721 1722 //------------------------------------------- 1723 // Display the histograms 1724 // 1725 TObject *c; 1726 1727 c = pliston.FindObject("MHHadronness")->DrawClone(); 1728 c->Print(); 1729 1730 c = pliston.FindObject("MHHillas")->DrawClone(); 1731 1732 c = pliston.FindObject("MHHillasSrc")->DrawClone(); 1733 1734 //pliston.FindObject("MHHillasExt")->DrawClone(); 1735 c = pliston.FindObject(nxt)->DrawClone(); 1736 1737 c = pliston.FindObject("MHStarMap")->DrawClone(); 1738 1739 1740 1741 cout << "Macro AnalyseCT1 : End of Job 5" << endl; 1742 cout << "===================================================" << endl; 1743 } 1744 1745 1746 //--------------------------------------------------------------------- 1747 // Job 6 1748 //====== 1749 1750 // - read in the matrices of look-alike events for gammas and hadrons 1751 1752 // then read ON1 data file 1753 // 1754 // - perform the g/h separation 1755 // - apply the final cuts 1756 // 1757 // - write a file of ON events (ON2) 1758 // (after the g/h separation and after the final cuts) 1759 // (to be used for the flux calculation) 1760 // 1761 // - do the flux calculation 1762 1763 if (Job6) 1764 { 1765 cout << "=====================================================" << endl; 1766 cout << "Macro AnalyseCT1 : Start of Job 6" << endl; 1767 1768 cout << "" << endl; 1769 cout << "Macro AnalyseCT1 : Job6, WON2 = " 1770 << Job6 << ", " << WON2 << endl; 1771 1772 TString typeInput = "ON"; 1773 cout << "typeInput = " << typeInput << endl; 1774 1775 TString typeMC = "MC"; 1776 cout << "typeMC = " << typeMC << endl; 1777 1778 // use hadron matrix from ON or OFF data 1779 TString typeData = "ON"; 1780 //TString typeData = "OFF"; 1781 cout << "typeData = " << typeData << endl; 1782 1783 1784 //************************************************************************* 1785 // read in matrices of look-alike events 1786 1787 const char* mtxName = "MatrixGammas"; 1788 1789 TString outName = "matrix_gammas_"; 1790 outName += typeMC; 1791 outName += "_"; 1792 outName += typeData; 1793 outName += ".root"; 1794 1795 1796 cout << "" << endl; 1797 cout << "========================================================" << endl; 1798 cout << "Get matrix for (gammas)" << endl; 1799 cout << "matrix name = " << mtxName << endl; 1800 cout << "name of root file = " << outName << endl; 1801 cout << "" << endl; 1802 1803 1804 // read in the object with the name 'mtxName' from file 'outName' 1805 // 1806 TFile *fileg = new TFile(outName); 1807 1808 MHMatrix matrixg; 1809 matrixg.Read(mtxName); 1810 pmatrixg = &matrixg; 1811 1812 delete fileg; 1813 1814 //***************************************************************** 1815 1816 const char* mtxName = "MatrixHadrons"; 1817 1818 TString outName = "matrix_hadrons_"; 1819 outName += typeData; 1820 outName += ".root"; 1821 1822 1823 cout << "" << endl; 1824 cout << "========================================================" << endl; 1825 cout << " Get matrix for (hadrons)" << endl; 1826 cout << "matrix name = " << mtxName << endl; 1827 cout << "name of root file = " << outName << endl; 1828 cout << "" << endl; 1829 1830 1831 // read in the object with the name mtxName 1832 // 1833 TFile *fileh = new TFile(outName); 1834 1835 MHMatrix matrixh; 1836 matrixh.Read(mtxName); 1837 pmatrixh = &matrixh; 1838 1839 delete fileh; 1840 1841 //************************************************************************* 1842 1843 1844 MTaskList tliston; 1845 MParList pliston; 1846 pliston.AddToList(&tliston); 1847 1848 pliston.AddToList(pmatrixg); 1849 pliston.AddToList(pmatrixh); 1850 1851 //:::::::::::::::::::::::::::::::::::::::::: 1852 1853 MBinning binssize("BinningSize"); 1854 binssize.SetEdgesLog(50, 10, 1.0e5); 1855 pliston.AddToList(&binssize); 1856 1857 MBinning binsdistc("BinningDist"); 1858 binsdistc.SetEdges(50, 0, 1.4); 1859 pliston.AddToList(&binsdistc); 1860 1861 MBinning binswidth("BinningWidth"); 1862 binswidth.SetEdges(50, 0, 1.0); 1863 pliston.AddToList(&binswidth); 1864 1865 MBinning binslength("BinningLength"); 1866 binslength.SetEdges(50, 0, 1.0); 1867 pliston.AddToList(&binslength); 1868 1869 MBinning binsalpha("BinningAlpha"); 1870 binsalpha.SetEdges(100, -100, 100); 1871 pliston.AddToList(&binsalpha); 1872 1873 MBinning binsasym("BinningAsym"); 1874 binsasym.SetEdges(50, -1.5, 1.5); 1875 pliston.AddToList(&binsasym); 1876 1877 MBinning binsht("BinningHeadTail"); 1878 binsht.SetEdges(50, -1.5, 1.5); 1879 pliston.AddToList(&binsht); 1880 1881 MBinning binsm3l("BinningM3Long"); 1882 binsm3l.SetEdges(50, -1.5, 1.5); 1883 pliston.AddToList(&binsm3l); 1884 1885 MBinning binsm3t("BinningM3Trans"); 1886 binsm3t.SetEdges(50, -1.5, 1.5); 1887 pliston.AddToList(&binsm3t); 1888 1889 1890 //..... 1891 MBinning binsb("BinningSigmabar"); 1892 binsb.SetEdges( 100, 0.0, 5.0); 1893 pliston.AddToList(&binsb); 1894 1895 MBinning binth("BinningTheta"); 1896 binth.SetEdges( 70, -0.5, 69.5); 1897 pliston.AddToList(&binth); 1898 1899 MBinning binsdiff("BinningDiffsigma2"); 1900 binsdiff.SetEdges(100, -5.0, 20.0); 1901 pliston.AddToList(&binsdiff); 1902 //:::::::::::::::::::::::::::::::::::::::::: 1903 1904 1905 //------------------------------------------- 1906 // create the tasks which should be executed 1907 // 1908 1909 TString filenameData = typeInput; 1910 filenameData += "1.root"; 1911 1912 readname = "ReadCT1ONdata"; 1913 MReadMarsFile read("Events", filenameData); 1914 read.DisableAutoScheme(); 1915 1916 1917 //....................................................................... 1918 // g/h separation 1919 1920 MMultiDimDistCalc ghsep; 1921 ghsep.SetUseNumRows(25); 1922 ghsep.SetUseKernelMethod(kFALSE); 1923 1924 1925 1926 //....................................................................... 1927 1928 // geometry is needed in MHHillas... classes 1929 MGeomCam *fGeom = 1930 (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); 1931 1932 TString fHilName = "Hillas"; 1933 TString fHilSrcName = "HillasSrc"; 1934 1935 Float_t hadcut = 0.2; 1936 MSelFinal selfinalgh(fHilName, fHilSrcName); 1937 selfinalgh.SetHadronnessCut(hadcut); 1938 selfinalgh.SetAlphaCut(100.0); 1939 1940 MFillH fillh("MHHadronness"); 1941 1942 MSelFinal selfinal(fHilName, fHilSrcName); 1943 selfinal.SetHadronnessCut(hadcut); 1944 selfinal.SetAlphaCut(20.0); 1945 1946 if (WON2) 1947 { 1948 TString outNameImage = typeInput; 1949 outNameImage += "_"; 1950 outNameImage += typeData; 1951 outNameImage += "2.root"; 1952 MWriteRootFile write(outNameImage); 1953 write.AddContainer("MHadronness", "Events"); 1954 write.AddContainer("MSigmabar", "Events"); 1955 write.AddContainer("Hillas", "Events"); 1956 write.AddContainer("MMcEvt", "Events"); 1957 write.AddContainer("HillasSrc", "Events"); 1958 write.AddContainer("MRawRunHeader", "RunHeaders"); 1959 //write.AddContainer("MMcRunHeader","RunHeaders"); 1960 write.AddContainer("MSrcPosCam", "RunHeaders"); 1961 } 1962 1963 1964 MFillH hfill1("MHHillas", fHilName); 1965 MFillH hfill2("MHStarMap", fHilName); 1966 1967 TString nxt = "HillasExt"; 1968 MHHillasExt hhilext(nxt, "", fHilName); 1969 pliston.AddToList(&hhilext); 1970 TString namext = nxt; 1971 namext += "[MHHillasExt]"; 1972 MFillH hfill3(namext, fHilSrcName); 1973 1974 MFillH hfill4("MHHillasSrc", fHilSrcName); 1975 1976 1977 1978 //***************************** 1979 // tasks to be executed 1980 1981 tliston.AddToList(&read); 1982 1983 tliston.AddToList(&ghsep); 1984 tliston.AddToList(&fillh); 1985 1986 tliston.AddToList(&selfinalgh); 1987 1988 tliston.AddToList(&hfill1); 1989 tliston.AddToList(&hfill2); 1990 tliston.AddToList(&hfill3); 1991 tliston.AddToList(&hfill4); 1992 1993 tliston.AddToList(&selfinal); 1994 if (WON2) 1995 tliston.AddToList(&write); 1996 1997 1998 1999 //***************************** 2000 2001 //------------------------------------------- 2002 // Execute event loop 2003 // 2004 MProgressBar bar; 2005 MEvtLoop evtloop; 2006 evtloop.SetParList(&pliston); 2007 evtloop.SetProgressBar(&bar); 2008 2009 Int_t maxevents = 1000000000; 2010 //Int_t maxevents = 1000; 2011 if ( !evtloop.Eventloop(maxevents) ) 2012 return; 2013 2014 tliston.PrintStatistics(0, kTRUE); 2015 2016 2017 //------------------------------------------- 2018 // Display the histograms 2019 // 2020 TObject *c; 2021 2022 c = pliston.FindObject("MHHadronness")->DrawClone(); 2023 c->Print(); 2024 2025 //c = pliston.FindObject("MHHillas")->DrawClone(); 2026 2027 c = pliston.FindObject("MHHillasSrc")->DrawClone(); 2028 2029 ////pliston.FindObject("MHHillasExt")->DrawClone(); 2030 //c = pliston.FindObject(nxt)->DrawClone(); 2031 2032 c = pliston.FindObject("MHStarMap")->DrawClone(); 2033 2034 2035 2036 cout << "Macro AnalyseCT1 : End of Job 6" << endl; 2037 cout << "=======================================================" << endl; 2038 } 2039 //--------------------------------------------------------------------- 718 2040 } 719 2041 … … 721 2043 722 2044 2045 2046 2047 2048 2049 2050 -
trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc
r1660 r1809 208 208 } 209 209 210 //fHadronness->SetHadronness(dg/(dg+dh)); 211 fHadronness->SetHadronness(exp(-dh/dg)); 210 Double_t arg; 211 212 if (dg+dh != 0.0) 213 arg = dg / (dg+dh); 214 else 215 arg = 1.e10; 216 //fHadronness->SetHadronness(arg); 217 218 if (dg != 0.0) 219 arg = exp(-dh/dg); 220 else 221 arg = 0.0; 222 fHadronness->SetHadronness(arg); 223 212 224 213 225 return kTRUE; -
trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc
r1772 r1809 1 2 1 /* ======================================================================== *\ 3 2 ! … … 102 101 fHDiffPixTheta = fHist3Diff; 103 102 103 fHSigmaTheta->SetDirectory(NULL); 104 fHSigmaPixTheta->SetDirectory(NULL); 105 fHDiffPixTheta->SetDirectory(NULL); 106 104 107 Print(); 105 108 } … … 239 242 fHDiffPixTh->SetZTitle("Sigma^2-Sigmabar^2"); 240 243 244 //-------------------------------------------------------------------- 245 246 memset(fErrors, 0, sizeof(fErrors)); 247 241 248 return kTRUE; 242 249 } … … 251 258 //*fLog << "Entry MPadSchweizer::Process();" << endl; 252 259 260 Int_t rc; 261 253 262 const UInt_t npix = fEvt->GetNumPixels(); 263 264 Double_t SigbarOld; 265 266 //*fLog << "before padding : " << endl; 267 //SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); 268 //fSigmabar->Print(""); 269 254 270 255 271 //$$$$$$$$$$$$$$$$$$$$$$$$$$ 256 272 // to simulate the situation that before the padding the NSB and 257 273 // electronic noise are zero : set Sigma = 0 for all pixels 258 for (UInt_t i=0; i<npix; i++)259 {260 MCerPhotPix &pix = fEvt->operator[](i);261 Int_t j = pix.GetPixId();262 263 MPedestalPix &ppix = fPed->operator[](j);264 ppix.SetMeanRms(0.0);265 }274 //for (UInt_t i=0; i<npix; i++) 275 //{ 276 // MCerPhotPix &pix = fEvt->operator[](i); 277 // Int_t j = pix.GetPixId(); 278 279 // MPedestalPix &ppix = fPed->operator[](j); 280 // ppix.SetMeanRms(0.0); 281 //} 266 282 //$$$$$$$$$$$$$$$$$$$$$$$$$$ 267 283 … … 269 285 // Calculate average sigma of the event 270 286 // 271 Double_tSigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt);287 SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); 272 288 //fSigmabar->Print(""); 273 289 274 //if (SigbarOld > 0.0)275 //{290 if (SigbarOld > 0.0) 291 { 276 292 //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : " 277 293 // << SigbarOld << ". Stop event loop " << endl; 278 294 // input data should have Sigmabar = 0; stop event loop 279 // return kFALSE; 280 //} 295 296 rc = 1; 297 fErrors[rc]++; 298 return kCONTINUE; 299 } 281 300 282 301 Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); … … 296 315 << Theta << ", " << binNumber << "; Skip event " << endl; 297 316 // event cannot be padded; skip event 317 318 rc = 2; 319 fErrors[rc]++; 298 320 return kCONTINUE; 299 321 } … … 308 330 // event cannot be padded; skip event 309 331 delete fHSigma; 332 333 rc = 3; 334 fErrors[rc]++; 310 335 return kCONTINUE; 311 336 } … … 331 356 *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : " 332 357 << Sigmabar << ", " << SigbarOld << ", Skip event" << endl; 358 359 rc = 4; 360 fErrors[rc]++; 333 361 return kCONTINUE; 334 362 } … … 426 454 << binTheta << " and pixel bin " << binPixel 427 455 << " has no entries; aborting " << endl; 428 return kERROR; 456 457 rc = 5; 458 fErrors[rc]++; 459 return kCONTINUE; 429 460 } 430 461 … … 465 496 << binTheta << " and pixel bin " << binPixel 466 497 << " has no entries; aborting " << endl; 467 return kERROR; 498 499 rc = 6; 500 fErrors[rc]++; 501 return kCONTINUE; 468 502 } 469 503 … … 553 587 554 588 // Calculate Sigmabar again and crosscheck 589 555 590 Double_t SigbarNew = fSigmabar->Calc(*fCam, *fPed, *fEvt); 591 //*fLog << "after padding : " << endl; 556 592 //fSigmabar->Print(""); 557 593 … … 585 621 //*fLog << "Exit MPadSchweizer::Process();" << endl; 586 622 623 rc = 0; 624 fErrors[rc]++; 625 587 626 return kTRUE; 588 627 … … 594 633 Bool_t MPadSchweizer::PostProcess() 595 634 { 635 if (GetNumExecutions() != 0) 636 { 637 638 *fLog << inf << endl; 639 *fLog << GetDescriptor() << " execution statistics:" << endl; 640 *fLog << dec << setfill(' '); 641 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 642 << (int)(fErrors[1]*100/GetNumExecutions()) 643 << "%) Evts skipped due to: Sigmabar_old > 0" << endl; 644 645 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 646 << (int)(fErrors[2]*100/GetNumExecutions()) 647 << "%) Evts skipped due to: Zenith angle out of range" << endl; 648 649 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) 650 << (int)(fErrors[3]*100/GetNumExecutions()) 651 << "%) Evts skipped due to: No data for generating Sigmabar" << endl; 652 653 *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) 654 << (int)(fErrors[4]*100/GetNumExecutions()) 655 << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; 656 657 *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) 658 << (int)(fErrors[5]*100/GetNumExecutions()) 659 << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; 660 661 *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) 662 << (int)(fErrors[6]*100/GetNumExecutions()) 663 << "%) Evts skipped due to: No data for generating Sigma" << endl; 664 665 *fLog << " " << fErrors[0] << " (" 666 << (int)(fErrors[0]*100/GetNumExecutions()) 667 << "%) Evts survived the padding!" << endl; 668 *fLog << endl; 669 670 } 671 672 //--------------------------------------------------------------- 596 673 TCanvas &c = *(MH::MakeDefCanvas("PadSchweizer", "", 900, 1200)); 597 674 c.Divide(3, 4); … … 600 677 601 678 c.cd(1); 679 fHSigmaTheta->SetDirectory(NULL); 602 680 fHSigmaTheta->SetTitle("2D : Sigmabar, \\Theta (reference sample)"); 603 fHSigmaTheta->DrawC lone();681 fHSigmaTheta->DrawCopy(); 604 682 fHSigmaTheta->SetBit(kCanDelete); 605 683 606 684 //c.cd(1); 607 //fHSigmaPixTheta->DrawC lone();685 //fHSigmaPixTheta->DrawCopy(); 608 686 //fHSigmaPixTheta->SetBit(kCanDelete); 609 687 610 688 c.cd(4); 611 fHSigbarTheta->DrawClone(); 689 fHSigbarTheta->SetDirectory(NULL); 690 fHSigbarTheta->DrawCopy(); 612 691 fHSigbarTheta->SetBit(kCanDelete); 613 692 614 693 615 694 c.cd(7); 616 fHNSB->DrawClone(); 695 fHNSB->SetDirectory(NULL); 696 fHNSB->DrawCopy(); 617 697 fHNSB->SetBit(kCanDelete); 618 698 … … 625 705 c.cd(2); 626 706 l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zx"); 627 707 l->SetDirectory(NULL); 628 708 l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); 629 709 l->SetXTitle("\\Theta [\\circ]"); 630 710 l->SetYTitle("Sigma^2-Sigmabar^2"); 631 711 632 *fLog << "before box" << endl; 633 634 l->Draw("box"); 712 l->DrawCopy("box"); 635 713 l->SetBit(kCanDelete);; 636 714 637 715 c.cd(5); 638 716 l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zy"); 717 l->SetDirectory(NULL); 639 718 l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); 640 719 l->SetXTitle("pixel"); 641 720 l->SetYTitle("Sigma^2-Sigmabar^2"); 642 721 643 l->Draw ("box");722 l->DrawCopy("box"); 644 723 l->SetBit(kCanDelete);; 645 724 … … 655 734 c.cd(3); 656 735 k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zx"); 736 k->SetDirectory(NULL); 657 737 k->SetTitle("Sigma vs. \\Theta (all pixels)"); 658 738 k->SetXTitle("\\Theta [\\circ]"); 659 739 k->SetYTitle("Sigma"); 660 740 661 k->Draw ("box");741 k->DrawCopy("box"); 662 742 k->SetBit(kCanDelete); 663 743 664 744 c.cd(6); 665 745 k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zy"); 746 k->SetDirectory(NULL); 666 747 k->SetTitle("Sigma vs. pixel number (all \\Theta)"); 667 748 k->SetXTitle("pixel"); 668 749 k->SetYTitle("Sigma"); 669 750 670 k->Draw ("box");751 k->DrawCopy("box"); 671 752 k->SetBit(kCanDelete);; 672 753 … … 677 758 678 759 c.cd(10); 679 fHSigmaPedestal->DrawClone(); 760 fHSigmaPedestal->SetDirectory(NULL); 761 fHSigmaPedestal->DrawCopy(); 680 762 fHSigmaPedestal->SetBit(kCanDelete); 681 763 682 764 c.cd(11); 683 fHPhotons->DrawClone(); 765 fHPhotons->SetDirectory(NULL); 766 fHPhotons->DrawCopy(); 684 767 fHPhotons->SetBit(kCanDelete); 685 768 -
trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h
r1771 r1809 39 39 Int_t fRunType; 40 40 Int_t fGroup; 41 42 Int_t fErrors[7]; 41 43 42 44 // plots used for the padding -
trunk/MagicSoft/Mars/manalysis/MSelBasic.cc
r1762 r1809 44 44 #include "MGeomCam.h" 45 45 #include "MPedestalCam.h" 46 #include "MPedestalPix.h" 46 47 #include "MGeomPix.h" 47 48 … … 119 120 Bool_t MSelBasic::Process() 120 121 { 122 /* 123 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 124 *fLog << "=========================================================" << endl; 125 *fLog << "" << endl; 126 *fLog << "MmcEvt data : " << endl; 127 *fLog << "" << endl; 128 fMcEvt->Print(); 129 *fLog << "" << endl; 130 *fLog << "PartId() = " << fMcEvt->GetPartId() << endl; 131 *fLog << "Energy() = " << fMcEvt->GetEnergy() << endl; 132 *fLog << "Theta() = " << fMcEvt->GetTheta() << endl; 133 134 *fLog << "Phi() = " << fMcEvt->GetPhi() << endl; 135 //*fLog << "CoreD() = " << fMcEvt->GetCoreD() << endl; 136 //*fLog << "CoreX() = " << fMcEvt->GetCoreX() << endl; 137 138 //*fLog << "CoreY() = " << fMcEvt->GetCoreY() << endl; 139 *fLog << "Impact() = " << fMcEvt->GetImpact() << endl; 140 //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl; 141 142 //*fLog << "PassPhotAtm() = " << fMcEvt->GetPassPhotAtm() << endl; 143 //*fLog << "PassPhotRef() = " << fMcEvt->GetPassPhotRef() << endl; 144 //*fLog << "PassPhotCone() = " << fMcEvt->GetPassPhotCone() << endl; 145 146 //*fLog << "PhotElfromShower() = " << fMcEvt->GetPhotElfromShower() << endl; 147 //*fLog << "PhotElinCamera() = " << fMcEvt->GetPhotElinCamera() << endl; 148 *fLog << "TelescopePhi() = " << fMcEvt->GetTelescopePhi() << endl; 149 150 *fLog << "TelescopeTheta() = " << fMcEvt->GetTelescopeTheta() << endl; 151 *fLog << "Impact() = " << fMcEvt->GetImpact() << endl; 152 //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl; 153 *fLog << "" << endl; 154 *fLog << "=========================================================" << endl; 155 156 *fLog << "=========================================================" << endl; 157 *fLog << "" << endl; 158 *fLog << "MPedestalPix data : " << endl; 159 *fLog << "" << endl; 160 161 Int_t ntot; 162 ntot = fPed->GetSize(); 163 *fLog << "MeanRms() :" << endl; 164 for (Int_t i=0; i<ntot; i++) 165 { 166 MPedestalPix &pix = (*fPed)[i]; 167 //*fLog << "Mean() = " << i << ", " << pix.GetMean() << endl; 168 //*fLog << "Sigma() = " << i << ", " << pix.GetSigma() << endl; 169 *fLog << pix.GetMeanRms() << " "; 170 //*fLog << "SigmaRms()= " << i << ", " << pix.GetSigmaRms() << endl; 171 } 172 *fLog << "" << endl; 173 174 *fLog << "" << endl; 175 ntot = fEvt->GetNumPixels(); 176 *fLog << "MCerPhotPix : pix.GetNumPhotons()" << endl; 177 for (Int_t i=0; i<ntot; i++) 178 { 179 MCerPhotPix &pix = (*fEvt)[i]; 180 *fLog << pix.GetNumPhotons() << " "; 181 } 182 *fLog << "" << endl; 183 184 *fLog << "" << endl; 185 ntot = fEvt->GetNumPixels(); 186 *fLog << "MCerPhotPix : pix.GetErrorPhot()" << endl; 187 for (Int_t i=0; i<ntot; i++) 188 { 189 MCerPhotPix &pix = (*fEvt)[i]; 190 *fLog << pix.GetErrorPhot() << " "; 191 } 192 *fLog << "" << endl; 193 194 //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 195 */ 196 197 121 198 Int_t rc = 0; 122 199 123 //if ( fRawRun->GetRunNumber() < 1 025)200 //if ( fRawRun->GetRunNumber() < 16279 ) 124 201 //{ 125 202 // rc = 1; … … 128 205 129 206 Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); 130 if (Theta > 45.0 || !SwTrigger() ) 131 { 132 //*fLog << "MSelBasic::Process; Theta = " << Theta << endl; 207 if ( Theta < 0.0 ) 208 { 209 *fLog << "MSelBasic::Process; Run, Event, Theta = " 210 << fRawRun->GetRunNumber()<< ", " 211 << fMcEvt->GetEvtNumber() << ", " << Theta << endl; 133 212 rc = 1; 213 } 214 215 else if ( Theta > 45.0 ) 216 { 217 rc = 2; 218 } 219 220 else if ( !SwTrigger() ) 221 { 222 //*fLog << "MSelBasic::Process; SwTrigger = " << SwTrigger << endl; 223 rc = 3; 134 224 } 135 225 … … 223 313 *fLog << GetDescriptor() << " execution statistics:" << endl; 224 314 *fLog << dec << setfill(' '); 225 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Basic selections are not fullfilled" << endl; 315 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle < 0" << endl; 316 317 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle too high" << endl; 318 319 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Software trigger not fullfilled" << endl; 226 320 227 321 *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Basic selections!" << endl; … … 230 324 return kTRUE; 231 325 } 326 -
trunk/MagicSoft/Mars/manalysis/MSelBasic.h
r1762 r1809 30 30 const MRawRunHeader *fRawRun; 31 31 32 Int_t fErrors[ 2];32 Int_t fErrors[4]; 33 33 34 34 public: -
trunk/MagicSoft/Mars/manalysis/MSelFinal.cc
r1781 r1809 41 41 42 42 #include "MHillas.h" 43 #include "MHillasExt.h" 43 44 #include "MHillasSrc.h" 44 45 #include "MCerPhotEvt.h" … … 57 58 // Default constructor. 58 59 // 59 MSelFinal::MSelFinal( MHillas *parhil, MHillasSrc *parhilsrc,60 MSelFinal::MSelFinal(const char *HilName, const char *HilSrcName, 60 61 const char *name, const char *title) 61 62 { … … 63 64 fTitle = title ? title : "Task to evaluate the Final Cuts"; 64 65 65 fHil = parhil; 66 fHilsrc = parhilsrc; 66 fHilName = HilName; 67 fHilSrcName = HilSrcName; 68 69 fHadronnessCut = 0.2; 70 fAlphaCut = 100.0; //degrees 67 71 } 68 72 … … 75 79 Bool_t MSelFinal::PreProcess(MParList *pList) 76 80 { 81 fHil = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt"); 82 if (!fHil) 83 { 84 *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl; 85 return kFALSE; 86 } 87 88 fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); 89 if (!fHilSrc) 90 { 91 *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl; 92 return kFALSE; 93 } 94 95 77 96 fHadronness = (MHadronness*)pList->FindObject("MHadronness"); 78 97 if (!fHadronness) … … 90 109 } 91 110 92 fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");93 if (!fEvt)94 {95 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;96 return kFALSE;97 }98 99 100 fCam = (MGeomCam*)pList->FindObject("MGeomCam");101 if (!fCam)102 {103 *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;104 return kFALSE;105 }106 fMm2Deg = fCam->GetConvMm2Deg();107 108 *fLog << "fMm2Deg = " << fMm2Deg << endl;109 110 111 memset(fErrors, 0, sizeof(fErrors)); 111 112 … … 122 123 Bool_t MSelFinal::Process() 123 124 { 125 //*fLog << "Entry MSelFinal; fHilSrc = " << fHilSrc << endl; 126 127 128 124 129 Int_t rc = 0; 125 130 126 Double_t alphacut = 20.0;131 Double_t modalpha = fabs( fHilSrc->GetAlpha() ); 127 132 128 Double_t modalpha = fabs( fHilsrc->GetAlpha() );129 133 Double_t h = fHadronness->GetHadronness(); 130 134 131 if ( h> 0.5 || modalpha > alphacut )135 if ( h>fHadronnessCut ) 132 136 { 133 137 //*fLog << "MSelFinal::Process; h, alpha = " << h << ", " 134 // << fHil src->GetAlpha() << endl;138 // << fHilSrc->GetAlpha() << endl; 135 139 rc = 1; 140 } 141 142 else if ( modalpha > fAlphaCut ) 143 { 144 //*fLog << "MSelFinal::Process; h, alpha = " << h << ", " 145 // << fHilSrc->GetAlpha() << endl; 146 rc = 2; 136 147 } 137 148 … … 153 164 *fLog << GetDescriptor() << " execution statistics:" << endl; 154 165 *fLog << dec << setfill(' '); 155 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Final selections are not fullfilled" << endl; 166 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 167 << (int)(fErrors[1]*100/GetNumExecutions()) 168 << "%) Evts skipped due to: g/h separation cut (" << fHadronnessCut 169 << ")" << endl; 156 170 157 *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Final selections!" << endl; 171 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 172 << (int)(fErrors[2]*100/GetNumExecutions()) 173 << "%) Evts skipped due to: cut in ALPHA (" << fAlphaCut 174 << " degrees)" << endl; 175 176 *fLog << " " << fErrors[0] << " (" 177 << (int)(fErrors[0]*100/GetNumExecutions()) 178 << "%) Evts survived Final selections!" << endl; 158 179 *fLog << endl; 159 180 160 181 return kTRUE; 161 182 } 183 184 185 -
trunk/MagicSoft/Mars/manalysis/MSelFinal.h
r1781 r1809 28 28 MMcEvt *fMcEvt; 29 29 MHillas *fHil; 30 MHillasSrc *fHil src;30 MHillasSrc *fHilSrc; 31 31 MHadronness *fHadronness; 32 32 33 33 Double_t fMm2Deg; // conversion mm to degrees in camera 34 Int_t fErrors[2]; 34 Int_t fErrors[3]; 35 TString fHilName; 36 TString fHilSrcName; 37 38 Float_t fHadronnessCut; 39 Float_t fAlphaCut; 35 40 36 41 public: 37 MSelFinal( MHillas *fHil, MHillasSrc *fHilsrc,42 MSelFinal(const char *HilName, const char *HilSrcName, 38 43 const char *name=NULL, const char *title=NULL); 39 44 … … 42 47 Bool_t PostProcess(); 43 48 49 void SetHadronnessCut(Float_t hadcut) { fHadronnessCut = hadcut; } 50 void SetAlphaCut(Float_t alpha) { fAlphaCut = alpha; } 44 51 45 52 ClassDef(MSelFinal, 0) // Task to evaluate final cuts -
trunk/MagicSoft/Mars/manalysis/MSelStandard.cc
r1762 r1809 1 1 2 /* ======================================================================== *\ 2 3 ! … … 39 40 40 41 #include "MHillas.h" 42 #include "MHillasExt.h" 41 43 #include "MHillasSrc.h" 42 44 #include "MCerPhotEvt.h" … … 54 56 // Default constructor. 55 57 // 56 MSelStandard::MSelStandard(const MHillas *parhil, const MHillasSrc *parhilsrc,58 MSelStandard::MSelStandard(const char *HilName, const char *HilSrcName, 57 59 const char *name, const char *title) 58 60 { … … 60 62 fTitle = title ? title : "Task to evaluate the Standard Cuts"; 61 63 62 fHil = parhil;63 fHil src = parhilsrc;64 fHilName = HilName; 65 fHilSrcName = HilSrcName; 64 66 } 65 67 … … 72 74 Bool_t MSelStandard::PreProcess(MParList *pList) 73 75 { 76 fHil = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt"); 77 if (!fHil) 78 { 79 *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl; 80 return kFALSE; 81 } 82 83 fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); 84 if (!fHilSrc) 85 { 86 *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl; 87 return kFALSE; 88 } 89 90 74 91 fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); 75 92 if (!fMcEvt) … … 95 112 fMm2Deg = fCam->GetConvMm2Deg(); 96 113 97 *fLog << "fMm2Deg = " << fMm2Deg << endl;114 //*fLog << "fMm2Deg = " << fMm2Deg << endl; 98 115 99 116 memset(fErrors, 0, sizeof(fErrors)); … … 114 131 Int_t rc = 0; 115 132 116 //Double_t fLength = fHil->GetLength() * fMm2Deg;117 //Double_t fWidth = fHil->GetWidth() * fMm2Deg;118 Double_t fDist = fHil src->GetDist()* fMm2Deg;133 Double_t fLength = fHil->GetLength() * fMm2Deg; 134 Double_t fWidth = fHil->GetWidth() * fMm2Deg; 135 Double_t fDist = fHilSrc->GetDist()* fMm2Deg; 119 136 //Double_t fDelta = fHil->GetDelta() * kRad2Deg; 120 137 Double_t fSize = fHil->GetSize(); … … 122 139 Int_t fNumCorePixels = fHil->GetNumCorePixels(); 123 140 124 if ( fSize <= 60.0 || fDist< 0.4 || fDist > 1.1 125 || fNumUsedPixels >= 92 || fNumCorePixels < 4) 141 if ( fNumUsedPixels >= 92 || fNumCorePixels < 4 ) 126 142 { 127 143 //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = " … … 129 145 // << fNumCorePixels << endl; 130 146 rc = 1; 147 } 148 149 else if ( fSize <= 60.0 || fDist< 0.4 || fDist > 1.1 ) 150 { 151 //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = " 152 // << fSize << ", " << fDist << ", " << fNumUsedPixels << ", " 153 // << fNumCorePixels << endl; 154 rc = 2; 155 } 156 157 else if ( fLength <= 0.0 || fWidth <= 0.0 ) 158 { 159 //*fLog << "MSelStandard::Process; fLength, fWidth = " 160 // << fLength << ", " << fWidth << endl; 161 rc = 3; 131 162 } 132 163 … … 149 180 *fLog << GetDescriptor() << " execution statistics:" << endl; 150 181 *fLog << dec << setfill(' '); 151 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Standard selections are not fullfilled" << endl; 182 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on no.of used or core pxels not fullfilled" << endl; 183 184 *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on SIZE or DIST not fullfilled" << endl; 185 186 *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Length or Width is <= 0" << endl; 152 187 153 188 *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Standard selections!" << endl; … … 156 191 return kTRUE; 157 192 } 193 194 -
trunk/MagicSoft/Mars/manalysis/MSelStandard.h
r1762 r1809 23 23 { 24 24 private: 25 constMGeomCam *fCam; // Camera Geometry26 constMCerPhotEvt *fEvt; // Cerenkov Photon Event27 constMMcEvt *fMcEvt;28 constMHillas *fHil;29 const MHillasSrc *fHilsrc;25 MGeomCam *fCam; // Camera Geometry 26 MCerPhotEvt *fEvt; // Cerenkov Photon Event 27 MMcEvt *fMcEvt; 28 MHillas *fHil; 29 MHillasSrc *fHilSrc; 30 30 31 Double_t fMm2Deg; // conversion mm to degrees in camera 32 Int_t fErrors[2]; 31 Double_t fMm2Deg; // conversion mm to degrees in camera 32 Int_t fErrors[4]; 33 TString fHilName; 34 TString fHilSrcName; 33 35 34 36 public: 35 MSelStandard(const MHillas *fHil, const MHillasSrc *fHilsrc,37 MSelStandard(const char *HilName, const char *HilSrcName, 36 38 const char *name=NULL, const char *title=NULL); 37 39 -
trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc
r1749 r1809 122 122 { 123 123 const char *name = gSystem->ExpandPathName(txt); 124 124 125 TString fname(name); 125 126 delete [] name; … … 336 337 // float frms_pedsig_phot[iMAXNUMPIX]; // standard deviation of the calibrated signals from the pedestal run */ 337 338 fPedest->InitSize(iMAXNUMPIX); 339 *fLog << "PedestalRMS : "; 338 340 for (Int_t i=0; i<iMAXNUMPIX; i++) 341 { 339 342 (*fPedest)[i].SetMeanRms(outpars.frms_pedsig_phot[i]); 343 *fLog << outpars.frms_pedsig_phot[i] << " "; 344 345 //$$$$$$$$$$$$$$$$$$$$$$$$$ 346 savePedRMS[i] = outpars.frms_pedsig_phot[i]; 347 //$$$$$$$$$$$$$$$$$$$$$$$$$ 348 } 349 *fLog << endl; 340 350 341 351 fPedest->SetReadyToSave(); … … 373 383 * be treated further. */ 374 384 385 *fLog << "outpars.bmontecarlo = " << outpars.bmontecarlo << endl; 386 *fLog << "outpars.imcparticle = " << outpars.imcparticle << endl; 387 *fLog << "outpars.dsourcera_hours = " << outpars.dsourcera_hours << endl; 388 *fLog << "outpars.dsourcedec_deg = " << outpars.dsourcedec_deg << endl; 389 390 //*fLog << "File is a "; 391 //*fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data"); 392 //*fLog << " file." << endl; 393 394 395 // Next statement commented out because bmontecarlo was set wrongly 396 //fIsMcFile = outpars.bmontecarlo==TRUE; 397 fIsMcFile = (outpars.dsourcera_hours == 0.0 && 398 outpars.dsourcedec_deg == 0.0 && 399 outpars.imcparticle != 0 ); 400 375 401 *fLog << "File is a "; 376 *fLog << ( outpars.bmontecarlo? "Monte Carlo" : "Real Data");402 *fLog << (fIsMcFile ? "Monte Carlo" : "Real Data"); 377 403 *fLog << " file." << endl; 378 379 fIsMcFile = outpars.bmontecarlo==TRUE; 380 404 *fLog << " " << endl; 381 405 } 382 406 … … 433 457 ProcessRunHeader(outpars); 434 458 459 435 460 //rwagner: ReInit whenever new run commences 436 461 // rc==-1 means: ReInit didn't work out 462 437 463 if (!fTaskList->ReInit(fParList)) 438 464 return -1; … … 453 479 TString m = cEND_EVENTS_TEMPLATE; 454 480 Int_t p = m.First('%'); 481 455 482 456 483 if (!s.BeginsWith(m(0,p))) … … 515 542 516 543 // 517 // Check for the exist ance of a next file to read544 // Check for the existence of a next file to read 518 545 // 519 546 TNamed *file = (TNamed*)fFileNames->First(); … … 538 565 539 566 if (!CheckHeader(fname)) 540 return kFALSE; 567 { 568 *fLog << "OpenNextFile : CheckHeader(fname) is FALSE" << endl; 569 return kFALSE; 570 } 541 571 542 572 fIn = new ifstream(fname); 543 573 544 574 *fLog << inf << "-----------------------------------------------------------------------" << endl; 575 545 576 546 577 switch (ReadRunHeader()) … … 553 584 return kFALSE; 554 585 default: 586 *fLog << "OpenNextFile : after ReadRunHeader, FnumEventsInRun = " 587 << fNumEventsInRun << endl; 555 588 return kTRUE; 556 589 } 590 591 592 557 593 } 558 594 … … 752 788 void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event) 753 789 { 790 if (fRawRunHeader->GetRunNumber() == 1) 791 { 792 *fLog << "eventrecord" << endl; 793 *fLog << "isecs_since_midday = " << event.isecs_since_midday << endl; 794 *fLog << "isecfrac_200ns = " << event.isecfrac_200ns << endl; 795 *fLog << "snot_ok_flags = " << event.snot_ok_flags << endl; 796 *fLog << "ialt_arcs = " << event.ialt_arcs << endl; 797 *fLog << "iaz_arcs = " << event.iaz_arcs << endl; 798 *fLog << "ipreproc_alt_arcs = " << event.ipreproc_alt_arcs << endl; 799 *fLog << "ipreproc_az_arcs = " << event.ipreproc_az_arcs << endl; 800 *fLog << "ifieldrot_arcs = " << event.ifieldrot_arcs << endl; 801 802 *fLog << "srate_millihz = " << event.srate_millihz << endl; 803 *fLog << "fhourangle = " << event.fhourangle << endl; 804 *fLog << "fmcenergy_tev = " << event.fmcenergy_tev << endl; 805 *fLog << "fmcsize_phel = " << event.fmcsize_phel << endl; 806 *fLog << "imcimpact_m = " << event.imcimpact_m << endl; 807 *fLog << "imcparticle = " << event.imcparticle << endl; 808 *fLog << "imctriggerflag = " << event.imctriggerflag << endl; 809 } 810 811 812 for (Int_t i=0; i<iMAXNUMPIX; i++) 813 { 814 (*fPedest)[i].SetMeanRms(savePedRMS[i]); 815 } 816 817 754 818 // int isecs_since_midday; // seconds passed since midday before sunset (JD of run start) 755 819 // int isecfrac_200ns; // fractional part of isecs_since_midday … … 788 852 // actual number of pixels (outputpars.inumpixels) is written out 789 853 // short spixsig_10thphot[iMAXNUMPIX]; 854 //*fLog << "spixsig_10thphot : " << endl; 790 855 for (Int_t i=0; i<iMAXNUMPIX; i++) 791 856 { 857 //*fLog << event.spixsig_10thphot[i] << " "; 858 792 859 if (event.spixsig_10thphot[i]==0) 793 860 continue; … … 796 863 (*fPedest)[i].GetMeanRms()); 797 864 } 865 //*fLog << "" << endl; 866 798 867 fNphot->SetReadyToSave(); 799 868 … … 801 870 // int ipreproc_az_arcs; // "should be" az according to preproc (arcseconds) 802 871 803 fMcEvt->Fill( 0, /*fEvtNum*/872 fMcEvt->Fill(event.isecs_since_midday, //0, /*fEvtNum*/ 804 873 fIsMcFile ? event.imcparticle : 0, /*corsika particle type*/ 805 874 fIsMcFile ? event.fmcenergy_tev*1000 : 0, … … 867 936 // an event 868 937 // 869 if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0]) 938 939 // "goto TONI" was introduced in order to check for a footer record 940 // after reading a run header; this is necessary for runs with 941 // zero events after the filter 942 TONI: 943 944 if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0]) 870 945 return kTRUE; 871 946 … … 874 949 // must be an event 875 950 // 951 876 952 switch (ReadRunFooter()) 877 953 { … … 906 982 *fLog << "-----------------------------------------------------------------------" << endl; 907 983 984 908 985 if (ReadRunHeader() < 0) 909 986 { … … 911 988 return kFALSE; 912 989 } 990 991 goto TONI; 992 913 993 return kTRUE; 914 994 } -
trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h
r1735 r1809 48 48 UInt_t fNumFilterEvts; // number of events mentioned in the runs footers 49 49 50 Float_t savePedRMS[127]; 51 52 50 53 Bool_t OpenNextFile(); 51 54 -
trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc
r1762 r1809 118 118 Bool_t MFEventSelector::PreProcess(MParList *plist) 119 119 { 120 memset(fErrors, 0, sizeof(fErrors)); 121 120 122 fNumSelectedEvts = 0; 121 123 if (fSelRatio>0) … … 153 155 Bool_t MFEventSelector::Process() 154 156 { 157 Int_t rc; 158 155 159 const Float_t evt = gRandom->Uniform(); 156 160 … … 161 165 162 166 if (fResult) 167 { 163 168 fNumSelectedEvts++; 164 169 170 rc = 0; 171 fErrors[rc]++; 172 return kTRUE; 173 } 174 175 rc = 1; 176 fErrors[rc]++; 177 165 178 return kTRUE; 166 179 } … … 172 185 Bool_t MFEventSelector::PostProcess() 173 186 { 187 //--------------------------------- 188 if (GetNumExecutions() != 0) 189 { 190 *fLog << inf << endl; 191 *fLog << GetDescriptor() << " execution statistics:" << endl; 192 *fLog << dec << setfill(' '); 193 *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 194 << (int)(fErrors[1]*100/GetNumExecutions()) 195 << "%) Events not selected" << endl; 196 197 *fLog << " " << fErrors[0] << " (" 198 << (int)(fErrors[0]*100/GetNumExecutions()) 199 << "%) Events selected!" << endl; 200 *fLog << endl; 201 } 202 203 //--------------------------------- 174 204 if (TestBit(kNumTotalFromFile)) 175 205 fNumTotalEvts = -1; … … 199 229 200 230 } 201 -
trunk/MagicSoft/Mars/mfilter/MFEventSelector.h
r1762 r1809 37 37 */ 38 38 39 Int_t fErrors[2]; 40 39 41 public: 40 42 // MFEventSelector(); -
trunk/MagicSoft/Mars/mhist/MHHadronness.cc
r1668 r1809 235 235 } 236 236 237 return val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5); 237 Float_t retValue; 238 if (val2x-val1x != 0.0) 239 retValue = val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5); 240 else 241 retValue = 0.0; 242 243 return retValue; 238 244 } 239 245 … … 252 258 fQfac->Set(n); 253 259 254 const Stat_t sumg = fGhness->Integral(1, n+1); 255 const Stat_t sump = fPhness->Integral(1, n+1); 260 Stat_t sumg; 261 Stat_t sump; 262 263 sumg = fGhness->Integral(1, n); 264 sump = fPhness->Integral(1, n); 265 266 if (sumg == 0.0 || sump == 0.0) 267 { 268 *fLog << "MHHadronness::Finalize; sumg or sump is zero; sumg, sump = " 269 << sumg << ", " << sump << ". Cannot calculate hadronness" 270 << endl; 271 } 272 273 274 // Normalize photon distribution 275 Stat_t con; 276 if (sumg > 0.0) 277 for (Int_t i=1; i<=n; i++) 278 { 279 con = (fGhness->GetBinContent(i)) / sumg; 280 fGhness->SetBinContent(i, con); 281 } 282 283 // Normalize hadron distribution 284 if (sump > 0.0) 285 for (Int_t i=1; i<=n; i++) 286 { 287 con = (fPhness->GetBinContent(i)) / sump; 288 fPhness->SetBinContent(i, con); 289 } 290 291 // Calculate acceptances 292 sumg = fGhness->Integral(1, n); 293 sump = fPhness->Integral(1, n); 294 295 *fLog << "MHHadronness::Finalize; sumg, sump = " << sumg << ", " 296 << sump << endl; 256 297 257 298 Float_t max=0; … … 259 300 for (Int_t i=1; i<=n; i++) 260 301 { 261 const Stat_t ip = fPhness->Integral(1, i)/sump; 262 const Stat_t ig = fGhness->Integral(1, i)/sumg; 302 Stat_t ip; 303 if (sump != 0.0) 304 ip = fPhness->Integral(1, i)/sump; 305 else 306 ip = 0; 307 308 Stat_t ig; 309 if (sumg != 0.0) 310 ig = fGhness->Integral(1, i)/sumg; 311 else 312 ig = 0; 263 313 264 314 fIntPhness->SetBinContent(i, ip); … … 409 459 410 460 c.cd(1); 411 gStyle->SetOptStat(10);461 //gStyle->SetOptStat(10); 412 462 Getghness()->DrawCopy(); 413 463 Getphness()->SetLineColor(kRed); 414 464 Getphness()->DrawCopy("same"); 415 465 c.cd(2); 416 gStyle->SetOptStat(0);466 //gStyle->SetOptStat(0); 417 467 Getighness()->DrawCopy(); 418 468 Getiphness()->SetLineColor(kRed); … … 483 533 484 534 gPad->cd(1); 485 gStyle->SetOptStat(10);535 //gStyle->SetOptStat(10); 486 536 Getghness()->Draw(); 487 537 Getphness()->SetLineColor(kRed); 488 538 Getphness()->Draw("same"); 489 539 gPad->cd(2); 490 gStyle->SetOptStat(0);540 //gStyle->SetOptStat(0); 491 541 Getighness()->Draw(); 492 542 Getiphness()->SetLineColor(kRed); -
trunk/MagicSoft/Mars/mhist/MHMatrix.cc
r1772 r1809 504 504 } 505 505 506 // -------------------------------------------------------------------------- 507 // 506 508 void MHMatrix::Reassign() 507 509 { … … 691 693 // 692 694 // Define the reference matrix 693 // refcolumn number of the column containing the variable, for which a694 // target distribution may be given;695 // refcolumn number of the column (starting at 1)containing the variable, 696 // for which a target distribution may be given; 695 697 // if refcolumn is negative the target distribution will be set 696 698 // equal to the real distribution; the events in the reference … … 748 750 //--------------------------------------------------------- 749 751 // 750 // if refcol < 0 : select reference events randomly 751 // i.e. set the normaliztion factotrs equal to 1.0 752 // if refcolumn < 0 : select reference events randomly 753 // i.e. set the normaliztion factotrs equal to 1.0 754 // refcol is the column number starting at 0; it is >= 0 752 755 753 756 if (refcolumn<0) 754 757 { 755 frefcol = -refcolumn ;758 frefcol = -refcolumn - 1; 756 759 } 757 760 else 758 761 { 759 frefcol = refcolumn;762 frefcol = refcolumn - 1; 760 763 } 761 764 … … 790 793 // << ", " << fM(j-1,frefcol) << endl; 791 794 } 792 793 // if refcolumn<0 set target distribution equal to the real distribution794 // in order to obtain normalization factors = 1.0795 //if (refcolumn<0)796 //{797 // for (Int_t j=1; j<=fnbins; j++)798 // {799 // float cont = fHth-> GetBinContent(j);800 // fHthsh->SetBinContent(j, cont);801 // }802 //}803 795 804 796 //--------------------------------------------------------- … … 816 808 817 809 // if refcolumn < 0 set the correction factors equal to 1 818 if ( refcolumn>=0 .0)810 if ( refcolumn>=0 ) 819 811 b = fHthsh->GetBinContent(j); 820 812 else … … 948 940 { 949 941 *fLog <<ir <<" "; 950 for (Int_t ic=0; ic<13;ic++) cout<<Mnew(ir,ic)<<" ";942 for (Int_t ic=0; ic<Mnew.GetNcols(); ic++) cout<<Mnew(ir,ic)<<" "; 951 943 *fLog <<endl; 952 944 } … … 957 949 { 958 950 float a = fHthaft->GetBinContent(j); 959 if (a>0)*fLog << j << " "<< a << " ";951 *fLog << j << " "<< a << " "; 960 952 } 961 953 *fLog <<endl; … … 967 959 968 960 th1->cd(1); 969 ((TH1F*)fHthsh)->Draw (); // target961 ((TH1F*)fHthsh)->DrawCopy(); // target 970 962 971 963 th1->cd(2); 972 ((TH1F*)fHth)->Draw (); // real histogram before964 ((TH1F*)fHth)->DrawCopy(); // real histogram before 973 965 974 966 th1->cd(3); 975 ((TH1F*)fHthd) -> Draw();// correction factors967 ((TH1F*)fHthd)->DrawCopy(); // correction factors 976 968 977 969 th1->cd(4); 978 ((TH1F*)fHthaft) -> Draw(); // histogram after 979 980 //--------------------------------------------------------- 981 // --- write onto output file 982 // 983 //TFile *outfile = new TFile(fileNameout, "RECREATE", ""); 984 //Mnew.Write(fileNameout); 985 //outfile->Close(); 970 ((TH1F*)fHthaft)->DrawCopy(); // histogram after 971 972 //--------------------------------------------------------- 986 973 987 974 return kTRUE; -
trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc
r1785 r1809 224 224 TObject *MHSigmaTheta::DrawClone(Option_t *opt) 225 225 { 226 TCanvas &c = *MakeDefCanvas("SigmaTheta ", "Sigmabar vs. Theta",226 TCanvas &c = *MakeDefCanvas("SigmaThetaPlot", "Sigmabar vs. Theta", 227 227 900, 900); 228 228 c.Divide(3, 3); … … 236 236 c.cd(1); 237 237 h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E"); 238 h->SetDirectory(NULL); 238 239 h->SetTitle("Distribution of \\Theta"); 239 240 h->SetXTitle("\\Theta [\\circ]"); 240 241 h->SetYTitle("No.of events"); 241 242 242 h->Draw (opt);243 h->DrawCopy(opt); 243 244 h->SetBit(kCanDelete);; 244 245 gPad->SetLogy(); … … 246 247 c.cd(4); 247 248 h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E"); 249 h->SetDirectory(NULL); 248 250 h->SetTitle("Distribution of Sigmabar"); 249 251 h->SetXTitle("Sigmabar"); 250 252 h->SetYTitle("No.of events"); 251 253 252 h->Draw (opt);254 h->DrawCopy(opt); 253 255 h->SetBit(kCanDelete);; 254 256 … … 263 265 c.cd(2); 264 266 l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zx"); 267 l->SetDirectory(NULL); 265 268 l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); 266 269 l->SetXTitle("\\Theta [\\circ]"); 267 270 l->SetYTitle("Sigma^2-Sigmabar^2"); 268 271 269 l->Draw ("box");272 l->DrawCopy("box"); 270 273 l->SetBit(kCanDelete);; 271 274 272 275 c.cd(5); 273 276 l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zy"); 277 l->SetDirectory(NULL); 274 278 l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); 275 279 l->SetXTitle("pixel"); 276 280 l->SetYTitle("Sigma^2-Sigmabar^2"); 277 281 278 l->Draw ("box");282 l->DrawCopy("box"); 279 283 l->SetBit(kCanDelete);; 280 284 … … 290 294 c.cd(3); 291 295 k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zx"); 296 k->SetDirectory(NULL); 292 297 k->SetTitle("Sigma vs. \\Theta (all pixels)"); 293 298 k->SetXTitle("\\Theta [\\circ]"); 294 299 k->SetYTitle("Sigma"); 295 300 296 k->Draw ("box");301 k->DrawCopy("box"); 297 302 k->SetBit(kCanDelete);; 298 303 299 304 c.cd(6); 300 305 k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zy"); 306 k->SetDirectory(NULL); 301 307 k->SetTitle("Sigma vs. pixel number (all \\Theta)"); 302 308 k->SetXTitle("pixel"); 303 309 k->SetYTitle("Sigma"); 304 310 305 k->Draw ("box");311 k->DrawCopy("box"); 306 312 k->SetBit(kCanDelete);; 307 313
Note:
See TracChangeset
for help on using the changeset viewer.