source: trunk/MagicSoft/Mars/mtemp/mpadova/macros/area.C@ 6068

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;

}