Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 1001)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 1078)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.28 $
+// $Revision: 1.29 $
 // $Author: blanch $ 
-// $Date: 2001-10-26 16:31:45 $
+// $Date: 2001-11-14 17:38:23 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -69,4 +69,5 @@
 #include "MMcEvt.hxx"
 #include "MMcTrig.hxx"
+#include "MMcRunHeader.hxx"
 #include "MMcTrigHeader.hxx"
 #include "MMcFadcHeader.hxx"
@@ -431,5 +432,5 @@
 
   int simulateNSB;            //@< Will we simulate NSB?
-  int nphe2NSB;               //@< From how many ohe will we simulate NSB?
+  int nphe2NSB;               //@< From how many phe will we simulate NSB?
   float meanNSB;              //@< diffuse NSB mean value (phe per ns per central pixel)
   float diffnsb_phepns[iMAXNUMPIX];  //@< diffuse NSB values for each pixel  
@@ -455,4 +456,5 @@
     EXTWAVEBAND5
   };                          
+  float zenfactor=1.0;  //  correction factor calculated from the extinction
 
   float qTailCut;             //@< Tail Cut value
@@ -473,5 +475,5 @@
   float fpixelthres[TRIGGER_PIXELS];         //@< Threshold values
 
-  TArrayC *fadcValues;  //@< the analog Fadc siganl for pixels
+  TArrayC *fadcValues;  //@< the analog Fadc signal for pixels
 
   float plateScale_cm2deg;    //@< plate scale (deg/cm)
@@ -479,4 +481,22 @@
 
   int still_in_loop = FALSE;
+
+  //@< Variables to fill the McRunHeader
+  Int_t sfRaH = 5;
+  Int_t sfRaM = 34;
+  Int_t sfRaS = 32;
+  Int_t sfDeD = 22;
+  Int_t sfDeM = 00;
+  Int_t sfDeS = 52;
+  Float_t shthetamax = 0.0;
+  Float_t shthetamin = 0.0;
+  Float_t shphimax = 0.0;
+  Float_t shphimin = 0.0;
+  Float_t telestheta = 0.0;
+  Float_t telesphi = 0.0;
+  Float_t sofftheta = 0.0;
+  Float_t soffphi = 0.0;
+  UInt_t corsika = 5200 ;
+
 
   struct camera cam; // structure holding the camera definition
@@ -743,7 +763,7 @@
   //  Initialise McTrig information class if we want to save trigger informtion
 
-  MMcTrig **McTrig; 
-  MMcTrigHeader **HeaderTrig; 
-  MMcFadcHeader **HeaderFadc;
+  MMcTrig **McTrig = NULL; 
+  MMcTrigHeader **HeaderTrig = NULL; 
+  MMcFadcHeader **HeaderFadc = NULL;
 
   if (Write_McTrig){
@@ -796,12 +816,10 @@
        // Header Tree
 
-  MRawRunHeader **RunHeader;
-
-  RunHeader = new MRawRunHeader * [1]; 
-  RunHeader[0] = new MRawRunHeader();
+  MRawRunHeader *RunHeader= new MRawRunHeader();
+  MMcRunHeader  *McRunHeader = new MMcRunHeader();
 
        // Header branch
 
-  MRawEvtHeader **EvtHeader;
+  MRawEvtHeader **EvtHeader = NULL;
 
   if (Write_RawEvt) {
@@ -815,5 +833,5 @@
        // Data branch
 
-  MRawEvtData **EvtData;    // Data branch
+  MRawEvtData **EvtData = NULL;    // Data branch
 
   if (Write_RawEvt) {
@@ -822,4 +840,6 @@
     for (i=0;i<icontrigger;i++) {
       EvtData[i] = new MRawEvtData();
+      EvtData[i]->Init(RunHeader);     //  We need thr RunHeader to read
+                                       //  number of pixels
     }
   }
@@ -833,6 +853,4 @@
   TFile outfile_temp ( rootname , "RECREATE" ); 
 
-  Int_t bsize=128000; Int_t split=1;
-
   //      create a Tree for the Header Event
   TTree HeaderTree("RunHeaders","Headers of Run");
@@ -843,10 +861,13 @@
 
   HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
-  		    &RunHeader[0], bsize, split);
+  		    &RunHeader);
+
+  HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
+  		    &McRunHeader);
 
   if(!Trigger_Loop && Write_McTrig){
     
     HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader", 
-		 &HeaderTrig[0], bsize, split);    
+		      &HeaderTrig[0]);
   }
   if (Trigger_Loop && Write_McTrig){
@@ -858,5 +879,5 @@
       strcat (branchname, ".");
       HeaderTree.Branch(branchname,"MMcTrigHeader", 
-		     &HeaderTrig[i], bsize, split);
+			&HeaderTrig[i]);
     }
   }  
