Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 6465) +++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 6486) @@ -19,5 +19,11 @@ -*-*- END OF LINE -*-*- - 2005/02/10 Javier Rico + 2005/02/15 Javier Rico + * library/MEffAreaAndCoeffCalc.[h,cc] + * macros/computeCoeff.C + - extend the energy bounds to compute weights + - add possibility to set the zenith angle binning from outside + + 2005/02/14 Javier Rico * library/MEffAreaAndCoeffCalc.[h,cc] * macros/computeCoeff.C Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (revision 6465) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (revision 6486) @@ -47,7 +47,11 @@ using namespace std; -const Int_t fNTbins = 2; -const Double_t fTbin[fNTbins+1] = {0,10,20}; - +const Int_t ntbins=1; // default number of theta bins +const Double_t tbin[ntbins+1] = {0,90}; // default theta bins bounds +const Int_t nebins = 10; // default number of energy bins +const Float_t emin = 10.; // default minimum energy value +const Float_t emax = 10000.; // default maximum energy value +const Int_t nsubbins = 20; // default number of subbins per bin +const Char_t* deff = "4.e9*pow(x,-2.6+1)"; // default spectrum function // ------------------------------------------------------------------------- @@ -56,12 +60,11 @@ // MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc() - : fSpec(NULL), fHorig(NULL), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20), - fCoeff(NULL), fEffA(NULL), fFile(NULL) + : fSpec(NULL), fEmin(emin), fEmax(emax), fEbins(nebins), fEsubbins(nsubbins), + fWeight(NULL), fCoeff(NULL), fEffA(NULL)//, fFile(NULL) { // set the default function - SetFunction("4.e9*pow(x,-2.6+1)"); + SetFunction(deff); // create the TChains - // OJO!! make sure they read the appropriate branch fCini = new TChain("OriginalMC"); fCcut = new TChain("Events"); @@ -76,7 +79,10 @@ fCcut->SetBranchAddress("MHillas.",&fHillas); fCcut->SetBranchAddress("MMcEvt.",&fMcEvt); + + // initial value of the zenith angle binning + SetThetaBinning(ntbins,tbin); // borra - fFile = new TFile("test.root","RECREATE"); + // fFile = new TFile("coeftest.root","RECREATE"); } @@ -90,4 +96,7 @@ delete fSpec; + if(fTbin) + delete [] fTbin; + if(fWeight) delete [] fWeight; @@ -99,9 +108,20 @@ delete fEffA; - if(fHorig) - delete fHorig; delete fCini; delete fCcut; +} + +// ------------------------------------------------------------------------- +// +// Set the binning in zenith angle +// +void MEffAreaAndCoeffCalc::SetThetaBinning(Int_t n, const Double_t* binlist) +{ + fNTbins=n; + if(fTbin) delete [] fTbin; + fTbin = new Double_t[n+1]; + for(Int_t i=0;iDraw("logenergy>>fhorig","","goff"); - // borra - fHorig->Write(); } @@ -161,37 +172,54 @@ void MEffAreaAndCoeffCalc::ComputeWeights() { - const Double_t logemin = TMath::Log10(fEmin); - const Double_t logemax = TMath::Log10(fEmax); - const Double_t de = (logemax-logemin)/fEbins; // bin size (in log) - const Double_t desub = de/fEsubbins; // subbin size (in log) - const Int_t esbins = fEbins*fEsubbins; // total number of subbins - + // OJO!! maybe this should be hard-coded somewhere else + const Float_t abslogmin=0.; + const Float_t abslogmax=5.; + + const Float_t logemin = TMath::Log10(fEmin); + const Float_t logemax = TMath::Log10(fEmax); + const Float_t de = (logemax-logemin)/fEbins; // bin size (in log) + const Float_t desub = de/fEsubbins; // subbin size (in log) + const Int_t esbins = fEbins*fEsubbins; // total number of subbins + + // compute the binning for weights histogram + const Int_t nmindist = (logemin>abslogmin)? Int_t((logemin-abslogmin)/desub) : 0; + const Int_t nmaxdist = (logemaxGetBinContent(i+1); // number of events in initial spectrum - const Double_t ew = TMath::Power(10,logemin+(i+0.5)*desub); // real energy - const Double_t numer = fSpec->Eval(ew); // number of events for the required spectrum - if(denom) + + TH1F* horigs = new TH1F("horigs","Original energy spectrum",fEWbins,fLogEWmin,fLogEWmax); + fCini->Draw("logenergy>>horigs","","goff"); + // borra + // horigs->Write(); + + // borra + // TH1F hw("hw","hw",fEWbins,fLogEWmin,fLogEWmax); + for(Int_t i=0;iGetBinContent(i+1); // number of events in initial spectrum + const Float_t ew = TMath::Power(10,fLogEWmin+(i+0.5)*desub); // real energy + const Float_t numer = fSpec->Eval(ew); // number of events for the required spectrum + if(denom>10) fWeight[i]=numer/denom; else { - cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 " << endl; + // cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 " << endl; fWeight[i]=-1; } // borra - hw.SetBinContent(i+1,fWeight[i]); - } - // borra - hw.Write(); + // hw.SetBinContent(i+1,fWeight[i]); + } + // borra + // hw.Write(); + + delete horigs; } @@ -208,9 +236,8 @@ } - const Double_t logemin = TMath::Log10(fEmin); - const Double_t logemax = TMath::Log10(fEmax); - const Double_t de = (logemax-logemin)/fEbins; // bin size (in log) - const Double_t desub = de/fEsubbins; // subbin size (in log) - const Int_t esbins = fEbins*fEsubbins; // total number of subbins + const Float_t logemin = TMath::Log10(fEmin); + const Float_t logemax = TMath::Log10(fEmax); + const Float_t de = (logemax-logemin)/fEbins; // bin size (in log) + const Float_t desub = de/fEsubbins; // subbin size (in log) const Int_t nentries = Int_t(fCcut->GetEntries()); @@ -220,7 +247,7 @@ // borra - TH1F* hest1 = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax); - TH1F* hmc1 = new TH1F("hmc1","MC energy",fEbins,logemin,logemax); - TH1F* coef1 = new TH1F("coef1","coefficients",fEbins,logemin,logemax); + // TH1F* hest1 = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax); + // TH1F* hmc1 = new TH1F("hmc1","MC energy",fEbins,logemin,logemax); + // TH1F* coef1 = new TH1F("coef1","coefficients",fEbins,logemin,logemax); // reset the coefficients @@ -238,23 +265,33 @@ Float_t emc = fMcEvt->GetEnergy(); Float_t estim = fHillas->GetSize()/0.18/15.; + + if((emcfEmax && estim>fEmax) || + (emcfEmax) || (emc>fEmax && estimGetTheta()*180./3.14159; // zenith angle (in degrees) // compute the bin where the weight is taken from - Int_t wbin = Int_t((TMath::Log10(emc)-logemin)/desub); + Int_t wbin = Int_t((TMath::Log10(emc)-fLogEWmin)/desub); // fill the histograms - if(wbinFill(TMath::Log10(estim),theta,fWeight[wbin]); - hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]); - // borra - if(theta0) { - hest1->Fill(TMath::Log10(estim),fWeight[wbin]); - hmc1->Fill(TMath::Log10(emc),fWeight[wbin]); + hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]); + hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]); + // borra +// if(thetaFill(TMath::Log10(estim),fWeight[wbin]); +// hmc1->Fill(TMath::Log10(emc),fWeight[wbin]); +// } } + else + cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has no computed weight (lack of MC statistics), skipping" << endl; } else - cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl; + cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl; } @@ -266,12 +303,12 @@ const Float_t den = hest->GetBinContent(i+1,j+1); //borra - const Float_t num1 = hmc1->GetBinContent(i+1); - const Float_t den1 = hest1->GetBinContent(i+1); +// const Float_t num1 = hmc1->GetBinContent(i+1); +// const Float_t den1 = hest1->GetBinContent(i+1); if(den) { fCoeff->SetBinContent(i+1,j+1,num/den); //borra - if(j==0 && den1) - coef1->SetBinContent(i+1,num1/den1); +// if(j==0 && den1) +// coef1->SetBinContent(i+1,num1/den1); } else @@ -284,7 +321,7 @@ //borra - hmc1->Write(); - hest1->Write(); - coef1->Write(); +// hmc1->Write(); +// hest1->Write(); +// coef1->Write(); delete hmc; @@ -305,12 +342,11 @@ //OJO!! radius should be read from somewhere! - const Double_t radius = 30000.; // generation radius (in cm) - // const Double_t Atot = 3.14159*radius*radius; //total area (in cm^2) - const Double_t Atot = 1.; //total area (in cm^2) - const Double_t logemin = TMath::Log10(fEmin); - const Double_t logemax = TMath::Log10(fEmax); + const Float_t radius = 30000.; // generation radius (in cm) + const Float_t Atot = 3.14159*radius*radius; //total area (in cm^2) + const Float_t logemin = TMath::Log10(fEmin); + const Float_t logemax = TMath::Log10(fEmax); const Int_t esbins = fEbins*fEsubbins; // total number of subbins - const Double_t de = (logemax-logemin)/fEbins; // bin size (in log) - const Double_t desub = de/fEsubbins; // subbin size (in log) + const Float_t de = (logemax-logemin)/fEbins; // bin size (in log) + const Float_t desub = de/fEsubbins; // subbin size (in log) // reset the effective areas @@ -319,5 +355,5 @@ fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin); //borra - TH1F eff("eff","Effective area",fEbins,logemin,logemax); +// TH1F eff("eff","Effective area",fEbins,logemin,logemax); // fill the spectrum of the survivors @@ -343,6 +379,6 @@ else { - const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub); - const Double_t ww = fSpec->Eval(ew); + const Float_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub); + const Float_t ww = fSpec->Eval(ew); effarea+=Atot*numerator/denominator*ww; wtot += ww; @@ -357,14 +393,15 @@ { fEffA->SetBinContent(i+1,j+1,effarea/wtot); - if(j==0) - { - cout << "*****" << i << ": " << effarea/wtot << endl; - eff.SetBinContent(i+1,effarea/wtot); - } + // borra +// if(j==0) +// { +// // cout << "*****" << i << ": " << effarea/wtot << endl; +// eff.SetBinContent(i+1,effarea/wtot); +// } } } // borra - eff.Write(); +// eff.Write(); delete hpass; Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h (revision 6465) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h (revision 6486) @@ -19,10 +19,16 @@ TF1* fSpec; // function used to parametrize the spectrum - TH1F* fHorig; // histogram with the original sample energy spectrum - Double_t fEmin; // Minimum energy in GeV - Double_t fEmax; // Maximum energy in GeV + Float_t fEmin; // Minimum energy in GeV + Float_t fEmax; // Maximum energy in GeV Int_t fEbins; // number of bins to build spectrum Int_t fEsubbins; // number of subbins per big bin (to compute weights, eff areas...) + + Float_t fLogEWmin; // Log Minimum energy for weights (in GeV) + Float_t fLogEWmax; // Log Maximum energy for weights (in GeV) + Int_t fEWbins; // Number of bins for weights + + Int_t fNTbins; // Number of bins in zenith angle + Double_t* fTbin; // array containing bin boundaries (size must be fNTbins+) Double_t* fWeight; // array containing weights @@ -58,4 +64,6 @@ void SetEmax(Float_t x) {fEmax=x;} + void SetThetaBinning(Int_t n, const Double_t* binlist); + void AddFile(const Char_t* name) {fCini->Add(name); fCcut->Add(name);} Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C (revision 6465) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C (revision 6486) @@ -10,14 +10,16 @@ /* now it's just a cross check with the same MC */ /************************************************/ - const Int_t ebins = 10; // number of bins to build spectrum + const Int_t ebins = 20; // number of bins to build spectrum const Double_t emin=10; - const Double_t emax=600; - const Double_t logemin = TMath::Log10(emin); - const Double_t logemax = TMath::Log10(emax); + const Double_t emax=200; + const Int_t tbins = 2; + const Double_t thetab[tbins+1] = {0,10,20}; calc.SetEbins(ebins); - // calc.SetEsubbins(1); - + 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)",emin,emax); + calc.SetFunction("4.e9*pow(x,-2.6+1)"); calc.ComputeAllFactors(); @@ -31,4 +33,7 @@ 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(); @@ -42,5 +47,5 @@ Float_t corrval; if(effa) - corrval = uncorrval*unfold/effa; + corrval = uncorrval*unfold/effa*1e9; else corrval = 0;