Index: trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc
===================================================================
--- trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc	(revision 2250)
+++ trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc	(revision 2258)
@@ -254,100 +254,144 @@
     {
       // This theta is not exactly the one of the MC events, just about 
-      // the same:
+      // the same (bins have been selected so):
+
       Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin);
 
-      Float_t emin1, emax1, emin2, emax2;
-      Float_t index, expo, k1, k2;
-      Float_t numevts1, numevts2;
-      Float_t r1, r2;        // Impact parameter range (on ground).
-
-      emin1 = 0;  emax1 = 0; emin2 = 0;  emax2 = 0;
-      expo = 0.; k1 = 0.; k2 = 0.; r1 = 0.; r2 = 0.;
-      numevts1 = 0; numevts2 = 0;
+      Float_t emin[4];
+      Float_t emax[4];
+      Float_t index[4];
+      Float_t numevts[4];
+      Float_t rmin[4];
+      Float_t rmax[4];     // Impact parameter range (on ground).
+
+      memset(emin,    0, 4*sizeof(Float_t));
+      memset(emax,    0, 4*sizeof(Float_t));
+      memset(index,   0, 4*sizeof(Float_t));
+      memset(numevts, 0, 4*sizeof(Float_t));
+      memset(rmin,    0, 4*sizeof(Float_t));
+      memset(rmax,    0, 4*sizeof(Float_t));
+
+      //
+      // rmin and rmax are the minimum and maximum values of the impact 
+      // parameter of the shower on the ground (horizontal plane).
+      //
+
+      Int_t num_MC_samples = 0;
 
       if (theta > 8 && theta < 9)   // 8.75 deg
 	{
-	  r1 = 0.;
-	  r2 = 250.;     //meters
-	  emin1 = 300.;
-	  emax1 = 400.;  // Energies in GeV.
-	  emin2 = 400.;
-	  emax2 = 30000.;
-	  numevts1 = 4000.;
-	  numevts2 = 25740.;
+	  rmin[0] = 0.;
+	  rmax[0] = 250.;     //meters
+	  emin[0] = 300.;
+	  emax[0] = 400.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 4000.;
+
+	  rmin[1] = 0.;
+	  rmax[1] = 250.;     //meters
+	  emin[1] = 400.;
+	  emax[1] = 30000.;
+	  index[1] = 1.5;
+	  numevts[1] = 25740.;
+
+	  num_MC_samples = 2;
 	}
       else if (theta > 20 && theta < 21)  // 20.5 deg 
 	{
-	  r1 = 0.;
-	  r2 = 263.;     //meters
-	  emin1 = 300.;
-	  emax1 = 400.;  // Energies in GeV.
-	  emin2 = 400.;
-	  emax2 = 30000.;
-	  numevts1 = 6611.;
-	  numevts2 = 24448.;
+	  rmin[0] = 0.;
+	  rmax[0] = 263.;     //meters
+	  emin[0] = 300.;
+	  emax[0] = 400.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 6611.;
+
+	  rmin[1] = 0.;
+	  rmax[1] = 263.;     //meters
+	  emin[1] = 400.;
+	  emax[1] = 30000.;
+	  index[1] = 1.5;
+	  numevts[1] = 24448.;
+
+	  num_MC_samples = 2;
 	}
       else if (theta > 26 && theta < 27)  // 26.5 degrees
 	{
-	  r1 = 0.;
-	  r2 = 290.;     //meters
-	  emin1 = 300.;
-	  emax1 = 400.;  // Energies in GeV.
-	  emax2 = emax1;	  emin2 = 400.;
-	  emax2 = 30000.;
-	  numevts1 = 4000.;
-	  numevts2 = 26316.;
+	  rmin[0] = 0.;
+	  rmax[0] = 290.;     //meters
+	  emin[0] = 300.;
+	  emax[0] = 400.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 4000.;
+
+	  rmin[1] = 0.;
+	  rmax[1] = 290.;     //meters
+	  emin[1] = 400.;
+	  emax[1] = 30000.;
+	  index[1] = 1.5;
+	  numevts[1] = 26316.;
+
+	  num_MC_samples = 2;
 	}
       else if (theta > 32 && theta < 33)  // 32.5 degrees
 	{
-	  r1 = 0.;
-	  r2 = 350.;     //meters
-	  emin1 = 300.;
-	  emax1 = 30000.;  // Energies in GeV.
-	  emax2 = emax1;
-	  numevts1 = 33646.;
+	  rmin[0] = 0.;
+	  rmax[0] = 350.;     //meters
+	  emin[0] = 300.;
+	  emax[0] = 30000.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 33646.;
+
+	  num_MC_samples = 1;
 	}
       else if (theta > 38 && theta < 39)  // 38.75 degrees
 	{
-	  r1 = 0.;
-	  r2 = 380.;     //meters
-	  emin1 = 300.;
-	  emax1 = 30000.;  // Energies in GeV.
-	  emax2 = emax1;
-	  numevts1 = 38415.;
+	  rmin[0] = 0.;
+	  rmax[0] = 380.;     //meters
+	  emin[0] = 300.;
+	  emax[0] = 30000.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 38415.;
+
+	  num_MC_samples = 1;
 	}
       else if (theta > 45 && theta < 47)  // 46 degrees
 	{
-	  r1 = 0.;
-	  r2 = 565.;     //meters
-	  emin1 = 300.;
-	  emax1 = 50000.;  // Energies in GeV.
-	  emax2 = emax1;
-	  numevts1 = 30197.;
-	}
-
-      index = 1.5;     // Differential spectral Index.
-      expo = 1.-index;
-      k1 = numevts1 / (pow(emax1,expo) - pow(emin1,expo));
-      k2 = numevts2 / (pow(emax2,expo) - pow(emin2,expo));
+	  rmin[0] = 0.;
+	  rmax[0] = 565.;     //meters
+	  emin[0] = 300.;
+	  emax[0] = 50000.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 30197.;
+
+	  num_MC_samples = 1;
+	}
+
+      //
+      // TO BE FIXED: Add the cases for 55 and 65 degrees zenith angle.
+      //
 
       for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
 	{
-	  const Float_t e1 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i));
-	  const Float_t e2 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i+1));
-
-	  if (e1 < emin1 || e2 > emax2)
-            continue;
-
-	  Float_t events;
-
-	  if (e2 <= emax1)
-	    events = k1 * (pow(e2, expo) - pow(e1, expo));
-	  else if (e1 >= emin2)
-	    events = k2 * (pow(e2, expo) - pow(e1, expo));
-	  else
-	    events = 
-	      k1 * (pow(emax1, expo) - pow(e1, expo))+
-	      k2 * (pow(e2, expo) - pow(emin2, expo));
+	  Float_t e1 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i));
+	  Float_t e2 = pow(10.,fHistAll->GetXaxis()->GetBinLowEdge(i+1));
+
+	  Float_t events = 0.;
+
+	  for (Int_t sample = 0; sample < num_MC_samples; sample++)
+	    {
+	      Float_t expo = 1.-index[sample];
+	      Float_t k = numevts[sample] / (pow(emax[sample],expo) - pow(emin[sample],expo));
+
+	      if (e2 < emin[sample] || e1 > emax[sample])
+		continue;
+
+	      if (emin[sample] > e1) 
+		e1 = emin[sample];
+
+	      if (emax[sample] < e2) 
+		e2 = emax[sample];
+
+	      events += k * (pow(e2, expo) - pow(e1, expo));
+	    }
 
 	  fHistAll->SetBinContent(i, thetabin, events);
@@ -357,5 +401,5 @@
       // -----------------------------------------------------------
 
-      const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
+      const Float_t dr = TMath::Pi() * (rmax[0]*rmax[0] - rmin[0]*rmin[0]);
 
       for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
@@ -392,5 +436,8 @@
 	  const Double_t efferr = sqrt((1.-eff)*Ns)/Na;
 
-
+	  //
+	  // Total area, perpendicular to the observation direction
+	  // in which the events were generated (correct for cos theta):
+	  //
 	  const Float_t area = dr * cos(theta*TMath::Pi()/180.);