Index: trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc
===================================================================
--- trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc	(revision 2258)
+++ trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc	(revision 2261)
@@ -250,4 +250,7 @@
   //
   //
+  //
+  // Only for the binning taken from D. Kranich :
+  //
 
   for (Int_t thetabin = 1; thetabin <= fHistAll->GetNbinsY(); thetabin++)
@@ -258,10 +261,13 @@
       Float_t theta = fHistAll->GetYaxis()->GetBinCenter(thetabin);
 
-      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).
+      Float_t emin[4];       // Minimum energy in MC sample
+      Float_t emax[4];       // Maximum energy in MC sample
+      Float_t index[4];      // Spectral index
+      Float_t numevts[4];    // Number of events
+      Float_t multfactor[4]; // Factor by which the original number of events in an MC 
+                             // sample has been multiplied to account for the differences
+                             // in the generation areas of the various samples.
+      Float_t rmax;          // Maximum impact parameter range (on ground up to 45 degrees, 
+                             // on a plane perpendicular to Shower axis for 55 and 65 deg).
 
       memset(emin,    0, 4*sizeof(Float_t));
@@ -269,6 +275,10 @@
       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));
+      rmax = 0.;
+
+      multfactor[0] = 1.;
+      multfactor[1] = 1.;
+      multfactor[2] = 1.;
+      multfactor[3] = 1.;
 
       //
@@ -281,6 +291,4 @@
       if (theta > 8 && theta < 9)   // 8.75 deg
 	{
-	  rmin[0] = 0.;
-	  rmax[0] = 250.;     //meters
 	  emin[0] = 300.;
 	  emax[0] = 400.;  // Energies in GeV.
@@ -288,6 +296,4 @@
 	  numevts[0] = 4000.;
 
-	  rmin[1] = 0.;
-	  rmax[1] = 250.;     //meters
 	  emin[1] = 400.;
 	  emax[1] = 30000.;
@@ -295,10 +301,9 @@
 	  numevts[1] = 25740.;
 
+	  rmax = 250.;     //meters
 	  num_MC_samples = 2;
 	}
       else if (theta > 20 && theta < 21)  // 20.5 deg 
 	{
-	  rmin[0] = 0.;
-	  rmax[0] = 263.;     //meters
 	  emin[0] = 300.;
 	  emax[0] = 400.;  // Energies in GeV.
@@ -306,6 +311,4 @@
 	  numevts[0] = 6611.;
 
-	  rmin[1] = 0.;
-	  rmax[1] = 263.;     //meters
 	  emin[1] = 400.;
 	  emax[1] = 30000.;
@@ -313,10 +316,9 @@
 	  numevts[1] = 24448.;
 
+	  rmax = 263.;
 	  num_MC_samples = 2;
 	}
       else if (theta > 26 && theta < 27)  // 26.5 degrees
 	{
-	  rmin[0] = 0.;
-	  rmax[0] = 290.;     //meters
 	  emin[0] = 300.;
 	  emax[0] = 400.;  // Energies in GeV.
@@ -324,6 +326,4 @@
 	  numevts[0] = 4000.;
 
-	  rmin[1] = 0.;
-	  rmax[1] = 290.;     //meters
 	  emin[1] = 400.;
 	  emax[1] = 30000.;
@@ -331,10 +331,9 @@
 	  numevts[1] = 26316.;
 
+	  rmax = 290.;     //meters
 	  num_MC_samples = 2;
 	}
       else if (theta > 32 && theta < 33)  // 32.5 degrees
 	{
-	  rmin[0] = 0.;
-	  rmax[0] = 350.;     //meters
 	  emin[0] = 300.;
 	  emax[0] = 30000.;  // Energies in GeV.
@@ -342,10 +341,9 @@
 	  numevts[0] = 33646.;
 
+	  rmax = 350.;     //meters
 	  num_MC_samples = 1;
 	}
       else if (theta > 38 && theta < 39)  // 38.75 degrees
 	{
-	  rmin[0] = 0.;
-	  rmax[0] = 380.;     //meters
 	  emin[0] = 300.;
 	  emax[0] = 30000.;  // Energies in GeV.
@@ -353,10 +351,9 @@
 	  numevts[0] = 38415.;
 
+	  rmax = 380.;     //meters
 	  num_MC_samples = 1;
 	}
       else if (theta > 45 && theta < 47)  // 46 degrees
 	{
-	  rmin[0] = 0.;
-	  rmax[0] = 565.;     //meters
 	  emin[0] = 300.;
 	  emax[0] = 50000.;  // Energies in GeV.
@@ -364,10 +361,75 @@
 	  numevts[0] = 30197.;
 
+	  rmax = 565.;     //meters
 	  num_MC_samples = 1;
 	}
