Index: /trunk/MagicSoft/Mars/macros/getThreshold.C =================================================================== --- /trunk/MagicSoft/Mars/macros/getThreshold.C (revision 860) +++ /trunk/MagicSoft/Mars/macros/getThreshold.C (revision 861) @@ -1,7 +1,7 @@ -void getThreshold( char* filename = "data/oscar_protons.root") +void getThreshold( char* filename = "/hdb1/users/jlopez/temp_loop.root") { // - // This macro fill the container MMcEnerThre using the task - // MMcEnerThreCalc and shows the results. + // This macro fill the container MHMcEnergies using the task + // MMcThresholdCalc and shows the results. // MParList parlist; @@ -12,5 +12,5 @@ // // Setup the parameter list - // - You need create the container MMcEnerThre. + // - You need create the container MHMcEnergies. // + You need to put the number of trigger conditions when // you declarete the MHMcEnergies @@ -18,5 +18,5 @@ // taking only the trigger information from MMcTrig // - const UInt_t numtriggerconditions = 1; + const UInt_t numtriggerconditions = 5; MHMcEnergies *hists = new MHMcEnergies(numtriggerconditions); @@ -49,5 +49,7 @@ // Now you can display the results // - for (UInt_t i=0; i0) sprintf(aux+10, ";%i", idx); - hLogEner = new TH1F(aux, "", 100, 0.5, 4.5); + hLogEner = new TH1F(aux, "", 40, 0.5, 4.5); hLogEner->SetXTitle("log(E) [GeV]"); hLogEner->SetYTitle("dN/dE"); - //hLogEner->SetBins(60); } +//------------------------------------------------------------------------- +// +// Defualt Destructor +// MHMcEnergy::~MHMcEnergy() { @@ -70,4 +82,8 @@ } +//-------------------------------------------------------------------------- +// +// Fill the histogram with the log10 of the energy for triggered events. +// void MHMcEnergy::Fill(Float_t log10E, Float_t w) { @@ -75,17 +91,21 @@ } +// ------------------------------------------------------------------------- +// +// Fitting function +// void MHMcEnergy::Fit(Axis_t xxmin, Axis_t xxmax) { // // 0: don't draw the function (it is drawn together with the histogram) - // +: add these function to the list of fits. Don't delete the last fit. + // Q: quiet mode // - // FIXME: R means: use the range specified in the function (xxmin, xxmax are ignored!) - // Q means: quiet (why?) - // - // - hLogEner->Fit(fLogEner->GetName(), "Q0+", "", xxmin, xxmax); + hLogEner->Fit(fLogEner->GetName(), "Q0", "", xxmin, xxmax); } +// ------------------------------------------------------------------------ +// +// Drawing function. It creates its own canvas. +// void MHMcEnergy::Draw(Option_t *option) { @@ -111,4 +131,16 @@ } +// -------------------------------------------------------------------------- +// +// Set the number of bins in the histogran. +// +void MHMcEnergy::SetBins(Int_t nbins) +{ + hLogEner->SetBins(nbins,0.5,4.5); +} +// -------------------------------------------------------------------------- +// +// Write the threshold and its error in the standard output +// void MHMcEnergy::Print(Option_t*) { @@ -116,4 +148,8 @@ } +// ------------------------------------------------------------------------- +// +// Return the threshold +// Float_t MHMcEnergy::GetThreshold() const { @@ -123,4 +159,8 @@ } +// ------------------------------------------------------------------------- +// +// Return the error of the threshold. +// Float_t MHMcEnergy::GetThresholdErr() const { @@ -129,7 +169,12 @@ const Float_t p1err = fLogEner->GetParError(1); + // The error has into accuont the error in the fit return pow(10, p1) * p1err * lg10; } +// ------------------------------------------------------------------------- +// +// Return the peak of the fitted gaussan function. +// Float_t MHMcEnergy::GetGaussPeak() const { @@ -137,4 +182,8 @@ } +// ------------------------------------------------------------------------- +// +// Return the sigma of the fitted gaussan function. +// Float_t MHMcEnergy::GetGaussSigma() const { @@ -142,2 +191,5 @@ } + + + Index: /trunk/MagicSoft/Mars/mhist/MHMcEnergy.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMcEnergy.h (revision 860) +++ /trunk/MagicSoft/Mars/mhist/MHMcEnergy.h (revision 861) @@ -32,9 +32,10 @@ void Fill(Float_t log10E, Float_t w); void Fit(Axis_t xxmin, Axis_t xxmax); + void SetBins(Int_t nbins = 100); void Draw(Option_t* option = ""); void Print(Option_t* option = NULL); - ClassDef(MHMcEnergy, 1) // Histogram container for montecarlo energy + ClassDef(MHMcEnergy, 1) // Histogram container for montecarlo energy threshold }; Index: /trunk/MagicSoft/Mars/mmontecarlo/MMcThresholdCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mmontecarlo/MMcThresholdCalc.cc (revision 860) +++ /trunk/MagicSoft/Mars/mmontecarlo/MMcThresholdCalc.cc (revision 861) @@ -24,14 +24,14 @@ \* ======================================================================== */ -///////////////////////////////////////////////////////////////////////////// -// // -// MMcThresholdCalc // -// // -// Input Containers: // -// MMcEvt // -// // -// Output Containers: // -// MMcTrig;*, MHMcEnergy;* // -// // +/////////////////////////////////////////////////////////////////////////// +// +// MMcThresholdCalc +// +// Input Containers: +// MMcEvt, MMcTrig;* +// +// Output Containers: +// MHMcEnergies +// ///////////////////////////////////////////////////////////////////////////// @@ -52,4 +52,8 @@ ClassImp(MMcThresholdCalc) +// -------------------------------------------------------------------------- +// +// Default Constructor. +// MMcThresholdCalc::MMcThresholdCalc(const UInt_t dim, const char* name, const char* title) : fDimension(dim) @@ -58,4 +62,6 @@ *fTitle = title ? title : "Task to calculate the energy threshold from Monte Carlo"; + // Arrays of MMcTrig* and MHMcEnergy* are created in order to be + // able to work with root files with several trigger conditions. fMcTrig = new MMcTrig*[fDimension]; fHMcEnergy = new MHMcEnergy*[fDimension]; @@ -66,4 +72,8 @@ } +// ------------------------------------------------------------------------- +// +// Default Destructor. +// MMcThresholdCalc::~MMcThresholdCalc() { @@ -81,4 +91,7 @@ // connect Monte Carlo data with this task + // This task has into accout if the root file has one trigger + // condition (MMcTrig) or severl of them (MMcTrig;#.) + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); if (!fMcEvt) @@ -93,5 +106,5 @@ { if (fDimension>1) - sprintf(auxname+7, ";%i", i+1); + sprintf(auxname+7, ";%i.", i+1); fMcTrig[i] = (MMcTrig*)pList->FindObject(auxname); @@ -110,6 +123,6 @@ for (unsigned int i=0; i1) - sprintf(auxname+10, ";%i", i+1); + if (fDimension>1&&i!=0) + sprintf(auxname+10, ";%i", i); fHMcEnergy[i] = (MHMcEnergy*)pList->FindObject(auxname); @@ -119,5 +132,5 @@ *fLog << dbginf << "'" << auxname << "' not found in list... creating." << endl; - fHMcEnergy[i] = new MHMcEnergy(fDimension>1 ? i+1 : 0); + fHMcEnergy[i] = new MHMcEnergy(fDimension>1&&i!=0 ? i : 0); fMustDelete[i] = kTRUE; pList->AddToList(fHMcEnergy[i]); @@ -129,4 +142,9 @@ Bool_t MMcThresholdCalc::Process() { + + // The histograms are filled with log10 of the energy for triggered + // events and weighted with 1/E because it is needed the dN/dE vs. logE + // distribution to get the energy threshold. + const Float_t energy = fMcEvt->GetEnergy(); const Float_t lg10 = log10(energy); @@ -147,4 +165,5 @@ { // fit the energy distribution to get the threshold + // Some iterations are done to be sure the fit parameters converge. const Float_t sqrt2 = sqrt(2); @@ -152,12 +171,16 @@ for (unsigned int i=0; iFit(1, 3); + + peak = fHMcEnergy[i]->GetGaussPeak(); + sigma = fHMcEnergy[i]->GetGaussSigma(); + fHMcEnergy[i]->Fit(peak - 2. *sigma, peak + 2. *sigma); + + peak = fHMcEnergy[i]->GetGaussPeak(); + sigma = fHMcEnergy[i]->GetGaussSigma(); + fHMcEnergy[i]->Fit(peak - sqrt2*sigma, peak + sqrt2*sigma); } return kTRUE; Index: /trunk/MagicSoft/Mars/mmontecarlo/MMcThresholdCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mmontecarlo/MMcThresholdCalc.h (revision 860) +++ /trunk/MagicSoft/Mars/mmontecarlo/MMcThresholdCalc.h (revision 861) @@ -39,5 +39,5 @@ Bool_t PostProcess(); - ClassDef(MMcThresholdCalc, 0) // Task to fill the energy threshold information + ClassDef(MMcThresholdCalc, 0) // Task to compute the energy threshold };