Index: trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/Changelog (revision 4136) +++ trunk/MagicSoft/Mars/mtemp/Changelog (revision 4137) @@ -18,4 +18,10 @@ -*-*- END OF LINE -*-*- + + 2004/05/23: Thomas Hengstebeck + * mberlin/macros/ScaleHillas.C , SrcCorrect.C , Convert2Matrix.C , + CutOptim.C , AlphaPlot.C , IOMkn421.h + - Added analysis macros + 2004/05/18: Daniel Mazin Index: trunk/MagicSoft/Mars/mtemp/mberlin/macros/AlphaPlot.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mberlin/macros/AlphaPlot.C (revision 4137) +++ trunk/MagicSoft/Mars/mtemp/mberlin/macros/AlphaPlot.C (revision 4137) @@ -0,0 +1,250 @@ +#include "TString.h" +#include "TChain.h" +#include "TCanvas.h" +#include "TH1F.h" +#include "TF1.h" +#include "TLegend.h" +#include "TArrayF.h" +#include "TGraphErrors.h" +#include "MRawRunHeader.h" +#include "TPad.h" +#include "TPaveText.h" +#include "iostream.h" +#include "TObjArray.h" +#include "TFile.h" + +//-------------------------------------------------------------------------------------- +// MAIN INPUT/OUTPUT + +#include "IOMkn421.h" +//-------------------------------------------------------------------------------------- + +Double_t signal(Double_t *x, Double_t *par) +{ + return par[0]/(TMath::Sqrt(2.*TMath::Pi())*par[1])*TMath::Exp(-(x[0]-par[2])*(x[0]-par[2])/2./par[1]/par[1]); +} + +Double_t background(Double_t *x, Double_t *par) +{ + return par[0]+par[1]*x[0]+par[2]*x[0]*x[0]+par[3]*x[0]*x[0]*x[0]; +} + +Double_t fitFunction(Double_t *x, Double_t *par) +{ + return signal(x,par) + background(x,&par[3]); +} + +void AlphaPlot() +{ + DistLo*=deg2mm; + DistUp*=deg2mm; + + Double_t acut = 10; // alpha region of excess + Int_t nbins = 36; + Int_t nbin_flat= nbins/3; + + // input file chain + TChain *fEvents=new TChain("Events"); + fEvents->Add(fileOn1.Data()); + fEvents->Add(fileOn2.Data()); + + TChain *fEventsOff=new TChain("Events"); + fEventsOff->Add(fileOff1.Data()); + fEventsOff->Add(fileOff2.Data()); + + + cout<<"Non="<GetEntries()<<" Noff="<GetEntries()<SetFillColor(0); + c1->SetBorderMode(0); + c1->cd(); + + TH1F *halpha=new TH1F("halpha","",nbins,0.,90.); + halpha->GetXaxis()->SetTitle("Alpha [deg]"); + halpha->GetYaxis()->SetTitle("Counts"); + + TH1F *halphaOff=new TH1F("halphaOff","",nbins,0.,90.); + halpha->SetStats(kFALSE); + halphaOff->SetStats(kFALSE); + + TString strcut=Form("MHillas.fSize>%f",SizeLo); + + strcut+=Form("&& MHillasSrc.fDist>%f",DistLo); + strcut+=Form("&& MHillasSrc.fDist<%f",DistUp); + + strcut+=Form("&& MHillas.fWidth>%f",WidthLo); + strcut+=Form("&& MHillas.fWidth<%f",WidthUp); + + strcut+=Form("&& MHillas.fLength>%f",LengthLo); + strcut+=Form("&& MHillas.fLength<%f",LengthUp); + + // fill histograms + fEvents->Project("halpha","TMath::Abs(MHillasSrc.fAlpha)",strcut.Data()); + halpha->SetLineWidth(2); + halpha->DrawCopy("e"); + + fEventsOff->Project("halphaOff","TMath::Abs(MHillasSrc.fAlpha)",strcut.Data()); + halphaOff->SetLineColor(2); + halphaOff->SetLineWidth(2); + + + // scale Off histogram to On data using bins nbin_flat(first bin after peak) - nbins + Double_t nsample_on=0; + Double_t nsample_off=0; + for(Int_t i=nbin_flat;iGetBinContent(i); + for(Int_t i=nbin_flat;iGetBinContent(i); + + Double_t scal=nsample_on/nsample_off; + halphaOff->Sumw2(); + halphaOff->Scale(scal); + halphaOff->Draw("e same"); + + gPad->SetGridx(); + gPad->SetGridy(); + + //-------------------------------------------------------------------------- + // fit alpha hist + TF1 *backFcn = new TF1("backFcn",background,20.,89.5,4); + backFcn->SetLineWidth(2); + backFcn->SetLineColor(kRed); + + TF1 *signalFcn = new TF1("signalFcn",signal,0.,89.5,3); + signalFcn->SetLineWidth(2); + signalFcn->SetLineColor(kBlue); + signalFcn->SetNpx(500); + + TF1 *fitFcn = new TF1("fitFcn",fitFunction,0.,89.5,7); + fitFcn->SetNpx(500); + fitFcn->SetLineWidth(2); + fitFcn->SetLineColor(kMagenta); + + // start-parameters + Double_t p3start=halpha->GetBinContent(halpha->GetNbinsX()-1); + Double_t p4start=0.; + Double_t p5start=0.; + Double_t p6start=0.; + + Double_t p0start=(halpha->GetEntries()-p3start*nbins)*halpha->GetBinWidth(1); + Double_t p1start=5.; + Double_t p2start=0.; + + cout<<"Start values for fit:"<GetParameters(parb); + + // total fit + Double_t par[7]; + fitFcn->SetParameters(p0start,p1start,p2start, + parb[0],parb[1],parb[2],parb[3]); + + fitFcn->FixParameter(2,0.); + fitFcn->FixParameter(3,parb[0]); + fitFcn->FixParameter(4,parb[1]); + fitFcn->FixParameter(5,parb[2]); + fitFcn->FixParameter(6,parb[3]); + halpha->Fit ("fitFcn","RN"); + + // draw fit results + legend + TLegend *legend=new TLegend(0.4,0.5,0.93,0.65); + legend->SetTextFont(40); + legend->SetTextSize(0.04); + legend->SetFillColor(19); + legend->AddEntry(halpha,"Data","lp"); + TString strpol(""); + strpol=Form("Polynom"); + legend->AddEntry(backFcn,strpol.Data(),"l"); + legend->AddEntry(fitFcn,"Global Fit","l"); + //legend->Draw(); + + fitFcn->GetParameters(par); + fitFcn->SetRange(-90.,90.); + fitFcn->SetLineWidth(2); + fitFcn->Draw("same"); + + signalFcn->SetParameters(par); + //signalFcn->SetLineWidth(2); + //signalFcn->Draw("same"); + + backFcn->SetParameters(&par[3]); + backFcn->SetRange(-90,90); + backFcn->SetLineWidth(2); + backFcn->Draw("same"); + + //----------------------------------------------------------------------- + // calculate significance + + // analytical method + Double_t nSig=signalFcn->Integral(0.,acut)/halpha->GetBinWidth(1); + Double_t nOffScaled=backFcn->Integral(0.,acut)/halpha->GetBinWidth(1); + + Double_t nOn=nSig+nOffScaled; + + // significance after Li/Ma + Double_t theta=1.; // should be = Ton/Toff ??? + Double_t nOff=nOffScaled/theta; + Double_t siglima = sqrt(2*(nOn*TMath::Log((1+theta)*nOn/(theta*(nOn+nOff))) + +nOff*TMath::Log((1+theta)*(nOff/(nOn+nOff))))); + + // "standard" significance + Double_t sigstd = nSig/sqrt(nOn+nOffScaled); + + TPaveText *pt = new TPaveText(0.25,0.75,0.56,0.98,"brNDC"); + pt->SetTextFont(60); + pt->SetTextSize(0.03); + pt->SetTextAlign(12); + + TString strcut1(Form("cuts:")); + TString strcut3(Form(" |alpha| < %4.2f ",acut)); + pt->AddText(strcut3.Data()); + TString str(""); + pt->AddText(str); + + TString strsig(Form("NExcess = %d",Int_t(nSig))); + pt->AddText(strsig.Data()); + TString strbak(Form("NOff = %d",Int_t(nOffScaled))); + pt->AddText(strbak.Data()); + pt->AddText(str); + + TString strsg(Form("Significance (Li/Ma) = %4.2f",siglima)); + pt->AddText(strsg.Data()); + strsg=Form("Significance (std) = %4.2f",sigstd); + pt->AddText(strsg.Data()); + pt->Draw(); + + cout<GetNrows();i++) + { + if( (*matOn)(i,isize)d1) ||((*matOn)(i,idist)w1) ||((*matOn)(i,iwidth)l1)||((*matOn)(i,ilength)GetNrows();i++) + { + if( (*matOff)(i,isize)d1) ||((*matOff)(i,idist)w1) ||((*matOff)(i,iwidth)l1)||((*matOff)(i,ilength)GetNrows(); + Int_t nrowOff = matOff->GetNrows(); + + Int_t ncolOn = matOn->GetNcols(); + Int_t ncolOff = matOff->GetNcols(); + + cout<<"matrix On: ncol="<SetFillColor(0); + + //----------------------------------------------------------------------- + // optimize cuts + + Int_t iflag; + Double_t val; + Double_t vstart[6]; + + Double_t *gin=NULL; + + Double_t arglist[10]; + Int_t ierflg = 0; + + TMinuit *gMinuit = new TMinuit(10); + TGraph* gr=NULL; + + if(!aplot) + c->Divide(2,1); + else + goto APlot; + + // start values + vstart[0]=WidthUp; + vstart[1]=WidthLo; + vstart[2]=LengthUp; + vstart[3]=LengthLo; + vstart[4]=DistUp; + vstart[5]=DistLo; + + gMinuit->SetFCN(fcn); + + arglist[0] = 1; + gMinuit->mnexcm("SET ERR", arglist ,1,ierflg); + + // Set starting values and step sizes for parameters + static Double_t step[6]; + const Double_t def_step=0.2; + step[0]=def_step; + step[1]=def_step; + step[2]=def_step; + step[3]=def_step; + step[4]=def_step; + step[5]=def_step; + + gMinuit->mnparm(0,"w1", vstart[0], step[0], 0.,0.,ierflg); + gMinuit->mnparm(1,"w2", vstart[1], step[1], 0.,0.,ierflg); + gMinuit->mnparm(2,"l1", vstart[2], step[2], 0.,0.,ierflg); + gMinuit->mnparm(3,"l2", vstart[3], step[3], 0.,0.,ierflg); + gMinuit->mnparm(4,"d1", vstart[4], step[4], 0.2,1.,ierflg); + gMinuit->mnparm(5,"d2", vstart[5], step[5], 0.5,1.5,ierflg); + + //gMinuit->FixParameter(0); // WidthUp + gMinuit->FixParameter(1); // WidthLo + //gMinuit->FixParameter(2); // LengthUp + gMinuit->FixParameter(3); // LengthLo + gMinuit->FixParameter(4); // DistUp + gMinuit->FixParameter(5); // DistLo + + gMinuit->Eval(6,gin,val,vstart,iflag); + + arglist[0] = 5000; + arglist[1] = 1.; + + //gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg); + gMinuit->mnexcm("SIM", arglist ,2,ierflg); + + // Print results + Double_t amin,edm,errdef; + Int_t nvpar,nparx,icstat; + gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat); + gMinuit->mnprin(3,amin); + + Double_t fpar[6]; + Double_t fparerr[6]; + for(Int_t i=0;i<6;i++) + gMinuit->GetParameter(i,fpar[i],fparerr[i]); + + WidthUp = fpar[0]; + WidthLo = fpar[1]; + LengthUp = fpar[2]; + LengthLo = fpar[3]; + DistUp = fpar[4]; + DistLo = fpar[5]; + + + // parameter's sigma plots + c->cd(2); + // 2 sigma contour + /* + gMinuit->SetErrorDef(4); + gr = (TGraph*) gMinuit->Contour(100,0,2); + gr->DrawClone("AL");*/ + + // 1 sigma contour + gMinuit->SetErrorDef(1); + gr = (TGraph*) gMinuit->Contour(100,0,2); + gr->DrawClone("AL"); + /* + // contour of min. function + c->cd(2); + arglist[0] = 1; // WidthUp + gMinuit->mnexcm("SCA", arglist ,1,ierflg); + TGraph* gr = (TGraph*)gMinuit->GetPlot(); + gr->DrawClone("al"); + + c->cd(4); + arglist[0] = 4; // LengthLo + gMinuit->mnexcm("SCA", arglist ,1,ierflg); + gr = (TGraph*)gMinuit->GetPlot(); + gr->DrawClone("al");//*/ + + c->cd(1); + +APlot: + //----------------------------------------------------------------------- + // plot alpha histogram with current cuts and calc significance + + // alpha histogram + TH1F *halpha=new TH1F("halpha","",18,0.,90.); + TH1F *halphaOff=new TH1F("halphaOff","",18,0.,90.); + + halpha->SetStats(kFALSE); + halphaOff->SetStats(kFALSE); + + for(Int_t i=0;iDistUp)) continue; + + if((mOn(i,iwidth) < WidthLo) ||(mOn(i,iwidth) > WidthUp)) continue; + if((mOn(i,ilength)< LengthLo)||(mOn(i,ilength)> LengthUp)) continue; + + halpha->Fill(TMath::Abs(mOn(i,ialpha))); + } + halpha->SetLineWidth(2); + halpha->DrawClone(); + + for(Int_t i=0;iDistUp)) continue; + + if((mOff(i,iwidth) < WidthLo) ||(mOff(i,iwidth) > WidthUp)) continue; + if((mOff(i,ilength)< LengthLo)||(mOff(i,ilength)> LengthUp)) continue; + + halphaOff->Fill(TMath::Abs(mOff(i,ialpha))); + + } + + // scale Off histogram to On data using bins 6-18 + Double_t nsample_on=0; + Double_t nsample_off=0; + for(Int_t i=6;i<18;i++) + nsample_on+=halpha->GetBinContent(i); + for(Int_t i=6;i<18;i++) + nsample_off+=halphaOff->GetBinContent(i); + + Double_t scal=nsample_on/nsample_off; + + halphaOff->SetLineColor(2); + halphaOff->Scale(scal); + halphaOff->Draw("same"); + + gPad->SetGridx(); + gPad->SetGridy(); + + // define functions + TF1 *backFcn = new TF1("backFcn",background,20.,89.5,4); + backFcn->SetLineWidth(2); + backFcn->SetLineColor(kRed); + + TF1 *signalFcn = new TF1("signalFcn",signal,0.,89.5,3); + signalFcn->SetLineWidth(2); + signalFcn->SetLineColor(kBlue); + signalFcn->SetNpx(500); + + TF1 *fitFcn = new TF1("fitFcn",fitFunction,0.,89.5,7); + fitFcn->SetNpx(500); + fitFcn->SetLineWidth(2); + fitFcn->SetLineColor(kMagenta); + + // start-parameters + Double_t p3start=halpha->GetBinContent(halpha->GetNbinsX()-1); + Double_t p4start=0.; + Double_t p5start=0.; + Double_t p6start=0.; + + Double_t p0start=(halpha->GetEntries()*halpha->GetBinWidth(1)-p3start)*0.05; + Double_t p1start=5.; + Double_t p2start=0.; + + cout<<"start values:"<Fit ("backFcn","RN"); + backFcn->GetParameters(parb); + + // total fit + Double_t par[7]; + fitFcn->SetParameters(p0start,p1start,p2start, + parb[0],parb[1],parb[2],parb[3]); + + fitFcn->FixParameter(2,0.); + fitFcn->FixParameter(3,parb[0]); + fitFcn->FixParameter(4,parb[1]); + fitFcn->FixParameter(5,parb[2]); + fitFcn->FixParameter(6,parb[3]); + halpha->Fit ("fitFcn","RN"); + + // draw the legend + TLegend *legend=new TLegend(0.4,0.5,0.93,0.65); + legend->SetTextFont(40); + legend->SetTextSize(0.04); + legend->SetFillColor(19); + legend->AddEntry(halpha,"Data","lp"); + TString strpol(""); + strpol=Form("Polynom"); + legend->AddEntry(backFcn,strpol.Data(),"l"); + legend->AddEntry(fitFcn,"Global Fit","l"); + //legend->Draw(); + + // drawing the functions + fitFcn->GetParameters(par); + fitFcn->SetRange(-90.,90.); + fitFcn->SetLineWidth(2); + fitFcn->Draw("same"); + + signalFcn->SetParameters(par); + signalFcn->SetLineWidth(2); + //signalFcn->Draw("same"); + + backFcn->SetParameters(&par[3]); + backFcn->SetRange(-90,90); + backFcn->SetLineWidth(2); + backFcn->Draw("same"); + + + //----------------------------------------------------------------------- + // calculate significance + + // analytical method + Double_t nSig=signalFcn->Integral(0.,acut)/halpha->GetBinWidth(1); + Double_t nOffScaled=backFcn->Integral(0.,acut)/halpha->GetBinWidth(1); + + Double_t nOn=nSig+nOffScaled; + + // significance after Li/Ma + Double_t theta=1.; // should be = Ton/Toff ??? + nSig=nOn-nOffScaled; + Double_t nOff=nOffScaled/theta; + + Double_t siglima = sqrt(2*(nOn*TMath::Log((1+theta)*nOn/(theta*(nOn+nOff))) + +nOff*TMath::Log((1+theta)*(nOff/(nOn+nOff))))); + + // "standard" significance + Double_t sigstd = nSig/sqrt(nOn+nOffScaled); + + + TPaveText *pt = new TPaveText(0.25,0.75,0.56,0.98,"brNDC"); + + pt->SetTextFont(60); + pt->SetTextSize(0.03); + pt->SetTextAlign(12); + + TString strcut1(Form("cuts:")); + TString strcut3(Form(" |alpha| < %4.2f ",acut)); + pt->AddText(strcut3.Data()); + TString str(""); + pt->AddText(str); + + TString strsig(Form("NExcess = %d",Int_t(nSig))); + pt->AddText(strsig.Data()); + TString strbak(Form("NOff = %d",Int_t(nOffScaled))); + pt->AddText(strbak.Data()); + pt->AddText(str); + + TString strsg(Form("Significance (Li/Ma) = %4.2f",siglima)); + pt->AddText(strsg.Data()); + + cout<<"results for events with alpha < "< 1300 +Double_t LengthLo = 0.318; +Double_t LengthUp = 0.809; +Double_t WidthLo = 0.3; +Double_t WidthUp = 0.705; +//*/ Index: trunk/MagicSoft/Mars/mtemp/mberlin/macros/ScaleHillas.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mberlin/macros/ScaleHillas.C (revision 4137) +++ trunk/MagicSoft/Mars/mtemp/mberlin/macros/ScaleHillas.C (revision 4137) @@ -0,0 +1,626 @@ +#include "TH1F.h" +#include "TH2F.h" +#include "TH3F.h" +#include "TLine.h" +#include "TTree.h" + +#include "TCanvas.h" +#include "TStyle.h" +#include "TString.h" +#include "TFile.h" +#include "TMath.h" + +#include "MReadTree.h" +#include "MSrcPosCam.h" +#include "MWriteRootFile.h" +#include "MParList.h" +#include "MTaskList.h" +#include "MEvtLoop.h" + +#include "MHillas.h" +#include "MHillasExt.h" +#include "MHillasExt.h" +#include "MHillasSrc.h" +#include "MNewImagePar.h" + +#include "MTaskInteractive.h" + +#include "iostream.h" +#include "stdlib.h" +#include "stdio.h" +#include "string.h" + +//-------------------------------------------------------------------------------------- +// MAIN INPUT/OUTPUT + +#include "IOMkn421.h" +//-------------------------------------------------------------------------------------- + +Double_t estposx = 0.; +Double_t estposy = 0.; + +//-------------------------------------------------------------------------------------- +// setting for hists +const Int_t nBinSize=200; +const Int_t nBinDist=1; +const Int_t nBinLW=1000; + +// limits for size and dist for histograms +Double_t hSizeLo=2.; +Double_t hSizeUp=6.; +Double_t hDistLo=0.5; +Double_t hDistUp=1.15; + +// limits for length, width for histograms +Double_t hLengthLo=0.; +Double_t hLengthUp=0.8; +Double_t hWidthLo =0.; +Double_t hWidthUp =0.6; + +// limits for scaled Hillas for histograms +Double_t hSLenLo =0.; +Double_t hSLenUp =15.; +Double_t hSWidLo =0.; +Double_t hSWidUp =7.; + +// special limits for hLength, hWidth (see ReadHillas()) +const Int_t nBinLWFine=1000; +Double_t hLenFineLo=0.; +Double_t hLenFineUp=1.; +Double_t hWidFineLo=0.; +Double_t hWidFineUp=1.; + +// hists for mean of unscaled Hillas +TH2F hLenMean("hLenMean","",nBinSize,hSizeLo,hSizeUp, nBinDist,hDistLo,hDistUp); +TH2F hWidMean("hWidMean","",nBinSize,hSizeLo,hSizeUp, nBinDist,hDistLo,hDistUp); + +TH1F hLenMeanSize("hLenMeanSize","",nBinSize,hSizeLo,hSizeUp); +TH1F hLenMeanDist("hLenMeanDist","",nBinDist,hDistLo,hDistUp); +TH1F hWidMeanSize("hWidMeanSize","",nBinSize,hSizeLo,hSizeUp); +TH1F hWidMeanDist("hWidMeanDist","",nBinDist,hDistLo,hDistUp); + +// hists for unscaled Hillas +TH2F hLenSize("hLenSize","",nBinSize,hSizeLo,hSizeUp, nBinLW,hLengthLo,hLengthUp); +TH2F hLenDist("hLenDist","",nBinDist,hDistLo,hDistUp, nBinLW,hLengthLo,hLengthUp); + +TH2F hWidSize("hWidSize","",nBinSize,hSizeLo,hSizeUp, nBinLW,hWidthLo,hWidthUp); +TH2F hWidDist("hWidDist","",nBinDist,hDistLo,hDistUp, nBinLW,hWidthLo,hWidthUp); + +// hists for mean-scaled Hillas +TH2F hSLenSize("hSLenSize","",nBinSize,hSizeLo,hSizeUp,nBinLW,hSLenLo,hSLenUp); +TH2F hSLenDist("hSLenDist","",nBinDist,hDistLo,hDistUp,nBinLW,hSLenLo,hSLenUp); + +TH2F hSWidSize("hSWidSize","",nBinSize,hSizeLo,hSizeUp,nBinLW,hSWidLo,hSWidUp); +TH2F hSWidDist("hSWidDist","",nBinDist,hDistLo,hDistUp,nBinLW,hSWidLo,hSWidUp); + + +//**************************************************************************************** +// first part: prepare for scaling, read Hillas to calc mean +//**************************************************************************************** + +void ReadHillas() +{ + // fill hillas into histograms to calc mean + TFile *file=new TFile(fileHillas.Data()); + TTree *fEvents=(TTree*)file->Get("Events"); + + // hists for Hillas + // in following 2 hists the binning in length/width should be + // very fine (like unbinned) and the limits should be generous + TH3F hLength("hLength","",nBinSize,hSizeLo,hSizeUp, nBinDist,hDistLo,hDistUp, nBinLWFine,hLenFineLo,hLenFineUp); + TH3F hWidth("hWidth", "",nBinSize,hSizeLo,hSizeUp, nBinDist,hDistLo,hDistUp, nBinLWFine,hWidFineLo,hWidFineUp); + + TString strlen = Form("MHillas.fLength*%f:MHillasSrc.fDist*%f:log10(MHillas.fSize)",mm2deg,mm2deg); + TString strwid = Form("MHillas.fWidth*%f :MHillasSrc.fDist*%f:log10(MHillas.fSize)",mm2deg,mm2deg); + + TString strcut="MHillas.fSize>0"; + fEvents->Project("hLength",strlen.Data(),strcut.Data()); + fEvents->Project("hWidth", strwid.Data(),strcut.Data()); + + file->Close(); + + // fill hists hLenMean, hWidMean (2-dim (size,dist)-hists) + for(Int_t i=1;i<=nBinSize;i++) + for(Int_t j=1;j<=nBinDist;j++) + { + Double_t lsum=0.; Double_t lsum2=0.; Double_t lcnt=0.; + Double_t wsum=0.; Double_t wsum2=0.; Double_t wcnt=0.; + + for(Int_t k=1;k<=nBinLWFine;k++) + { + Double_t nl = hLength.GetBinContent(i,j,k); + Double_t xl = hLength.GetZaxis()->GetBinCenter(k); + + lsum +=nl*xl; + lsum2+=nl*xl*xl; + lcnt +=nl; + + Double_t nw = hWidth.GetBinContent(i,j,k); + Double_t xw = hWidth.GetZaxis()->GetBinCenter(k); + + wsum +=nw*xw; + wsum2+=nw*xw*xw; + wcnt +=nw; + } + + Double_t lmean=0.; + Double_t lmeanerr=0.; + if(lcnt>1)lmean=lsum/lcnt; + if(lcnt>2)lmeanerr=(lsum2/lcnt-lmean*lmean)/(lcnt-1.); + lmeanerr=TMath::Sqrt(lmeanerr); + + hLenMean.SetBinContent(i,j,lmean); + hLenMean.SetBinError(i,j,lmeanerr); + + Double_t wmean=0.; + Double_t wmeanerr=0.; + if(wcnt>1)wmean=wsum/wcnt; + if(wcnt>2)wmeanerr=(wsum2/wcnt-wmean*wmean)/(wcnt-1.); + wmeanerr=TMath::Sqrt(wmeanerr); + + hWidMean.SetBinContent(i,j,wmean); + hWidMean.SetBinError(i,j,wmeanerr); + } + + // fill hists hLenMeanSize, hWidMeanSize (1-dim hists, x-axis = size) + for(Int_t i=1;i<=nBinSize;i++) + { + Double_t lsum=0.; Double_t lsum2=0.; Double_t lcnt=0.; + Double_t wsum=0.; Double_t wsum2=0.; Double_t wcnt=0.; + + for(Int_t j=1;j<=nBinDist;j++) + for(Int_t k=1;k<=nBinLWFine;k++) + { + Double_t nl = hLength.GetBinContent(i,j,k); + Double_t xl = hLength.GetZaxis()->GetBinCenter(k); + + lsum +=nl*xl; + lsum2+=nl*xl*xl; + lcnt +=nl; + + Double_t nw = hWidth.GetBinContent(i,j,k); + Double_t xw = hWidth.GetZaxis()->GetBinCenter(k); + + wsum +=nw*xw; + wsum2+=nw*xw*xw; + wcnt +=nw; + } + + Double_t lmean=0.; + Double_t lmeanerr=0.; + if(lcnt>1)lmean=lsum/lcnt; + if(lcnt>2)lmeanerr=(lsum2/lcnt-lmean*lmean)/(lcnt-1.); + lmeanerr=TMath::Sqrt(lmeanerr); + + hLenMeanSize.SetBinContent(i,lmean); + hLenMeanSize.SetBinError(i,lmeanerr); + + Double_t wmean=0.; + Double_t wmeanerr=0.; + if(wcnt>1)wmean=wsum/wcnt; + if(wcnt>2)wmeanerr=(wsum2/wcnt-wmean*wmean)/(wcnt-1.); + wmeanerr=TMath::Sqrt(wmeanerr); + + hWidMeanSize.SetBinContent(i,wmean); + hWidMeanSize.SetBinError(i,wmeanerr); + } + + // fill hists hLenMeanDist, hWidMeanDist (1-dim hists, x-axis = dist) + for(Int_t j=1;j<=nBinDist;j++) + { + Double_t lsum=0.; Double_t lsum2=0.; Double_t lcnt=0.; + Double_t wsum=0.; Double_t wsum2=0.; Double_t wcnt=0.; + + for(Int_t i=1;i<=nBinSize;i++) + for(Int_t k=1;k<=nBinLWFine;k++) + { + Double_t nl = hLength.GetBinContent(i,j,k); + Double_t xl = hLength.GetZaxis()->GetBinCenter(k); + + lsum +=nl*xl; + lsum2+=nl*xl*xl; + lcnt +=nl; + + Double_t nw = hWidth.GetBinContent(i,j,k); + Double_t xw = hWidth.GetZaxis()->GetBinCenter(k); + + wsum +=nw*xw; + wsum2+=nw*xw*xw; + wcnt +=nw; + } + + Double_t lmean=0.; + Double_t lmeanerr=0.; + if(lcnt>1)lmean=lsum/lcnt; + if(lcnt>2)lmeanerr=(lsum2/lcnt-lmean*lmean)/(lcnt-1.); + lmeanerr=TMath::Sqrt(lmeanerr); + + hLenMeanDist.SetBinContent(j,lmean); + hLenMeanDist.SetBinError(j,lmeanerr); + + Double_t wmean=0.; + Double_t wmeanerr=0.; + if(wcnt>1)wmean=wsum/wcnt; + if(wcnt>2)wmeanerr=(wsum2/wcnt-wmean*wmean)/(wcnt-1.); + wmeanerr=TMath::Sqrt(wmeanerr); + + hWidMeanDist.SetBinContent(j,wmean); + hWidMeanDist.SetBinError(j,wmeanerr); + } + + return; +} + + +//**************************************************************************************** +// second part: do scaling +//**************************************************************************************** + +// pointer to Hillas ParContainer +MHillas *fHillas=NULL; +MHillasSrc *fHillasSrc=NULL; + + +Int_t PreProcess(MParList *plist) +{ + fHillas = (MHillas*) plist->FindObject("MHillas"); + fHillasSrc = (MHillasSrc*) plist->FindObject("MHillasSrc"); + + return (fHillas && fHillasSrc) ? kTRUE : kFALSE; +} + +Int_t Process() +{ + // estimated source pos. + MSrcPosCam srcpos; + srcpos.SetXY(estposx*deg2mm,estposy*deg2mm); + fHillasSrc->SetSrcPos(&srcpos); + fHillasSrc->Calc(fHillas); + + Float_t width = fHillas->GetWidth(); + Float_t length = fHillas->GetLength(); + Float_t size = fHillas->GetSize(); + Float_t meanx = fHillas->GetMeanX(); + Float_t meany = fHillas->GetMeanY(); + Float_t logsize = TMath::Log10(fHillas->GetSize()); + Float_t delta = fHillas->GetDelta(); + Float_t dist = fHillasSrc->GetDist(); + + width *=mm2deg; + length*=mm2deg; + dist *=mm2deg; + + hLenSize.Fill(logsize,length); + hLenDist.Fill(dist,length); + hWidSize.Fill(logsize,width); + hWidDist.Fill(dist,width); + + // (log(size), dist) - bin + Int_t ix=hLenMean.GetXaxis()->FindBin(logsize); + Int_t iy=hLenMean.GetYaxis()->FindBin(dist); + + // check for over-/underflow + /*if(ix*iy==0 || ix>nBinSize || iy>nBinDist) + { + cout<<"Warning, over-/underflow: Log(size)-bin="<SetPalette(1); + + + TCanvas *c1=new TCanvas("c1","",800,400); + if(nBinDist>1) + c1->Divide(2,2); + else + c1->Divide(2,1); + + c1->cd(1); + gPad->SetGrid(); + hSLenSize.SetStats(kFALSE); + hSLenSize.GetXaxis()->SetTitle("log(size [#phot])"); + hSLenSize.GetYaxis()->SetTitle("scaled length"); + hSLenSize.SetLineWidth(2); + hSLenSize.DrawCopy("colz"); + hSLenSize.Write(); + + TLine lLen(2.7,1.,hSizeUp,1.); + TLine lLenUp(2.7,LengthUp,hSizeUp,LengthUp); + TLine lLenLo(2.7,LengthLo,hSizeUp,LengthLo); + + lLen.SetLineColor(2); + lLen.SetLineWidth(2); + lLenUp.SetLineWidth(2); + lLenLo.SetLineWidth(2); + lLen.DrawClone(); + lLenUp.DrawClone(); + lLenLo.DrawClone(); + + + c1->cd(2); + gPad->SetGrid(); + hLenSize.SetStats(kFALSE); + hLenSize.SetTitle(""); + hLenSize.GetXaxis()->SetTitle("log(size [#phot])"); + hLenSize.GetYaxis()->SetTitle("length [deg]"); + hLenSize.SetLineWidth(2); + hLenSize.DrawCopy("colz"); + hLenSize.Write(); + hLenMeanSize.SetLineColor(2); + hLenMeanSize.SetLineWidth(2); + hLenMeanSize.DrawCopy("e same"); + hLenMeanSize.Write(); + + TH1F hLenCutLo("hLenCutLo","",nBinSize,hSizeLo,hSizeUp); + TH1F hLenCutUp("hLenCutUp","",nBinSize,hSizeLo,hSizeUp); + for(Int_t i=1;i<=nBinSize;i++) + { + Float_t meanval=hLenMeanSize.GetBinContent(i); + hLenCutLo.SetBinContent(i,meanval*LengthLo); + hLenCutUp.SetBinContent(i,meanval*LengthUp); + } + hLenCutLo.SetLineWidth(2); + hLenCutUp.SetLineWidth(2); + hLenCutLo.DrawCopy("][ same"); + hLenCutUp.DrawCopy("][ same"); + + + if(nBinDist>1) + { + c1->cd(3); + gPad->SetGrid(); + hSLenDist.SetStats(kFALSE); + hSLenDist.GetXaxis()->SetTitle("dist [deg]"); + hSLenDist.GetYaxis()->SetTitle("scaled length"); + hSLenDist.SetLineWidth(2); + hSLenDist.DrawCopy("colz"); + hSLenDist.Write(); + + c1->cd(4); + gPad->SetGrid(); + hLenDist.SetStats(kFALSE); + hLenDist.SetTitle(""); + hLenDist.GetXaxis()->SetTitle("dist [deg]"); + hLenDist.GetYaxis()->SetTitle("length"); + hLenDist.SetLineWidth(2); + hLenDist.DrawCopy("colz"); + hLenDist.Write(); + hLenMeanDist.SetLineColor(2); + hLenMeanDist.SetLineWidth(2); + hLenMeanDist.DrawCopy("e same"); + hLenMeanDist.Write(); + } + + TCanvas *c2=new TCanvas("c2","",800,400); + if(nBinDist>1) + c2->Divide(2,2); + else + c2->Divide(2,1); + + c2->cd(1); + gPad->SetGrid(); + hSWidSize.SetStats(kFALSE); + hSWidSize.GetXaxis()->SetTitle("log(size [#phot])"); + hSWidSize.GetYaxis()->SetTitle("scaled width"); + hSWidSize.SetLineWidth(2); + hSWidSize.DrawCopy("colz"); + hSWidSize.Write(); + + TLine lWid(2.7,1.,hSizeUp,1.); + TLine lWidUp(2.7,WidthUp,hSizeUp,WidthUp); + TLine lWidLo(2.7,WidthLo,hSizeUp,WidthLo); + + lWid.SetLineColor(2); + lWid.SetLineWidth(2); + lWidUp.SetLineWidth(2); + lWidLo.SetLineWidth(2); + lWid.DrawClone(); + lWidUp.DrawClone(); + lWidLo.DrawClone(); + + c2->cd(2); + gPad->SetGrid(); + hWidSize.SetStats(kFALSE); + hWidSize.SetTitle(""); + hWidSize.GetXaxis()->SetTitle("log(size [#phot])"); + hWidSize.GetYaxis()->SetTitle("width [deg]"); + hWidSize.SetLineWidth(2); + hWidSize.DrawCopy("colz"); + hWidSize.Write(); + hWidMeanSize.SetLineColor(2); + hWidMeanSize.SetLineWidth(2); + hWidMeanSize.DrawCopy("e same"); + hWidMeanSize.Write(); + + TH1F hWidCutLo("hWidCutLo","",nBinSize,hSizeLo,hSizeUp); + TH1F hWidCutUp("hWidCutUp","",nBinSize,hSizeLo,hSizeUp); + for(Int_t i=1;i<=nBinSize;i++) + { + Float_t meanval=hWidMeanSize.GetBinContent(i); + hWidCutLo.SetBinContent(i,meanval*WidthLo); + hWidCutUp.SetBinContent(i,meanval*WidthUp); + } + hWidCutLo.SetLineWidth(2); + hWidCutUp.SetLineWidth(2); + hWidCutLo.DrawCopy("][ same"); + hWidCutUp.DrawCopy("][ same"); + + if(nBinDist>1) + { + c2->cd(3); + gPad->SetGrid(); + hSWidDist.SetStats(kFALSE); + hSWidDist.GetXaxis()->SetTitle("dist [deg]"); + hSWidDist.GetYaxis()->SetTitle("scaled width"); + hSWidDist.SetLineWidth(2); + hSWidDist.DrawCopy("colz"); + hSWidDist.Write(); + + c2->cd(4); + gPad->SetGrid(); + hWidDist.SetStats(kFALSE); + hWidDist.SetTitle(""); + hWidDist.GetXaxis()->SetTitle("dist [deg]"); + hWidDist.GetYaxis()->SetTitle("width"); + hWidDist.SetLineWidth(2); + hWidDist.DrawCopy("colz"); + hWidDist.Write(); + hWidMeanDist.SetLineColor(2); + hWidMeanDist.SetLineWidth(2); + hWidMeanDist.DrawCopy("e same"); + hWidMeanDist.Write(); + } + + // mean values + if(nBinDist>1) + { + TCanvas *c4=new TCanvas("c4","",800,400); + c4->Divide(1,2); + + c4->cd(1); + gPad->SetGrid(); + hLenMean.SetStats(kFALSE); + hLenMean.SetTitle("mean_length(log(size),dist)"); + hLenMean.GetXaxis()->SetTitle("log(size [#phot])"); + hLenMean.GetYaxis()->SetTitle("dist [deg])"); + hLenMean.DrawCopy("colz"); + hLenMean.Write(); + + c4->cd(2); + gPad->SetGrid(); + hWidMean.SetStats(kFALSE); + hWidMean.SetTitle("mean_width(log(size),dist)"); + hWidMean.GetXaxis()->SetTitle("log(size [#phot])"); ; + hWidMean.GetYaxis()->SetTitle("dist [deg])"); + hWidMean.DrawCopy("colz"); + hWidMean.Write("colz"); + } + + fileHist.Close(); + + return; +} + +//**************************************************************************************** +// main +//**************************************************************************************** + +void RunScaling(Int_t nent=-1) +{ + printf("Prepare for scaling..."); + ReadHillas(); + printf("done.\n"); + + ScaleHillas(nent); + + CheckScaling(); + + return; +} + Index: trunk/MagicSoft/Mars/mtemp/mberlin/macros/SrcCorrect.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mberlin/macros/SrcCorrect.C (revision 4137) +++ trunk/MagicSoft/Mars/mtemp/mberlin/macros/SrcCorrect.C (revision 4137) @@ -0,0 +1,768 @@ +#include "MParList.h" +#include "MTaskList.h" +#include "MEvtLoop.h" + +#include "MReadTree.h" + +#include "MHillasSrc.h" +#include "MSrcPosCam.h" + +#include "MHillasCalc.h" +#include "MHillasSrcCalc.h" +#include "MTaskInteractive.h" +#include "MContinue.h" +#include "MRawRunHeader.h" +#include "MWriteRootFile.h" + +#include "TH2F.h" +#include "TF2.h" +#include "TArrayF.h" +#include "TObjArray.h" +#include "TStyle.h" +#include "TFile.h" +#include "TMath.h" +#include "TString.h" +#include "TCanvas.h" +#include "TGraphErrors.h" +#include "TMarker.h" + +#include "iostream.h" + +//---------------------------------------------------------------- +// MAIN INPUT/OUTPUT + +#include "IOMkn421.h" +//---------------------------------------------------------------- + +Double_t rmean=0.3; +Double_t dtlimit=10; // time slice width [min] +//Double_t hlimit =0.85; +Double_t hlimit =0.9; + +//---------------------------------------------------------------- +// false source hist --------------------------------------------- +//---------------------------------------------------------------- +const Int_t nx = 100; // no. of points in x-grid +const Int_t ny = 100; // no. of points in y-grid + +Double_t alimit = 8.; // signal region in alpha-plot [degree] + +Double_t xmin = -1.; // units in degree !! +Double_t xmax = 1.; +Double_t ymin = -1.; +Double_t ymax = 1.; + +Double_t dx=(xmax-xmin)/Double_t(nx); +Double_t dy=(ymax-ymin)/Double_t(ny); + +TH2F FSrcHist("FSrcHist","",nx,xmin-dx,xmax+dx,ny,ymin-dy,ymax+dy); +TObjArray FSrcArray(1); + +TH1F APlot("APlot","",36,0,90); +TObjArray APlotArray(1); + +TArrayF xpeak(100); +TArrayF ypeak(100); + +TArrayF xpeakerr(100); +TArrayF ypeakerr(100); + +Int_t runno = -1; +Int_t runno_old = -1; +Int_t runno_lo = -1; +Int_t runno_up = -1; + +Double_t mjdend_old = -1.; +Double_t mjdstart_lo = -1.; +Double_t mjdstart_up = -1.; + +Double_t mjd_lo = -1.; +Double_t mjd_up = -1.; + +Double_t dt = 0.; +Int_t islice = 0; + +TF2 *fpeak = NULL; +FILE *fsrcpos = NULL; + +//**************************************************************************************** +// first part: track the source in time-intervalls of duration dtlimit +//**************************************************************************************** + +MHillas *hil=NULL; +MRawRunHeader *run=NULL; + +Double_t gaus2d(Double_t *x,Double_t *par) +{ + Double_t N=par[0]; //integral + Double_t mx=par[1]; //mean_x + Double_t my=par[2]; //mean_y + Double_t sx=par[3]; //sigma_x + Double_t sy=par[4]; //sigma_y + Double_t rho=par[5];//correlation factor + Double_t b=par[6]; //constant background + + Double_t z=(x[0]-mx)*(x[0]-mx)/sx/sx; + z+=(x[1]-my)*(x[1]-my)/sy/sy; + z-=2.*rho*(x[0]-mx)*(x[1]-my)/sx/sy; + + Double_t val=N/(2.*TMath::Pi()*sx*sy*sqrt(1-rho*rho)); + val*=exp( -z/(2.*(1.-rho*rho)) ); + val+=b; + return val; +} + +Bool_t CalcPeakXY(Double_t *xmean, Double_t *ymean, + Double_t *xmeanerr, Double_t *ymeanerr) +{ + Int_t mbin=FSrcHist.GetMaximumBin(); + Int_t biny = mbin/(nx+2); + Int_t binx = mbin%(nx+2); + + if(mbin!=FSrcHist.GetBin(binx,biny)) + { + cout<<"Error extracting binx, biny from global bin!"<GetParameter(1); + ybar=fpeak->GetParameter(2); + xbarerr=fpeak->GetParameter(3); + ybarerr=fpeak->GetParameter(4); + */ + + *xmean=xbar; + *ymean=ybar; + *xmeanerr=xbarerr; + *ymeanerr=ybarerr; + + return kTRUE; +} + +void AddFSrcHist(Int_t i,Double_t delta) +{ + Int_t n=FSrcArray.GetSize(); + if(i>=n) FSrcArray.Expand(n+1); + + TH2F *htemp=(TH2F*)FSrcHist.Clone(); + TString name=Form("%d",i); + TString title=Form("TimeIntervall = %f min",delta); + + htemp->SetName(name.Data()); + htemp->SetTitle(title.Data()); + htemp->SetDrawOption("colz"); + + FSrcArray[i] = htemp; + + return; +} + +void AddAPlot(Int_t i,Double_t delta) +{ + Int_t n=APlotArray.GetSize(); + if(i>=n) APlotArray.Expand(n+1); + + TH1F *htemp=(TH1F*)APlot.Clone(); + TString name=Form("%d",i); + TString title=Form("TimeIntervall = %f min",delta); + + htemp->SetName(name.Data()); + htemp->SetTitle(title.Data()); + htemp->SetFillColor(16); + + APlotArray[i] = htemp; + + return; +} + +//---------------------------------------------------------------- +// Source Position Grid ------------------------------------------ +//---------------------------------------------------------------- + +Int_t STrackPreProc(MParList *plist) +{ + hil = (MHillas*) plist->FindObject("MHillas"); + run = (MRawRunHeader*) plist->FindObject("MRawRunHeader"); + + return (hil && run) ? kTRUE : kFALSE; +} + +Int_t STrackProc() +{ + runno = run->GetRunNumber(); + if(runno==0) return kTRUE; + + Double_t mjdstart = (run->GetRunStart()).GetMjd(); + Double_t mjdend = (run->GetRunEnd()).GetMjd(); + Double_t t = (mjdend-mjdstart)*24.*60.; + + Double_t dtrun; + if(runno!=runno_old) + { + if(dt<0.001) {runno_lo=runno;mjd_lo=mjdstart;mjdend_old=mjdstart;} // very first run + + dtrun=(mjdend_old-mjd_lo)*24.*60.; + //if(dt>dtlimit) + if(dtrun>dtlimit) + { + runno_up=runno_old; + mjd_up=mjdend_old; + + Double_t x; Double_t xerr; + Double_t y; Double_t yerr; + + if(!CalcPeakXY(&x,&y,&xerr,&yerr)) return kFALSE; + AddFSrcHist(islice,dt); + + xpeak[islice]=x; + ypeak[islice]=y; + xpeakerr[islice]=xerr; + ypeakerr[islice]=yerr; + + printf("\n TimeSlice %d RunNo %d-%d Mjd %f-%f Dt=%f \n",islice,runno_lo,runno_up,mjd_lo,mjd_up,dt); + printf(" X=%f+-%f Y=%f+-%f \n\n",x,xerr,y,yerr); + fprintf(fsrcpos,"%d %d %d %f %f %f %f %f %f %f\n",islice,runno_lo,runno_up,x,y,xerr,yerr,mjd_lo,mjd_up,dt); + + FSrcHist.Reset(); + dt=0.; + + runno_lo=runno; + mjd_lo=mjdstart; + + islice++; + } + + dt+=t; + printf("Runno=%d MjdStart=%f MjdEnde=%f dt=%f \n",runno,mjdstart,mjdend,dt); + runno_old=runno; + mjdend_old=mjdend; + } + + + MHillasSrc hillas_src; + MSrcPosCam srcpos; + + Double_t length=hil->GetLength(); + Double_t width=hil->GetWidth(); + + // size cut + if(hil->GetSize()LengthUp) || (lengthWidthUp) || (widthDistUp || dist0.001) + { + runno_up=runno_old; + mjd_up=mjdend_old; + + Double_t x; Double_t xerr; + Double_t y; Double_t yerr; + + if(!CalcPeakXY(&x,&y,&xerr,&yerr)) return kFALSE; + AddFSrcHist(islice,dt); + + xpeak[islice]=x; + ypeak[islice]=y; + xpeakerr[islice]=xerr; + ypeakerr[islice]=yerr; + + printf("\n TimeSlice %d RunNo %d-%d Mjd %f-%f Dt=%f \n",islice,runno_lo,runno_up,mjd_lo,mjd_up,dt); + printf(" X=%f+-%f Y=%f+-%f \n\n",x,xerr,y,yerr); + fprintf(fsrcpos,"%d %d %d %f %f %f %f %f %f %f\n",islice,runno_lo,runno_up,x,y,xerr,yerr,mjd_lo,mjd_up,dt); + + islice++; + } + + return kTRUE; +} + +//---------------------------------------------------------------- +//---------------------------------------------------------------- + +void SourceTrack(Int_t nent) +{ + + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + MReadTree read("Events",fileHScaled); + read.DisableAutoScheme(); + + //TString runcuts("(MRawRunHeader.fRunNumber<17427)"); + //runcuts+="|| (MRawRunHeader.fRunNumber>17428)"; + //MContinue runsel(runcuts.Data()); + + MTaskInteractive strack; + strack.SetPreProcess(STrackPreProc); + strack.SetProcess(STrackProc); + strack.SetPostProcess(STrackPostProc); + + tlist.AddToList(&read); + //tlist.AddToList(&runsel); + tlist.AddToList(&strack); + + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + //------------------------------------------ + // initializations + + FSrcHist.SetDrawOption("colz"); + FSrcHist.SetStats(kFALSE); + + fpeak=new TF2("fpeak",gaus2d,-1,1,-1,1,6); + fsrcpos=fopen(fileSrcPos.Data(),"w"); + + //------------------------------------------ + // standard eventloop + + if (!evtloop.Eventloop(nent)) + return; + + tlist.PrintStatistics(); + + fclose(fsrcpos); + + TFile file(fileFSrcPlot.Data(),"recreate"); + FSrcArray.Write(); + FSrcHist.Write(); + file.Close(); + + const Int_t dim=islice; + TGraphErrors g(dim,xpeak.GetArray(),ypeak.GetArray(), + xpeakerr.GetArray(),ypeakerr.GetArray()); + + //------------------------------------------ + // canvas with complete hist + TCanvas *c1=new TCanvas("c1","",400,800); + gStyle->SetPalette(1); + + c1->Divide(2,1); + c1->cd(1); + gPad->Update(); + + FSrcHist.Reset(); + for(Int_t i=0;i%f ",SizeLo); + label+=Form("%fSetTitle("X"); + FSrcHist.GetYaxis()->SetTitle("Y"); + + FSrcHist.DrawClone("colz"); + g.DrawClone("pl"); + + c1->cd(2); + gPad->Update(); + g.GetXaxis()->SetRangeUser(-1.,1.); + g.GetYaxis()->SetRangeUser(-1.,1.); + g.DrawClone("apl"); + + //------------------------------------------ + // canvas with slides + TCanvas *c2=new TCanvas("c2","",800,800); + gStyle->SetPalette(1); + + Int_t n=Int_t(TMath::Ceil(sqrt(double(dim)))); + + c2->Divide(n,n); + for(Int_t i=0;icd(row*n+col+1); + + gPad->Update(); + TH2F *hptr = (TH2F*) (FSrcArray[i]); + + TMarker m; + m.SetMarkerStyle(2); + m.SetMarkerSize(4); + m.SetMarkerColor(1); + m.SetX(xpeak[i]); + m.SetY(ypeak[i]); + + hptr->DrawClone("colz"); + m.DrawClone(); + } + + return; +} + +//**************************************************************************************** +// second part: correct for moving source position +//**************************************************************************************** + +// initilizations + +Double_t dtref=0.; +Int_t islice_old; +Double_t dtref_old=0; + +Double_t xref = 0.; Double_t xreferr = 0.; +Double_t yref = 0.; Double_t yreferr = 0.; + +MSrcPosCam *src=NULL; + +//---------------------------------------------------------------- +// Source Position Grid ------------------------------------------ +//---------------------------------------------------------------- +Int_t SCorrectPreProc(MParList *plist) +{ + hil = (MHillas*) plist->FindObject("MHillas"); + src = (MSrcPosCam*) plist->FindCreateObj("MSrcPosCam"); + run = (MRawRunHeader*) plist->FindObject("MRawRunHeader"); + + return (hil && run && src) ? kTRUE : kFALSE; +} + +Int_t SCorrectProc() +{ + if(hil->GetSize()<=0.) return kTRUE; + + Int_t runno = run->GetRunNumber(); + if(runno==0) return kTRUE; + + Double_t mjdstart = (run->GetRunStart()).GetMjd(); + Double_t mjdend = (run->GetRunEnd()).GetMjd(); + Double_t t=(mjdend-mjdstart)*24.*60.; + + if(runno!=runno_old) + { + if((runnorunno_up)) + { + dt=t; + + fscanf(fsrcpos,"%d %d %d %lf %lf %lf %lf %lf %lf %lf\n",&islice,&runno_lo,&runno_up,&xref,&yref,&xreferr,&yreferr,&mjd_lo,&mjd_up,&dtref); + printf("\nApplying correction for runs %d - %d (MJD %f - %f), DtRef=%f :\n",runno_lo,runno_up,mjd_lo,mjd_up,dtref); + printf("Source pos set to X=%f Y=%f \n\n",xref,yref); + + if(islice>0) + { + AddAPlot(islice_old,dtref_old); + APlot.Reset(); + } + islice_old=islice; + dtref_old=dtref; + } + else + { + dt+=t; + } + printf("Runno=%d MjdStart=%f MjdEnde=%f Dt=%f \n",runno,mjdstart,mjdend,dt); + runno_old=runno; + + } + + src->SetXY(xref*deg2mm,yref*deg2mm); + src->SetReadyToSave(); + // false source grid + + MHillasSrc hillas_src; + MSrcPosCam srcpos; + + Double_t length=hil->GetLength(); + Double_t width =hil->GetWidth(); + + // size cut + if(hil->GetSize()LengthUp) || (lengthWidthUp) || (widthDistUp || distSetPalette(1); + TFile file(fileFSrcPlotC.Data(),"recreate"); + FSrcHist.Write(); + APlotArray.Write(); + file.Close(); + + //------------------------------------------ + // canvas with complete hist + TCanvas *c=new TCanvas("c","",800,400); + c->cd(); + + TString label=Form("Size>%f ",SizeLo); + label+=Form("%fSetTitle("X"); + FSrcHist.GetYaxis()->SetTitle("Y"); + + FSrcHist.DrawCopy("colz"); + + //------------------------------------------ + // canvas with alpha hists + TCanvas *c2=new TCanvas("c2","",800,800); + gStyle->SetPalette(1); + + const Int_t dim=islice+1; + Int_t n=Int_t(TMath::Ceil(sqrt(double(dim)))); + + c2->Divide(n,n); + for(Int_t i=0;icd(row*n+col+1); + + gPad->Update(); + TH1F *hptr = (TH1F*) (APlotArray[i]); + + hptr->DrawClone(); + } + + return; +} + +//**************************************************************************************** +// main +//**************************************************************************************** + +void SrcCorrect(Int_t icorrect,Int_t nent=-1) +{ + + if(!icorrect) + { + // (re-) initilizations + islice=0; + runno_old=-1; + runno_lo=-1; + runno_up=-1; + + dt=0; + + mjd_lo=-1; + mjd_up=-1; + + hil=NULL; run=NULL; + + cout<<"Starting source track reconstruction."<