source: trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C@ 7005

void computeCoeff()
{    
  MEffAreaAndCoeffCalc calc;
  calc.AddFile("/data/star_gamma_test.root");

  /************************************************/
  /* Build spectrum                               */
  /* now it's just a cross check with the same MC */
  /************************************************/
  const Int_t ebins = 20;              // number of bins to build spectrum
  const Double_t emin=10;
  const Double_t emax=200;
  const Int_t tbins = 2;
  const Double_t thetab[tbins+1] = {0,10,20};
  calc.SetEbins(ebins);
  calc.SetEmin(emin);
  calc.SetEmax(emax);
  calc.SetThetaBinning(tbins,thetab);
  
  // define the funtion of the desired spectrum
  calc.SetFunction("4.e9*pow(x,-2.6+1)");
  calc.ComputeAllFactors();
  
  const UInt_t ncutfiles=1;
  Char_t* cutName[ncutfiles]={"/data/star_gamma_test.root"};

  TChain* ccut = new TChain("Events");
  for(Int_t i = 0; i < ncutfiles; i++)
    ccut->Add(cutName[i]);
  ccut->SetAlias("logestenergy","log10(MHillas.fSize/0.18/15.)");
  ccut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");


  const Double_t logemin = TMath::Log10(emin);
  const Double_t logemax = TMath::Log10(emax);
  TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax);
  hspec->Sumw2();
  ccut->Draw("logestenergy>>hspec","theta<10","goff");

  for(Int_t i=0;i<ebins;i++)
    {
      const Float_t uncorrval = hspec->GetBinContent(i+1);
      const Float_t effa  =  calc.GetEffectiveAreaHisto()->GetBinContent(i+1,1);
      const Float_t unfold = calc.GetCoefficientHisto()->GetBinContent(i+1,1);
      Float_t corrval;
      if(effa)
        corrval = uncorrval*unfold/effa*1e9;
      else
        corrval = 0;
      hspec->SetBinContent(i+1,corrval);
    }

  // SAVE RESULTS
  TFile file("prueba.root","RECREATE");
  hspec->Write();
  calc.GetEffectiveAreaHisto()->Write();
  calc.GetCoefficientHisto()->Write();
  return;
}