Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 8881) +++ /trunk/MagicSoft/Mars/Changelog (revision 8882) @@ -19,8 +19,44 @@ -*-*- END OF LINE -*-*- + 2008/03/18 Thomas Bretz + + * sponde.cc: + - added new option "--force-runtime" + + * mbase/MEnv.h: + - added WriteFile to context menu + + * mjobs/MJSpectrum.[h,cc]: + - added a new option to force using the runtime instead of the + effective observation time (this might bw wrong for very + short datasets) + - added a check if the effective observation time is out of + the histogram range... print a warning if so and include + the overflow bins into the eff. obs time + - added an estimated sensitivity curve for high and low za + to the spectrum plots + - added description text for 1553 and crab spectrum + - write out the MC events after cuts including their weights + - do not fit at 1TeV but 500GeV instead + + * mjobs/MJob.cc: + - check in WriteContainer whether the file is already open + + * mpointing/MPointingDevCalc.cc: + - added some more comments + + + 2008/03/14 Daniel Hoehne * datacenter/macros/filldotrun.C: - inserted new arehucas version + + + + 2008/03/04 Thomas Bretz + + * condor/program.submit, condor/macro.submit, condor/script.submit: + - added Index: /trunk/MagicSoft/Mars/mbase/MEnv.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MEnv.h (revision 8881) +++ /trunk/MagicSoft/Mars/mbase/MEnv.h (revision 8882) @@ -73,4 +73,7 @@ Int_t ReadFile(const char *fname, EEnvLevel level); + Int_t WriteFile(const char *filename, EEnvLevel level) { return TEnv::WriteFile(filename, level); } + Int_t WriteFile(const char *filename) { return WriteFile(filename, kEnvLocal); } //*MENU* + void PrintEnv(EEnvLevel level = kEnvAll) const; void Print(Option_t *option) const { TEnv::Print(option); } Index: /trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8881) +++ /trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8882) @@ -38,4 +38,5 @@ #include #include +#include #include #include @@ -56,4 +57,5 @@ #include "MHn.h" #include "MBinning.h" +#include "MParameters.h" #include "MDataSet.h" #include "MMcCorsikaRunHeader.h" @@ -79,4 +81,5 @@ #include "MFDataMember.h" #include "MContinue.h" +#include "MWriteRootFile.h" ClassImp(MJSpectrum); @@ -239,7 +242,8 @@ } - TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime"); - TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); - if (!vstime || !size) + TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime"); + TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); + TH1D *time = (TH1D*)arr.FindObjectInCanvas("ExcessTime", "TH1D", "Hist"); + if (!vstime || !size || !time) return -1; @@ -280,5 +284,24 @@ return -1; - return vstime->Integral(); + if (vstime->GetBinContent(0)>0) + { + *fLog << err << "ERROR - Undeflow bin of OnTime histogram not empty as it ought to be." << endl; + return -1; + } + + const Double_t ofl = vstime->GetBinContent(vstime->GetNbinsX()); + const Double_t eff = vstime->Integral()+ofl; + if (ofl>0) + *fLog << warn << "WARNING - " << Form("%.1f%% (%.0fs)", 100*ofl/eff, ofl) << " of the eff. observation time is out of histogram range." << endl; + + const Double_t obs = time->GetXaxis()->GetXmax()-time->GetXaxis()->GetXmin(); + if (fForceRunTime) + { + *fLog << inf; + *fLog << "Total run time: " << obs/60 << "min" << endl; + *fLog << "Eff. obs. time: " << eff/60 << "min (" << Form("%.1f%%", 100*eff/obs) << ")" << endl; + } + + return fForceRunTime ? obs : eff; } @@ -479,7 +502,5 @@ // ----- Complete incomplete energy ranges ----- // ----- and apply production area weights ----- - weight.CompleteEnergySpectrum(*hfill, Emin); - - hfill->Scale(scale*scale); + weight.CompleteEnergySpectrum(*hfill, Emin, scale*scale); // Add weighted data from file to final histogram @@ -537,5 +558,5 @@ // Calculate the Probability temp1.Divide(&temp2); - temp1.Scale(1./temp1.Integral()); + temp1.Scale(1./temp1.Integral(1, temp1.GetNbinsX())); // Some cosmetics: Name, Axis, etc. @@ -849,5 +870,5 @@ { const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1); - const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0); + const TString fmt1 = Form("(\\frac{E}{500GeV})^{%%%sf#pm%%%sf}", f0, f0); str = Form(fmt0.Data(), p1/exp, e1/exp, np); @@ -966,12 +987,9 @@ spectrum.SetBit(TH1::kNoStats); - for (int i=0; iSetGridx(); gPad->SetGridy(); + + TGraph gsensl;//("Sensitivity LZA", ""); + TGraph gsensh;//("Sensitivity HZA", ""); + + const Double_t sqrton = TMath::Sqrt(ontime/3600.); + + for (int i=0; iSetSelectedPad(0); + + gsensl.SetLineStyle(5); + gsensl.SetLineColor(14); + gsensl.SetLineWidth(2); + gsensl.DrawClone("C")->SetBit(kCanDelete); + + gsensh.SetLineStyle(5); + gsensh.SetLineColor(14); + gsensh.SetLineWidth(2); + gsensh.DrawClone("C")->SetBit(kCanDelete); // Display dN/dE*E^2 for conveinience @@ -1015,4 +1073,12 @@ spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e); spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) *e*e); + } + + for (int i=0; i(fc).Copy(fa2); fa2.DrawCopy("same"); + + gsensl.DrawClone("C")->SetBit(kCanDelete); + gsensh.DrawClone("C")->SetBit(kCanDelete); // Fit Spectrum @@ -1317,4 +1386,45 @@ hist.InitTitle(";Slope;Disp-Dist [\\circ];"); hist.SetDrawOption("colz profx"); +} + +// -------------------------------------------------------------------------- +// +// Setup write to write: +// container tree optional? +// -------------- ---------- ----------- +// "MHillas" to "Events" +// "MHillasSrc" to "Events" +// "Hadronness" to "Events" yes +// "MEnergyEst" to "Events" yes +// "DataType" to "Events" +// +void MJSpectrum::SetupWriter(MWriteRootFile *write/*, const char *name*/) const +{ + if (!write) + return; + + //write->SetName(name); + write->AddContainer("MHillas", "Events"); + write->AddContainer("MHillasSrc", "Events"); + write->AddContainer("MHillasExt", "Events"); + //write->AddContainer("MPointingPos", "Events"); + write->AddContainer("MHillasSrcAnti", "Events", kFALSE); + write->AddContainer("MImagePar", "Events", kFALSE); + write->AddContainer("MNewImagePar", "Events", kFALSE); + write->AddContainer("MNewImagePar2", "Events", kFALSE); + write->AddContainer("Hadronness", "Events", kFALSE); + write->AddContainer("MSrcPosCam", "Events", kFALSE); + write->AddContainer("MSrcPosAnti", "Events", kFALSE); + write->AddContainer("ThetaSquared", "Events", kFALSE); + write->AddContainer("OpticalAxis", "Events", kFALSE); + write->AddContainer("Disp", "Events", kFALSE); + write->AddContainer("Ghostbuster", "Events", kFALSE); + write->AddContainer("MEnergyEst", "Events", kFALSE); + write->AddContainer("MTime", "Events", kFALSE); + write->AddContainer("MMcEvt", "Events", kFALSE); + write->AddContainer("MWeight", "Events"); + write->AddContainer("DataType", "Events"); + write->AddContainer("RunNumber", "Events"); + write->AddContainer("EvtNumber", "Events"); } @@ -1598,4 +1708,11 @@ taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); + MWriteRootFile write(GetPathOut()); + SetupWriter(&write); + + MParameterI par("DataType"); + plist.AddToList(&par); + par.SetVal(2); + tlist2.AddToList(&read); // If no weighting should be applied but the events should @@ -1621,4 +1738,6 @@ tlist2.AddToList(fCut3); tlist2.AddToList(&taskenv2); + if (!GetPathOut().IsNull()) + tlist2.AddToList(&write); tlist2.AddToList(&fill31); tlist2.AddToList(&fill4); @@ -1670,5 +1789,5 @@ // Spectrum fitted (convert res[1] from TeV to GeV) - TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0])); + TF1 flx("flx", Form("%e*pow(x/500, %f)", res[1]/500, res[0])); // Number of events this spectrum would produce per s and m^2 @@ -1701,4 +1820,5 @@ cont.Add((TObject*)GetEnv()); // const_cast cont.Add((TObject*)&set); // const_cast + cont.Add(plist.FindObject("MAlphaFitter")); cont.Add(&area0); cont.Add(&area1); @@ -1708,5 +1828,5 @@ cont.Add(fDisplay); - if (!WriteContainer(cont, "", "RECREATE")) + if (!WriteContainer(cont, "", GetPathOut().IsNull()?"RECREATE":"UPDATE")) { *fLog << err << GetDescriptor() << ": Writing result failed." << endl; Index: /trunk/MagicSoft/Mars/mjobs/MJSpectrum.h =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 8881) +++ /trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 8882) @@ -19,4 +19,5 @@ class MAlphaFitter; class MStatusArray; +class MWriteRootFile; class MHCollectionArea; class MMcSpectrumWeight; @@ -36,4 +37,5 @@ Bool_t fForceTheta; + Bool_t fForceRunTime; // Setup Histograms @@ -50,4 +52,7 @@ Bool_t Refill(MParList &plist, TH1D &h) /*const*/; Bool_t InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const; + + // Write output + void SetupWriter(MWriteRootFile *write/*, const char *name*/) const; // Display Output @@ -67,5 +72,6 @@ Bool_t Process(const MDataSet &set); - void ForceTheta(Bool_t b=kTRUE) { fForceTheta=b; } + void ForceTheta(Bool_t b=kTRUE) { fForceTheta=b; } + void ForceRunTime(Bool_t b=kTRUE) { fForceRunTime=b; } void SetEnergyEstimator(const MTask *task); Index: /trunk/MagicSoft/Mars/mjobs/MJob.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJob.cc (revision 8881) +++ /trunk/MagicSoft/Mars/mjobs/MJob.cc (revision 8882) @@ -431,6 +431,20 @@ title += fName; + // In case the update-option is selected check whether + // the file is already open + if (TString(option).Contains("update", TString::kIgnoreCase)) + { + TFile *file = dynamic_cast(gROOT->GetListOfFiles()->FindObject(oname)); + if (file && file->IsOpen() && !file->IsZombie()) + { + *fLog << inf << "Open file found." << endl; + file->cd(); + return WriteContainer(cont); + } + } + + // Open a new file with the defined option for writing TFile file(oname, option, title, compr); - if (!file.IsOpen()) + if (!file.IsOpen() || file.IsZombie()) { *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl; Index: /trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc (revision 8881) +++ /trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc (revision 8882) @@ -144,5 +144,6 @@ // New pointing models have been installed (if the pointing model // is different, than the previous one it might mean, that also -// the starguider calibration is different.) +// the starguider calibration is different.) The date only means +// day-time (noon) at which the model has been changed. // // 29. Apr. 2004 ~25800 @@ -156,5 +157,5 @@ // 17. Oct. 2006 ~103130 // 17. Jun. 2007 ~248193 -// 18. Oct. 2007 +// 18. Oct. 2007 ~291039(?) // // From 2.2.2006 beginnig of the night (21:05h, run >=81855) to 24.2.2006 Index: /trunk/MagicSoft/Mars/sponde.cc =================================================================== --- /trunk/MagicSoft/Mars/sponde.cc (revision 8881) +++ /trunk/MagicSoft/Mars/sponde.cc (revision 8882) @@ -51,4 +51,5 @@ gLog << " Operation Mode:" << endl; gLog << " --force-theta Force execution even with non-fitting theta distributions." << endl; + gLog << " --force-runtime Force usage of runtime instead of eff. observation time." << endl; gLog << endl; gLog << " Options:" << endl; @@ -120,4 +121,5 @@ const Bool_t kForceTheta = arg.HasOnlyAndRemove("--force-theta"); + const Bool_t kForceRunTime = arg.HasOnlyAndRemove("--force-runtime"); // @@ -244,4 +246,5 @@ job.ForceTheta(kForceTheta); + job.ForceRunTime(kForceRunTime); if (!job.Process(seq))