Index: trunk/MagicSoft/Mars/.rootrc =================================================================== --- trunk/MagicSoft/Mars/.rootrc (revision 8679) +++ trunk/MagicSoft/Mars/.rootrc (revision 8680) @@ -45,2 +45,16 @@ Root.Html.SearchEngine: http://magic.astro.uni-wuerzburg.de/mars/search.html Root.Html.ViewCVS: http://www.astro.uni-wuerzburg.de/cgi-bin/viewcvs.cgi/ + +############################################################################# +# # +# Some default changes # +# # +############################################################################# + +#Canvas.HighLightColor: 3 +#Canvas.MoveOpaque: 0 +#Canvas.ResizeOpaque: 0 +Canvas.ShowEventStatus: Yes +#Canvas.ShowToolBar: Yes +#Canvas.ShowEditor: Yes +#Canvas.AutoExec: No Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 8679) +++ trunk/MagicSoft/Mars/Changelog (revision 8680) @@ -18,4 +18,48 @@ -*-*- END OF LINE -*-*- + + 2007/08/20 Thomas Bretz + + * .rootrc: + - added some comments about defaults + - set the ShowEventStatus for the canvases to yes as default + + * sponde.cc: + - removed the refill option (it was just a dummy) + - removed the accurate option. It didn't give more accurate + results at all + + * mbase/MStatusDisplay.[h,cc]: + - added an update option to SetProgressBarPosition + + * mhflux/MMcSpectrumWeight.[h,cc]: + - allow to give a container name to GetFormula* + - changed default for spectral slope from -9 to -99 + - allow to give integral range to GetSpec*Integral + - added a new member function CompeleteEnergySpectrum which completes + a simulated spectrum starting at an energy fEnergyMin down to + an energy emin. + - added two options ("new" and "old") to Print function + - do not stop anymore if lower energy boundary changes + + * mjobs/MJSpectrum.[h,cc]: + - cleaned the code in general + - removed fRefill (was not used in the code at all) + - added MJSpectrum to global ListOfCleanups to handle + the display more properly + - removed reading of the first (it was the second!) + MMcCorsikaRunHeader. It is now read for each file individually + - The read monte carlo events are now weighted with the mc + production area (events per area) + - incomplete (to lower energies) spectra are completed + - removed accurate mode, it was not more accurate + - we fit the spectrum now from the first to the last bin + - for comparison crab and 1553 are plotted + - changed the processing such that first the MCs are processed + and then the spectrum is refilled + - now the MC distribution from OriginalMC is read only once + - added new tab showing the basic event distribution + + 2007/08/19 Thomas Bretz Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 8679) +++ trunk/MagicSoft/Mars/NEWS (revision 8680) @@ -223,4 +223,10 @@ give any improvement in accuracy, only decreased execution speed. + - sponde: the so called "simple"-mode has been removed. It didn't + give any improvement in simple. + + - sponde: the so called "refill"-mode has been removed. It was anyhow + not implemented. + - sponde: now checks whether the theta distribution of your on-data and the theta-distribution of your Monte Carlo sample (after @@ -236,4 +242,10 @@ different maximum impact parameters. + - sponde: If MC files with different lower energy limits are used + the primary MC spectrum is artificially completed down to the + lowest energy used at all. WARNING: that this gives correct + collection areas ONLY if none of the events in this region would + survive your cuts at all. + - sponde: the output file now contains more information about the spectrum (eg. the full 2D collection area histogram). @@ -248,4 +260,11 @@ The same information you get from the error bars of the weighted histograms, but this is less intuitive. + + - sponde: The spectrum plots now show the crab- and 1553-spectrum + for comparison. It is not meant to show these curves in + publications, they are only for production. + + - sponde: The OriginalMC tree with the events produced by corsika + is now processed only once Index: trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc (revision 8679) +++ trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc (revision 8680) @@ -602,7 +602,9 @@ // is assumed to be (0,1) // -void MStatusDisplay::SetProgressBarPosition(Float_t p) +void MStatusDisplay::SetProgressBarPosition(Float_t p, Bool_t upd) { fBar->SetPosition(p); + if (upd) + gClient->ProcessEventsFor(fBar); } @@ -685,4 +687,8 @@ // p==NULL means: Take gClient->GetRoot() if not in batch mode // see TGWindow::TGWindow() + + // Make sure that the display is removed via RecursiveRemove + // from whereever possible. + SetBit(kMustCleanup); // Index: trunk/MagicSoft/Mars/mbase/MStatusDisplay.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MStatusDisplay.h (revision 8679) +++ trunk/MagicSoft/Mars/mbase/MStatusDisplay.h (revision 8680) @@ -168,5 +168,5 @@ void SetUpdateTime(Long_t t); - void SetProgressBarPosition(Float_t p); + void SetProgressBarPosition(Float_t p, Bool_t upd=kFALSE); TGProgressBar *GetBar() const { return (TGProgressBar*)fBar; } Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 8679) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 8680) @@ -19,5 +19,5 @@ ! Author(s): Marcos Lopez 10/2003 ! -! Copyright: MAGIC Software Development, 2000-2006 +! Copyright: MAGIC Software Development, 2000-2007 ! ! @@ -40,5 +40,5 @@ // Method: // ------- -// +// // - Corsika spectrun: dN/dE = A * E^(a) // with a = fOldSlope, and A = N/integral{E*de} from ELowLim to EUppLim @@ -72,4 +72,6 @@ // MMcEvt // MMcCorsikaRunHeader +// [MPointingPos] +// [MHillas] // // Output Container: @@ -81,4 +83,6 @@ #include #include +#include +#include #include @@ -109,6 +113,6 @@ fNameMcEvt = "MMcEvt"; - fNewSlope = -9; - fOldSlope = -9; + fNewSlope = -99; + fOldSlope = -99; fEnergyMin = -1; @@ -210,7 +214,7 @@ // The slope is returned as %.3f // -TString MMcSpectrumWeight::GetFormulaSpecOld() const -{ - return Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fOldSlope); +TString MMcSpectrumWeight::GetFormulaSpecOld(const char *name) const +{ + return Form("pow(%s.fEnergy, %.3f)", name, fOldSlope); } @@ -225,9 +229,9 @@ // unchanged. // -TString MMcSpectrumWeight::GetFormulaSpecNew() const -{ - TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fNewSlope) : fFormula.Data(); +TString MMcSpectrumWeight::GetFormulaSpecNew(const char *name) const +{ + TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", name, fNewSlope) : fFormula.Data(); if (!fFormula.IsNull()) - str.ReplaceAll("X", Form("(%s.fEnergy)", fNameMcEvt.Data())); + str.ReplaceAll("X", Form("(%s.fEnergy)", name)); return str; @@ -258,5 +262,5 @@ // Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600)) // -TString MMcSpectrumWeight::GetFormulaWeights() const +TString MMcSpectrumWeight::GetFormulaWeights(const char *name) const { if (GetFormulaSpecOld()==GetFormulaSpecNew()) @@ -268,5 +272,5 @@ const Double_t norm = fNorm*iold/inew; - return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data()); + return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew(name).Data(), GetFormulaSpecOld(name).Data()); } @@ -276,8 +280,8 @@ // GetFormulaSpecNewX() describing the destination spectrum // -Double_t MMcSpectrumWeight::GetSpecNewIntegral() const +Double_t MMcSpectrumWeight::GetSpecNewIntegral(Double_t emin, Double_t emax) const { TF1 funcnew("Dummy", GetFormulaSpecNewX()); - return funcnew.Integral(fEnergyMin, fEnergyMax); + return funcnew.Integral(emin, emax); } @@ -287,8 +291,8 @@ // GetFormulaSpecOldX() describing the simulated spectrum // -Double_t MMcSpectrumWeight::GetSpecOldIntegral() const +Double_t MMcSpectrumWeight::GetSpecOldIntegral(Double_t emin, Double_t emax) const { TF1 funcold("Dummy", GetFormulaSpecOldX()); - return funcold.Integral(fEnergyMin, fEnergyMax); + return funcold.Integral(emin, emax); } @@ -357,4 +361,12 @@ if (fEnergyMax>fEnergyMin && !fAllowChange) { + if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10) + { + *fLog << err; + *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change ("; + *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl; + return kFALSE; + } + if (TMath::Abs(fOldSlope-rh.GetSlopeSpec())>1e-10) { @@ -364,19 +376,14 @@ return kFALSE; } - if (TMath::Abs(fEnergyMin-rh.GetELowLim())>1e-10) - { - *fLog << err; - *fLog << "ERROR - The minimum simulated Monte Carlo energy is not allowed to change "; - *fLog << "(" << fEnergyMin << " --> " << rh.GetELowLim() << ")... abort." << endl; - return kFALSE; - } - if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10) - { - *fLog << err; - *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change ("; - *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl; - return kFALSE; - } - return kTRUE; + + // No change happened + if (TMath::Abs(fEnergyMin-rh.GetELowLim())<=1e-10) + return kTRUE; + + // The lower energy limit has changed + *fLog << warn; + *fLog << "The minimum simulated Monte Carlo energy has changed from "; + *fLog << fEnergyMin << "GeV to " << rh.GetELowLim() << "GeV." << endl; + fEnergyMin = rh.GetELowLim(); } @@ -385,7 +392,7 @@ fEnergyMax = rh.GetEUppLim(); - if (fNewSlope==-9) - { - *fLog << inf << "The new slope of the power law is undefined (-9)... no weighting applied." << endl; + if (fNewSlope==-99) + { + *fLog << inf << "A new slope for the power law has not been defined... no weighting applied." << endl; fNewSlope = fOldSlope; } @@ -403,18 +410,101 @@ // --------------------------------------------------------------------------- // +// completes a simulated spectrum starting at an energy fEnergyMin down to +// an energy emin. +// +// It is assumed that the contents of MMcSpectrumWeight for the new spectrum +// correctly describe the spectrum within the histogram, and fEnergyMin +// and fEnergyMax correctly describe the range. +// +// In the 1D case it is assumed that the x-axis is a zenith angle binning. +// In the 2D case the x-axis is assumed to be zenith angle, the y-axis +// to be energy. +// +void MMcSpectrumWeight::CompleteEnergySpectrum(TH1 &h, Double_t emin) const +{ + if (h.InheritsFrom(TH3::Class())) + { + *fLog << "ERROR - MMcSpctrumWeight::CompleteEnergySpectrum doesn't support TH3." << endl; + return; + } + + if (fEnergyMin <= emin) + return; + + const Double_t norm = GetSpecNewIntegral(); + + if (!h.InheritsFrom(TH2::Class())) + { + h.Scale(GetSpecNewIntegral(emin, fEnergyMax)/norm); + return; + } + + const TAxis &axey = *h.GetYaxis(); + + const Int_t first = axey.FindFixBin(emin); + const Int_t last = axey.FindFixBin(fEnergyMin); + + for (int x=1; x<=h.GetNbinsX(); x++) + { + const Double_t f = h.Integral(x, x, -1, 9999)/norm; + + for (int y=first; y<=last; y++) + { + const Double_t lo = axey.GetBinLowEdge(y) fEnergyMin ? fEnergyMin : axey.GetBinLowEdge(y+1); + + h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi)); + } + } +} + +// --------------------------------------------------------------------------- +// // The current contants are printed // void MMcSpectrumWeight::Print(Option_t *o) const { + const TString opt(o); + + const Bool_t hasnew = opt.Contains("new") || opt.IsNull(); + const Bool_t hasold = opt.Contains("old") || opt.IsNull(); + *fLog << all << GetDescriptor() << endl; - *fLog << " Simulated energy range: " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl; - *fLog << " Simulated spectral slope: " << fOldSlope << endl; - *fLog << " New spectral slope: " << fNewSlope << endl; + + if (hasold) + { + *fLog << " Simulated energy range: " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl; + *fLog << " Simulated spectral slope: "; + if (fOldSlope==-99) + *fLog << "undefined" << endl; + else + *fLog << fOldSlope << endl; + } + if (hasnew) + { + *fLog << " New spectral slope: "; + if (fNewSlope==-99) + *fLog << "undefined" << endl; + else + *fLog << fNewSlope << endl; + } *fLog << " Additional user norm.: " << fNorm << endl; *fLog << " Spectra are normalized: " << (fNormEnergy<0?"by integral":Form("at %.1fGeV", fNormEnergy)) << endl; - *fLog << " Old Spectrum: " << GetFormulaSpecOldX() << " (I=" << GetSpecOldIntegral() << ")" << endl; - *fLog << " New Spectrum: " << GetFormulaSpecNewX() << " (I=" << GetSpecNewIntegral() << ")" << endl; + if (hasold) + { + *fLog << " Old Spectrum: " << GetFormulaSpecOldX(); + if (fEnergyMin>=0 && fEnergyMax>0) + *fLog << " (I=" << GetSpecOldIntegral() << ")"; + *fLog << endl; + } + if (hasnew) + { + *fLog << " New Spectrum: " << GetFormulaSpecNewX(); + if (fEnergyMin>=0 && fEnergyMax>0) + *fLog << " (I=" << GetSpecNewIntegral() << ")"; + *fLog << endl; + } if (fFunc) - *fLog << " Weight function: " << fFunc->GetTitle() << endl; + *fLog << " Weight func: " << fFunc->GetTitle() << endl; } @@ -436,4 +526,31 @@ } +/* + * This could be used to improve the Zd-weighting within a bin. + * Another option is to use more bins, or both. + * Note that it seems unnecessary, because the shape within the + * theta-bins should be similar in data and Monte Carlo... hopefully. + * +void MMcSpectrumWeight::InitZdWeights() +{ + TH2D w(*fWeightsZd); + + for (int i=1; i<=w.GetNbinsX(); i++) + { + const Double_t tmin = w.GetBinLowEdge(i) *TMath::DegToRad(); + const Double_t tmax = w.GetBinLowEdge(i+1)*TMath::DegToRad(); + + const Double_t wdth = tmax-tmin; + const Double_t integ = cos(tmin)-cos(tmax); + + w.SetBinContent(i, w.GetBinContent(i)*wdth/integ); + } + + // const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd()); + // const Double_t theta = fPointing->GetZd()*TMath::DegToRad(); + // w = sin(theta)*w.GetBinContent(i); +} +*/ + // ---------------------------------------------------------------------------- // @@ -447,18 +564,4 @@ 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); @@ -472,4 +575,5 @@ const Double_t e = fMcEvt->GetEnergy(); + fWeight->SetVal(fFunc->Eval(e)*w); Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h (revision 8679) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h (revision 8680) @@ -77,7 +77,10 @@ // Getter - TString GetFormulaSpecOld() const; - TString GetFormulaSpecNew() const; - TString GetFormulaWeights() const; + TString GetFormulaSpecOld(const char *name) const; + TString GetFormulaSpecNew(const char *name) const; + TString GetFormulaWeights(const char *name) const; + TString GetFormulaSpecOld() const { return GetFormulaSpecOld(fNameMcEvt); } + TString GetFormulaSpecNew() const { return GetFormulaSpecNew(fNameMcEvt); } + TString GetFormulaWeights() const { return GetFormulaWeights(fNameMcEvt); } TString GetFormulaSpecOldX() const { return ReplaceX(GetFormulaSpecOld()); } @@ -85,6 +88,9 @@ TString GetFormulaWeightsX() const { return ReplaceX(GetFormulaWeights()); } - Double_t GetSpecNewIntegral() const; - Double_t GetSpecOldIntegral() const; + Double_t GetSpecNewIntegral(Double_t emin, Double_t emax) const; + Double_t GetSpecOldIntegral(Double_t emin, Double_t emax) const; + + Double_t GetSpecNewIntegral() const { return GetSpecNewIntegral(fEnergyMin, fEnergyMax); } + Double_t GetSpecOldIntegral() const { return GetSpecOldIntegral(fEnergyMin, fEnergyMax); } Double_t CalcSpecNew(Double_t e) const; @@ -93,4 +99,7 @@ Double_t GetEnergyMin() const { return fEnergyMin; } Double_t GetEnergyMax() const { return fEnergyMax; } + + // Functions + void CompleteEnergySpectrum(TH1 &h, Double_t emin) const; // TObject Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8679) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8680) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 4/2005 ! -! Copyright: MAGIC Software Development, 2000-2005 +! Copyright: MAGIC Software Development, 2000-2007 ! ! @@ -62,5 +62,4 @@ #include "../mhflux/MHAlpha.h" #include "../mhflux/MHCollectionArea.h" -//#include "../mhflux/MHThreshold.h" #include "../mhflux/MHEnergyEst.h" #include "../mhflux/MMcSpectrumWeight.h" @@ -80,5 +79,4 @@ #include "MFillH.h" #include "MHillasCalc.h" -//#include "MSrcPosCalc.h" #include "MContinue.h" @@ -89,12 +87,18 @@ MJSpectrum::MJSpectrum(const char *name, const char *title) : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0), - fEstimateEnergy(0), fCalcHadronness(0), fRefill(kFALSE), - /*fSimpleMode(kTRUE),*/ fForceTheta(kFALSE), - fRawMc(kFALSE), fNoThetaWeights(kFALSE) + fEstimateEnergy(0), fCalcHadronness(0), fForceTheta(kFALSE) { fName = name ? name : "MJSpectrum"; fTitle = title ? title : "Standard program to calculate spectrum"; -} - + + // Make sure that fDisplay is maintained properly + // (i.e. removed via RecursiveRemove if closed) + gROOT->GetListOfCleanups()->Add(this); +} + +// -------------------------------------------------------------------------- +// +// Delete all the fCut* data members and fCalcHadronness/fEstimateEnergy +// MJSpectrum::~MJSpectrum() { @@ -128,6 +132,14 @@ } +// -------------------------------------------------------------------------- +// +// Read a MTask named name into task from the open file. If this task is +// required mustexist can be set. Depending on success kTRUE or kFALSE is +// returned. If the task is a MContinue SetAllowEmpty is called. +// The name of the task is set to name. +// Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const { + // If a task is set delete task if (task) { @@ -136,4 +148,5 @@ } + // Try to get task from file task = (MTask*)gFile->Get(name); if (!task) @@ -144,4 +157,6 @@ return kFALSE; } + + // Check if task inherits from MTask if (!task->InheritsFrom(MTask::Class())) { @@ -151,6 +166,8 @@ } + // Set name of task task->SetName(name); + // SetAllowEmpty for MContinue tasks if (dynamic_cast(task)) dynamic_cast(task)->SetAllowEmpty(); @@ -159,4 +176,8 @@ } +// -------------------------------------------------------------------------- +// +// Print setup used for the MC processing, namely MAlphaFitter ans all fCut*. +// void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const { @@ -194,28 +215,6 @@ } - // Read the first corsika RunHeader from the MC files - TChain chain("RunHeaders"); - if (!set.AddFilesOn(chain)) - return kFALSE; - - MMcCorsikaRunHeader *h=0; - chain.SetBranchAddress("MMcCorsikaRunHeader.", &h); - chain.GetEntry(1); - - if (!h) - { - *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl; - return kFALSE; - } - - // Propagate the run header to MMcSpectrumWeight - if (!w.Set(*h)) - { - *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl; - return kFALSE; - } - - // Print new setup of MMcSpectrumWeight - w.Print(); + w.Print("new"); + return kTRUE; } @@ -281,4 +280,51 @@ } +// -------------------------------------------------------------------------- +// +// return maximum value of MMcRunHeader.fImpactMax stored in the RunHeaders +// of the files from the MC dataset +// +Bool_t MJSpectrum::AnalyzeMC(const MDataSet &set, Float_t &impactmax, Float_t &emin/*, Float_t emax*/) const +{ + if (fDisplay) + fDisplay->SetStatusLine1("Analyzing Monte Carlo headers..."); + + // ----- Create chain ----- + *fLog << inf << "Getting files... " << flush; + TChain chain("RunHeaders"); + if (!set.AddFilesOn(chain)) + return kFALSE; + *fLog << "done. " << endl; + + *fLog << all; + *fLog << "Searching for maximum impact... " << flush; + impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax"); + *fLog << "found " << impactmax/100 << "m" << endl; + + *fLog << "Searching for minimum lower energy limit... " << flush; + emin = chain.GetMinimum("MMcCorsikaRunHeader.fELowLim"); + *fLog << "found " << emin << "GeV" << endl; + + *fLog << "Searching for maximum lower energy limit... " << flush; + const Float_t emin2 = chain.GetMaximum("MMcCorsikaRunHeader.fELowLim"); + *fLog << "found " << emin2 << "GeV" << endl; + + if (emin2>emin) + { + *fLog << warn; + *fLog << "WARNING - You are using files with different lower limits for the simulated" << endl; + *fLog << " energy, i.e. that the collection area and your correction" << endl; + *fLog << " coefficients between " << emin << "GeV and "; + *fLog << emin2 << "GeV might be wrong." << endl; + } + +/* + *fLog << "Getting maximum energy... " << flush; + emax = chain.GetMaximum("MMcCorsikaRunHeader.fEUppLim"); + *fLog << "found " << emax << "GeV" << endl; + */ + return kTRUE; +} + Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const { @@ -286,33 +332,15 @@ fLog->Separator("Compiling original MC distribution"); - // The name of the input container is MMcEvtBasic - weight.SetNameMcEvt("MMcEvtBasic"); - // Get the corresponding rule for the weights - const TString w(weight.GetFormulaWeights()); - // Set back to MMcEvt - weight.SetNameMcEvt(); - - *fLog << inf << "Using weights: " << w << endl; - *fLog << "Please stand by, this may take a while..." << endl; - - if (fDisplay) - fDisplay->SetStatusLine1("Getting maximum impact..."); - - // ----- Create chain ----- - *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; + Float_t impactmax=0, Emin=0; + if (!AnalyzeMC(set, impactmax, Emin)) + return kFALSE; + + *fLog << inf << "Getting files... " << flush; MDirIter iter; if (!set.AddFilesOn(iter)) return kFALSE; *fLog << "done. " << endl; + + const Int_t tot = iter.GetNumEntries(); // Prepare histogram @@ -339,16 +367,30 @@ if (fDisplay) - fDisplay->SetStatusLine1("Reading OriginalMC"); + fDisplay->SetStatusLine1("Reading OriginalMC distribution..."); TH1 *hfill = (TH1*)h.Clone(h.InheritsFrom(TH2::Class())?"ThetaEMC":"ThetaMC"); hfill->SetDirectory(0); + *fLog << all << "Compiling simulated source spectrum... please stand by, this may take a while." << endl; + + Double_t evts = 0; + Int_t num = 0; + + // Reading this with a eventloop is five times slower :( TString fname; - while (1) - { + while (fDisplay) + { + if (fDisplay) + fDisplay->SetProgressBarPosition(Float_t(num++)/tot); + + // Get next filename fname = iter.Next(); if (fname.IsNull()) break; + if (fDisplay) + fDisplay->SetStatusLine2(fname); + + // open file TFile file(fname); if (file.IsZombie()) @@ -358,4 +400,18 @@ } + // Get tree OriginalMC + TTree *t = dynamic_cast(file.Get("OriginalMC")); + if (!t) + { + *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl; + return kFALSE; + } + if (t->GetEntries()==0) + { + *fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl; + continue; + } + + // Get tree RunHeaders TTree *rh = dynamic_cast(file.Get("RunHeaders")); if (!rh) @@ -364,5 +420,4 @@ return kFALSE; } - if (rh->GetEntries()!=1) { @@ -371,27 +426,35 @@ } - TTree *t = dynamic_cast(file.Get("OriginalMC")); - if (!t) + // Get corsika run header + MMcCorsikaRunHeader *head=0; + rh->SetBranchAddress("MMcCorsikaRunHeader.", &head); + rh->GetEntry(0); + if (!head) { - *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl; + *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from " << fname << "." << endl; return kFALSE; } + // Get the maximum impact parameter of this file. Due to different + // production areas an additional scale-factor is applied. + // + // Because it is assumed that the efficiency outside the production + // area is nearly zero no additional weight has to be applied to the + // events after cuts. + const Double_t impact = rh->GetMaximum("MMcRunHeader.fImpactMax"); + const Double_t scale = impactmax/impact; + + // Propagate the run header to MMcSpectrumWeight + if (!weight.Set(*head)) + { + *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl; + return kFALSE; + } + + // Get the corresponding rule for the weights + const TString weights(weight.GetFormulaWeights("MMcEvtBasic")); + + // No we found everything... go on reading contents *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 @@ -403,22 +466,33 @@ } hfill->SetDirectory(0); + + evts += hfill->GetEntries(); + + // ----- Complete incomplete energy ranges ----- + // ----- and apply production area weights ----- + weight.CompleteEnergySpectrum(*hfill, Emin); + + hfill->Scale(scale*scale); + + // Add weighted data from file to final histogram h.Add(hfill); } delete hfill; - *fLog << "Total number of events from file: " << h.GetEntries() << endl; + *fLog << all << "Total number of events from files in Monte Carlo dataset: " << evts << endl; if (fDisplay) fDisplay->SetStatusLine2("done."); - if (h.GetEntries()>0) + if (!fDisplay || h.GetEntries()>0) return kTRUE; - *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl; - + *fLog << err << "ERROR - Histogram with distribution from OriginalMC empty..." << endl; return kFALSE; } -Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const -{ +void MJSpectrum::GetThetaDistribution(TH1D &temp1, TH2D &hist2) const +{ + TH1D &temp2 = *hist2.ProjectionX(); + // Display some stuff if (fDisplay) @@ -430,4 +504,6 @@ c.cd(1); gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); temp1.DrawCopy(); @@ -435,4 +511,6 @@ c.cd(2); gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); temp2.SetName("NVsTheta"); temp2.DrawCopy("hist"); @@ -440,12 +518,16 @@ c.cd(4); gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); c.cd(3); gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); } // Calculate the Probability temp1.Divide(&temp2); - temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral()); + temp1.Scale(1./temp1.Integral()); // Some cosmetics: Name, Axis, etc. @@ -456,5 +538,5 @@ temp1.DrawCopy("hist"); - return kTRUE; + delete &temp2; } @@ -492,4 +574,5 @@ } + // Check whether histogram is empty if (proj.GetMaximum()==0) { @@ -497,12 +580,15 @@ *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; + *fLog << " Maybe the energy binning is undefined or wrong (from "; + *fLog << h2.GetYaxis()->GetXmin() << "GeV to " << h2.GetYaxis()->GetXmax() << "GeV)" << endl; theta->SetLineColor(kRed); return kFALSE;; } + // scale histogram and set new maximum for display proj.Scale(theta->GetMaximum()/proj.GetMaximum()); theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum())); + // draw project proj.Draw("same"); @@ -688,5 +774,5 @@ fill0.SetFilter(&sel1); - if (!fRawMc) + //if (!fRawMc) tlist1.AddToList(&sel1); tlist1.AddToList(&fill0); @@ -767,9 +853,13 @@ TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const { - TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4); - f.SetParameter(0, -2.5); - f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000))); + Axis_t lo, hi; + MH::GetRangeUser(spectrum, lo, hi); + + TF1 f("f", "[1]*(x/500)^[0]", lo, hi); + f.SetParameter(0, -3.0); + f.SetParameter(1, spectrum.GetMaximum()); f.SetLineColor(kBlue); - spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped + f.SetLineWidth(2); + spectrum.Fit(&f, "NIR"); // M skipped f.DrawCopy("same"); @@ -881,12 +971,13 @@ spec = (TH1D*)spectrum.DrawCopy(); - // Calculate Spectrum * E^2 - for (int i=0; i(fc).Copy(fa); + fa.DrawCopy("same"); // Display dN/dE*E^2 for conveinience @@ -896,10 +987,27 @@ gPad->SetGrid(); + // Calculate Spectrum * E^2 + for (int i=0; i(fc).Copy(fc2); + fc2.DrawCopy("same"); + + TF1 fa2("PG1553*E^2", "x^2*PG1553*1e-6", 100, 6000); + static_cast(fc).Copy(fa2); + fa2.DrawCopy("same"); // Fit Spectrum @@ -943,37 +1051,4 @@ return res; */ -/* - // Plot other spectra from Whipple - f.SetParameter(0, -2.45); - f.SetParameter(1, 3.3e-7); - f.SetRange(300, 8000); - f.SetLineColor(kBlack); - f.SetLineStyle(kDashed); - f.DrawCopy("same"); - - // Plot other spectra from Cangaroo - f.SetParameter(0, -2.53); - f.SetParameter(1, 2.0e-7); - f.SetRange(7000, 50000); - f.SetLineColor(kBlack); - f.SetLineStyle(kDashed); - f.DrawCopy("same"); - - // Plot other spectra from Robert - f.SetParameter(0, -2.59); - f.SetParameter(1, 2.58e-7); - f.SetRange(150, 1500); - f.SetLineColor(kBlack); - f.SetLineStyle(kDashed); - f.DrawCopy("same"); - - // Plot other spectra from HEGRA - f.SetParameter(0, -2.61); - f.SetParameter(1, 2.7e-7); - f.SetRange(1000, 20000); - f.SetLineColor(kBlack); - f.SetLineStyle(kDashed); - f.DrawCopy("same"); - */ } @@ -1179,4 +1254,5 @@ *fLog << endl; + // Setup everything which is read from the ganymed file MBinning bins1("BinningAlpha"); MBinning bins2("BinningEnergyEst"); @@ -1191,5 +1267,4 @@ MBinning binsa("BinningAsym"); MBinning binsb("BinningConc1"); - MAlphaFitter fit; @@ -1209,6 +1284,7 @@ plist.AddToList(&fit); - TH1D temp1, size; - const Float_t ontime = ReadInput(plist, temp1, size); + // Read from the ganymed file + TH1D htheta, size; + const Float_t ontime = ReadInput(plist, htheta, size); if (ontime<0) { @@ -1217,28 +1293,65 @@ } - plist.AddToList(&bins2); - - // Initialize weighting to a new spectrum - // as defined in the resource file + // Set Zenith angle binning to binning from the ganymed-file + bins3.SetEdges(htheta, 'x'); + + // Read energy binning from resource file + if (!CheckEnv(bins2)) + { + *fLog << err << "ERROR - Reading energy binning from resources failed." << endl; + return kFALSE; + } + plist.AddToList(&bins2); // For later use in MC processing + + // Initialize weighting to a new spectrum as defined in the resource file MMcSpectrumWeight weight; if (!InitWeighting(set, weight)) return kFALSE; - bins3.SetEdges(temp1, 'x'); - - // Read the MC distribution as produced by corsika into - // temp2 (1D) and apply the weights previously determined - TH1D temp2(temp1); - if (!ReadOrigMCDistribution(set, temp2, weight)) - return kFALSE; - - // Calculate the weights according to the zenith angle distribution - // of the raw-data which have to be applied to the MC events - if (!GetThetaDistribution(temp1, temp2)) - return kFALSE; - - // Tell the weighting task about the ZA-weights - if (!fNoThetaWeights) - weight.SetWeightsZd(&temp1); + // Print Theta and energy binning for cross-checks + *fLog << all << endl; + bins2.Print(); + bins3.Print(); + + // Now we read the MC distribution as produced by corsika + // vs zenith angle and energy. + // Weight for the new energy spectrum defined in MMcSpectumWeight + // are applied. + // Also correction for different lower energy bounds and + // different production areas (impact parameters) are applied. + TH2D hist; + hist.UseCurrentStyle(); + MH::SetBinning(&hist, &bins3, &bins2); + if (!ReadOrigMCDistribution(set, hist, weight)) + return kFALSE; + + // Check if user has closed the display + if (!fDisplay) + return kTRUE; + + // Display histograms and calculate za-weights into htheta + GetThetaDistribution(htheta, hist); + + // Give the zenoith angle weights to the weighting task + weight.SetWeightsZd(&htheta); + + // No we apply the the zenith-angle-weights to the corsika produced + // MC distribution. Unfortunately this must be done manually + // because we are multiplying column by column + for (int y=0; y