Ignore:
Timestamp:
07/04/03 08:00:18 (21 years ago)
Author:
wittek
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/macros
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/macros/CT1Analysis.C

    r2255 r2262  
    210210
    211211
    212     // Job E_EST_UP :
    213     //  - read MC1.root file
    214     //  - select g/h separation method XX
    215     //  - optimize energy estimation for events passing the final cuts
    216     //  - write parameters of energy estimator onto file
    217     //  - update ON1.root and MC1.root files with estimated energy
    218     //    (ON_XX1.root and MC_XX1.root)
    219 
    220     Bool_t JobE_EST_UP  = kFALSE; 
    221     Bool_t WESTUP       = kFALSE;  // update root files ?
    222 
    223 
    224 
    225     // Job F_XX : 
     212    // Job E_XX : 
    226213    //  - select g/h separation method XX
    227214    //  - read MC_XX2.root file
    … …  
    232219    //  - write root file for ON data after final cuts (ON3.root))
    233220
    234     Bool_t JobF_XX  = kTRUE; 
     221    Bool_t JobE_XX  = kTRUE; 
    235222    Bool_t WXX      = kFALSE;  // write out root file ON3.root ?
    236223
    … …  
    23612348
    23622349
     2350
     2351
    23632352  //---------------------------------------------------------------------
    2364   // Job E_EST_UP
    2365   //============
    2366 
    2367     //  - read MC2.root file
     2353  // Job E_XX
     2354  //=========
     2355
    23682356    //  - select g/h separation method XX
    2369     //  - optimize energy estimator for events passing final cuts
    2370     //  - write parameters of energy estimator onto file "energyest_XX.root"
    2371     //
    2372     //  - read ON2.root and MC2.root files
    2373     //  - update input root file with the estimated energy
    2374     //    (ON_XX2.root, MC_XX2.root)
    2375 
    2376 
    2377  if (JobE_EST_UP)
     2357    //  - read MC_XX2.root file
     2358    //  - calculate eff. collection area
     2359    //  - read ON_XX2.root file
     2360    //  - apply final cuts
     2361    //  - calculate flux
     2362    //  - write root file for ON data after final cuts (ON_XX3.root))
     2363
     2364
     2365 if (JobE_XX)
    23782366 {
    23792367    gLog << "=====================================================" << endl;
    2380     gLog << "Macro CT1Analysis : Start of Job E_EST_UP" << endl;
     2368    gLog << "Macro CT1Analysis : Start of Job E_XX" << endl;
    23812369
    23822370    gLog << "" << endl;
    2383     gLog << "Macro CT1Analysis : JobE_EST_UP, WESTUP = "
    2384          << JobE_EST_UP  << ",  " << WESTUP << endl;
    2385 
    2386 
    2387     TString typeON = "ON";
    2388     TString typeMC = "MC";
    2389     TString ext    = "3.root";
    2390     TString extout = "4.root";
     2371    gLog << "Macro CT1Analysis : JobE_XX, WXX = "
     2372         << JobE_XX  << ",  " << WXX << endl;
     2373
     2374    // type of data to be analysed
     2375    TString typeData = "ON";
     2376    //TString typeData = "OFF";
     2377    //TString typeData = "MC";
     2378    gLog << "typeData = " << typeData << endl;
     2379
     2380    TString typeMC   = "MC";
     2381    TString ext      = "3.root";
     2382    TString extout   = "4.root";
     2383
     2384    //------------------------------
     2385    // selection of g/h separation method
     2386    // and definition of final selections
     2387
     2388    //TString XX("NN");
     2389    //TString XX("SC");
     2390    TString XX("RF");
     2391    TString fhadronnessName("Had");
     2392    fhadronnessName += XX;
     2393    gLog << "fhadronnessName = " << fhadronnessName << endl;
     2394
     2395    // maximum values of the hadronness, |ALPHA| and DIST
     2396    Float_t maxhadronness   = 0.40;
     2397    Float_t maxalpha        = 20.0;
     2398    Float_t maxdist         = 10.0;
     2399    gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
     2400         << maxhadronness << ",  " << maxalpha << ",  "
     2401         << maxdist << endl;
    23912402
    23922403    //------------------------------
    … …  
    23982409
    23992410    //------------------------------
    2400     // selection of g/h separation method
    2401     // and definition of final selections
    2402 
    2403     TString XX("NN");
    2404     //TString XX("SC");
    2405     //TString XX("RF");
    2406     TString fhadronnessName("Had");
    2407     fhadronnessName += XX;
    2408     gLog << "fhadronnessName = " << fhadronnessName << endl;
    2409 
    2410     // maximum values of the hadronness, |alpha| and dist
    2411     Float_t maxhadronness   = 0.40;
    2412     Float_t maxalpha        = 20.0;
    2413     Float_t maxdist         = 10.0;
    2414     gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
    2415          << maxhadronness << ",  " << maxalpha << ",  "
    2416          << maxdist << endl;
    2417 
    24182411    // name of file containing the parameters of the energy estimator
    24192412    TString energyParName(outPath);
    24202413    energyParName += "energyest_";
    2421     energyParName += XX;
     2414    //energyParName += XX;
    24222415    energyParName += ".root";
    24232416    gLog << "energyParName = " << energyParName << endl;
    24242417
    2425 
    24262418    //------------------------------
    2427     // name of ON file to be updated
    2428     TString filenameON(outPath);
    2429     filenameON += typeON;
    2430     filenameON += ext;
    2431     gLog << "filenameON = " << filenameON << endl;
    2432 
    2433     // name of MC file to be updated
    2434     TString filenameMC(outPath);
    2435     filenameMC += typeMC;
    2436     filenameMC += ext;
    2437     gLog << "filenameMC = " << filenameMC << endl;
     2419    // name of MC file to be used for calculating the eff. collection areas
     2420    TString filenameArea(outPath);
     2421    filenameArea += typeMC;
     2422    //filenameArea += "_";
     2423    //filenameArea += XX;
     2424    filenameArea += ext;
     2425    gLog << "filenameArea = " << filenameArea << endl;
    24382426
    24392427    //------------------------------
    2440     // name of updated ON file
    2441     TString filenameONup(outPath);
    2442     filenameONup += typeON;
    2443     filenameONup += "_";
    2444     filenameONup += XX;
    2445     filenameONup += extout;
    2446     gLog << "filenameONup = " << filenameONup << endl;
    2447 
    2448     // name of updated MC file
    2449     TString filenameMCup(outPath);
    2450     filenameMCup += typeMC;
    2451     filenameMCup += "_";
    2452     filenameMCup += XX;
    2453     filenameMCup += extout;
    2454     gLog << "filenameMCup = " << filenameMCup << endl;
    2455 
    2456     //-----------------------------------------------------------
     2428    // name of file containing the eff. collection areas
     2429    TString collareaName(outPath);
     2430    collareaName += "area_";
     2431    //collareaName += XX;
     2432    collareaName += ".root";
     2433    gLog << "collareaName = " << collareaName << endl;
     2434
     2435    //------------------------------
     2436    // name of data file to be analysed
     2437    TString filenameData(outPath);
     2438    filenameData += typeData;
     2439    //filenameData += "_";
     2440    //filenameData += XX;
     2441    filenameData += ext;
     2442    gLog << "filenameData = " << filenameData << endl;
     2443
     2444    //------------------------------
     2445    // name of output data file (after the final cuts)
     2446    TString filenameDataout(outPath);
     2447    filenameDataout += typeData;
     2448    //filenameDataout += "_";
     2449    //filenameDataout += XX;
     2450    filenameDataout += extout;
     2451    gLog << "filenameDataout = " << filenameDataout << endl;
     2452
     2453
     2454    //====================================================================
     2455    gLog << "-----------------------------------------------" << endl;
     2456    gLog << "Start calculation of effective collection areas" << endl;
     2457    MParList  parlist;
     2458    MTaskList tasklist;
     2459
     2460    //---------------------------------------
     2461    // Setup the tasks to be executed
     2462    //
     2463    MReadMarsFile reader("Events", filenameArea);
     2464    reader.DisableAutoScheme();
     2465
     2466    MFCT1SelFinal cuthadrons;
     2467    cuthadrons.SetHadronnessName(fhadronnessName);
     2468    cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
     2469
     2470    MContinue conthadrons(&cuthadrons);
     2471
     2472    //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
     2473    //MHMcCT1CollectionArea* collarea;
     2474
     2475    MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
     2476    filler.SetName("CollectionArea");
     2477
     2478    //********************************
     2479    // entries in MParList
     2480
     2481    parlist.AddToList(&tasklist);
     2482    InitBinnings(&parlist);
     2483    //parlist.AddToList(collarea);
     2484
     2485    //********************************
     2486    // entries in MTaskList
     2487
     2488    tasklist.AddToList(&reader);   
     2489    tasklist.AddToList(&conthadrons);
     2490    tasklist.AddToList(&filler);
     2491
     2492    //********************************
     2493
     2494    //-----------------------------------------
     2495    // Execute event loop
     2496    //
     2497    MEvtLoop magic;
     2498    magic.SetParList(&parlist);
     2499
     2500    MProgressBar bar;
     2501    magic.SetProgressBar(&bar);
     2502    if (!magic.Eventloop())
     2503        return;
     2504
     2505    tasklist.PrintStatistics(0, kTRUE);
     2506
     2507    // Calculate effective collection areas
     2508    // and display the histograms
     2509    //
     2510    MHMcCT1CollectionArea *collarea =
     2511         (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
     2512    collarea->CalcEfficiency();
     2513    collarea->DrawClone("lego");
     2514
     2515    // save binnings for call to CT1EEst
     2516    MBinning *binsE     = (MBinning*)parlist.FindObject("BinningE");
     2517    if (!binsE)
     2518        {
     2519          gLog << "Object 'BinningE' not found in MParList" << endl;
     2520          return;
     2521        }
     2522    MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
     2523    if (!binsTheta)
     2524        {
     2525          gLog << "Object 'BinningTheta' not found in MParList" << endl;
     2526          return;
     2527        }
     2528
     2529
     2530    //---------------------------------------------
     2531    // Write histograms to a file
     2532    //
     2533
     2534    TFile f(collareaName, "RECREATE");
     2535    collarea->GetHist()->Write();
     2536    collarea->GetHAll()->Write();
     2537    collarea->GetHSel()->Write();
     2538    f.Close();
     2539
     2540
     2541    gLog << "Calculation of effective collection areas done" << endl;
     2542    gLog << "-----------------------------------------------" << endl;   
     2543    //------------------------------------------------------------------
     2544
    24572545
    24582546    TString fHilName    = "MHillas";
    … …  
    24602548    TString fHilNameSrc = "MHillasSrc";
    24612549    TString fImgParName = "MNewImagePar";
     2550
    24622551
    24632552    //===========================================================
    … …  
    24712560    CT1EEst(inpath,   filenameOpt,   outpath, energyParName,
    24722561            fHilName, fHilNameSrc,   fhadronnessName,
    2473             howMany,  maxhadronness, maxalpha, maxdist);
     2562            howMany,  maxhadronness, maxalpha, maxdist,
     2563            binsE, binsTheta);
    24742564
    24752565    //-----------------------------------------------------------
    24762566    //
    2477     // Read in parameters of energy estimator
     2567    // Read in parameters of energy estimator ("MMcEnergyEst")
     2568    //                   and migration matrix ("MHMcEnergyMigration")
    24782569    //
    24792570    gLog << "================================================================"
    24802571         << endl;
    2481     gLog << "Macro CT1Analysis.C : read parameters of energy estimator from file '"
     2572    gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
    24822573         << energyParName << "'" << endl;
    24832574    TFile enparam(energyParName);
    2484     MMcEnergyEst mcest("MMcEnergyEst");
    2485     mcest.Read("MMcEnergyEst");
    2486     enparam.Close();
     2575    enparam.ls();
     2576
     2577    //MMcEnergyEst mcest("MMcEnergyEst");
     2578    //mcest.Read("MMcEnergyEst");
     2579    MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
     2580
     2581    gLog << "Parameters of energy estimator were read in" << endl;
    24872582
    24882583    TArrayD parA(5);
    … …  
    24932588      parB[i] = mcest.GetCoeff( i+parA.GetSize() );
    24942589
    2495 
    2496    if (WESTUP)
    2497    {
    2498     //==========   start update   ============================================
    2499     //
    2500     // Update ON and MC root files with the estimated energy
    2501 
    2502     //---------------------------------------------------
    2503     // Update ON data
     2590    gLog << "Read in Migration matrix" << endl;   
     2591
     2592    //MHMcEnergyMigration mighiston("MHMcEnergyMigration");
     2593    //mighiston.Read("MHMcEnergyMigration");
     2594    MHMcEnergyMigration &mighiston =
     2595          *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
     2596
     2597    gLog << "Migration matrix was read in" << endl;
     2598
     2599    //*************************************************************************
     2600    //
     2601    // Analyse the data
    25042602    //
    25052603    gLog << "============================================================"
    25062604         << endl;
    2507     gLog << "Macro CT1Analysis.C : update file '" << filenameON
    2508          << "'" << endl;
     2605    gLog << "Analyse the data" << endl;
    25092606
    25102607    MTaskList tliston;
    25112608    MParList pliston;
    2512 
    25132609
    25142610    // geometry is needed in  MHHillas... classes
    … …  
    25162612             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    25172613
     2614
    25182615    //-------------------------------------------
    25192616    // create the tasks which should be executed
    25202617    //
    25212618
    2522     MReadMarsFile read("Events", filenameON);
     2619    MReadMarsFile read("Events", filenameData);
    25232620    read.DisableAutoScheme();
    25242621
    2525     //---------------------------
    2526     // calculate estimated energy
    2527 
    2528     MEnergyEstParam eest2(fHilName);
    2529     eest2.Add(fHilNameSrc);
    2530 
    2531     eest2.SetCoeffA(parA);
    2532     eest2.SetCoeffB(parB);
    2533 
    25342622    //.......................................................................
    25352623
    2536       MWriteRootFile write(filenameONup);
     2624
     2625      MWriteRootFile write(filenameDataout);
    25372626
    25382627      write.AddContainer("MRawRunHeader", "RunHeaders");
    … …  
    25532642      write.AddContainer("MEnergyEst",    "Events");
    25542643
    2555     //-----------------------------------------------------------------
    2556 
    2557     MFCT1SelFinal selfinal(fHilNameSrc);
    2558     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2559     selfinal.SetHadronnessName(fhadronnessName);
    2560     MContinue contfinal(&selfinal);
    2561 
    2562 
    2563     //*****************************
    2564     // entries in MParList
    2565 
    2566     pliston.AddToList(&tliston);
    2567     InitBinnings(&pliston);
    2568 
    2569 
    2570     //*****************************
    2571     // entries in MTaskList
    2572    
    2573     tliston.AddToList(&read);
    2574     tliston.AddToList(&eest2);
    2575     tliston.AddToList(&write);
    2576     tliston.AddToList(&contfinal);
    2577 
    2578     //*****************************
    2579 
    2580     //-------------------------------------------
    2581     // Execute event loop
    2582     //
    2583     MProgressBar bar;
    2584     MEvtLoop evtloop;
    2585     evtloop.SetParList(&pliston);
    2586     evtloop.SetProgressBar(&bar);
    2587 
    2588     Int_t maxevents = -1;
    2589     //Int_t maxevents = 1000;
    2590     if ( !evtloop.Eventloop(maxevents) )
    2591         return;
    2592 
    2593     tliston.PrintStatistics(0, kTRUE);
    2594     DeleteBinnings(&pliston);
    2595 
    2596     //---------------------------------------------------
    2597     //---------------------------------------------------
    2598     // Update MC data
    2599     //
    2600     gLog << "============================================================"
    2601          << endl;
    2602     gLog << "Macro CT1Analysis.C : update file '" << filenameMC
    2603          << "'" << endl;
    2604 
    2605     MTaskList tlistmc;
    2606     MParList plistmc;
    2607 
    2608     //-------------------------------------------
    2609     // create the tasks which should be executed
    2610     //
    2611 
    2612     MReadMarsFile read("Events", filenameMC);
    2613     read.DisableAutoScheme();
    2614 
    2615     //---------------------------
    2616     // calculate estimated energy
    2617 
    2618     MEnergyEstParam eest2(fHilName);
    2619     eest2.Add(fHilNameSrc);
    2620 
    2621     eest2.SetCoeffA(parA);
    2622     eest2.SetCoeffB(parB);
    2623 
    2624     //.......................................................................
    2625 
    2626       MWriteRootFile write(filenameMCup);
    2627 
    2628       write.AddContainer("MRawRunHeader", "RunHeaders");
    2629       write.AddContainer("MTime",         "Events");
    2630       write.AddContainer("MMcEvt",        "Events");
    2631       write.AddContainer("ThetaOrig",     "Events");
    2632       write.AddContainer("MSrcPosCam",    "Events");
    2633       write.AddContainer("MSigmabar",     "Events");
    2634       write.AddContainer("MHillas",       "Events");
    2635       write.AddContainer("MHillasExt",    "Events");
    2636       write.AddContainer("MHillasSrc",    "Events");
    2637       write.AddContainer("MNewImagePar",  "Events");
    2638 
    2639       write.AddContainer("HadNN",         "Events");
    2640       write.AddContainer("HadSC",         "Events");
    2641       write.AddContainer("HadRF",         "Events");
    2642 
    2643       write.AddContainer("MEnergyEst",    "Events");
    2644 
    2645     //-----------------------------------------------------------------
    2646 
    2647     MFCT1SelFinal selfinal(fHilNameSrc);
    2648     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    2649     selfinal.SetHadronnessName(fhadronnessName);
    2650     MContinue contfinal(&selfinal);
    2651 
    2652 
    2653     //*****************************
    2654     // entries in MParList
    2655 
    2656     plistmc.AddToList(&tlistmc);
    2657     InitBinnings(&plistmc);
    2658 
    2659 
    2660     //*****************************
    2661     // entries in MTaskList
    2662    
    2663     tlistmc.AddToList(&read);
    2664     tlistmc.AddToList(&eest2);
    2665     tlistmc.AddToList(&write);
    2666     tlistmc.AddToList(&contfinal);
    2667 
    2668     //*****************************
    2669 
    2670     //-------------------------------------------
    2671     // Execute event loop
    2672     //
    2673     MProgressBar bar;
    2674     MEvtLoop evtloop;
    2675     evtloop.SetParList(&plistmc);
    2676     evtloop.SetProgressBar(&bar);
    2677 
    2678     Int_t maxevents = -1;
    2679     //Int_t maxevents = 1000;
    2680     if ( !evtloop.Eventloop(maxevents) )
    2681         return;
    2682 
    2683     tlistmc.PrintStatistics(0, kTRUE);
    2684     DeleteBinnings(&plistmc);
    2685 
    2686 
    2687     //==========   end update   ============================================
    2688    }
    2689    
    2690     enparam.Close();
    2691 
    2692     gLog << "Macro CT1Analysis : End of Job E_EST_UP" << endl;
    2693     gLog << "=======================================================" << endl;
    2694  }
    2695   //---------------------------------------------------------------------
    2696 
    2697 
    2698   //---------------------------------------------------------------------
    2699   // Job F_XX
    2700   //=========
    2701 
    2702     //  - select g/h separation method XX
    2703     //  - read MC_XX2.root file
    2704     //  - calculate eff. collection area
    2705     //  - read ON_XX2.root file
    2706     //  - apply final cuts
    2707     //  - calculate flux
    2708     //  - write root file for ON data after final cuts (ON_XX3.root))
    2709 
    2710 
    2711  if (JobF_XX)
    2712  {
    2713     gLog << "=====================================================" << endl;
    2714     gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;
    2715 
    2716     gLog << "" << endl;
    2717     gLog << "Macro CT1Analysis : JobF_XX, WXX = "
    2718          << JobF_XX  << ",  " << WXX << endl;
    2719 
    2720     // type of data to be analysed
    2721     TString typeData = "ON";
    2722     //TString typeData = "OFF";
    2723     //TString typeData = "MC";
    2724     gLog << "typeData = " << typeData << endl;
    2725 
    2726     TString typeMC   = "MC";
    2727     TString ext      = "3.root";
    2728     TString extout   = "4.root";
    2729 
    2730     //------------------------------
    2731     // selection of g/h separation method
    2732     // and definition of final selections
    2733 
    2734     //TString XX("NN");
    2735     //TString XX("SC");
    2736     TString XX("RF");
    2737     TString fhadronnessName("Had");
    2738     fhadronnessName += XX;
    2739     gLog << "fhadronnessName = " << fhadronnessName << endl;
    2740 
    2741     // maximum values of the hadronness, |ALPHA| and DIST
    2742     Float_t maxhadronness   = 0.40;
    2743     Float_t maxalpha        = 20.0;
    2744     Float_t maxdist         = 10.0;
    2745     gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
    2746          << maxhadronness << ",  " << maxalpha << ",  "
    2747          << maxdist << endl;
    2748 
    2749     //------------------------------
    2750     // name of MC file to be used for optimizing the energy estimator
    2751     TString filenameOpt(outPath);
    2752     filenameOpt += typeMC;
    2753     filenameOpt += ext;
    2754     gLog << "filenameOpt = " << filenameOpt << endl;
    2755 
    2756     //------------------------------
    2757     // name of file containing the parameters of the energy estimator
    2758     TString energyParName(outPath);
    2759     energyParName += "energyest_";
    2760     //energyParName += XX;
    2761     energyParName += ".root";
    2762     gLog << "energyParName = " << energyParName << endl;
    2763 
    2764     //------------------------------
    2765     // name of MC file to be used for calculating the eff. collection areas
    2766     TString filenameArea(outPath);
    2767     filenameArea += typeMC;
    2768     //filenameArea += "_";
    2769     //filenameArea += XX;
    2770     filenameArea += ext;
    2771     gLog << "filenameArea = " << filenameArea << endl;
    2772 
    2773     //------------------------------
    2774     // name of file containing the eff. collection areas
    2775     TString collareaName(outPath);
    2776     collareaName += "area_";
    2777     //collareaName += XX;
    2778     collareaName += ".root";
    2779     gLog << "collareaName = " << collareaName << endl;
    2780 
    2781     //------------------------------
    2782     // name of data file to be analysed
    2783     TString filenameData(outPath);
    2784     filenameData += typeData;
    2785     //filenameData += "_";
    2786     //filenameData += XX;
    2787     filenameData += ext;
    2788     gLog << "filenameData = " << filenameData << endl;
    2789 
    2790     //------------------------------
    2791     // name of output data file (after the final cuts)
    2792     TString filenameDataout(outPath);
    2793     filenameDataout += typeData;
    2794     //filenameDataout += "_";
    2795     //filenameDataout += XX;
    2796     filenameDataout += extout;
    2797     gLog << "filenameDataout = " << filenameDataout << endl;
    2798 
    2799 
    2800     //====================================================================
    2801     gLog << "-----------------------------------------------" << endl;
    2802     gLog << "Start calculation of effective collection areas" << endl;
    2803     MParList  parlist;
    2804     MTaskList tasklist;
    2805 
    2806     //---------------------------------------
    2807     // Setup the tasks to be executed
    2808     //
    2809     MReadMarsFile reader("Events", filenameArea);
    2810     reader.DisableAutoScheme();
    2811 
    2812     MFCT1SelFinal cuthadrons;
    2813     cuthadrons.SetHadronnessName(fhadronnessName);
    2814     cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
    2815 
    2816     MContinue conthadrons(&cuthadrons);
    2817 
    2818     //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
    2819     //MHMcCT1CollectionArea* collarea;
    2820 
    2821     MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
    2822     filler.SetName("CollectionArea");
    2823 
    2824     //********************************
    2825     // entries in MParList
    2826 
    2827     parlist.AddToList(&tasklist);
    2828     InitBinnings(&parlist);
    2829     //parlist.AddToList(collarea);
    2830 
    2831     //********************************
    2832     // entries in MTaskList
    2833 
    2834     tasklist.AddToList(&reader);   
    2835     tasklist.AddToList(&conthadrons);
    2836     tasklist.AddToList(&filler);
    2837 
    2838     //********************************
    2839 
    2840     //-----------------------------------------
    2841     // Execute event loop
    2842     //
    2843     MEvtLoop magic;
    2844     magic.SetParList(&parlist);
    2845 
    2846     MProgressBar bar;
    2847     magic.SetProgressBar(&bar);
    2848     if (!magic.Eventloop())
    2849         return;
    2850 
    2851     tasklist.PrintStatistics(0, kTRUE);
    2852 
    2853     // Calculate effective collection areas
    2854     // and display the histograms
    2855     //
    2856     MHMcCT1CollectionArea *collarea =
    2857          (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
    2858     collarea->CalcEfficiency();
    2859     collarea->DrawClone("lego");
    2860 
    2861     //---------------------------------------------
    2862     // Write histograms to a file
    2863     //
    2864 
    2865     TFile f(collareaName, "RECREATE");
    2866     collarea->GetHist()->Write();
    2867     collarea->GetHAll()->Write();
    2868     collarea->GetHSel()->Write();
    2869     f.Close();
    2870 
    2871 
    2872     gLog << "Calculation of effective collection areas done" << endl;
    2873     gLog << "-----------------------------------------------" << endl;   
    2874     //------------------------------------------------------------------
    2875 
    2876 
    2877     TString fHilName    = "MHillas";
    2878     TString fHilNameExt = "MHillasExt";
    2879     TString fHilNameSrc = "MHillasSrc";
    2880     TString fImgParName = "MNewImagePar";
    2881 
    2882 
    2883     //===========================================================
    2884     //
    2885     // Optimization of energy estimator
    2886     //
    2887 
    2888     TString inpath("");
    2889     TString outpath("");
    2890     Int_t howMany = 2000;
    2891     CT1EEst(inpath,   filenameOpt,   outpath, energyParName,
    2892             fHilName, fHilNameSrc,   fhadronnessName,
    2893             howMany,  maxhadronness, maxalpha, maxdist);
    2894 
    2895     //-----------------------------------------------------------
    2896     //
    2897     // Read in parameters of energy estimator ("MMcEnergyEst")
    2898     //                   and migration matrix ("MHMcEnergyMigration")
    2899     //
    2900     gLog << "================================================================"
    2901          << endl;
    2902     gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
    2903          << energyParName << "'" << endl;
    2904     TFile enparam(energyParName);
    2905     enparam.ls();
    2906 
    2907     //MMcEnergyEst mcest("MMcEnergyEst");
    2908     //mcest.Read("MMcEnergyEst");
    2909     MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
    2910 
    2911     gLog << "Parameters of energy estimator were read in" << endl;
    2912 
    2913     TArrayD parA(5);
    2914     TArrayD parB(7);
    2915     for (Int_t i=0; i<parA.GetSize(); i++)
    2916       parA[i] = mcest.GetCoeff(i);
    2917     for (Int_t i=0; i<parB.GetSize(); i++)
    2918       parB[i] = mcest.GetCoeff( i+parA.GetSize() );
    2919 
    2920     gLog << "Read in Migration matrix" << endl;   
    2921 
    2922     //MHMcEnergyMigration mighiston("MHMcEnergyMigration");
    2923     //mighiston.Read("MHMcEnergyMigration");
    2924     MHMcEnergyMigration &mighiston =
    2925           *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
    2926 
    2927     gLog << "Migration matrix was read in" << endl;
    2928 
    2929     //*************************************************************************
    2930     //
    2931     // Analyse the data
    2932     //
    2933     gLog << "============================================================"
    2934          << endl;
    2935     gLog << "Analyse the data" << endl;
    2936 
    2937     MTaskList tliston;
    2938     MParList pliston;
    2939 
    2940     // geometry is needed in  MHHillas... classes
    2941     MGeomCam *fGeom =
    2942              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    2943 
    2944 
    2945     //-------------------------------------------
    2946     // create the tasks which should be executed
    2947     //
    2948 
    2949     MReadMarsFile read("Events", filenameData);
    2950     read.DisableAutoScheme();
    2951 
    2952     //.......................................................................
    2953 
    2954 
    2955       MWriteRootFile write(filenameDataout);
    2956 
    2957       write.AddContainer("MRawRunHeader", "RunHeaders");
    2958       write.AddContainer("MTime",         "Events");
    2959       write.AddContainer("MMcEvt",        "Events");
    2960       write.AddContainer("ThetaOrig",     "Events");
    2961       write.AddContainer("MSrcPosCam",    "Events");
    2962       write.AddContainer("MSigmabar",     "Events");
    2963       write.AddContainer("MHillas",       "Events");
    2964       write.AddContainer("MHillasExt",    "Events");
    2965       write.AddContainer("MHillasSrc",    "Events");
    2966       write.AddContainer("MNewImagePar",  "Events");
    2967 
    2968       write.AddContainer("HadNN",         "Events");
    2969       write.AddContainer("HadSC",         "Events");
    2970       write.AddContainer("HadRF",         "Events");
    2971 
    2972       write.AddContainer("MEnergyEst",    "Events");
    2973 
    29742644
    29752645    //-----------------------------------------------------------------
    … …  
    30912761    DeleteBinnings(&pliston);
    30922762
    3093     gLog << "Macro CT1Analysis : End of Job F_XX" << endl;
     2763    gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
    30942764    gLog << "=======================================================" << endl;
    30952765 }

  • trunk/MagicSoft/Mars/macros/CT1EgyEst.C

    r2243 r2262  
    299299  out2.Close();
    300300
     301  TCanvas *c = new TCanvas;
     302  c->cd();
    301303  mighist.Draw();
    302304
    303   cout << "Quality histograms were added onto the file '" << paramout << endl;
     305  cout << "Quality histograms and migration matrix were added onto the file '" << paramout << endl;
    304306  cout << endl;
    305307  cout << "End of energy estimation part" << endl;
Note: See TracChangeset for help on using the changeset viewer.