Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 1694)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 1706)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.50 $
+// $Revision: 1.51 $
 // $Author: blanch $ 
-// $Date: 2003-01-07 16:33:31 $
+// $Date: 2003-01-14 13:40:17 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -443,4 +443,7 @@
   int addElecNoise;           //@< Will we add ElecNoise?
 
+  float riseDiskThres;
+  float secureDiskThres;
+
   int simulateNSB;            //@< Will we simulate NSB?
   int nphe2NSB;               //@< From how many phe will we simulate NSB?
@@ -472,5 +475,8 @@
 
   int anaPixels;
-    
+
+  int nstoredevents = 0;
+  int flagstoring = 0;
+
   float fCorrection;          //@< Factor to apply to pixel values (def. 1.)
 
@@ -609,4 +615,6 @@
   Individual_Thres_Pixel = get_indi_thres_pixel();
 
+  get_secure_threhold(&riseDiskThres,&secureDiskThres);
+
   get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm);
 
@@ -811,4 +819,11 @@
   Trigger.CheckThreshold(&qThreshold[0]);
 
+  // Set flag in pixel 0 (not used for trigger) that indicates if secure pixel
+  // is active: secureDiskThres*10000+riseDiskThres
+
+  if(riseDiskThres>0.0)
+    qThreshold[0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres;
+  else
+    qThreshold[0]=0.0;
 
   //  Initialise McTrig information class if we want to save trigger informtion
@@ -963,5 +978,5 @@
 
   RunHeader->SetMagicNumber(kMagicNumber);
-  RunHeader->SetFormatVersion(3);
+  RunHeader->SetFormatVersion(4);
   RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
   RunHeader->SetRunType(256);
@@ -1288,11 +1303,7 @@
 	if (reflector_file_version<6){
 	  fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
-	  thetashw = mcevth.get_theta();
-	  phishw = mcevth.get_phi();
 	}
 	else{
 	  fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile );
-	  thetashw = mcevth_2.get_theta();
-	  phishw = mcevth_2.get_phi();
 	}
 
@@ -1341,48 +1352,78 @@
 	m1 = sin(thetashw)*sin(phishw);
 	n1 = cos(thetashw);
-	
-	// read the deviation of the telescope with respect to the shower
-	
-	if (reflector_file_version<6)
+
+	if (reflector_file_version<6){
+
+	  mcevth.get_core(&coreX, &coreY);	
+	  
+	  // read the deviation of the telescope with respect to the shower
 	  mcevth.get_deviations ( &thetaCT, &phiCT );
-	else
-	  mcevth_2.get_deviations ( &thetaCT, &phiCT );
-
-	if ( (thetaCT == 0.) && (phiCT == 0.) ) {
+
+	  if ( (thetaCT == 0.) && (phiCT == 0.) ) {
 	  
-	  // CT was looking to the source (both lines are parallel)
-	  // therefore, we calculate the impact parameter as the distance 
-	  // between the CT axis and the core position
-
-	  if (reflector_file_version<6)
-	    mcevth.get_core(&coreX, &coreY);
-	  else
-	    mcevth_2.get_core(&coreX, &coreY);
-
-	  impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
+	    // CT was looking to the source (both lines are parallel)
+	    // therefore, we calculate the impact parameter as the distance 
+	    // between the CT axis and the core position
+
+	    impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
+	    thetaCT += thetashw;
+	    phiCT   += phishw;
+	    
+	  } else {
+	    
+	    // the shower comes off-axis
+	    
+	    // obtain with this the final direction of the CT
+	    
+	    thetaCT += thetashw;
+	    phiCT   += phishw;
+	    
+	    // calculate vector for telescope
+	    
+	    l2 = sin(thetaCT)*cos(phiCT);
+	    m2 = sin(thetaCT)*sin(phiCT);
+	    n2 = cos(thetaCT);
+	    
+	    num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
+	    den = (SQR(l1*m2 - l2*m1) + 
+		   SQR(m1*n2 - m2*n1) + 
+		   SQR(n1*l2 - n2*l1));
+	    den = sqrt(den);
+	    
+	    impactD = fabs(num)/den;
+	    
+	  }
+	}
+	else{
+	  mcevth_2.get_core(&coreX, &coreY);	
+	  thetaCT=mcevth_2.get_theta_CT();
+	  phiCT=mcevth_2.get_phi_CT();
+	  if ( (thetaCT == thetashw) && (phiCT == phishw) ) {
 	  
-	} else {
-	  
-	  // the shower comes off-axis
-	  
-	  // obtain with this the final direction of the CT
-	  
-	  thetaCT += thetashw;
-	  phiCT   += phishw;
-	  
-	  // calculate vector for telescope
-	  
-	  l2 = sin(thetaCT)*cos(phiCT);
-	  m2 = sin(thetaCT)*sin(phiCT);
-	  n2 = cos(thetaCT);
-	  
-	  num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
-	  den = (SQR(l1*m2 - l2*m1) + 
-		 SQR(m1*n2 - m2*n1) + 
-		 SQR(n1*l2 - n2*l1));
-	  den = sqrt(den);
-	  
-	  impactD = fabs(num)/den;
-	  
+	    // CT was looking to the source (both lines are parallel)
+	    // therefore, we calculate the impact parameter as the distance 
+	    // between the CT axis and the core position
+
+	    impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
+	    
+	  } else {
+	    
+	    // the shower comes off-axis
+	    
+	    // calculate vector for telescope
+	    
+	    l2 = sin(thetaCT)*cos(phiCT);
+	    m2 = sin(thetaCT)*sin(phiCT);
+	    n2 = cos(thetaCT);
+	    
+	    num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
+	    den = (SQR(l1*m2 - l2*m1) + 
+		   SQR(m1*n2 - m2*n1) + 
+		   SQR(n1*l2 - n2*l1));
+	    den = sqrt(den);
+	    
+	    impactD = fabs(num)/den;
+	    
+	  }
 	}
 
@@ -1465,5 +1506,4 @@
 			     rho);
 	    