@@ -865,5 +886,5 @@
     
     HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader", 
-		 &HeaderFadc[0], bsize, split);    
+		      &HeaderFadc[0]);
   }
   if (Trigger_Loop && Write_McFADC){
@@ -875,5 +896,5 @@
       strcat (branchname, ".");
       HeaderTree.Branch(branchname,"MMcFadcHeader", 
-		     &HeaderFadc[i], bsize, split);
+			&HeaderFadc[i]);
     }
   }  
@@ -881,10 +902,12 @@
   //  Fill branches for MRawRunHeader
 
-  RunHeader[0]->SetMagicNumber(kMagicNumber);
-  RunHeader[0]->SetFormatVersion(2);
-  RunHeader[0]->SetSoftVersion((UShort_t) (VERSION*10));
-  RunHeader[0]->SetRunType(256);
-  RunHeader[0]->SetRunNumber(0);
-  RunHeader[0]->SetNumSamples(0,FADC_SLICES);
+  RunHeader->SetMagicNumber(kMagicNumber);
+  RunHeader->SetFormatVersion(2);
+  RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
+  RunHeader->SetRunType(256);
+  RunHeader->SetRunNumber(0);
+  RunHeader->SetNumSamples(0,FADC_SLICES);
+  RunHeader->SetNumCrates(1);
+  RunHeader->SetNumPixInCrate(CAMERA_PIXELS);
 
 
@@ -956,5 +979,5 @@
 
   //  Fill the Header Tree with the current leaves of each branch
-  HeaderTree.Fill() ;
+  // HeaderTree.Fill() ;
 	    
 
@@ -965,5 +988,5 @@
 
     EvtTree.Branch("MMcEvt","MMcEvt", 
-		   &McEvt, bsize, split);  
+		   &McEvt);
   }
 
@@ -972,11 +995,11 @@
     if (Write_RawEvt){
       EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", 
-      	     &EvtHeader[0], bsize, split);
+		     &EvtHeader[0]);
       EvtTree.Branch("MRawEvtData","MRawEvtData", 
-		     &EvtData[0], bsize, split);
+		     &EvtData[0]);
     }
     if (Write_McTrig){
       EvtTree.Branch("MMcTrig","MMcTrig", 
-		     &McTrig[0], bsize, split);
+		     &McTrig[0]);
     }    
   }
@@ -990,5 +1013,5 @@
 	strcat (branchname, ".");
 	EvtTree.Branch(branchname,"MMcTrig", 
-		       &McTrig[i], bsize, split);
+		       &McTrig[i]);
       }
     }
@@ -1003,5 +1026,5 @@
       strcat (branchname, ".");
       EvtTree.Branch(branchname,"MRawEvtHeader", 
-		     &EvtHeader[i], bsize, split);
+		     &EvtHeader[i]);
     }
     for(char branchname[15],i=0;i<icontrigger;i++){
@@ -1012,24 +1035,23 @@
       strcat (branchname, ".");
       EvtTree.Branch(branchname,"MRawEvtData", 
-		     &EvtData[i], bsize, split);
-    }
-  }
-
-  unsigned short ulli = 0 ; 
+		     &EvtData[i]);
+    }
+  }
 
   TApplication theAppTrigger("App", &argc, argv);
 
-  if (gROOT->IsBatch()) {
-    fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
-    //    return 1;
-  }
-
   if(FADC_Scan){
-  TApplication theAppFadc("App", &argc, argv);
-
-  if (gROOT->IsBatch()) {
-    fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
-    //    return 1;
-  }
+      if (gROOT->IsBatch()) {
+	  fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
+	  //    return 1;
+      }
+  }
+  if(FADC_Scan){
+      //TApplication theAppFadc("App", &argc, argv);
+
+      if (gROOT->IsBatch()) {
+	  fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
+	  //    return 1;
+      }
   }
   
