Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5775) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 5776) @@ -227,6 +227,6 @@ return kFALSE; - fChiSqSignal = fit.GetChiSqSignal(); - fChiSqBg = fit.GetChiSqBg(); + fChiSqSignal = fit.GetChiSqSignal(); + fChiSqBg = fit.GetChiSqBg(); fCoefficients = fit.GetCoefficients(); @@ -393,4 +393,5 @@ const TString name1(Form("TempAlpha%06d_on", gRandom->Integer(1000000))); const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000))); + TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E"); TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E"); @@ -486,3 +487,5 @@ if (scaleof!=0) of.Scale(scaleon/scaleof); -} + else + of.Reset(); +} Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5775) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h (revision 5776) @@ -55,5 +55,5 @@ public: // Implementing the function yourself is only about 5% faster - MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kEntries), fScaleUser(1) + MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kOffRegion), fScaleUser(1) { fName = name ? name : "MAlphaFitter"; Index: trunk/MagicSoft/Mars/mhflux/MHAlpha.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5775) +++ trunk/MagicSoft/Mars/mhflux/MHAlpha.cc (revision 5776) @@ -185,4 +185,6 @@ const Int_t n = fHAlpha.GetNbinsY(); + fHEnergy.SetEntries(0); + for (int i=1; i<=n; i++) { @@ -201,4 +203,6 @@ const Int_t n = fHAlpha.GetNbinsX(); + + fHTheta.SetEntries(0); for (int i=1; i<=n; i++) @@ -218,7 +222,7 @@ fHTheta.Reset(); - if (fName!=(TString)"MHAlphaOff") - { - MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff"); + if (fName!=(TString)"MHAlphaOff" && fOffData==NULL) + { + MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff", "MHAlpha"); if (!hoff) *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl; @@ -400,4 +404,5 @@ // Reset Histogram fHAlphaTime.Reset(); + fHAlphaTime.SetEntries(0); // Get number of bins @@ -646,5 +651,5 @@ h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E"); - h->SetBit(TH1::kNoTitle|TH1::kNoStats); + h->SetBit(TH1::kNoTitle); h->SetXTitle("\\alpha [\\circ]"); h->SetYTitle("Counts"); Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 5775) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (revision 5776) @@ -147,5 +147,5 @@ using namespace std; -class MHillasExt; +//class MHillasExt; // -------------------------------------------------------------------------- @@ -155,5 +155,6 @@ MHFalseSource::MHFalseSource(const char *name, const char *title) : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5), - fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinDW(-1), fMaxDW(-1) + fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinDW(-1), fMaxDW(-1), + fHistOff(0) { // @@ -174,4 +175,7 @@ void MHFalseSource::MakeSymmetric(TH1 *h) { + h->SetMinimum(); + h->SetMaximum(); + const Float_t min = TMath::Abs(h->GetMinimum()); const Float_t max = TMath::Abs(h->GetMaximum()); @@ -248,5 +252,5 @@ // Also search for MTime, MObservatory and MPointingPos // -MHillasExt *ext=0; +//MHillasExt *ext=0; Bool_t MHFalseSource::SetupFill(const MParList *plist) { @@ -257,4 +261,5 @@ return kFALSE; } + /* ext = (MHillasExt*)plist->FindObject("MHillasExt"); if (!ext) @@ -263,21 +268,38 @@ return kFALSE; } - + */ fMm2Deg = geom->GetConvMm2Deg(); - MBinning binsa; - binsa.SetEdges(18, 0, 90); - - const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource"); - if (!bins) + if (fName!=(TString)"MHFalseSourceOff" && fHistOff==NULL) { - const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg; - - MBinning b; - b.SetEdges(20, -r, r); - SetBinning(&fHist, &b, &b, &binsa); + MHFalseSource *hoff = (MHFalseSource*)plist->FindObject("MHFalseSourceOff", "MHFalseSource"); + if (!hoff) + *fLog << inf << "No MHFalseSourceOff [MHFalseSource] found... using current data only!" << endl; + else + { + *fLog << inf << "MHFalseSource [MHFalseSource] found... using on-off mode!" << endl; + SetOffData(*hoff); + } } + + if (fHistOff) + MH::SetBinning(&fHist, fHistOff); else - SetBinning(&fHist, bins, bins, &binsa); + { + MBinning binsa; + binsa.SetEdges(18, 0, 90); + + const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource"); + if (!bins) + { + const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg; + + MBinning b; + b.SetEdges(20, -r, r); + SetBinning(&fHist, &b, &b, &binsa); + } + else + SetBinning(&fHist, bins, bins, &binsa); + } fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos")); @@ -399,7 +421,7 @@ // the same number of bins than for on-data // -void MHFalseSource::ProjectOff(TH2D *h2, TH2D *all) -{ - TAxis &axe = *fHist.GetZaxis(); +void MHFalseSource::ProjectOff(const TH3D &src, TH2D *h2, TH2D *all) +{ + TAxis &axe = *src.GetZaxis(); // Find range to cut (left edge and width) @@ -413,5 +435,5 @@ // Get projection for range - TH2D *p = (TH2D*)fHist.Project3D("yx_off"); + TH2D *p = (TH2D*)src.Project3D("yx_off"); // Reset range @@ -435,7 +457,7 @@ // range (0, fAlphaCut) // -void MHFalseSource::ProjectOn(TH2D *h3, TH2D *all) -{ - TAxis &axe = *fHist.GetZaxis(); +void MHFalseSource::ProjectOn(const TH3D &src, TH2D *h3, TH2D *all) +{ + TAxis &axe = *src.GetZaxis(); // Find range to cut @@ -445,5 +467,5 @@ // Get projection for range - TH2D *p = (TH2D*)fHist.Project3D("yx_on"); + TH2D *p = (TH2D*)src.Project3D("yx_on"); // Reset range @@ -482,4 +504,23 @@ // Set Minimum as minimum value Greater Than 0 h3->SetMinimum(GetMinimumGT(*h3)); +} + +void MHFalseSource::ProjectOnOff(TH2D *h2, TH2D *h0) +{ + ProjectOn(*fHistOff, h2, h0); + + TH2D h; + MH::SetBinning(&h, h2); + + // Divide by number of entries in off region (of off-data) + ProjectOff(*fHistOff, &h, h0); + h2->Divide(&h); + + // Multiply by number of entries in off region (of on-data) + ProjectOff(fHist, &h, h0); + h2->Multiply(&h); + + // Recalculate Minimum + h2->SetMinimum(GetMinimumGT(*h2)); } @@ -518,10 +559,15 @@ padsave->GetPad(1)->cd(1); if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on"))) - ProjectOn(h3, h0); + ProjectOn(fHist, h3, h0); // Update projection of off-events padsave->GetPad(1)->cd(3); if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off"))) - ProjectOff(h2, h0); + { + if (!fHistOff) + ProjectOff(fHist, h2, h0); + else + ProjectOnOff(h2, h0); + } padsave->GetPad(2)->cd(2); @@ -579,4 +625,16 @@ TH1 *h = fHist.ProjectionZ("Alpha_z", maxx, maxx, maxy, maxy); h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s)); + + TH1D *h0=0; + if ((h0 = (TH1D*)gPad->FindObject("AlphaOff_z"))) + { + fHistOff->ProjectionZ("AlphaOff_z", maxx, maxx, maxy, maxy); + + const Int_t f = h0->GetXaxis()->FindFixBin(fBgMean-fAlphaCut/2); + const Int_t l = h0->GetXaxis()->FindFixBin(fAlphaCut)+f-1; + h0->Scale(h1->Integral(f, l)/h0->Integral(f, l)); + //h0->Scale(h1->GetEntries()/h0->GetEntries()); + + } } } @@ -620,10 +678,13 @@ pad->Divide(1, 2, 0, 0.03); - TObject *catalog = GetCatalog(); +// TObject *catalog = GetCatalog(); // Initialize upper part pad->cd(1); // Make sure that the catalog is deleted only once - gROOT->GetListOfCleanups()->Add(gPad); + // Normally this is not done by root, because it is not necessary... + // but in this case it is necessary, because the catalog is + // deleted by the first pad and the second one tries to do the same! +// gROOT->GetListOfCleanups()->Add(gPad); gPad->SetBorderMode(0); gPad->Divide(3, 1); @@ -640,5 +701,5 @@ h3->Draw("colz"); h3->SetBit(kCanDelete); - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); // PAD #2 @@ -654,11 +715,22 @@ h4->Draw("colz"); h4->SetBit(kCanDelete); - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); // PAD #3 pad->GetPad(1)->cd(3); gPad->SetBorderMode(0); - fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); - TH1 *h2 = fHist.Project3D("yx_off"); + TH1 *h2 = 0; + if (fHistOff) + { + fHistOff->GetZaxis()->SetRangeUser(0,fAlphaCut); + h2 = fHistOff->Project3D("yx_off"); + fHistOff->GetZaxis()->SetRange(0,9999); + } + else + { + fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); + h2 = fHist.Project3D("yx_off"); + fHist.GetZaxis()->SetRange(0,9999); + } h2->SetDirectory(NULL); h2->SetXTitle(fHist.GetXaxis()->GetTitle()); @@ -666,10 +738,10 @@ h2->Draw("colz"); h2->SetBit(kCanDelete); - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); // Initialize lower part pad->cd(2); // Make sure that the catalog is deleted only once - gROOT->GetListOfCleanups()->Add(gPad); +// gROOT->GetListOfCleanups()->Add(gPad); gPad->SetBorderMode(0); gPad->Divide(3, 1); @@ -685,4 +757,17 @@ h1->Draw(); h1->SetBit(kCanDelete); + if (fHistOff) + { + h1->SetLineColor(kGreen); + + h1 = fHistOff->ProjectionZ("AlphaOff_z"); + h1->SetDirectory(NULL); + h1->SetTitle("Distribution of \\alpha"); + h1->SetXTitle(fHistOff->GetZaxis()->GetTitle()); + h1->SetYTitle("Counts"); + h1->Draw("same"); + h1->SetBit(kCanDelete); + h1->SetLineColor(kRed); + } // PAD #5 @@ -695,5 +780,5 @@ h5->Draw("colz"); h5->SetBit(kCanDelete); - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); // PAD #6 @@ -706,5 +791,5 @@ h0->Draw("colz"); h0->SetBit(kCanDelete); - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); } @@ -785,5 +870,5 @@ void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom) { - TObject *catalog = GetCatalog(); +// TObject *catalog = GetCatalog(); TH1D h0a("A", "", 50, 0, 4000); @@ -868,5 +953,5 @@ fit.SetPolynomOrder(polynom); - TH1D *h=0; + TH1D *h=0, *hoff=0; for (int ix=1; ix<=nx; ix++) for (int iy=1; iy<=ny; iy++) @@ -879,8 +964,13 @@ continue; - h->Scale(entries/h->GetEntries()); - if (!fit.Fit(*h)) + if (fHistOff) + { + hoff = fHistOff->ProjectionZ("AlphaFitOff", ix, ix, iy, iy); + hoff->Scale(entries/hoff->GetEntries()); + } + + if (!fit.Fit(*h, hoff)) continue; @@ -894,5 +984,5 @@ h1.Fill(fit.GetGausMu()); // mu h2.Fill(fit.GetGausSigma()); // sigma-gaus - if (polynom>1) + if (polynom>1 && !fHistOff) h3.Fill(fit.GetCoefficient(5)); h4b.Fill(fit.GetChiSqSignal()); @@ -945,5 +1035,5 @@ hists->Draw("colz"); hists->SetBit(kCanDelete); - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); c->cd(2); gPad->SetBorderMode(0); @@ -976,10 +1066,10 @@ - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); c->cd(3); gPad->SetBorderMode(0); histb->Draw("colz"); histb->SetBit(kCanDelete); - catalog->Draw("mirror same *"); +// catalog->Draw("mirror same *"); c->cd(4); gPad->Divide(1,3, 0, 0); @@ -1018,16 +1108,48 @@ hist->GetYaxis()->GetBinCenter(maxy), maxs); - TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy); - result->Scale(entries/h->GetEntries()); - - result->SetDirectory(NULL); - result->SetNameTitle("Result \\alpha", title); - result->SetBit(kCanDelete); - result->SetXTitle("\\alpha [\\circ]"); - result->SetYTitle("Counts"); - result->Draw(); - - TF1 f1("f1", Form("gaus(0) + pol%d(3)", polynom), 0, 90); - TF1 f2("f2", Form("gaus(0) + pol%d(3)", polynom), 0, 90); + h = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy); + h->Scale(entries/h->GetEntries()); + + h->SetDirectory(NULL); + h->SetNameTitle("Result \\alpha", title); + h->SetBit(kCanDelete); + h->SetXTitle("\\alpha [\\circ]"); + h->SetYTitle("Counts"); + h->Draw(); + + if (fHistOff) + { + h->SetLineColor(kGreen); + + TH1D *hof=fHistOff->ProjectionZ("AlphaFitOff", maxx, maxx, maxy, maxy); + hof->Scale(entries/hof->GetEntries()); + + fit.Scale(*(TH1D*)hof, *(TH1D*)h); + + hof->SetLineColor(kRed); + hof->SetDirectory(NULL); + hof->SetNameTitle("Result \\alpha", title); + hof->SetBit(kCanDelete); + hof->SetXTitle("\\alpha [\\circ]"); + hof->SetYTitle("Counts"); + hof->Draw("same"); + + TH1D *diff = (TH1D*)h->Clone("AlphaFitOnOff"); + diff->Add(hof, -1); + diff->SetLineColor(kBlue); + diff->SetNameTitle("Result On-Off \\alpha", title); + diff->SetBit(kCanDelete); + diff->SetXTitle("\\alpha [\\circ]"); + diff->SetYTitle("Counts"); + diff->Draw("same"); + + h->SetMinimum(diff->GetMinimum()<0 ? diff->GetMinimum()*1.2 : 0); + + TLine l; + l.DrawLine(0, 0, 90, 0); + } + + TF1 f1("f1", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90); + TF1 f2("f2", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90); f1.SetParameters(maxpar.GetArray()); f2.SetParameters(maxpar.GetArray()); @@ -1053,5 +1175,5 @@ leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11); if (polynom>1) - leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11); + leg->AddText(0.60, 0.14, Form("c=%.2f", !fHistOff?maxpar[5]:0))->SetTextAlign(11); leg->SetBit(kCanDelete); leg->Draw(); Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 5775) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 5776) @@ -39,4 +39,6 @@ TH3D fHist; // Alpha vs. x and y + const TH3D *fHistOff; + Double_t fRa; Double_t fDec; @@ -46,6 +48,7 @@ void ProjectAll(TH2D *h); - void ProjectOff(TH2D *h, TH2D *all); - void ProjectOn(TH2D *h, TH2D *all); + void ProjectOff(const TH3D &src, TH2D *h, TH2D *all); + void ProjectOn(const TH3D &src, TH2D *h, TH2D *all); + void ProjectOnOff(TH2D *h, TH2D *all); TObject *GetCatalog(); @@ -77,4 +80,14 @@ void SetBgMeanMinus5() { SetBgMean(fBgMean-5); } //*MENU* + void SetOffData(const MHFalseSource &fs) { + fHistOff = &fs.fHist; + fMinDist = fs.fMinDist; + fMaxDist = fs.fMaxDist; + fMaxDW = fs.fMaxDW; + fMinDW = fs.fMinDW; + fAlphaCut = fs.fAlphaCut; + fBgMean = fs.fBgMean; + } + void Paint(Option_t *opt=""); void Draw(Option_t *option=""); Index: trunk/MagicSoft/Mars/mraw/MRawFileRead.cc =================================================================== --- trunk/MagicSoft/Mars/mraw/MRawFileRead.cc (revision 5775) +++ trunk/MagicSoft/Mars/mraw/MRawFileRead.cc (revision 5776) @@ -165,5 +165,5 @@ // This opens the next file in the list and deletes its name from the list. // -Bool_t MRawFileRead::OpenNextFile() +Int_t MRawFileRead::OpenNextFile() { // @@ -209,10 +209,10 @@ // if (!fRawRunHeader->ReadEvt(*fIn)) - return kFALSE; + return kERROR; if (!(*fIn)) { *fLog << err << "Error: Accessing file '" << name << "'" << endl; - return kFALSE; + return kERROR; } @@ -242,5 +242,5 @@ { *fLog << err << dbginf << "MTaskList not found... abort." << endl; - return kFALSE; + return kERROR; } @@ -270,5 +270,5 @@ fNumFile=0; fNumEvents=0; - return OpenNextFile(); + return OpenNextFile()==kTRUE; } @@ -348,6 +348,7 @@ // If an event could not be read from file try to open new file // - if (!OpenNextFile()) - return kFALSE; + const Int_t rc = OpenNextFile(); + if (rc!=kTRUE) + return rc; } return kTRUE; Index: trunk/MagicSoft/Mars/mraw/MRawFileRead.h =================================================================== --- trunk/MagicSoft/Mars/mraw/MRawFileRead.h (revision 5775) +++ trunk/MagicSoft/Mars/mraw/MRawFileRead.h (revision 5776) @@ -22,5 +22,5 @@ UInt_t fInterleave; - Bool_t OpenNextFile(); + Int_t OpenNextFile(); Int_t PreProcess(MParList *pList); Index: trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc =================================================================== --- trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc (revision 5775) +++ trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc (revision 5776) @@ -83,5 +83,5 @@ const UShort_t MRawRunHeader::kMagicNumber = 0xc0c0; -const Byte_t MRawRunHeader::kMaxFormatVersion = 5; +const Byte_t MRawRunHeader::kMaxFormatVersion = 2; // -------------------------------------------------------------------------- @@ -139,5 +139,5 @@ if (fMagicNumber != kMagicNumber && fMagicNumber != kMagicNumber+1) { - *fLog << err << "Error: Wrong Magic Number (0x" << hex << fMagicNumber << "): Not a Magic File!" << endl; + *fLog << err << "ERROR - Wrong Magic Number (0x" << hex << fMagicNumber << "): Not a Magic File!" << endl; return kFALSE; } @@ -151,5 +151,5 @@ if (fFormatVersion>kMaxFormatVersion) { - *fLog << err << "WARNING - File vormat V" << fFormatVersion << " not implemented!" << endl; + *fLog << err << "ERROR - File vormat V" << fFormatVersion << " not implemented!" << endl; return kFALSE; }