Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 6464) +++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 6465) @@ -19,4 +19,11 @@ -*-*- END OF LINE -*-*- + 2005/02/10 Javier Rico + * library/MEffAreaAndCoeffCalc.[h,cc] + * macros/computeCoeff.C + - compute the coefficients in the big bins + - migrate to new data format + + 2005/02/10 Javier Rico * library/MEffAreaAndCoeffCalc.[h,cc] Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (revision 6464) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (revision 6465) @@ -38,4 +38,5 @@ #include "MMcEvt.hxx" #include "TH2F.h" +#include "TFile.h" #include "MLog.h" @@ -55,29 +56,27 @@ // MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc() - : fSpec(NULL), fHorig(NULL), fEbins(10), fEsubbins(20), fEmin(10.), fEmax(10000.), - fCoeff(NULL), fEffA(NULL) + : fSpec(NULL), fHorig(NULL), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20), + fCoeff(NULL), fEffA(NULL), fFile(NULL) { // set the default function - // OJO!! this should be changed SetFunction("4.e9*pow(x,-2.6+1)"); // create the TChains // OJO!! make sure they read the appropriate branch - fCini = new TChain("Events"); + fCini = new TChain("OriginalMC"); fCcut = new TChain("Events"); - // define some usefull aliases - fCini->SetAlias("logenergy","log10(MMcEvt.fEnergy)"); - fCini->SetAlias("theta","MMcEvt.fTheta*180./3.14159"); - - // OJO!! the estimated energy must be changed for a real one + // define some useful aliases + fCini->SetAlias("logenergy","log10(MMcEvtBasic.fEnergy)"); + fCini->SetAlias("theta","MMcEvtBasic.fTelescopeTheta*180./3.14159"); + fCcut->SetAlias("logenergy","log10(MMcEvt.fEnergy)"); - fCcut->SetAlias("theta","MMcEvt.fTheta*180./3.14159"); - - fCcut->SetBranchStatus("*",0); - fCcut->SetBranchStatus("MHillas.*",1); - fCcut->SetBranchStatus("MMcEvt.*",1); + fCcut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159"); + fCcut->SetBranchAddress("MHillas.",&fHillas); fCcut->SetBranchAddress("MMcEvt.",&fMcEvt); + + // borra + fFile = new TFile("test.root","RECREATE"); } @@ -119,4 +118,9 @@ emin = fEmin; emax = fEmax; + } + else + { + fEmin = emin; + fEmax = emax; } fSpec = new TF1("fspec",chfunc,emin,emax); @@ -147,4 +151,6 @@ fHorig = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax); fCini->Draw("logenergy>>fhorig","","goff"); + // borra + fHorig->Write(); } @@ -169,4 +175,6 @@ FillOriginalSpectrum(); + // borra + TH1F hw("hw","hw",esbins,fEmin,fEmax); for(Int_t i=0;iGetEnergy(); - Float_t estim = fHillas->GetSize()/15.; + Float_t estim = fHillas->GetSize()/0.18/15.; Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees) @@ -231,7 +248,13 @@ 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 energy out of bounds, skpping" << endl; + cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl; } @@ -242,6 +265,14 @@ const Float_t num = hmc->GetBinContent(i+1,j+1); 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); if(den) - fCoeff->SetBinContent(i+1,j+1,num/den); + { + fCoeff->SetBinContent(i+1,j+1,num/den); + //borra + if(j==0 && den1) + coef1->SetBinContent(i+1,num1/den1); + } else { @@ -252,4 +283,9 @@ + //borra + hmc1->Write(); + hest1->Write(); + coef1->Write(); + delete hmc; delete hest; @@ -270,14 +306,19 @@ //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 = 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 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) // reset the effective areas if(fEffA) delete fEffA; - fEffA = new TH2F("feffa","Effective area",esbins,logemin,logemax,fNTbins,fTbin); - + fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin); + //borra + TH1F eff("eff","Effective area",fEbins,logemin,logemax); + // fill the spectrum of the survivors TH2F* hpass= new TH2F("hpass","Survivors",esbins,logemin,logemax,fNTbins,fTbin); @@ -289,16 +330,41 @@ // compute the effective areas for(Int_t j=0;jGetBinContent(i+1,j+1)<=0) + Float_t effarea =0; + Float_t wtot = 0; + for(Int_t k=0;kGetBinContent(i*fEsubbins+k+1,j+1); + Float_t denominator = horig->GetBinContent(i*fEsubbins+k+1,j+1); + + if(denominator<=0) + cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy subbin " << i*fEsubbins+k <<", theta bin " << j << " contains no events" << endl; + else + { + const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub); + const Double_t ww = fSpec->Eval(ew); + effarea+=Atot*numerator/denominator*ww; + wtot += ww; + } + } + if(!wtot) { cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i <<", theta bin " << j << " contains no events setting effective area to 0" << endl; fEffA->SetBinContent(i+1,j+1,0); + } + else + { + fEffA->SetBinContent(i+1,j+1,effarea/wtot); + if(j==0) + { + cout << "*****" << i << ": " << effarea/wtot << endl; + eff.SetBinContent(i+1,effarea/wtot); + } } - else - fEffA->SetBinContent(i+1,j+1,Atot*hpass->GetBinContent(i+1,j+1)/horig->GetBinContent(i+1,j+1)); } - - + + // borra + eff.Write(); delete hpass; Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h (revision 6464) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h (revision 6465) @@ -18,11 +18,11 @@ private: - TF1* fSpec; // function used to parametrized the spectrum - TH1F* fHorig; // histogram with the original energy spectrum + TF1* fSpec; // function used to parametrize the spectrum + TH1F* fHorig; // histogram with the original sample energy spectrum - Int_t fEbins; // number of bins to build spectrum - Int_t fEsubbins; // number of subbins per big bin Double_t fEmin; // Minimum energy in GeV Double_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...) Double_t* fWeight; // array containing weights @@ -30,9 +30,11 @@ TH2F* fEffA; // histogram containing effective areas - TChain* fCini; // chain for initial MC events (even those not triggering) + TChain* fCini; // chain for initial MC files (before trigger) TChain* fCcut; // chain for surviving MC events (after cuts) MHillas* fHillas; // pointer to the MHillas Branch MMcEvt* fMcEvt; // pointer to the MMcEvt Branch + + TFile* fFile; // output file (for debugging only) protected: @@ -44,4 +46,5 @@ public: + MEffAreaAndCoeffCalc(); @@ -55,6 +58,5 @@ void SetEmax(Float_t x) {fEmax=x;} - void AddFileToInitialMC(const Char_t* name) {fCini->Add(name);} - void AddFileToFinalMC(const Char_t* name) {fCcut->Add(name);} + void AddFile(const Char_t* name) {fCini->Add(name); fCcut->Add(name);} TH2F* GetEffectiveAreaHisto() {return fEffA;} Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C (revision 6464) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C (revision 6465) @@ -4,27 +4,6 @@ { MEffAreaAndCoeffCalc calc; + calc.AddFile("/data/star_gamma_test.root"); - // initial MC files - calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1000to1009_w0.root"); - calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1130to1139_w0.root"); - calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1520to1529_w0.root"); - calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1260to1269_w0.root"); - calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1390to1399_w0.root"); - calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1650to1659_w0.root"); - - - // files with remaining events after cuts (or where to apply cuts) see warning above - calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root"); - calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root"); - calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root"); - calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root"); - calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root"); - calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root"); - - - // define the funtion of the desired spectrum - calc.SetFunction("4.e9*pow(x,-2.6+1)",10.,10000.); - calc.ComputeAllFactors(); - /************************************************/ /* Build spectrum */ @@ -32,57 +11,40 @@ /************************************************/ const Int_t ebins = 10; // number of bins to build spectrum - const Int_t esubbins = 20; // number of subbins per big bin const Double_t emin=10; - const Double_t emax=1000; + const Double_t emax=600; const Double_t logemin = TMath::Log10(emin); const Double_t logemax = TMath::Log10(emax); - const Double_t de = (logemax-logemin)/ebins; // bin size (in log) - const Double_t desub = de/esubbins; // subbin size (in log) - const Int_t esbins = ebins*esubbins; // total number of subbins + calc.SetEbins(ebins); + // calc.SetEsubbins(1); - // remaining events after cuts - const UInt_t ncutfiles=6; - Char_t* cutName[ncutfiles]={ - "/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root", - "/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root", - "/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root", - "/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root", - "/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root", - "/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root" - }; + // define the funtion of the desired spectrum + calc.SetFunction("4.e9*pow(x,-2.6+1)",emin,emax); + 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/15.)"); - const Int_t nentries = Int_t(ccut->GetEntries()); - - MHillas* hillas; - MMcEvt* mcevt; - ccut->SetBranchStatus("*",0); - ccut->SetBranchStatus("MHillas.*",1); - ccut->SetBranchStatus("MMcEvt.*",1); - ccut->SetBranchAddress("MHillas.",&hillas); - ccut->SetBranchAddress("MMcEvt.",&mcevt); + ccut->SetAlias("logestenergy","log10(MHillas.fSize/0.18/15.)"); + ccut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159"); TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax); hspec->Sumw2(); + ccut->Draw("logestenergy>>hspec","theta<10","goff"); - for(Int_t i=0;iGetEntry(i); - - // OJO!! Estimated energy will be taken directly from some input container - Float_t estim = hillas->GetSize()/15.; - - UInt_t effabin = UInt_t((TMath::Log10(estim)-logemin)/desub); - UInt_t coefbin = UInt_t((TMath::Log10(estim)-logemin)/de); - - Float_t effa = calc.GetEffectiveAreaHisto()->GetBinContent(effabin+1,1); - Float_t unfold = calc.GetCoefficientHisto()->GetBinContent(coefbin+1,1); - + 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) - hspec->Fill(TMath::Log10(estim),unfold/effa*1e9); - } + corrval = uncorrval*unfold/effa; + else + corrval = 0; + hspec->SetBinContent(i+1,corrval); + } // SAVE RESULTS