@@ -1060,5 +1082,5 @@
     k = produce_nsbrates( starfieldname,
 			  &cam,
-			  nsbrate_phepns );
+			  nsbrate_phepns);
 
     if (k != 0){
@@ -1075,8 +1097,13 @@
 
     // calculate nsb rate including diffuse and starlight
-    for(i=0; i<cam.inumpixels;i++){
-      nsb_phepns[i]= diffnsb_phepns[i];
-      for(j=0;j<iNUMWAVEBANDS;j++)
-	nsb_phepns[i]=nsb_phepns[i]+ nsbrate_phepns[i][j];
+    // we also include the extinction effect
+    for(j=0;j<iNUMWAVEBANDS;j++){
+	// calculate the effect of the atmospheric extinction 
+    
+	zenfactor = pow(10., -0.4 * ext[i] );
+	
+	for(i=0; i<cam.inumpixels;i++){
+	    nsb_phepns[i]+=diffnsb_phepns[i]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[i][j];
+	}
     }
   }
@@ -1357,4 +1384,10 @@
 	fadc.ElecNoise() ;
 
+	//  We should simulate the AC coupling behaviour:
+        //  For the FADC it is only done for the NSB while producinr
+        //  the StarResponse database.
+	//  For the trigger is done in the Trigger.Diskriminate(), which
+	//  is called later (it should be seperated to speed up the program
+
 	//  now a shift in the fadc signal due to the pedestlas is 
 	//  intorduced
@@ -1406,5 +1439,6 @@
 		  //   Start the First Level Trigger simulation
 		  //
-		  Lev1=Trigger.FirstLevel();
+		  if(Lev0!=0)
+		      Lev1=Trigger.FirstLevel();
 		  if (Write_McTrig)
 		    McTrig[iconcount]->SetFirstLevel (Lev1);
@@ -1700,4 +1734,44 @@
     } // end if else found start of run
   } // end big while loop
+
+  //<@ Finally we should fill th McRunHeader
+
+
+  get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
+  get_teles_axis(&telestheta,&telesphi);
+  get_source_off(&sofftheta,&soffphi);
+  mcevth.get_theta_range(&shthetamin, &shthetamax);
+  mcevth.get_phi_range(&shphimin,&shphimax);
+  corsika=get_corsika_ver();
+  if(!Trigger_Loop) icontrigger=0;
+  McRunHeader->Fill(icontrigger,
+		   !Write_All_Images,
+		   Write_McEvt,
+		   Write_McTrig,
+		   Write_McFADC,
+		   CAMERA_PIXELS,
+		   (UInt_t)ntshow,
+		   (UInt_t)ntrigger,
+		   sfRaH,
+		   sfRaM,
+		   sfRaS,
+		   sfDeD,
+		   sfDeM,
+		   sfDeS,
+		   meanNSB,
+		   telestheta,
+		   telesphi,
+		   sofftheta,
+		   soffphi,
+		   shthetamax,
+		   shthetamin,
+		   shphimax,
+		   shphimin,
+		   corsika,
+		   (UInt_t)(REFL_VERSION*100),
+		   (UInt_t)(VERSION*100));
+
+  //  Fill the Header Tree with the current leaves of each branch
+  HeaderTree.Fill() ;
   
   //++
@@ -2621,11 +2695,11 @@
 
 /******** bpoint_is_in_pix() ***************************************/
-
+/*
 int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
-  /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
-  /* the pixel is assumed to be a "closed set" */
-
-  double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
-  double c, xx, yy; /* auxiliary variable */
+  // return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system)
+  // the pixel is assumed to be a "closed set"
+
+  double a, b; // a = length of one of the edges of one pixel, b = half the width of one pixel
+  double c, xx, yy; // auxiliary variable
 
   b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
