Ignore:
Timestamp:
09/29/03 07:09:57 (21 years ago)
Author:
wittek
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2357 r2364  
    252252    //  - update the input files with the hadroness (ON1.root or MC1.root)
    253253
    254     Bool_t JobB_SC_UP  = kFALSE;
     254    Bool_t JobB_SC_UP  = kTRUE;
    255255    Bool_t RMatrix     = kFALSE;  // read training and test matrices from file 
    256     Bool_t WParSC      = kFALSE;  // do optimization and write supercuts
    257                                  // parameter values onto a file
    258     Bool_t WSC         = kFALSE; // update input root file ?
     256    Bool_t WOptimize   = kFALSE;  // do optimization using the training sample
     257                                  // and write supercuts parameter values
     258                                  // onto the file parSCfile
     259    Bool_t RTest       = kFALSE;  // test the supercuts using the test matrix
     260    Bool_t WSC         = kTRUE;  // update input root file ?
    259261
    260262
    … …  
    281283
    282284
    283     // Job E_XX : 
    284     //  - select g/h separation method XX
    285     //  - read MC root file
    286     //  - calculate eff. collection area
    287     //  - optimize energy estimator
    288     //  - read ON root file
    289     //  - apply final cuts
    290     //  - write root file for ON data after final cuts
    291 
    292     Bool_t JobE_XX  = kFALSE; 
    293     Bool_t WXX      = kFALSE;  // write out root file ?
    294 
    295 
    296     // Job F_XX : extended version of E_XX (including flux plots) 
     285
     286    // Job E_XX : extended version of E_XX (including flux plots) 
    297287    //  - select g/h separation method XX
    298288    //  - read MC root file
    … …  
    304294    //  - write root file for ON data after final cuts
    305295
    306     gLog << "before setting switches for JobF_XX" << endl;
    307 
    308     Bool_t JobF_XX  = kTRUE; 
    309     Bool_t WFX      = kTRUE;  // write out root file  ?
    310     Bool_t WRobert  = kTRUE;  // write out Robert's file  ?
     296
     297    Bool_t JobE_XX  = kFALSE; 
     298    Bool_t WEX      = kFALSE;  // write out root file  ?
     299    Bool_t WRobert  = kFALSE;  // write out Robert's file  ?
    311300
    312301
    … …  
    15501539
    15511540    gLog << "" << endl;
    1552     gLog << "Macro CT1Analysis : JobB_SC_UP, WParSC, WSC = "
     1541    gLog << "Macro CT1Analysis : JobB_SC_UP, RMatrix, WOptimize, RTest, WSC = "
    15531542         << (JobB_SC_UP ? "kTRUE" : "kFALSE")  << ",  "
    1554          << (WParSC     ? "kTRUE" : "kFALSE")  << ",  "
     1543         << (RMatrix    ? "kTRUE" : "kFALSE")  << ",  "
     1544         << (WOptimize  ? "kTRUE" : "kFALSE")  << ",  "
     1545         << (RTest      ? "kTRUE" : "kFALSE")  << ",  "
    15551546         << (WSC        ? "kTRUE" : "kFALSE")  << endl;
    1556 
    15571547
    15581548
    … …  
    15801570    // file to be updated (either ON or MC)
    15811571
    1582     //TString typeInput = "ON";
    1583     TString typeInput = "MC";
     1572    TString typeInput = "ON";
     1573    //TString typeInput = "MC";
    15841574    gLog << "typeInput = " << typeInput << endl;
    15851575
    … …  
    16521642    //--------------------------
    16531643    // create matrices and write them onto files
    1654     if (!RMatrix)
     1644    if (!RMatrix  &&  (WOptimize || RTest)  )
    16551645    {
    16561646      MH3 &mh3 = *(new MH3("MHillas.fSize"));
    … …  
    17051695
    17061696
    1707 if (WParSC)
     1697if (WOptimize)
    17081698  {
     1699    gLog << "CT1Analysis.C : optimize the supercuts using the training matrix"
     1700         << endl;
     1701
    17091702    TArrayD params(0);
    17101703    TArrayD steps(0);
    … …  
    18061799       return;
    18071800    }
     1801  }
    18081802
    18091803    //--------------------------------------
    18101804    // test the supercuts on the test sample
    18111805    //   
     1806
     1807 if (RTest)
     1808 {
     1809    gLog << "CT1Analysis.C : test the supercuts on the test matrix" << endl;
    18121810    Bool_t rt = findsuper.TestParams();
    1813   }
     1811    if (!rt)
     1812    {
     1813       gLog << "CT1Analysis.C : test of supercuts on the test matrix failed"
     1814            << endl;
     1815    }
     1816
     1817 }
    18141818
    18151819
    … …  
    19061910      //MWriteRootFile write(outNameImage, "UPDATE");
    19071911      //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
     1912
     1913    /*
    19081914      MWriteRootFile write(outNameImage, "RECREATE");
    19091915
    … …  
    19211927      write.AddContainer("HadRF",         "Events");
    19221928      write.AddContainer(hadSCName,       "Events");
    1923 
     1929    */
    19241930
    19251931    //-----------------------------------------------------------------
    … …  
    19952001    tliston.AddToList(&fillhadsc);
    19962002
    1997     tliston.AddToList(&write);
     2003    //tliston.AddToList(&write);
    19982004    tliston.AddToList(&contfinalgh);
    19992005
    … …  
    24622468
    24632469
     2470
    24642471  //---------------------------------------------------------------------
    24652472  // Job E_XX
    … …  
    24812488
    24822489    gLog << "" << endl;
    2483     gLog << "Macro CT1Analysis : JobE_XX, WXX = "
     2490    gLog << "Macro CT1Analysis : JobE_XX, WEX = "
    24842491         << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  "
    2485          << (WXX     ? "kTRUE" : "kFALSE")  << endl;
     2492         << (WEX     ? "kTRUE" : "kFALSE")  << endl;
    24862493
    24872494
    … …  
    25612568    gLog << "filenameDataout = " << filenameDataout << endl;
    25622569
     2570    //------------------------------
     2571    // name of file containing histograms for flux calculastion
     2572    TString filenameResults(outPath);
     2573    filenameResults += typeData;
     2574    filenameResults += "Results_";
     2575    filenameResults += XX;
     2576    filenameResults += extout;
     2577    gLog << "filenameResults = " << filenameResults << endl;
     2578
    25632579
    25642580    //====================================================================
    … …  
    25802596    MContinue conthadrons(&cuthadrons);
    25812597
    2582     //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
    2583     //MHMcCT1CollectionArea* collarea;
     2598    MHMcCT1CollectionArea collarea;
     2599    collarea.SetEaxis(MHMcCT1CollectionArea::kLinear);
    25842600
    25852601    MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
    … …  
    25912607    parlist.AddToList(&tasklist);
    25922608    InitBinnings(&parlist);
    2593     //parlist.AddToList(collarea);
     2609    parlist.AddToList(&collarea);
    25942610
    25952611    //********************************
    … …  
    26182634    // and display the histograms
    26192635    //
    2620     MHMcCT1CollectionArea *collarea =
    2621          (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
    2622     collarea->CalcEfficiency();
    2623     collarea->DrawClone("lego");
     2636    //MHMcCT1CollectionArea *collarea =
     2637    //     (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
     2638    collarea.CalcEfficiency();
     2639    collarea.DrawClone();
    26242640
    26252641    // save binnings for call to CT1EEst
    … …  
    26432659
    26442660    TFile f(collareaName, "RECREATE");
    2645     collarea->GetHist()->Write();
    2646     collarea->GetHAll()->Write();
    2647     collarea->GetHSel()->Write();
     2661    collarea.GetHist()->Write();
     2662    collarea.GetHAll()->Write();
     2663    collarea.GetHSel()->Write();
    26482664    f.Close();
    26492665
    … …  
    26602676    TString fImgParName = "MNewImagePar";
    26612677
    2662 
     2678 
    26632679    //===========================================================
    26642680    //
    … …  
    26772693
    26782694
    2679   if (WXX)
     2695  if (WEX)
    26802696  {
    26812697    //-----------------------------------------------------------
    … …  
    26972713    gLog << "Parameters of energy estimator were read in" << endl;
    26982714
    2699     TArrayD parA(5);
    2700     TArrayD parB(7);
     2715    TArrayD parA(mcest.GetNumCoeffA());
     2716    TArrayD parB(mcest.GetNumCoeffB());
    27012717    for (Int_t i=0; i<parA.GetSize(); i++)
    27022718      parA[i] = mcest.GetCoeff(i);
    … …  
    27382754    //.......................................................................
    27392755
    2740 
     2756      gLog << "CT1Analysis.C : write root file '" << filenameDataout
     2757           << "'" << endl;
     2758   
    27412759      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
     2760
    27422761      MWriteRootFile write(filenameDataout, "RECREATE");
    27432762
    … …  
    27532772      write.AddContainer("MNewImagePar",  "Events");
    27542773
    2755       write.AddContainer("HadNN",         "Events");
     2774      //write.AddContainer("HadNN",         "Events");
    27562775      write.AddContainer("HadSC",         "Events");
    27572776      write.AddContainer("HadRF",         "Events");
    … …  
    27722791    contfinalgh.SetName("ContSelFinalgh");
    27732792
    2774     MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
    2775     fillhadnn.SetName("HhadNN");
     2793    //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
     2794    //fillhadnn.SetName("HhadNN");
    27762795    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
    27772796    fillhadsc.SetName("HhadSC");
    … …  
    27912810    // calculate estimated energy using Daniel's parameters
    27922811
    2793     //MEnergyEstParam eeston(fHilName);
     2812    //MEnergyEstParamDanielMkn421 eeston(fHilName);
    27942813    //eeston.Add(fHilNameSrc);
     2814    //eeston.SetCoeffA(parA);
     2815    //eeston.SetCoeffB(parB);
     2816
    27952817
    27962818    //---------------------------
    … …  
    28112833    MFillH hfill5("MHNewImagePar", fImgParName);
    28122834    hfill5.SetName("HNewImagePar");   
     2835
     2836    //---------------------------
     2837    // new from Robert
     2838
     2839    MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
     2840    hfill6.SetName("HTimeDiffTheta");
     2841
     2842    MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
     2843    hfill6a.SetName("HTimeDiffTime");
     2844
     2845    MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
     2846    hfill7.SetName("HAlphaEnergyTheta");
     2847
     2848    MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
     2849    hfill7a.SetName("HAlphaEnergyTime");
     2850
     2851    MFillH hfill7b("MHThetabarTime", fHilNameSrc);
     2852    hfill7b.SetName("HThetabarTime");
     2853
     2854    MFillH hfill7c("MHEnergyTime", "MMcEvt");
     2855    hfill7c.SetName("HEnergyTime");
     2856
     2857
     2858    //---------------------------
    28132859
    28142860    MFCT1SelFinal selfinal(fHilNameSrc);
    … …  
    28322878    tliston.AddToList(&read);
    28332879
    2834     tliston.AddToList(&contfinalgh);
    2835 
    2836     if (WXX)
    2837       tliston.AddToList(&write);
    2838 
    2839     tliston.AddToList(&eeston);
    2840 
    2841     tliston.AddToList(&fillhadnn);
    2842     tliston.AddToList(&fillhadsc);
    2843     tliston.AddToList(&fillhadrf);
    2844 
    2845     tliston.AddToList(&hfill1);
    2846     tliston.AddToList(&hfill2);
    2847     tliston.AddToList(&hfill3);
    2848     tliston.AddToList(&hfill4);
    2849     tliston.AddToList(&hfill5);
    2850 
    2851     tliston.AddToList(&contfinal);
    2852 
    2853     //*****************************
    2854 
    2855     //-------------------------------------------
    2856     // Execute event loop
    2857     //
    2858     MProgressBar bar;
    2859     MEvtLoop evtloop;
    2860     evtloop.SetParList(&pliston);
    2861     evtloop.SetProgressBar(&bar);
    2862 
    2863     //Int_t maxevents = -1;
    2864     Int_t maxevents = 10000;
    2865     if ( !evtloop.Eventloop(maxevents) )
    2866         return;
    2867 
    2868     tliston.PrintStatistics(0, kTRUE);
    2869 
    2870 
    2871     //-------------------------------------------
    2872     // Display the histograms
    2873     //
    2874 
    2875     pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
    2876     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    2877     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    2878 
    2879     pliston.FindObject("MHHillas")->DrawClone();
    2880     pliston.FindObject("MHHillasExt")->DrawClone();
    2881     pliston.FindObject("MHHillasSrc")->DrawClone();
    2882     pliston.FindObject("MHNewImagePar")->DrawClone();
    2883     pliston.FindObject("MHStarMap")->DrawClone();
    2884 
    2885     DeleteBinnings(&pliston);
    2886 
    2887   }
    2888 
    2889     gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
    2890     gLog << "=======================================================" << endl;
    2891  }
    2892   //---------------------------------------------------------------------
    2893 
    2894 
    2895   //---------------------------------------------------------------------
    2896   // Job F_XX
    2897   //=========
    2898 
    2899     //  - select g/h separation method XX
    2900     //  - read MC_XX2.root file
    2901     //  - calculate eff. collection area
    2902     //  - read ON_XX2.root file
    2903     //  - apply final cuts
    2904     //  - calculate flux
    2905     //  - write root file for ON data after final cuts (ON_XX3.root))
    2906 
    2907 
    2908  if (JobF_XX)
    2909  {
    2910     gLog << "=====================================================" << endl;
    2911     gLog << "Macro CT1Analysis : Start of Job F_XX" << endl;
    2912 
    2913     gLog << "" << endl;
    2914     gLog << "Macro CT1Analysis : JobF_XX, WFX = "
    2915          << (JobF_XX ? "kTRUE" : "kFALSE")  << ",  "
    2916          << (WFX     ? "kTRUE" : "kFALSE")  << endl;
    2917 
    2918 
    2919     // type of data to be analysed
    2920     TString typeData = "ON";
    2921     //TString typeData = "OFF";
    2922     //TString typeData = "MC";
    2923     gLog << "typeData = " << typeData << endl;
    2924 
    2925     TString typeMC   = "MC";
    2926     TString ext      = "3.root";
    2927     TString extout   = "4.root";
    2928 
    2929     //------------------------------
    2930     // selection of g/h separation method
    2931     // and definition of final selections
    2932 
    2933     //TString XX("NN");
    2934     TString XX("SC");
    2935     //TString XX("RF");
    2936     TString fhadronnessName("Had");
    2937     fhadronnessName += XX;
    2938     gLog << "fhadronnessName = " << fhadronnessName << endl;
    2939 
    2940     // maximum values of the hadronness, |ALPHA| and DIST
    2941     Float_t maxhadronness   = 0.40;
    2942     Float_t maxalpha        = 20.0;
    2943     Float_t maxdist         = 10.0;
    2944     gLog << "Maximum values of hadronness, |ALPHA| and DIST = "
    2945          << maxhadronness << ",  " << maxalpha << ",  "
    2946          << maxdist << endl;
    2947 
    2948     //------------------------------
    2949     // name of MC file to be used for optimizing the energy estimator
    2950     TString filenameOpt(outPath);
    2951     filenameOpt += typeMC;
    2952     filenameOpt += ext;
    2953     gLog << "filenameOpt = " << filenameOpt << endl;
    2954 
    2955     //------------------------------
    2956     // name of file containing the parameters of the energy estimator
    2957     TString energyParName(outPath);
    2958     energyParName += "energyest_";
    2959     energyParName += XX;
    2960     energyParName += ".root";
    2961     gLog << "energyParName = " << energyParName << endl;
    2962 
    2963     //------------------------------
    2964     // name of MC file to be used for calculating the eff. collection areas
    2965     TString filenameArea(outPath);
    2966     filenameArea += typeMC;
    2967     filenameArea += ext;
    2968     gLog << "filenameArea = " << filenameArea << endl;
    2969 
    2970     //------------------------------
    2971     // name of file containing the eff. collection areas
    2972     TString collareaName(outPath);
    2973     collareaName += "area_";
    2974     collareaName += XX;
    2975     collareaName += ".root";
    2976     gLog << "collareaName = " << collareaName << endl;
    2977 
    2978     //------------------------------
    2979     // name of data file to be analysed
    2980     TString filenameData(outPath);
    2981     filenameData += typeData;
    2982     filenameData += ext;
    2983     gLog << "filenameData = " << filenameData << endl;
    2984 
    2985     //------------------------------
    2986     // name of output data file (after the final cuts)
    2987     TString filenameDataout(outPath);
    2988     filenameDataout += typeData;
    2989     filenameDataout += "_";
    2990     filenameDataout += XX;
    2991     filenameDataout += extout;
    2992     gLog << "filenameDataout = " << filenameDataout << endl;
    2993 
    2994     //------------------------------
    2995     // name of file containing histograms for flux calculastion
    2996     TString filenameResults(outPath);
    2997     filenameResults += typeData;
    2998     filenameResults += "Results_";
    2999     filenameResults += XX;
    3000     filenameResults += extout;
    3001     gLog << "filenameResults = " << filenameResults << endl;
    3002 
    3003 
    3004     //====================================================================
    3005     gLog << "-----------------------------------------------" << endl;
    3006     gLog << "Start calculation of effective collection areas" << endl;
    3007     MParList  parlist;
    3008     MTaskList tasklist;
    3009 
    3010     //---------------------------------------
    3011     // Setup the tasks to be executed
    3012     //
    3013     MReadMarsFile reader("Events", filenameArea);
    3014     reader.DisableAutoScheme();
    3015 
    3016     MFCT1SelFinal cuthadrons;
    3017     cuthadrons.SetHadronnessName(fhadronnessName);
    3018     cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist);
    3019 
    3020     MContinue conthadrons(&cuthadrons);
    3021 
    3022     MHMcCT1CollectionArea collarea;
    3023     collarea.SetEaxis(MHMcCT1CollectionArea::kLinear);
    3024 
    3025     MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
    3026     filler.SetName("CollectionArea");
    3027 
    3028     //********************************
    3029     // entries in MParList
    3030 
    3031     parlist.AddToList(&tasklist);
    3032     InitBinnings(&parlist);
    3033     parlist.AddToList(&collarea);
    3034 
    3035     //********************************
    3036     // entries in MTaskList
    3037 
    3038     tasklist.AddToList(&reader);   
    3039     tasklist.AddToList(&conthadrons);
    3040     tasklist.AddToList(&filler);
    3041 
    3042     //********************************
    3043 
    3044     //-----------------------------------------
    3045     // Execute event loop
    3046     //
    3047     MEvtLoop magic;
    3048     magic.SetParList(&parlist);
    3049 
    3050     MProgressBar bar;
    3051     magic.SetProgressBar(&bar);
    3052     if (!magic.Eventloop())
    3053         return;
    3054 
    3055     tasklist.PrintStatistics(0, kTRUE);
    3056 
    3057     // Calculate effective collection areas
    3058     // and display the histograms
    3059     //
    3060     //MHMcCT1CollectionArea *collarea =
    3061     //     (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
    3062     collarea.CalcEfficiency();
    3063     collarea.DrawClone();
    3064 
    3065     // save binnings for call to CT1EEst
    3066     MBinning *binsE     = (MBinning*)parlist.FindObject("BinningE");
    3067     if (!binsE)
    3068         {
    3069           gLog << "Object 'BinningE' not found in MParList" << endl;
    3070           return;
    3071         }
    3072     MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta");
    3073     if (!binsTheta)
    3074         {
    3075           gLog << "Object 'BinningTheta' not found in MParList" << endl;
    3076           return;
    3077         }
    3078 
    3079 
    3080     //---------------------------------------------
    3081     // Write histograms to a file
    3082     //
    3083 
    3084     TFile f(collareaName, "RECREATE");
    3085     collarea.GetHist()->Write();
    3086     collarea.GetHAll()->Write();
    3087     collarea.GetHSel()->Write();
    3088     f.Close();
    3089 
    3090     gLog << "Collection area plots written onto file " << collareaName << endl;
    3091 
    3092     gLog << "Calculation of effective collection areas done" << endl;
    3093     gLog << "-----------------------------------------------" << endl;   
    3094     //------------------------------------------------------------------
    3095 
    3096 
    3097     TString fHilName    = "MHillas";
    3098     TString fHilNameExt = "MHillasExt";
    3099     TString fHilNameSrc = "MHillasSrc";
    3100     TString fImgParName = "MNewImagePar";
    3101 
    3102  
    3103     //===========================================================
    3104     //
    3105     // Optimization of energy estimator
    3106     //
    3107     gLog << "Macro CT1Analysis.C : calling CT1EEst" << endl;
    3108 
    3109     TString inpath("");
    3110     TString outpath("");
    3111     Int_t howMany = 2000;
    3112     CT1EEst(inpath,   filenameOpt,   outpath, energyParName,
    3113             fHilName, fHilNameSrc,   fhadronnessName,
    3114             howMany,  maxhadronness, maxalpha, maxdist,
    3115             binsE, binsTheta);
    3116     gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
    3117 
    3118 
    3119   if (WFX)
    3120   {
    3121     //-----------------------------------------------------------
    3122     //
    3123     // Read in parameters of energy estimator ("MMcEnergyEst")
    3124     //                   and migration matrix ("MHMcEnergyMigration")
    3125     //
    3126     gLog << "================================================================"
    3127          << endl;
    3128     gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
    3129          << energyParName << "'" << endl;
    3130     TFile enparam(energyParName);
    3131     enparam.ls();
    3132 
    3133     //MMcEnergyEst mcest("MMcEnergyEst");
    3134     //mcest.Read("MMcEnergyEst");
    3135     MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
    3136 
    3137     gLog << "Parameters of energy estimator were read in" << endl;
    3138 
    3139     TArrayD parA(mcest.GetNumCoeffA());
    3140     TArrayD parB(mcest.GetNumCoeffB());
    3141     for (Int_t i=0; i<parA.GetSize(); i++)
    3142       parA[i] = mcest.GetCoeff(i);
    3143     for (Int_t i=0; i<parB.GetSize(); i++)
    3144       parB[i] = mcest.GetCoeff( i+parA.GetSize() );
    3145 
    3146     gLog << "Read in Migration matrix" << endl;   
    3147 
    3148     //MHMcEnergyMigration mighiston("MHMcEnergyMigration");
    3149     //mighiston.Read("MHMcEnergyMigration");
    3150     MHMcEnergyMigration &mighiston =
    3151           *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
    3152 
    3153     gLog << "Migration matrix was read in" << endl;
    3154 
    3155     //*************************************************************************
    3156     //
    3157     // Analyse the data
    3158     //
    3159     gLog << "============================================================"
    3160          << endl;
    3161     gLog << "Analyse the data" << endl;
    3162 
    3163     MTaskList tliston;
    3164     MParList pliston;
    3165 
    3166     // geometry is needed in  MHHillas... classes
    3167     MGeomCam *fGeom =
    3168              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    3169 
    3170 
    3171     //-------------------------------------------
    3172     // create the tasks which should be executed
    3173     //
    3174 
    3175     MReadMarsFile read("Events", filenameData);
    3176     read.DisableAutoScheme();
    3177 
    3178     //.......................................................................
    3179 
    3180       gLog << "CT1Analysis.C : write root file '" << filenameDataout
    3181            << "'" << endl;
    3182    
    3183       //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
    3184       /*
    3185       MWriteRootFile write(filenameDataout, "RECREATE");
    3186 
    3187       write.AddContainer("MRawRunHeader", "RunHeaders");
    3188       write.AddContainer("MTime",         "Events");
    3189       write.AddContainer("MMcEvt",        "Events");
    3190       write.AddContainer("ThetaOrig",     "Events");
    3191       write.AddContainer("MSrcPosCam",    "Events");
    3192       write.AddContainer("MSigmabar",     "Events");
    3193       write.AddContainer("MHillas",       "Events");
    3194       write.AddContainer("MHillasExt",    "Events");
    3195       write.AddContainer("MHillasSrc",    "Events");
    3196       write.AddContainer("MNewImagePar",  "Events");
    3197 
    3198       //write.AddContainer("HadNN",         "Events");
    3199       write.AddContainer("HadSC",         "Events");
    3200       write.AddContainer("HadRF",         "Events");
    3201 
    3202       write.AddContainer("MEnergyEst",    "Events");
    3203       */
    3204 
    3205     //-----------------------------------------------------------------
    3206     // geometry is needed in  MHHillas... classes
    3207     MGeomCam *fGeom =
    3208              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    3209 
    3210     MFCT1SelFinal selfinalgh(fHilNameSrc);
    3211     selfinalgh.SetCuts(maxhadronness, 100.0, maxdist);
    3212     selfinalgh.SetHadronnessName(fhadronnessName);
    3213     selfinalgh.SetName("SelFinalgh");
    3214     MContinue contfinalgh(&selfinalgh);
    3215     contfinalgh.SetName("ContSelFinalgh");
    3216 
    3217     //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
    3218     //fillhadnn.SetName("HhadNN");
    3219     MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
    3220     fillhadsc.SetName("HhadSC");
    3221     MFillH fillhadrf("hadRF[MHHadronness]", "HadRF");
    3222     fillhadrf.SetName("HhadRF");
    3223 
    3224     //---------------------------
    3225     // calculate estimated energy
    3226 
    3227     MEnergyEstParam eeston(fHilName);
    3228     eeston.Add(fHilNameSrc);
    3229 
    3230     eeston.SetCoeffA(parA);
    3231     eeston.SetCoeffB(parB);
    3232 
    3233     //---------------------------
    3234     // calculate estimated energy using Daniel's parameters
    3235 
    3236     //MEnergyEstParamDanielMkn421 eeston(fHilName);
    3237     //eeston.Add(fHilNameSrc);
    3238     //eeston.SetCoeffA(parA);
    3239     //eeston.SetCoeffB(parB);
    3240 
    3241 
    3242     //---------------------------
    3243 
    3244 
    3245     MFillH hfill1("MHHillas",    fHilName);
    3246     hfill1.SetName("HHillas");
    3247 
    3248     MFillH hfill2("MHStarMap",   fHilName);
    3249     hfill2.SetName("HStarMap");
    3250 
    3251     MFillH hfill3("MHHillasExt",   fHilNameSrc);
    3252     hfill3.SetName("HHillasExt");   
    3253 
    3254     MFillH hfill4("MHHillasSrc",   fHilNameSrc);
    3255     hfill4.SetName("HHillasSrc");   
    3256 
    3257     MFillH hfill5("MHNewImagePar", fImgParName);
    3258     hfill5.SetName("HNewImagePar");   
    3259 
    3260     //---------------------------
    3261     // new from Robert
    3262 
    3263     MFillH hfill6("MHTimeDiffTheta", "MMcEvt");
    3264     hfill6.SetName("HTimeDiffTheta");
    3265 
    3266     MFillH hfill6a("MHTimeDiffTime", "MMcEvt");
    3267     hfill6a.SetName("HTimeDiffTime");
    3268 
    3269     MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc);
    3270     hfill7.SetName("HAlphaEnergyTheta");
    3271 
    3272     MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc);
    3273     hfill7a.SetName("HAlphaEnergyTime");
    3274 
    3275     MFillH hfill7b("MHThetabarTime", fHilNameSrc);
    3276     hfill7b.SetName("HThetabarTime");
    3277 
    3278     MFillH hfill7c("MHEnergyTime", "MMcEvt");
    3279     hfill7c.SetName("HEnergyTime");
    3280 
    3281 
    3282     //---------------------------
    3283 
    3284     MFCT1SelFinal selfinal(fHilNameSrc);
    3285     selfinal.SetCuts(maxhadronness, maxalpha, maxdist);
    3286     selfinal.SetHadronnessName(fhadronnessName);
    3287     selfinal.SetName("SelFinal");
    3288     MContinue contfinal(&selfinal);
    3289     contfinal.SetName("ContSelFinal");
    3290 
    3291 
    3292     //*****************************
    3293     // entries in MParList
    3294 
    3295     pliston.AddToList(&tliston);
    3296     InitBinnings(&pliston);
    3297 
    3298 
    3299     //*****************************
    3300     // entries in MTaskList
    3301    
    3302     tliston.AddToList(&read);
    3303 
    33042880    // robert     
    33052881    tliston.AddToList(&hfill6);   //timediff
    … …  
    33092885    tliston.AddToList(&eeston);
    33102886
    3311     //tliston.AddToList(&write);
     2887    tliston.AddToList(&write);
    33122888
    33132889    //tliston.AddToList(&fillhadnn);
Note: See TracChangeset for help on using the changeset viewer.