-
-      //
-      // TO BE FIXED: Add the cases for 55 and 65 degrees zenith angle.
-      //
+      else if (theta > 54 && theta < 56)  // 55 degrees
+	{
+	  //
+	  // The value of numevts in the first sample (below) has been
+	  // changed to simplify calculations. We have multiplied it
+	  // times 1.2808997 to convert it to the number it would be if 
+	  // the generation area was equal to that of the other samples 
+	  // at 55 degrees (pi*600**2 m2). This has to be taken into account 
+	  // in the error in the number of events.
+	  //
+
+	  emin[0] = 500.;
+	  emax[0] = 50000.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 3298.;
+	  multfactor[0] = 1.2808997;
+
+	  emin[1] = 1500.;
+	  emax[1] = 50000.;  // Energies in GeV.
+	  index[1] = 1.5;
+	  numevts[1] = 22229.;
+
+	  emin[2] = 1500.;
+	  emax[2] = 50000.;  // Energies in GeV.
+	  index[2] = 1.7;
+	  numevts[2] = 7553.;
+
+	  rmax = 600;     //meters
+	  num_MC_samples = 3;
+	}
+
+      else if (theta > 64 && theta < 66)  // 65 degrees
+	{
+	  emin[0] = 2000.;
+	  emax[0] = 50000.;  // Energies in GeV.
+	  index[0] = 1.5;
+	  numevts[0] = 16310.;
+
+	  emin[1] = 2000.;
+	  emax[1] = 50000.;  // Energies in GeV.
+	  index[1] = 1.7;
+	  numevts[1] = 3000.;
+
+	  //
+	  // The value of numevts in the next two samples (below) has been
+	  // changed to simplify calculations. We have converted them to the
+	  // number it would be if the generation area was equal to that of 
+	  // the first two samples at 65 degrees (pi*800**2 m2) (four times 
+	  // as many, since the original maximum impact parameter was 400
+	  // instead of 800. This is taken into account in the error too.
+	  //
+
+	  emin[2] = 5000.;
+	  emax[2] = 50000.;  // Energies in GeV.
+	  index[2] = 1.5;
+	  numevts[2] = 56584.;
+	  multfactor[2] = 4;
+
+	  emin[3] = 5000.;
+	  emax[3] = 50000.;  // Energies in GeV.
+	  index[3] = 1.7;
+	  numevts[3] = 11464;
+	  multfactor[3] = 4;
+
+	  rmax = 800;     // meters
+	  num_MC_samples = 4;
+	}
+
 
       for (Int_t i=1; i <= fHistAll->GetNbinsX(); i++)
@@ -377,4 +439,5 @@
 
 	  Float_t events = 0.;
+	  Float_t errevents = 0.;
 
 	  for (Int_t sample = 0; sample < num_MC_samples; sample++)
@@ -393,13 +456,16 @@
 
 	      events += k * (pow(e2, expo) - pow(e1, expo));
+	      errevents += multfactor[sample] * events;
 	    }
 
+	  errevents= sqrt(errevents);
+
 	  fHistAll->SetBinContent(i, thetabin, events);
-	  fHistAll->SetBinError(i, thetabin, sqrt(events));
+	  fHistAll->SetBinError(i, thetabin, errevents);
 	}
 
       // -----------------------------------------------------------
 
-      const Float_t dr = TMath::Pi() * (rmax[0]*rmax[0] - rmin[0]*rmin[0]);
+      const Float_t dr = TMath::Pi() * rmax * rmax;
 
       for (Int_t ix = 1; ix <= fHistAll->GetNbinsX(); ix++)
@@ -415,5 +481,9 @@
 	      // below 1E4, it means that no events triggered -> eff area = 0
 	      //
-
+	      // NOW DISABLED: because collection area after analysis does not
+	      // saturate at high E!
+	      //
+
+	      /*
 	      if (fHistSel->GetXaxis()->GetBinLowEdge(ix) > 4.)
 		{
@@ -421,4 +491,5 @@
 		  fHistCol->SetBinError(ix, thetabin, fHistCol->GetBinError(ix-1, thetabin));
 		}
+	      */
 	      continue;
 	    }
@@ -437,8 +508,20 @@
 
 	  //
-	  // Total area, perpendicular to the observation direction
+	  // Now we get the 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.);
+
+	  Float_t area;
+
+	  if (theta < 50)
+	    area = dr * cos(theta*TMath::Pi()/180.);
+	  else
+	    area = dr;
+
+	  // Above 50 degrees MC was generated with Corsika 6.xx, and the cores
+	  // were distributed on a circle perpendicular to the observation direction, 
+	  // and not on ground, hence the correction for cos(theta) is not necessary.
+	  //
+
 
 	  fHistCol->SetBinContent(ix, thetabin, eff*area);