Index: /trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestPD.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestPD.C (revision 5122) +++ /trunk/MagicSoft/Mars/mtemp/mpadova/macros/RanForestPD.C (revision 5122) @@ -0,0 +1,302 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Hengstebeck 3/2003 ! +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +#include "MMcEvt.hxx" + +void RanForestPD(Int_t minsize=0) +{ + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + MReadMarsFile read("Events", "star_gamma_train.root"); + Float_t numgammas = read.GetEntries(); + + read.AddFile("star_proton_train.root"); + read.AddFile("star_helium_train.root"); + + Float_t numhadrons = read.GetEntries() - numgammas; + + // Fraction of gammas to be used in training: + // (less than total sample, to speed up execution) + // + Float_t gamma_frac = 1.5*(numhadrons / numgammas); + + + read.DisableAutoScheme(); + + tlist.AddToList(&read); + + MFParticleId fgamma("MMcEvt", '=', MMcEvt::kGAMMA); + tlist.AddToList(&fgamma); + + MFParticleId fhadrons("MMcEvt", '!', MMcEvt::kGAMMA); + tlist.AddToList(&fhadrons); + + + MHMatrix matrixg("MatrixGammas"); + + // matrixg.AddColumn("floor(0.5+(cos(MMcEvt.fTelescopeTheta)/0.01))"); + + // matrixg.AddColumn("MMcEvt.fTelescopePhi"); + // matrixg.AddColumn("MSigmabar.fSigmabar"); + + matrixg.AddColumn("MHillas.fSize"); + matrixg.AddColumn("MHillasSrc.fDist"); + + // matrixg.AddColumn("MHillas.fWidth"); + // matrixg.AddColumn("MHillas.fLength"); + + // Scaled Width and Length: + // + matrixg.AddColumn("MHillas.fWidth/(-13.1618+10.0492*log10(MHillas.fSize))"); + matrixg.AddColumn("MHillas.fLength/(-57.3784+32.6131*log10(MHillas.fSize))"); + + + // matrixg.AddColumn("(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); + matrixg.AddColumn("MNewImagePar.fConc"); + + // matrixg.AddColumn("MNewImagePar.fConc1"); + // matrixg.AddColumn("MNewImagePar.fAsym"); + // matrixg.AddColumn("MHillasExt.fM3Trans"); + // matrixg.AddColumn("abs(MHillasSrc.fHeadTail)"); + // matrixg.AddColumn("MNewImagePar.fLeakage1"); + + // + // Third moment, with sign referred to source position: + // + matrixg.AddColumn("MHillasExt.fM3Long*MHillasSrc.fCosDeltaAlpha/abs(MHillasSrc.fCosDeltaAlpha)"); + matrixg.AddColumn("MHillasSrc.fAlpha"); + + + // + // SIZE cut: + // + TString dummy("(MHillas.fSize<"); + dummy += minsize; + + // AM, useful FOR High Energy only: + // dummy += ") || (MImagePar.fNumSatPixelsLG>0"; + dummy += ")"; + + MContinue SizeCut(dummy); + tlist.AddToList(&SizeCut); + + plist.AddToList(&matrixg); + + MHMatrix matrixh("MatrixHadrons"); + matrixh.AddColumns(matrixg.GetColumns()); + plist.AddToList(&matrixh); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&fhadrons); + tlist.AddToList(&fillmath); + + + MFEventSelector reduce_training_gammas; + reduce_training_gammas.SetSelectionRatio(gamma_frac); + tlist.AddToList(&reduce_training_gammas); + + MFillH fillmatg("MatrixGammas"); + fillmatg.SetFilter(&fgamma); + fillmatg.SetFilter(&reduce_training_gammas); + tlist.AddToList(&fillmatg); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + MProgressBar bar; + evtloop.SetProgressBar(&bar); + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // second event loop: the trees of the random forest are grown, + // the event loop is now actually a tree loop (loop of training + // process) + // --------------------------------------------------------------- + // --------------------------------------------------------------- + + MTaskList tlist2; + plist.Replace(&tlist2); + + // + // Just a check: + // matrixh.ShuffleRows(9999); + // matrixg.ShuffleRows(1111); + // + + MRanForestGrow rfgrow2; + rfgrow2.SetNumTrees(100); + rfgrow2.SetNumTry(3); + rfgrow2.SetNdSize(10); + + MWriteRootFile rfwrite2("RF.root"); + rfwrite2.AddContainer("MRanTree","Tree"); //save all trees + MFillH fillh2("MHRanForestGini"); + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + tlist2.AddToList(&fillh2); + + // gRandom is accessed from MRanForest (-> bootstrap aggregating) + // and MRanTree (-> random split selection) and should be initialized + // here if you want to set a certain random number generator + if(gRandom) + delete gRandom; + // gRandom = new TRandom3(0); // Takes seed from the computer clock + gRandom = new TRandom3(); // Uses always same seed + // + // Execute tree growing (now the eventloop is actually a treeloop) + // + + evtloop.SetProgressBar(&bar); + if (!evtloop.Eventloop()) + return; + + tlist2.PrintStatistics(); + + plist.FindObject("MHRanForestGini")->DrawClone(); + + // --------------------------------------------------------------- + // --------------------------------------------------------------- + // third event loop: the control sample (star2.root) is processed + // through the previously grown random forest, + // + // the histograms MHHadronness (quality of g/h-separation) and + // MHRanForest are created and displayed. + // MHRanForest shows the convergence of the classification error + // as function of the number of grown (and combined) trees + // and tells the user how many trees are actually needed in later + // classification tasks. + // --------------------------------------------------------------- + // --------------------------------------------------------------- + + MTaskList tlist3; + + plist.Replace(&tlist3); + + MReadMarsFile read3("Events", "star_gamma_test.root"); + + read3.AddFile("star_proton_test.root"); + read3.AddFile("star_helium_test.root"); + + + /* + MReadMarsFile read3("Events", "star_gammas_test_extended.root"); + read3.AddFile("star_protons_test_extended.root"); + read3.AddFile("star_helium_test_extended.root"); + */ + + read3.DisableAutoScheme(); + tlist3.AddToList(&read3); + tlist3.AddToList(&SizeCut); + + MRanForestCalc calc; + tlist3.AddToList(&calc); + + // parameter containers you want to save + // + + TString outfname = "star_S"; + outfname += minsize; + outfname += "_all.root"; + MWriteRootFile write(outfname, "recreate"); + + + // MWriteRootFile write("star_surviving_hadrons.root", "recreate"); + // MWriteRootFile write("star_surviving_gammas.root", "recreate"); + + // write.AddContainer("MSigmabar", "Events"); + + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MHadronness", "Events"); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + write.AddContainer("MMcRunHeader", "RunHeaders"); + + /* + write.AddContainer("MMcTrig", "Events"); + write.AddContainer("MRawEvtData", "Events"); + write.AddContainer("MRawEvtHeader", "Events"); + */ + + + MFillH fillh3a("MHHadronness"); + MFillH fillh3b("MHRanForest"); + + tlist3.AddToList(&fillh3a); + tlist3.AddToList(&fillh3b); + + + tlist3.AddToList(&write); + + evtloop.SetProgressBar(&bar); + + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist3.PrintStatistics(); + + plist.FindObject("MHRanForest")->DrawClone(); + + TCanvas *had = new TCanvas; + had->cd(); + plist.FindObject("MHHadronness")->Draw(); + + TString gifname = "had_S"; + gifname += minsize; + gifname += ".gif"; + + had->SaveAs(gifname); + had->Close(); + + plist.FindObject("MHHadronness")->DrawClone(); + plist.FindObject("MHHadronness")->Print();//*/ + + return; +} Index: /trunk/MagicSoft/Mars/mtemp/mpadova/macros/area.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpadova/macros/area.C (revision 5122) +++ /trunk/MagicSoft/Mars/mtemp/mpadova/macros/area.C (revision 5122) @@ -0,0 +1,113 @@ +void area(Char_t* infilename = "star_S500_all.root", Float_t maxhadronness = 0.1) +{ + // + // N_gamma_original: Number of generated Corsika showers in the MC sample + // from which the star file has been produced. VERY IMPORTANT! + + const Double_t N_gamma_original = 3.9e6; + const Double_t MC_diff_spectrum = -2.6; + const Double_t E_MC_min = 10.; + const Double_t E_MC_max = 30000.; + const Double_t sizecut = 0.; + const Double_t mincorepixels = 5; + const Double_t leakagemax = 1.; + + MHMcCollectionArea collarea; + + // + // Calculate approximately the original number of events in each + // energy bin: + // + + const Double_t expo = 1 + MC_diff_spectrum; + const Double_t k = N_gamma_original / + (pow(E_MC_max,expo) - pow(E_MC_min,expo)) ; + + TH2D* hall = collarea.GetHistAll(); + + const Int_t nbinx = hall.GetNbinsX(); + + TAxis &axe = *hall.GetXaxis(); + for (Int_t i = 1; i <= nbinx; i++) + { + const Float_t e1 = axe.GetBinLowEdge(i); + const Float_t e2 = axe.GetBinLowEdge(i+1); + + if (e1 < E_MC_min || e2 > E_MC_max) + continue; + + const Float_t events = k * (pow(e2, expo) - pow(e1, expo)); + // + // We fill the i-th energy bin, with the total number of events + // Second argument of Fill would be impact parameter of each + // event, but we don't really need it for the collection area, + // so we just put a dummy value (1.) + // + + const Float_t energy = (e1+e2)/2.; + hall.Fill(energy, 1., events); + } + + // + // Now open input file, loop over the Events tree and fill the energy + // histogram for events passing all required cuts. + // + + TChain chain("Events"); + chain.Add(infilename); + + MMcEvt *mcevt; + MHillas *mhillas; + MNewImagePar *mnew; + MHadronness *mhadronness; + + chain.SetBranchAddress("MMcEvt.", &mcevt); + chain.SetBranchAddress("MHillas.", &mhillas); + chain.SetBranchAddress("MNewImagePar.", &mnew); + chain.SetBranchAddress("MHadronness.", &mhadronness); + + chain.SetBranchStatus("MMcEvt.", 1); + chain.SetBranchStatus("MHillas.", 1); + chain.SetBranchStatus("MNewImagePar.", 1); + chain.SetBranchStatus("MHadronness.", 1); + + for (Int_t ievt = 0; ievt < chain.GetEntries(); ievt++) + { + chain.GetEvent(ievt); + + if (mcevt->GetPartId() > 1) + continue; + + if (mhillas->GetSize() < sizecut) + continue; + + if (mnew->GetNumCorePixels() < mincorepixels) + continue; + + if (mnew->GetLeakage1() > leakagemax) + continue; + + if (mhadronness->GetHadronness() > maxhadronness) + continue; + + const Float_t energy = mcevt->GetEnergy(); + + collarea.FillSel(energy,1.); + // Second argument is again just a dummy value + } + + + gStyle->SetOptStat(0); + gStyle->SetOptLogy(); + collarea.CalcEfficiency2(); + collarea.DrawClone(); + + // infile->Close(); + + TFile *out = new TFile("area.root", "recreate"); + collarea.Write(); + out->Close(); + + return; + +} Index: /trunk/MagicSoft/Mars/mtemp/mpadova/macros/gammarate.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpadova/macros/gammarate.C (revision 5122) +++ /trunk/MagicSoft/Mars/mtemp/mpadova/macros/gammarate.C (revision 5122) @@ -0,0 +1,69 @@ +void gammarate() +{ + + gROOT->Reset(); + c1 = new TCanvas("c1","Gamma rate",700,300,500,500); + c1->SetGridx(); + c1->SetGridy(); + c1->SetLogy(1); + c1->SetLogx(1); + + TFile f("area.root"); + MHMcCollectionArea* mhmc = (MHMcCollectionArea*)f.Get("MHMcCollectionArea"); + TH1D* h1 = (TH1D*)mhmc->GetHist(); + + + h1->GetXaxis(); + Float_t xmin = h1->GetXaxis()->GetXmin(); + Float_t xmax = h1->GetXaxis()->GetXmax(); + Int_t nbin = h1->GetNbinsX(); + + // gamma spectrum: + fun1 = new TF1("fun1","10000*2.706e-11*pow(x,-2.47-0.11*log(x))",xmin/1000.,xmax/1000.); + + TH1F *h3f = (TH1F*)h1->Clone("h3f"); + h3f->SetTitle("Integrale per bin"); + h3f->Reset(); + h3f->SetMinimum(1.e-13); + h3f->SetMaximum(1.e-7); + + TH1F *h4f = (TH1F*)h1->Clone("h4f"); + + h4f->SetTitle("Rate"); + h4f->Reset(); + h4f->SetMinimum(1.e-4); + h4f->SetMaximum(1.e-1); + + + Float_t Error = 0.; + + cout << "Bin\tArea\tCenter\tFlux\t\tMolt\tInt\tRate"<GetXaxis()->FindBin(i) + Float_t valArea = h1->GetBinContent(i); + Float_t bincenter = h1->GetBinCenter(i); + Float_t valFlusso = fun1->Eval(bincenter/1000., 0., 0.); + Float_t Integral = fun1->Integral(h1->GetBinLowEdge(i)/1000.,h1->GetBinLowEdge(i+1)/1000.); + h3f->SetBinContent(i,Integral); + Float_t Rate = Integral * valArea; + Float_t ErArea = h1->GetBinError(i); + Float_t ErRate = ErArea * Integral; + Error += ErRate*ErRate; + h4f->SetBinContent(i,Rate); + h4f->SetBinError(i,ErRate); + cout << i << "\t" << bincenter << "\t" << Rate << "\t"<< endl; + } + //h1->DrawCopy(); + //h3f->DrawCopy(); + h4f->SetStats(kFALSE); + h4f->SetTitle("Gamma rate vs energy"); + h4f->GetYaxis()->SetTitle("Rate"); + h4f->DrawCopy("pe"); + + cout << "Rate (Hz):" << h4f->Integral(0,100)<< " +- " <Delete(); + h3f->Delete(); + h1->Delete(); +}