- 
 	  }
 
@@ -1551,11 +1591,20 @@
 	  // Set to zero the flag to know if some conditon has triggered
 	  btrigger=0;
+	  flagstoring = 0;
 	  //  Loop over trigger threshold
 	  int iconcount;
 	  for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
-	    for (i=0;i<CAMERA_PIXELS;i++)
+	    for (i=0;i<CAMERA_PIXELS;i++){
 	      fpixelthres[i]=
 		((Float_t)(fthrescount)>=qThreshold[i])?
 		(Float_t)(fthrescount):qThreshold[i];
+
+	      // Rise the discrimnator threshold to avoid huge rates
+
+	      if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
+		for(int ii=0;ii<CAMERA_PIXELS;ii++)
+		  if( nsb_phepns_rotated[ii]>riseDiskThres)
+		    fpixelthres[ii]=secureDiskThres;
+	    }
 	    Trigger.SetThreshold(fpixelthres);
 
@@ -1665,4 +1714,8 @@
 	    //
 
+	    if(!flagstoring)
+	      nstoredevents++;
+	    flagstoring = 1;
+
 	    if (Write_McEvt) {
 	      Float_t ftime, ltime;
@@ -1681,6 +1734,6 @@
 			     mcevth.get_coreY(),
 			     impactD,
-			     mcevth_2.get_theta_CT(),
-			     mcevth_2.get_phi_CT(),
+			     phiCT,
+			     thetaCT,
 			     ftime,
 			     ltime,
@@ -1701,5 +1754,5 @@
 	      else{
 		Float_t Nmax, t0, tmax, a, b, c, chi2;
-		mcevth.get_times(&ftime, &ltime);
+		mcevth_2.get_times(&ftime, &ltime);
 		chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
 		McEvt->Fill( mcevth_2.get_evt_number(),
@@ -1715,6 +1768,6 @@
 			     mcevth_2.get_coreY(),
 			     impactD,
+			     mcevth_2.get_phi_CT(),
 			     mcevth_2.get_theta_CT(),
-			     mcevth_2.get_phi_CT(),
 			     ftime,
 			     ltime,
@@ -1761,4 +1814,12 @@
 	  for (int i=0;i<CAMERA_PIXELS;i++)
 	    fpixelthres[i]=qThreshold[i];
+
+	  // Rise the discrimnator threshold to avoid huge rates
+	  if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
+	    for(int ii=0;ii<CAMERA_PIXELS;ii++){
+	      if( nsb_phepns_rotated[ii]>riseDiskThres)
+		fpixelthres[ii]=secureDiskThres;
+	    }
+	  
 	  Trigger.SetThreshold(fpixelthres);
 				 
@@ -1802,4 +1863,6 @@
 	    //
 	    fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
+
+	    nstoredevents++;
 
 	    if (Write_McTrig){
@@ -1851,6 +1914,6 @@
 			     mcevth.get_coreY(),
 			     impactD,
-			     mcevth_2.get_theta_CT(),
-			     mcevth_2.get_phi_CT(),
+			     phiCT,
+			     thetaCT,
 			     ftime,
 			     ltime,
@@ -1871,5 +1934,5 @@
 	      else{
 		Float_t Nmax, t0, tmax, a, b, c, chi2;
-		mcevth.get_times(&ftime, &ltime);
+		mcevth_2.get_times(&ftime, &ltime);
 		chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
 		McEvt->Fill( mcevth_2.get_evt_number(),
@@ -1885,6 +1948,6 @@
 			     mcevth_2.get_coreY(),
 			     impactD,
+			     mcevth_2.get_phi_CT(),
 			     mcevth_2.get_theta_CT(),
-			     mcevth_2.get_phi_CT(),
 			     ftime,
 			     ltime,
@@ -2022,17 +2085,33 @@
 
   get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
-  get_teles_axis(&telestheta,&telesphi);
+  if (reflector_file_version<6){
+    get_teles_axis(&telestheta,&telesphi);
+    mcevth.get_theta_range(&shthetamin, &shthetamax);
+    mcevth.get_phi_range(&shphimin,&shphimax);
+    mcevth.get_theta_range(&shthetamin, &shthetamax);
+    mcevth.get_phi_range(&shphimin,&shphimax);
+    corsika=UInt_t(mcevth.get_VersionPGM()*1000);
+    for (int i=0; i< 10;i++)
+      heights[i]=mcevth.get_HeightLev (i);
+    rnum=mcevth.get_RunNumber();
+  }
+  else{
+    telestheta=mcevth_2.get_theta_CT();
+    telesphi=mcevth_2.get_phi_CT();
+    mcevth_2.get_theta_range(&shthetamin, &shthetamax);
+    mcevth_2.get_phi_range(&shphimin,&shphimax);
+    corsika=UInt_t(mcevth_2.get_VersionPGM()*1000);
+    for (int i=0; i< 10;i++)
+      heights[i]=mcevth_2.get_HeightLev (i);
+    rnum=mcevth_2.get_RunNumber();
+  }
+
   get_source_off(&sofftheta,&soffphi);
-  mcevth.get_theta_range(&shthetamin, &shthetamax);
-  mcevth.get_phi_range(&shphimin,&shphimax);
-  corsika=UInt_t(mcevth.get_VersionPGM()*1000);
   if(!Trigger_Loop) icontrigger=0;
   time (&ltime);
   ftime = ((Float_t)ltime)/1000;
-  for (int i=0; i< 10;i++)
-    heights[i]=mcevth.get_HeightLev (i);
-  rnum=mcevth.get_RunNumber();
-
-  McRunHeader->Fill(rnum,
+
+  if (reflector_file_version<6)
+    McRunHeader->Fill(rnum,
 		    (UInt_t) 0,
 		    mcevth.get_DateRun(),
@@ -2047,5 +2126,5 @@
 		    CAMERA_PIXELS,
 		    (UInt_t)ntshow,
-		    (UInt_t)ntrigger,
+		    (UInt_t)nstoredevents,
 		    0,
 		    sfRaH,
@@ -2070,7 +2149,52 @@
 		    heights,
 		    corsika,
-		    (UInt_t)(REFL_VERSION_A*100),
+		    (UInt_t)(reflector_file_version*100),
 		    (UInt_t)(VERSION*100),
 		    0);
+  else
+    McRunHeader->Fill(rnum,
+		    (UInt_t) 0,
+		    mcevth_2.get_DateRun(),
+		    ftime,
+		    icontrigger,
+		    !Write_All_Images,
+		    Write_McEvt,
+		    Write_McTrig,
+		    Write_McFADC,
+		    Write_RawEvt,
+		    addElecNoise,
+		    CAMERA_PIXELS,
+		    (UInt_t)ntshow,
+		    (UInt_t)nstoredevents,
+		    0,
+		    sfRaH,
+		    sfRaM,
+		    sfRaS,
+		    sfDeD,
+		    sfDeM,
+		    sfDeS,
+		    meanNSB,
+		    telestheta,
+		    telesphi,
+		    sofftheta,
+		    soffphi,
+		    shthetamax,
+		    shthetamin,
+		    shphimax,
+		    shphimin,
+		    mcevth_2.get_CWaveLower(),
+		    mcevth_2.get_CWaveUpper(),
+		    mcevth_2.get_slope(),
+		    1,
+		    heights,
+		    corsika,
+		    (UInt_t)(reflector_file_version*100),
+		    (UInt_t)(VERSION*100),
+		    0);
+  // Fill some missing values for MRawRunHeader
+
+  RunHeader->SetRunNumber(rnum);
+  RunHeader->SetNumEvents(nstoredevents);
+
 
   //  Store qe for each PMT in output file
@@ -4028,4 +4152,10 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.50  2003/01/07 16:33:31  blanch
+// Star Field Rotation has been implemented by Raul Orduna. Now there is a
+// rotation for each shower. It is done by a non enter pixel rotation assuming
+// a circular symetry of the camera. It is not exact but it is accurate enough and
+// much faster.
+//
 // Revision 1.49  2002/12/13 10:04:07  blanch
 // *** empty log message ***
@@ -4179,4 +4309,10 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.50  2003/01/07 16:33:31  blanch
+// Star Field Rotation has been implemented by Raul Orduna. Now there is a
+// rotation for each shower. It is done by a non enter pixel rotation assuming
+// a circular symetry of the camera. It is not exact but it is accurate enough and
+// much faster.
+//
 // Revision 1.49  2002/12/13 10:04:07  blanch
 // *** empty log message ***
