Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 8673) +++ /trunk/MagicSoft/Mars/Changelog (revision 8674) @@ -18,4 +18,44 @@ -*-*- END OF LINE -*-*- + + 2007/08/18 Thomas Bretz + + * sponde.cc: + - added new option --force-theta + + * mbase/MEnv.[h,cc], mjobs/MSequence.[h,cc], mjobs/MDataSet.[h,cc]: + - GetName noe returns only the filename not the whole path. The + old behaviour made it impossible to access the container from + the file. + - GetRcName now returns the whole path/name. + - Print now outputs also path and file-name + + * mfileio/MReadTree.cc, mfileio/MWriteRootFile.cc: + - fixed typos in comments + + * mhflux/MHCollectionArea.cc: + - reset fCorsikaVersion to 0 in PreProcess + - print old and new Cosika version if mismatch is found + + * mhflux/MMcSpectrumWeight.cc: + - replaced the %.16f by %.16e. This is more accurate in cases + with high exponents + - added some sample/test code to weight the Zenith Angle + according to the sin-distribution produced by Corsika. + Currently not in use + + * mjobs/MJSpectrum.cc: + - removed the simple/accurate mode. There is now reason why + the previous "accurate"-mode should be more accurate at all. + It is only slower + - Reading the OriginalMC tree now is done such that the + events are properly weighted by the production area. This + allowes to use different impact paramters from dfferent files. + - a check has been implemented which compared the zenith angle + distribution of the data and the resulting monte carlo data. + Execution of the program can be forced with a new option. + - write more information to output file. + + 2007/08/17 Thomas Bretz Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 8673) +++ /trunk/MagicSoft/Mars/NEWS (revision 8674) @@ -221,4 +221,27 @@ contents where exchanged.... fixed. + - sponde: the so called "accurate"-mode has been removed. It didn't + give any improvement in accuracy, only decreased execution speed. + + - sponde: now checks whether the theta distribution of your on-data + and the theta-distribution of your Monte Carlo sample (after + weighting) fits. If it doesn't fit properly (eg. the Monte Carlo + sample is incomplete) execution is stopped. Execution can be forced + using the new option --force-theta. Use this option with care! + + - sponde: Proper collection areas can now be constructed also from + Monte Carlo samples generated with different maximum impact + parameters. Note that in previous version you neither got + a warning or failure, nor was there any obvious sign that the + collection area was overestimated due to usage of files with + different maximum impact parameters. + + - sponde: the output file now contains more information about + the spectrum (eg. the full 2D collection area histogram). + Note, that this information can only be written to the file + if it is stored automatically via command line argument. + If you only store the status display from within the display + the information is lost. + Index: /trunk/MagicSoft/Mars/mbase/MEnv.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MEnv.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mbase/MEnv.cc (revision 8674) @@ -71,4 +71,14 @@ if (GetEntries()<=0 && !fname.IsNull() && fname!=name) ReadFile(fname, kEnvLocal);; +} + +//--------------------------------------------------------------------------- +// +// Make sure that the name used for writing doesn't contain a full path +// +const char *MEnv::GetName() const +{ + const char *pos = strrchr(GetRcName(), '/'); + return pos>0 ? pos+1 : GetRcName(); } @@ -847,4 +857,15 @@ //--------------------------------------------------------------------------- // +// Add name and full path to output +// +void MEnv::PrintEnv(EEnvLevel level) const +{ + cout << "# Path: " << GetRcName() << endl; + cout << "# Name: " << GetName() << endl; + TEnv::PrintEnv(level); +} + +//--------------------------------------------------------------------------- +// // Print resources which have never been touched (for debugging) // Index: /trunk/MagicSoft/Mars/mbase/MEnv.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MEnv.h (revision 8673) +++ /trunk/MagicSoft/Mars/mbase/MEnv.h (revision 8674) @@ -38,5 +38,5 @@ const char *GetValue(const char *name, const char *dflt); - const char *GetName() const { return GetRcName(); } + const char *GetName() const { return strrchr(GetRcName(), '/')?strrchr(GetRcName(), '/')+1:GetRcName(); } Int_t GetColor(const char *name, Int_t dftl); @@ -67,6 +67,7 @@ void AddEnv(const TEnv &env, Bool_t overwrite=kTRUE); + void PrintEnv(EEnvLevel level = kEnvAll) const; void Print(Option_t *option) const { TEnv::Print(option); } - void Print() const { TEnv::PrintEnv(kEnvLocal); } //*MENU* + void Print() const { PrintEnv(kEnvLocal); } //*MENU* void PrintUntouched() const; Index: /trunk/MagicSoft/Mars/mfileio/MReadTree.cc =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 8674) @@ -267,5 +267,5 @@ // list gets a new address. Therefor we reset all the autodel flags. // -// MAKE SURE THAT ALL YOPUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY +// MAKE SURE THAT ALL YOUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY // Bool_t MReadTree::Notify() @@ -295,5 +295,5 @@ // list gets a new address. Therefor we reset all the autodel flags. // - // MAKE SURE THAT ALL YOPUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY + // MAKE SURE THAT ALL YOUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY // TIter NextB(fTree->GetListOfBranches()); Index: /trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc (revision 8674) @@ -35,5 +35,5 @@ // There is a special mode of operation which opens a new file for each new // file read by the reading task (opening the new file is initiated by -// ReInit()) For more details se the corresponding constructor. +// ReInit()) For more details see the corresponding constructor. // // Memory based trees Index: /trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 8674) @@ -183,5 +183,6 @@ MH::SetBinning(&fHistAll, &binst, &binse); - fMcAreaRadius = -1; + fMcAreaRadius = -1; + fCorsikaVersion = 0; return kTRUE; @@ -205,5 +206,9 @@ if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion()) - *fLog << warn << "Warning - Read files have different Corsika versions..." << endl; + { + *fLog << warn; + *fLog << "Warning - Read files have different Corsika versions..." << endl; + *fLog << " Last file=" << fCorsikaVersion << " New file=" << runheader->GetCorsikaVersion() << endl; + } fCorsikaVersion = runheader->GetCorsikaVersion(); @@ -285,5 +290,5 @@ // GetCollectionAreaAbs(), fMcAreaRadius); const TString txt = Form("r_{max}=%.0fm --> A_{max}=%.0fm^{2}", - GetCollectionAreaAbs(), fMcAreaRadius); + fMcAreaRadius, GetCollectionAreaAbs()); TLatex text(0.31, 0.95, txt); Index: /trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc =================================================================== --- /trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 8674) @@ -254,5 +254,5 @@ // are equal only fNorm is returned. // -// The normalization constant is returned as %.16f +// The normalization constant is returned as %.16e // // Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600)) @@ -261,5 +261,5 @@ { if (GetFormulaSpecOld()==GetFormulaSpecNew()) - return Form("%.16f", fNorm); + return Form("%.16e", fNorm); const Double_t iold = fNormEnergy<0 ? GetSpecOldIntegral() : CalcSpecOld(fNormEnergy); @@ -268,5 +268,5 @@ const Double_t norm = fNorm*iold/inew; - return Form("%.16f*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data()); + return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data()); } @@ -447,4 +447,18 @@ if (fWeightsZd) { + /* + const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd()); + + Double_t tmin = fWeightsZd->GetXaxis()->GetBinLowEdge(i) *TMath::DegToRad(); + Double_t tmax = fWeightsZd->GetXaxis()->GetBinLowEdge(i+1)*TMath::DegToRad(); + Double_t wdth = fWeightsZd->GetXaxis()->GetBinWidth(i) *TMath::DegToRad(); + + Double_t cont = fWeightsZd->GetBinContent(i); + + Double_t integ = cos(tmin)-cos(tmax); + + w = sin(tmax)*cont/integ*wdth; + */ + const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd()); w = fWeightsZd->GetBinContent(i); Index: /trunk/MagicSoft/Mars/mjobs/MDataSet.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MDataSet.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mjobs/MDataSet.cc (revision 8674) @@ -274,4 +274,15 @@ } +//--------------------------------------------------------------------------- +// +// Make sure that the name used for writing doesn't contain a full path +// +const char *MDataSet::GetName() const +{ + const char *pos = strrchr(GetRcName(), '/'); + return pos>0 ? pos+1 : GetRcName(); +} + + // -------------------------------------------------------------------------- // @@ -299,4 +310,7 @@ return; } + gLog << "# Path: " << GetRcName() << endl; + gLog << "# Name: " << GetName() << endl; + gLog << endl; gLog << "AnalysisNumber: " << fNumAnalysis << endl << endl; Index: /trunk/MagicSoft/Mars/mjobs/MDataSet.h =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MDataSet.h (revision 8673) +++ /trunk/MagicSoft/Mars/mjobs/MDataSet.h (revision 8674) @@ -67,4 +67,7 @@ MDataSet(Int_t num=(UInt_t)-1) : fNumAnalysis(num) { } MDataSet(const char *fname, TString sequences="", TString data=""); + + const char *GetName() const; + const char *GetRcName() const { return fName; } void Copy(TObject &obj) const @@ -153,5 +156,6 @@ // TObject - void Print(Option_t *o="") const; //*MENU* + void Print(Option_t *o) const; + void Print() const { Print(); } //*MENU* ClassDef(MDataSet, 2) Index: /trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8674) @@ -47,4 +47,6 @@ #include "MLogManip.h" +#include "MDirIter.h" + #include "MStatusArray.h" #include "MStatusDisplay.h" @@ -88,5 +90,6 @@ : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0), fEstimateEnergy(0), fCalcHadronness(0), fRefill(kFALSE), - fSimpleMode(kTRUE), fRawMc(kFALSE), fNoThetaWeights(kFALSE) + /*fSimpleMode(kTRUE),*/ fForceTheta(kFALSE), + fRawMc(kFALSE), fNoThetaWeights(kFALSE) { fName = name ? name : "MJSpectrum"; @@ -291,20 +294,29 @@ *fLog << inf << "Using weights: " << w << endl; - *fLog << "Please stand by, this may take a while..." << flush; + *fLog << "Please stand by, this may take a while..." << endl; if (fDisplay) - fDisplay->SetStatusLine1("Compiling MC distribution..."); + fDisplay->SetStatusLine1("Getting maximum impact..."); // Create chain - TChain chain("OriginalMC"); + *fLog << "Getting files... " << flush; + TChain chain("RunHeaders"); if (!set.AddFilesOn(chain)) return kFALSE; + *fLog << "done. " << endl; + + *fLog << "Getting maximum impact... " << flush; + const Double_t impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax"); + *fLog << "found " << impactmax/100 << "m" << endl; + + *fLog << "Getting files... " << flush; + MDirIter iter; + if (!set.AddFilesOn(iter)) + return kFALSE; + *fLog << "done. " << endl; // Prepare histogram h.Reset(); h.Sumw2(); - - // Fill histogram from chain - h.SetDirectory(gROOT); if (h.InheritsFrom(TH2::Class())) { @@ -312,6 +324,5 @@ h.SetXTitle("\\Theta [\\circ]"); h.SetYTitle("E [GeV]"); - h.SetZTitle("Counts"); - chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff"); + h.SetZTitle(Form("Counts normalized to r=%.0fm", impactmax/100)); } else @@ -319,10 +330,78 @@ h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)"); h.SetXTitle("\\Theta [\\circ]"); - h.SetYTitle("Counts"); - chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff"); + h.SetYTitle(Form("Counts normalized to r=%.0fm", impactmax/100)); } h.SetDirectory(0); - *fLog << "done." << endl; + const TString cont = h.InheritsFrom(TH2::Class()) ? + "MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC" : + "MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC"; + + if (fDisplay) + fDisplay->SetStatusLine1("Reading OriginalMC"); + + TH1 *hfill = (TH1*)h.Clone(h.InheritsFrom(TH2::Class())?"ThetaEMC":"ThetaMC"); + hfill->SetDirectory(0); + + TString fname; + while (1) + { + fname = iter.Next(); + if (fname.IsNull()) + break; + + TFile file(fname); + if (file.IsZombie()) + { + *fLog << err << "ERROR - Couldn't open file " << fname << endl; + return kFALSE; + } + + TTree *rh = dynamic_cast(file.Get("RunHeaders")); + if (!rh) + { + *fLog << err << "ERROR - File " << fname << " doesn't contain tree RunHeaders." << endl; + return kFALSE; + } + + if (rh->GetEntries()!=1) + { + *fLog << err << "ERROR - RunHeaders of " << fname << " doesn't contain exactly one entry." << endl; + return kFALSE; + } + + TTree *t = dynamic_cast(file.Get("OriginalMC")); + if (!t) + { + *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl; + return kFALSE; + } + + *fLog << inf << "Reading OriginalMC of " << fname << endl; + + if (t->GetEntries()==0) + { + *fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl; + continue; + } + + // Why does it crash if closed within loop (no update?) + //if (fDisplay) + // fDisplay->SetStatusLine2(fname); + + // Normalize by deviding through production area + const Double_t impact = rh->GetMaximum("MMcRunHeader.fImpactMax"); + const Double_t scale = impactmax/impact; + const TString weights = Form("%.16e*(%s)", scale*scale, w.Data()); + + // Fill histogram from tree + hfill->SetDirectory(&file); + t->Draw(cont, weights, "goff"); + hfill->SetDirectory(0); + h.Add(hfill); + } + delete hfill; + + *fLog << "Total number of events from file: " << h.GetEntries() << endl; if (fDisplay) fDisplay->SetStatusLine2("done."); @@ -333,5 +412,5 @@ *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl; - return h.GetEntries()>0; + return kFALSE; } @@ -412,6 +491,7 @@ { *fLog << err; - *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap"; - *fLog << " with the the Zenith Angle distribution of your observation."; + *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap" << endl; + *fLog << " with the Zenith Angle distribution of your observation." << endl; + *fLog << " Maybe the energy binning is not defined or wrong." << endl; theta->SetLineColor(kRed); return kFALSE;; @@ -424,5 +504,6 @@ // Compare both histograms - const Double_t prob = proj.Chi2Test(theta, ""); + *fLog << inf << "Comparing theta distributions for data and MCs." << endl; + const Double_t prob = proj.Chi2Test(theta, "P"); if (prob==1) return kTRUE; @@ -431,25 +512,26 @@ { *fLog << inf; - *fLog << "The Zenith Angle distribution of your Monte Carlos fits well"; - *fLog << "with the the Zenith Angle distribution of your observation."; + *fLog << "The Zenith Angle distribution of your Monte Carlos fits well" << endl; + *fLog << "with the Zenith Angle distribution of your observation." << endl; + *fLog << "The probability for identical Theta distributions is " << prob << endl; return kTRUE; } - // <0.01: passen gar nicht - // ==1: passen perfekt - // >0.99: passer sehr gut - if (prob<0.01) { *fLog << err; - *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit"; - *fLog << " with the the Zenith Angle distribution of your observation."; + *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit" << endl; + *fLog << " with the Zenith Angle distribution of your observation." << endl; + *fLog << " The probability for identical Theta distributions is " << prob << endl; + if (!fForceTheta) + *fLog << " To force processing use --force-theta (with care!)" << endl; theta->SetLineColor(kRed); - return kFALSE; + return fForceTheta; } *fLog << warn; - *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well"; - *fLog << " with the the Zenith Angle distribution of your observation."; + *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well" << endl; + *fLog << " with the Zenith Angle distribution of your observation." << endl; + *fLog << " The probability for identical Theta distributions is " << prob << endl; return kTRUE; } @@ -1167,5 +1249,5 @@ TH2D hist; MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy"); - if (fSimpleMode) + //if (fSimpleMode) { // Now we read the MC distribution as produced by corsika again @@ -1189,5 +1271,5 @@ } } - } + }/* else { @@ -1198,7 +1280,7 @@ return kFALSE; weight.SetNameMcEvt(); - } - - if (!DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist())) + }*/ + + if (!DisplayResult(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/)) return kFALSE; @@ -1245,6 +1327,6 @@ MHCollectionArea area0("TriggerArea"); MHCollectionArea area1; - area0.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist()); - area1.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist()); + area0.SetHistAll(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/); + area1.SetHistAll(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/); MHEnergyEst hest; @@ -1400,4 +1482,8 @@ TObjArray cont; cont.Add((TObject*)GetEnv()); + cont.Add(&area0); + cont.Add(&area1); + cont.Add(&hest); + cont.Add(const_cast(GetEnv())); if (fDisplay) Index: /trunk/MagicSoft/Mars/mjobs/MSequence.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MSequence.cc (revision 8673) +++ /trunk/MagicSoft/Mars/mjobs/MSequence.cc (revision 8674) @@ -478,4 +478,15 @@ } +//--------------------------------------------------------------------------- +// +// Make sure that the name used for writing doesn't contain a full path +// +const char *MSequence::GetName() const +{ + const char *pos = strrchr(GetRcName(), '/'); + return pos>0 ? pos+1 : GetRcName(); +} + + // -------------------------------------------------------------------------- // @@ -490,4 +501,7 @@ return; } + gLog << "# Path: " << GetRcName() << endl; + gLog << "# Name: " << GetName() << endl; + gLog << endl; gLog << "Sequence: " << fSequence << endl; if (fMonteCarlo) Index: /trunk/MagicSoft/Mars/mjobs/MSequence.h =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MSequence.h (revision 8673) +++ /trunk/MagicSoft/Mars/mjobs/MSequence.h (revision 8674) @@ -78,5 +78,9 @@ // TObject - void Print(Option_t *o="") const; //*MENU* + void Print(Option_t *o) const; + void Print() const { Print(); } //*MENU* + + const char *GetName() const; + const char *GetRcName() const { return fName; } // Genaral interface Index: /trunk/MagicSoft/Mars/sponde.cc =================================================================== --- /trunk/MagicSoft/Mars/sponde.cc (revision 8673) +++ /trunk/MagicSoft/Mars/sponde.cc (revision 8674) @@ -68,4 +68,5 @@ gLog << " --print-files Print Files taken from Sequences ('+' found, '-' missing)" << endl; gLog << " --config=sponde.rc Resource file [default=sponde.rc]" << endl; + gLog << " --force-theta Force execution even with non-fitting theta distributions." << endl; gLog << endl; gLog << " --version, -V Show startup message with version number" << endl; @@ -122,4 +123,5 @@ const Bool_t kRefill = arg.HasOnlyAndRemove("--refill"); // const Bool_t kSimple = !arg.HasOnlyAndRemove("--accurate"); + const Bool_t kForceTheta = arg.HasOnlyAndRemove("--force-theta"); const Bool_t kRawMc = arg.HasOnlyAndRemove("--raw-mc"); @@ -223,4 +225,5 @@ job.SetPathIn(kInfile); + job.ForceTheta(kForceTheta); job.EnableRefilling(kRefill); // job.EnableSimpleMode(kSimple);