Index: /trunk/MagicSoft/Mars/Changelog
===================================================================
--- /trunk/MagicSoft/Mars/Changelog	(revision 2260)
+++ /trunk/MagicSoft/Mars/Changelog	(revision 2261)
@@ -4,7 +4,10 @@
  2003/07/03: Abelardo Moralejo
 
-  * mhistmc/MHMcCT1CollectionArea.cc
-    - Added code to allow the calculation of CT1 collection areas 
-      at 55 and 65 degrees (from the events in DK's MC library)
+   * mhistmc/MHMcCT1CollectionArea.cc
+     - Added code to allow the calculation of CT1 collection areas 
+       at 55 and 65 degrees (from the events in DK's MC library)
+
+   * macros/CT1collarea.C
+     - Changed binning in theta to include high ZAs
 
 
Index: /trunk/MagicSoft/Mars/macros/CT1collarea.C
===================================================================
--- /trunk/MagicSoft/Mars/macros/CT1collarea.C	(revision 2260)
+++ /trunk/MagicSoft/Mars/macros/CT1collarea.C	(revision 2261)
@@ -25,5 +25,5 @@
 
 
-void CT1collarea(TString filename="MC1.root", TString outname="")
+void CT1collarea(TString filename="MC_SC4.root", TString outname="")
 { 
     //
@@ -46,16 +46,20 @@
 
     MBinning binsx("BinningE");
+
     //    binsx.SetEdges(30,2.,5.);
+    //    Double_t xedge[10] =
+    //    {2.50515, 2.69897, 2.89763, 3.09691, 3.30103, 3.49969, 3.69984, 3.89982, 4.10003, 4.29994};
 
-    Double_t xedge[10] =
-    {2.50515, 2.69897, 2.89763, 3.09691, 3.30103, 3.49969, 3.69984, 3.89982, 4.10003, 4.29994};
+
+    Double_t xedge[15] = {2.47712, 2.64345, 2.82607, 3., 3.17609, 3.35218, 3.52892, 3.70415, 3.88024, 4.05652, 4.23274, 4.40875, 4.58478, 4.76080, 4.90309};
+
     const TArrayD xed;
-    xed.Set(10,xedge);
+    xed.Set(15,xedge);
     binsx.SetEdges(xed);
 
     MBinning binsy("BinningTheta");
-    const Double_t yedge[7] = {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50.};
+    const Double_t yedge[9] = {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.};
     const TArrayD yed;
-    yed.Set(7,yedge);
+    yed.Set(9,yedge);
     binsy.SetEdges(yed);
 
@@ -66,11 +70,10 @@
     tasklist.AddToList(&reader);   
 
-    /*
-    MF filterhadrons("HadNN.fHadronness<0.25");
+    MF filterhadrons("HadSC.fHadronness<0.3 && {abs(MHillasSrc.fAlpha)} < 13.1");
     tasklist.AddToList(&filterhadrons);
-    */
+
 
     MFillH filler("MHMcCT1CollectionArea","MMcEvt");
-    //    filler.SetFilter(&filterhadrons);
+    filler.SetFilter(&filterhadrons);
     tasklist.AddToList(&filler);
 
Index: /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc
===================================================================
--- /trunk/MagicSoft/Mars/mhistmc/MHMcCT1CollectionArea.cc	(revision 2260)
+++ /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);