@@ -2643,9 +2717,50 @@
   if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) ||
      ((0. <  xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a)))   ){
-    return(TRUE); /* point is inside */
+    return(TRUE); // point is inside
   }
   else{
-    return(FALSE); /* point is outside */
-  }
+    return(FALSE); // point is outside
+  }
+}
+*/
+
+#define sqrt13 0.577350269 // = 1./sqrt(3.)
+#define sqrt3  1.732050807 // = sqrt(3.)
+
+int bpoint_is_in_pix(double dx, double dy, struct camera *pcam)
+{
+    /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
+    /* the pixel is assumed to be a "closed set" */
+
+    /*
+     a = length of one of the edges of one pixel,
+     b = half the width of one pixel
+     */
+
+    const int    num  = pcam->inumpixels;
+    const double diam = pcam->dpixdiameter_cm;
+
+    const double diam2 = diam/2;
+    const double diam3 = diam/sqrt3;
+
+    for (int i=0; i<num; i++)
+    {
+        const double size = pcam->dpixsizefactor[i];
+
+        const double b = diam2 * size;
+        const double a = diam3 * size;
+
+        const double xx = dx - pcam->dxc[i];
+        const double yy = dy - pcam->dyc[i];
+
+        if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
+           ((0. <  xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a)))   ){
+
+            return i; // inside i
+        }
+
+	//   return -1; // outside
+    }
+    return -1; // outside
 }
 
@@ -2789,9 +2904,13 @@
   // read the photons data
   
-  fread ( flag, SIZE_OF_FLAGS, 1, sp );
   
   // loop over the photons
   
-  while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
+  while (1) {
+
+    fread ( flag, SIZE_OF_FLAGS, 1, sp );
+
+    if (isA( flag, FLAG_END_OF_EVENT ))
+        break;
 	  
     memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
@@ -2808,11 +2927,5 @@
 
     if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
- 
-     //  read next Photon
-	  
-      fread( flag, SIZE_OF_FLAGS, 1, sp );
-      
-      // go to beginning of loop, the photon is lost
-      continue;
+        continue;
     }
   	  
@@ -2833,45 +2946,16 @@
 	  
     if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){ 
-	    
-      // cout << " lost first check\n";
-
-      // read next CPhoton
-
-      fread ( flag, SIZE_OF_FLAGS, 1, sp );
-	    
-      // go to beginning of loop, the photon is lost
       continue;
 	    
     }
 
-    ipixnum = -1;
-
-    for(i=0; i<cam->inumpixels; i++){
-      if( bpoint_is_in_pix( cx, cy, i, cam) ){
-	ipixnum = i;
-	break;
-      }
-    }
-
-    if(ipixnum==-1){// the photon is in none of the pixels
-
-      // cout << " lost pixlization\n";
-
-      // read next CPhoton
-
-      fread ( flag, SIZE_OF_FLAGS, 1, sp );
-      
-      // go to beginning of loop, the photon is lost
-      continue;
-    }
-
-    if(ipixnum==0) {// the phton is in the central pixel, which is not used for trigger
-      // read next CPhoton
-
-      fread ( flag, SIZE_OF_FLAGS, 1, sp );
-      
-      // go to beginning of loop, the photon is lost
-      continue;
-    }
+    ipixnum = bpoint_is_in_pix(cx, cy, cam);
+
+    // -1 = the photon is in none of the pixels
+    //  0 = the phton is in the central pixel, which is not used for trigger
+    if (ipixnum==-1 || ipixnum==0) {
+        continue;
+    }
+
     //+++
     // QE simulation
@@ -2885,12 +2969,4 @@
     
     if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
-
-      // cout << " lost\n"; 
-      
-      // read next Photon
-
-      fread ( flag, SIZE_OF_FLAGS, 1, sp );
-	    
-      // go to beginning of loop
       continue;
 	    
@@ -2911,14 +2987,5 @@
 
     if ( (RandomNumber) > qe ) {
-
-      // cout << " lost\n"; 
-
-      // read next Photon
-
-      fread ( flag, SIZE_OF_FLAGS, 1, sp );
-	    
-      // go to beginning of loop
       continue;
-	    
     }
 	  
@@ -2938,10 +3005,5 @@
     
     *itotnphe += 1;
-    
-    // read next Photon
-	  
-    fread( flag, SIZE_OF_FLAGS, 1, sp );
-	  
-  }  // end while, i.e. found end of event
+  }
 	  
   return(0);
@@ -3248,4 +3310,7 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.28  2001/10/26 16:31:45  blanch
+// Removing several warnings.
+//
 // Revision 1.27  2001/09/05 10:04:33  blanch
 // *** empty log message ***
@@ -3317,4 +3382,7 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.28  2001/10/26 16:31:45  blanch
+// Removing several warnings.
+//
 // Revision 1.27  2001/09/05 10:04:33  blanch
 // *** empty log message ***
