Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnerThreCalc.cc
===================================================================
--- trunk/MagicSoft/Mars/mmontecarlo/MMcEnerThreCalc.cc	(revision 838)
+++ trunk/MagicSoft/Mars/mmontecarlo/MMcEnerThreCalc.cc	(revision 838)
@@ -0,0 +1,159 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appear in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Javier Lopez 05/2001 (jlopez@ifae.es)
+!
+!   Copyright: MAGIC Software Development, 2000-2001
+!
+!
+\* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+//  MMcEnerThreCalc
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+#include "MMcEnerThreCalc.h"
+
+#include "MParList.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MMcEvt.hxx"
+#include "MMcTrig.hxx"
+#include "MMcEnerThre.h"
+#include "MMcEnerHisto.h"
+
+#include <math.h>
+
+ClassImp(MMcEnerThreCalc)
+
+MMcEnerThreCalc::MMcEnerThreCalc (const int dim,
+				  const char* name, const char* title)
+{
+  *fName  = name  ? name  : "MMcEnerThreCalc";
+  *fTitle = title ? title : "Task to calculate the energy threshold from Monte Carlo";
+
+  fDimension = dim;
+
+  if (fDimension > 0)
+    fMMcTrig = new MMcTrig * [fDimension];
+  else
+    fMMcTrig = new MMcTrig * [1];
+}
+
+Bool_t MMcEnerThreCalc::PreProcess(MParList* pList)
+{
+  // connect Monte Carlo data with this task
+
+  char auxname[15]; // string to write container names
+
+  fMMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
+  if (!fMMcEvt)
+  {
+    *fLog << dbginf << "MMcEvt not found... aborting." << endl;
+    return kFALSE;
+  }
+
+  if (fDimension == 0)
+    {
+      fMMcTrig[0] = (MMcTrig*)pList->FindObject("MMcTrig");
+      if (!fMMcTrig[0])
+	{
+	  *fLog << dbginf << "MMcTrig not found... aborting." << endl;
+	  return kFALSE;
+	}
+    }
+  else
+    {
+      for (int i=0; i<fDimension; i++)
+	{
+	  sprintf(auxname,"MMcTrig;%i.",i+1);
+	  fMMcTrig[i] = (MMcTrig*)pList->FindObject(auxname);
+	  if (!fMMcTrig[i])
+	    {
+	      *fLog << dbginf << "MMcTrig;" << i+1 << "  not found... aborting." << endl;
+	      return kFALSE;
+	    } 
+	}
+    }
+
+  fMMcEnerThre = (MMcEnerThre*)pList->FindObject("MMcEnerThre");
+  if (!fMMcEnerThre)
+  {
+    *fLog << dbginf << "MMcEnerThre not found... aborting." << endl;
+    return kFALSE; 
+  }
+    
+  return kTRUE;
+}
+
+Bool_t MMcEnerThreCalc::Process()
+{
+  const Float_t energy = fMMcEvt->GetEnergy();
+  if (fDimension <= 0)
+    {
+      Int_t trigger[1]; 
+      trigger[0] = fMMcTrig[0]->GetFirstLevel();
+      if (trigger[0]>0) (*fMMcEnerThre)[0]->Fill(log10(energy),1/energy);
+    }
+  else
+    {
+      Int_t *trigger;
+      trigger=new Int_t [fDimension];
+
+      for (int i=0; i<fDimension; i++)
+	{
+	  trigger[i] = fMMcTrig[i]->GetFirstLevel() ;
+	  if (trigger[i]>0) (*fMMcEnerThre)[i]->Fill(log10(energy),1/energy);
+	}
+    }
+  return kTRUE;
+}
+
+Bool_t MMcEnerThreCalc::PostProcess()
+{
+  // fit the energy distribution to get the threshold
+
+  char auxname[15];
+  
+  if (fDimension <= 0) 
+    {
+      (*fMMcEnerThre)[0]->Fit("fLogEner0","RQ","",1.,3.);
+      (*fMMcEnerThre)[0]->Fit("fLogEner0","RQ","",(*fMMcEnerThre)[0]->GetPeakAtLogE()-2*(*fMMcEnerThre)[0]->GetSigmaAtLogE(),(*fMMcEnerThre)[0]->GetPeakAtLogE()+2*(*fMMcEnerThre)[0]->GetSigmaAtLogE());
+      (*fMMcEnerThre)[0]->Fit("fLogEner0","RQ","",(*fMMcEnerThre)[0]->GetPeakAtLogE()-sqrt(2)*(*fMMcEnerThre)[0]->GetSigmaAtLogE(),(*fMMcEnerThre)[0]->GetPeakAtLogE()+sqrt(2)*(*fMMcEnerThre)[0]->GetSigmaAtLogE());
+ 
+    }
+  else
+    {
+      for (int i=0; i<fDimension; i++)
+	{
+	  sprintf (auxname,"fLogEner%i",i);
+	  
+	  (*fMMcEnerThre)[i]->Fit(auxname,"RQ","",1.,3.);
+	  (*fMMcEnerThre)[i]->Fit(auxname,"RQ","",(*fMMcEnerThre)[i]->GetPeakAtLogE()-2*(*fMMcEnerThre)[i]->GetSigmaAtLogE(),(*fMMcEnerThre)[i]->GetPeakAtLogE()+2*(*fMMcEnerThre)[i]->GetSigmaAtLogE());
+	  (*fMMcEnerThre)[i]->Fit(auxname,"RQ","",(*fMMcEnerThre)[i]->GetPeakAtLogE()-sqrt(2)*(*fMMcEnerThre)[i]->GetSigmaAtLogE(),(*fMMcEnerThre)[i]->GetPeakAtLogE()+sqrt(2)*(*fMMcEnerThre)[i]->GetSigmaAtLogE());
+	}
+    }
+  return kTRUE;
+}
+
+
+
+