Ignore:
Timestamp:
10/27/03 08:58:48 (21 years ago)
Author:
wittek
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2368 r2436  
    210210      //-----------------------------------------------
    211211
    212       TEnv env("macros/CT1env.rc");
    213       Bool_t printEnv = kFALSE;
     212      //TEnv env("macros/CT1env.rc");
     213      //Bool_t printEnv = kFALSE;
    214214
    215215    //************************************************************************
    … …  
    219219    //  - write root file for ON data (ON1.root);
    220220
    221     Bool_t JobA_ON = kTRUE; 
    222     Bool_t WHistON = kTRUE;   // write out histogram sigmabar vs. Theta ?
    223     Bool_t WON1    = kTRUE;   // write out root file ON1.root ?
     221    Bool_t JobA_ON = kFALSE; 
     222    Bool_t WHistON = kFALSE;   // write out histogram sigmabar vs. Theta ?
     223    Bool_t WON1    = kFALSE;   // write out root file ON1.root ?
    224224
    225225
    … …  
    234234
    235235    // Job B_RF_UP : read ON1.root (or MC1.root) file
    236     //  - depending on RLook : create (or read in) hadron and gamma matrix
    237     //  - depending on RTree : create (or read in) trees
     236    //  - if RTrainRF = TRUE : read in training matrices for hadrons and gammas
     237    //  - if CTrainRF = TRUE : create  training matrices for hadrons and gammas
     238    //  - if RTree    = TRUE : read in trees, otherwise create trees
    238239    //  - calculate hadroness for method of RANDOM FOREST
    239     //  - update the input files with the hadroness (ON1.root or MC1.root)
     240    //  - update the input files with the hadroness (==>ON2.root or MC2.root)
    240241
    241242    Bool_t JobB_RF_UP  = kFALSE; 
    242     Bool_t RLookRF     = kFALSE;  // read in look-alike events
     243    Bool_t RTrainRF    = kFALSE;  // read in matrices of training events
     244    Bool_t CTrainRF    = kFALSE;  // create matrices of training events
    243245    Bool_t RTree       = kFALSE;  // read in trees
    244246    Bool_t WRF         = kFALSE;  // update input root file ?
    … …  
    247249
    248250
    249     // Job B_SC_UP : read ON1.root (or MC1.root) file
     251    // Job B_SC_UP : read ON2.root (or MC2.root) file
    250252    //  - depending on WParSC : create (or read in) supercuts parameter values
    251253    //  - calculate hadroness for the SUPERCUTS
    252     //  - update the input files with the hadroness (ON1.root or MC1.root)
     254    //  - update the input files with the hadroness (==>ON3.root or MC3.root)
    253255
    254256    Bool_t JobB_SC_UP  = kFALSE;
    … …  
    263265
    264266    // Job C:
    265     //  - read ON1.root and MC1.root files
     267    //  - read ON3.root and MC3.root files
    266268    //    which should have been updated to contain the hadronnesses 
    267269    //    for the method of
    268     //              NEAREST NEIGHBORS 
     270    //              RF
    269271    //              SUPERCUTS and
    270     //              RF
    271272    //  - produce Neyman-Pearson plots
    272273
    … …  
    276277    // Job D : 
    277278    //  - select g/h separation method XX
    278     //  - read ON2 (or MC2) root file
     279    //  - read ON3 (or MC3) root file
    279280    //  - apply cuts in hadronness
    280281    //  - make plots
    … …  
    296297
    297298    Bool_t JobE_XX  = kFALSE; 
    298     Bool_t WEX      = kFALSE;  // write out root file  ?
     299    Bool_t OEEst    = kFALSE;  // optimize energy estimator
     300    Bool_t WEX      = kFALSE;  // update root file  ?
    299301    Bool_t WRobert  = kFALSE;  // write out Robert's file  ?
    300302
    … …  
    492494    MEvtLoop evtloop;
    493495    evtloop.SetParList(&pliston);
    494     evtloop.ReadEnv(env, "", printEnv);
     496    //evtloop.ReadEnv(env, "", printEnv);
    495497    evtloop.SetProgressBar(&bar);
    496498    //if (WON1)
    … …  
    822824    MEvtLoop evtloop;
    823825    evtloop.SetParList(&plist);
    824     evtloop.ReadEnv(env, "", printEnv);
     826    //evtloop.ReadEnv(env, "", printEnv);
    825827    evtloop.SetProgressBar(&bar);
    826828    //if (WMC1)   
    … …  
    863865
    864866
    865     //  - create (or read in) the matrices of look-alike events for gammas
     867    //  - create (or read in) the matrices of training events for gammas
    866868    //    and hadrons
    867869    //  - create (or read in) the trees
    … …  
    877879
    878880    gLog << "" << endl;
    879     gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = "
     881    gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
    880882         << (JobB_RF_UP ? "kTRUE" : "kFALSE")  << ",  "
    881          << (RLookRF    ? "kTRUE" : "kFALSE")  << ",  "
     883         << (RTrainRF    ? "kTRUE" : "kFALSE")  << ",  "
     884         << (CTrainRF    ? "kTRUE" : "kFALSE")  << ",  "
    882885         << (RTree      ? "kTRUE" : "kFALSE")  << ",  "
    883886         << (WRF        ? "kTRUE" : "kFALSE")  << endl;
    … …  
    890893    Int_t NdSize   =   1;
    891894
     895    // cut in RF hadronness
     896    TString hadRFName = "HadRF";
     897    Float_t maxhadronness =  0.23;
     898    Float_t maxalpha      =  20.0;
     899    Float_t maxdist       =  10.0;
     900
     901    TString fHilName    = "MHillas";
     902    TString fHilNameExt = "MHillasExt";
     903    TString fHilNameSrc = "MHillasSrc";
     904    TString fImgParName = "MNewImagePar";
     905
    892906
    893907    //--------------------------------------------
    894908    // file to be updated (either ON or MC)
    895909
    896     TString typeInput = "ON";
    897     //TString typeInput = "MC";
     910    //TString typeInput = "ON";
     911    TString typeInput = "MC";
    898912    gLog << "typeInput = " << typeInput << endl;
    899913
    … …  
    913927
    914928    //--------------------------------------------
    915     // files to be read for generating the look-alike events
     929    // files to be read for generating the matrices of training events
    916930    // "hadrons" :
    917931    TString filenameON = outPath;
    … …  
    932946   
    933947    //--------------------------------------------
    934     // files of look-alike events
     948    // files of training events
    935949
    936950    TString outNameGammas = outPath;
    … …  
    951965    matrixg.EnableGraphicalOutput();
    952966
     967    matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
     968    matrixg.AddColumn("MSigmabar.fSigmabar");
     969    matrixg.AddColumn("log10(MHillas.fSize)");
     970    matrixg.AddColumn("MHillasSrc.fDist");
     971    matrixg.AddColumn("MHillas.fWidth");
     972    matrixg.AddColumn("MHillas.fLength");
     973    matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
     974    matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
     975    matrixg.AddColumn("MNewImagePar.fConc");
     976    matrixg.AddColumn("MNewImagePar.fLeakage1");
     977
    953978    MHMatrix matrixh("MatrixHadrons");
    954979    matrixh.EnableGraphicalOutput();
    955980
     981    matrixh.AddColumns(matrixg.GetColumns());
     982
    956983    //--------------------------------------------
    957984    // file of trees of the random forest
    … …  
    962989
    963990   //*************************************************************************
    964    // read in matrices of look-alike events
    965 if (RLookRF)
     991   // read in matrices of training events
     992if (RTrainRF)
    966993  {
    967994    const char* mtxName = "MatrixGammas";
    … …  
    10051032
    10061033   //*************************************************************************
    1007    // create matrices of look-alike events
    1008 if (!RLookRF)
     1034   // create matrices of training events
     1035if (CTrainRF)
    10091036  {
    10101037    gLog << "" << endl;
    10111038    gLog << "========================================================" << endl;
    1012     gLog << " Create matrices of look-alike events" << endl;
     1039    gLog << " Create matrices of training events" << endl;
    10131040    gLog << " Gammas :" << endl;
    10141041
    … …  
    10851112    MEvtLoop evtloopg;
    10861113    evtloopg.SetParList(&plistg);
    1087     evtloopg.ReadEnv(env, "", printEnv);
     1114    //evtloopg.ReadEnv(env, "", printEnv);
    10881115    evtloopg.SetProgressBar(&matrixbar);
    10891116
    … …  
    11411168    MEvtLoop evtlooph;
    11421169    evtlooph.SetParList(&plisth);
    1143     evtlooph.ReadEnv(env, "", printEnv);
     1170    //evtlooph.ReadEnv(env, "", printEnv);
    11441171    evtlooph.SetProgressBar(&matrixbar);
    11451172
    … …  
    11551182
    11561183
    1157     // write out look-alike events
     1184    // write out training events
    11581185
    11591186    gLog << "" << endl;
    11601187    gLog << "========================================================" << endl;
    1161     gLog << "Write out look=alike events" << endl;
     1188    gLog << "Write out training events" << endl;
    11621189
    11631190
    … …  
    11711198
    11721199      gLog << "" << endl;
    1173       gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
     1200      gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file "
    11741201           << outNameGammas << endl;
    11751202
    … …  
    11831210
    11841211      gLog << "" << endl;
    1185       gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
     1212      gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file "
    11861213           << outNameHadrons << endl;
    11871214
    11881215  }
    1189    //**********   end of creating matrices of look-alike events   ***********
     1216   //**********   end of creating matrices of training events   ***********
    11901217
    11911218
    … …  
    13011328
    13021329
    1303 
    13041330    //-----------------------------------------------------------------
    13051331    // Update the input files with the RF hadronness
    13061332    //
    1307 
    13081333 if (WRF)
    13091334  {
    … …  
    13281353    read.DisableAutoScheme();
    13291354
    1330     TString fHilName    = "MHillas";
    1331     TString fHilNameExt = "MHillasExt";
    1332     TString fHilNameSrc = "MHillasSrc";
    1333     TString fImgParName = "MNewImagePar";
    13341355
    13351356    //.......................................................................
    13361357    // calculate hadronnes for method of RANDOM FOREST
    13371358
    1338     TString hadRFName = "HadRF";
     1359
    13391360    MRanForestCalc rfcalc;
    13401361    rfcalc.SetHadronnessName(hadRFName);
    … …  
    13601381
    13611382    //-----------------------------------------------------------------
    1362     // geometry is needed in  MHHillas... classes
    1363     MGeomCam *fGeom =
    1364              (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
    1365 
    1366 
    1367     Float_t maxhadronness =  0.40;
    1368     Float_t maxalpha      =  20.0;
    1369     Float_t maxdist       =  10.0;
     1383
    13701384
    13711385    MFCT1SelFinal selfinalgh(fHilNameSrc);
    … …  
    14821496    TH1  &alphaHist = absalpha->GetHist();
    14831497    alphaHist.SetXTitle("|alpha|  [\\circ]");
     1498    alphaHist.SetName("alpha-macro");
    14841499
    14851500    Double_t alphasig = 13.1;
    14861501    Double_t alphamin = 30.0;
    14871502    Double_t alphamax = 90.0;
    1488     Int_t    degree   =    4;
     1503    Int_t    degree   =    2;
    14891504    Double_t significance = -99.0;
    14901505    Bool_t   drawpoly  = kTRUE;
    … …  
    15631578
    15641579    TString parSCfile = outPath;
    1565     parSCfile += "parSC_2209a";
     1580    parSCfile += "parSC_2310a";
    15661581
    15671582    gLog << "parSCfile = " << parSCfile << endl;
    … …  
    15701585    // file to be updated (either ON or MC)
    15711586
    1572     TString typeInput = "ON";
    1573     //TString typeInput = "MC";
     1587    //TString typeInput = "ON";
     1588    TString typeInput = "MC";
    15741589    gLog << "typeInput = " << typeInput << endl;
    15751590
    … …  
    16121627
    16131628    //--------------------------------------------
    1614     // files to contain the matrices (genertated from filenameTrain and
    1615     //                                                filenameTest)
     1629    // files to contain the matrices (generated from filenameTrain and
     1630    //                                               filenameTest)
    16161631    //
    16171632    // for the training
    … …  
    16421657    //--------------------------
    16431658    // create matrices and write them onto files
    1644     if (!RMatrix  &&  (WOptimize || RTest)  )
     1659    if (!RMatrix)
    16451660    {
    16461661      MH3 &mh3 = *(new MH3("MHillas.fSize"));
    … …  
    19111926      //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
    19121927
    1913     /*
     1928   
    19141929      MWriteRootFile write(outNameImage, "RECREATE");
    19151930
    … …  
    19271942      write.AddContainer("HadRF",         "Events");
    19281943      write.AddContainer(hadSCName,       "Events");
    1929     */
     1944   
    19301945
    19311946    //-----------------------------------------------------------------
    … …  
    20012016    tliston.AddToList(&fillhadsc);
    20022017
    2003     //tliston.AddToList(&write);
     2018    tliston.AddToList(&write);
    20042019    tliston.AddToList(&contfinalgh);
    20052020
    … …  
    22922307    // and definition of final selections
    22932308
    2294     //TString XX("NN");
    2295     TString XX("SC");
    2296     //TString XX("RF");
     2309    //TString XX("SC");
     2310    TString XX("RF");
    22972311    TString fhadronnessName("Had");
    22982312    fhadronnessName += XX;
    … …  
    23002314
    23012315    // maximum values of the hadronness, |ALPHA| and DIST
    2302     Float_t maxhadronness   = 0.30;
     2316    Float_t maxhadronness   = 0.233;
    23032317    Float_t maxalpha        = 20.0;
    23042318    Float_t maxdist         = 10.0;
    … …  
    23552369    contfinalgh.SetName("ContSelFinalgh");
    23562370
    2357     MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
    2358     fillhadnn.SetName("HhadNN");
    23592371    MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
    23602372    fillhadsc.SetName("HhadSC");
    … …  
    23622374    fillhadrf.SetName("HhadRF");
    23632375
     2376    TString mh3name = "abs(Alpha)";
     2377    MBinning binsalphaabs("Binning"+mh3name);
     2378    binsalphaabs.SetEdges(50, -2.0, 98.0);
     2379
     2380    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
     2381    alphaabs.SetName(mh3name);
     2382
     2383    TH1  &alphahist = alphaabs->GetHist();
     2384
     2385    MFillH alpha(&alphaabs);
     2386    alpha.SetName("FillAlphaAbs");
    23642387
    23652388    MFillH hfill1("MHHillas",    fHilName);
    … …  
    23912414    pliston.AddToList(&tliston);
    23922415    InitBinnings(&pliston);
    2393 
     2416    pliston.AddToList(&binsalphaabs);
     2417    pliston.AddToList(&alphaabs);
    23942418
    23952419    //*****************************
    … …  
    24002424    tliston.AddToList(&contfinalgh);
    24012425
    2402     tliston.AddToList(&fillhadnn);
    24032426    tliston.AddToList(&fillhadsc);
    24042427    tliston.AddToList(&fillhadrf);
    24052428
     2429    tliston.AddToList(&alpha);
    24062430    tliston.AddToList(&hfill1);
    24072431    tliston.AddToList(&hfill2);
    … …  
    24342458    //
    24352459
    2436     pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
    24372460    pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
    24382461    pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
    … …  
    24462469
    24472470    //-------------------------------------------
    2448     // fit alpha distribution to get the number of excess events
    2449     //
    2450 
    2451     MHHillasSrc* hillasSrc =
    2452       (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
    2453     TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
     2471
     2472    // fit alpha distribution to get the number of excess events and
     2473    // calculate significance of gamma signal in the alpha plot
    24542474 
    2455     MHOnSubtraction onsub;
    2456     onsub.Calc(alphaHist, kTRUE, 13.1);
     2475    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
     2476    TH1  &alphaHist = absalpha->GetHist();
     2477    alphaHist.SetXTitle("|alpha|  [\\circ]");
     2478    alphaHist.SetName("alpha-JobD");
     2479
     2480    Double_t alphasig = 13.1;
     2481    Double_t alphamin = 30.0;
     2482    Double_t alphamax = 90.0;
     2483    Int_t    degree   =    2;
     2484    Double_t significance = -99.0;
     2485    Bool_t   drawpoly  = kTRUE;
     2486    Bool_t   fitgauss  = kTRUE;
     2487    Bool_t   print     = kTRUE;
     2488
     2489    MHFindSignificance findsig;
     2490    findsig.SetRebin(kTRUE);
     2491    findsig.SetReduceDegree(kFALSE);
     2492
     2493    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree,
     2494                        alphasig, drawpoly, fitgauss, print);
     2495    significance = findsig.GetSignificance();
     2496    Float_t alphasi = findsig.GetAlphasi();
     2497
     2498    gLog << "For file '" << filenameData << "' : " << endl;
     2499    gLog << "Significance of gamma signal after supercuts : "
     2500         << significance << " (for |alpha| < " << alphasi << " degrees)"
     2501         << endl;
     2502
     2503    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
     2504
    24572505    //-------------------------------------------
    24582506
    … …  
    24882536
    24892537    gLog << "" << endl;
    2490     gLog << "Macro CT1Analysis : JobE_XX, WEX = "
     2538    gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = "
    24912539         << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  "
     2540         << (OEEst ?   "kTRUE" : "kFALSE")  << ",  "
    24922541         << (WEX     ? "kTRUE" : "kFALSE")  << endl;
    24932542
    … …  
    25072556    // and definition of final selections
    25082557
    2509     //TString XX("NN");
    25102558    TString XX("SC");
    25112559    //TString XX("RF");
    … …  
    25152563
    25162564    // maximum values of the hadronness, |ALPHA| and DIST
    2517     Float_t maxhadronness   = 0.40;
     2565    Float_t maxhadronness   = 0.4;
    25182566    Float_t maxalpha        = 20.0;
    25192567    Float_t maxdist         = 10.0;
    … …  
    26762724    TString fImgParName = "MNewImagePar";
    26772725
    2678  
    2679     //===========================================================
     2726
     2727 if (OEEst)
     2728 {
     2729   //===========================================================
    26802730    //
    26812731    // Optimization of energy estimator
    … …  
    26912741            binsE, binsTheta);
    26922742    gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
    2693 
    2694 
    2695   if (WEX)
    2696   {
     2743 }
     2744
     2745 if (WEX)
     2746 {
    26972747    //-----------------------------------------------------------
    26982748    //
    … …  
    27062756    TFile enparam(energyParName);
    27072757    enparam.ls();
    2708 
    2709     //MMcEnergyEst mcest("MMcEnergyEst");
    2710     //mcest.Read("MMcEnergyEst");
    2711     MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
    2712 
     2758    MMcEnergyEst mcest("MMcEnergyEst");
     2759    mcest.Read("MMcEnergyEst");
     2760
     2761    //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
    27132762    gLog << "Parameters of energy estimator were read in" << endl;
     2763
     2764
     2765    gLog << "Read in Migration matrix" << endl;   
     2766
     2767    MHMcEnergyMigration mighiston("MHMcEnergyMigration");
     2768    mighiston.Read("MHMcEnergyMigration");
     2769    //MHMcEnergyMigration &mighiston =
     2770    //      *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
     2771
     2772    gLog << "Migration matrix was read in" << endl;
     2773
    27142774
    27152775    TArrayD parA(mcest.GetNumCoeffA());
    … …  
    27202780      parB[i] = mcest.GetCoeff( i+parA.GetSize() );
    27212781
    2722     gLog << "Read in Migration matrix" << endl;   
    2723 
    2724     //MHMcEnergyMigration mighiston("MHMcEnergyMigration");
    2725     //mighiston.Read("MHMcEnergyMigration");
    2726     MHMcEnergyMigration &mighiston =
    2727           *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
    2728 
    2729     gLog << "Migration matrix was read in" << endl;
    2730 
    27312782    //*************************************************************************
    27322783    //
    … …  
    27582809   
    27592810      //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
     2811
    27602812
    27612813      MWriteRootFile write(filenameDataout, "RECREATE");
    … …  
    29162968
    29172969    Int_t maxevents = -1;
    2918     //Int_t maxevents = 10000;
    29192970    if ( !evtloop.Eventloop(maxevents) )
    29202971        return;
    … …  
    29523003    gLog << "before DeleteBinnings" << endl;
    29533004
    2954     //DeleteBinnings(&pliston);
     3005    DeleteBinnings(&pliston);
    29553006
    29563007    gLog << "before Robert's code" << endl;
    … …  
    30363087  }
    30373088
    3038     gLog << "Macro CT1Analysis : End of Job F_XX" << endl;
     3089    gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
    30393090    gLog << "=======================================================" << endl;
    30403091 }
Note: See TracChangeset for help on using the changeset viewer.