Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 5247)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 5268)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.72 $
+// $Revision: 1.73 $
 // $Author: moralejo $ 
-// $Date: 2004-10-12 13:39:34 $
+// $Date: 2004-10-13 14:28:59 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -238,4 +238,6 @@
 static  float Spot_y=0.0;
 static  float Spotsigma=0.0;
+
+static  int CalibrationRun = 0;
 
   
@@ -715,13 +717,15 @@
   // get filenames
 
-  for(int ict=0;ict<ct_Number;ict++)
-    {
-      strcpy( inname_CT[ict], get_input_filename(ict) );
-      if (strlen(inname_CT[ict]) == 0)
-	{
-	  printf("\nError: missing input file name for CT id = %d. Exiting.\n\n", ict);
-	  exit(1);
-	}
-    }
+  CalibrationRun = 1;
+  if (! CalibrationRun)
+    for(int ict=0;ict<ct_Number;ict++)
+      {
+	strcpy( inname_CT[ict], get_input_filename(ict) );
+	if (strlen(inname_CT[ict]) == 0)
+	  {
+	    printf("\nError: missing input file name for CT id = %d. Exiting.\n\n", ict);
+	    exit(1);
+	  }
+      }
 
   strcpy( starfieldname, get_starfield_filename() );
@@ -1548,8 +1552,10 @@
 	zenfactor = pow(10., -0.4 * ext[j] );
 	
-	for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
-	  nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
-	  nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
-	}
+	for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++)
+	  {
+	    nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + 
+	      zenfactor * nsbrate_phepns[ui][j];
+	    nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
+	  }
       }
     }
@@ -1675,4 +1681,21 @@
   }
 
+
+  if (CalibrationRun)
+    {
+      DoCalibration(Fadc_CT, Trigger_CT, camgeom, nsb_trigresp, nsb_fadcresp,
+		    &lons, &lons_outer, nsb_phepns, addElecNoise, 
+		    &EvtTree, EvtHeader, McEvt, EvtData);
+
+      HeaderTree.Fill() ;
+
+      outfile_temp.Write();
+      outfile_temp.Close(); 
+
+      cout << endl << "Calibration run finished!" << endl << endl;
+
+      return (0);
+    }
+
   //
   // Read the reflector file with the Cherenkov data
@@ -1681,28 +1704,26 @@
   // select input file 
 
-  if ( Data_From_STDIN ) {
-
+  if ( Data_From_STDIN )
     inputfile[0] = stdin;
 
-  }
-  else{
-    for(int ict=0;ict<ct_Number;ict++){
-
-      log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
-      inputfile[ict] = fopen( inname_CT[ict], "r" );
-
-      if ( inputfile[ict] == NULL ) 
-	error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
-    }
-    
-  }
+  else
+    {
+      for(int ict = 0; ict < ct_Number; ict++)
+	{
+	  log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
+	  inputfile[ict] = fopen( inname_CT[ict], "r" );
+	  
+	  if ( inputfile[ict] == NULL ) 
+	    error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
+	}
+    }
   
   // get signature, and check it
 
-  for(int ict=0;ict<ct_Number;ict++){
-    if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE){
-      exit(1);
-    }
-  }  
+  for(int ict=0;ict<ct_Number;ict++)
+    {
+      if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE)
+	exit(1);
+    }  
 
   // open data file
@@ -2463,5 +2484,5 @@
 		//  Fill the header of this event 
 		//
-		
+
 		EvtHeader[ict]
 		  ->FillHeader ( (UInt_t) (ntshow + nshow) ,  20 ) ; 
@@ -2594,5 +2615,5 @@
 	      }
 	    }
-	    //   We don not count photons out of the camera.	
+	    //   We do not count photons out of the camera.	
 	      
 	    
@@ -2985,5 +3006,6 @@
        <<  SIGNATURE << '\n' << '\n'
        << "Processor of the reflector output\n"
-       << "J C Gonzalez, Jun 1998\n"
+       << "J. C. Gonzalez, Jun 1998\n"
+       << "O. Blanch, A. Moralejo, 2004\n"
        << "##################################################\n\n"
        << flush ;
@@ -4417,4 +4439,364 @@
 
 
+
+//-------------------------------------------------------------
+//
+// Function DoCalibration. A. Moralejo October 2004
+//
+//-------------------------------------------------------------
+
+int DoCalibration(MFadc **Fadc_CT, MTrigger **Trigger_CT,
+		  TObjArray camgeom,
+		  float *nsb_trigresp, float *nsb_fadcresp,
+		  MLons *lons, MLons *lons_outer,
+		  float **nsb_phepns, int addElecNoise,
+		  TTree *EvtTree, MRawEvtHeader **EvtHeader,
+		  MMcEvt **McEvt, MRawEvtData** EvtData)
+{
+
+  int ntshow = 0;
+  int lons_return;
+
+  int *ntotphe = new int[ct_Number];
+  int *ntotphot = new int[ct_Number];
+  float *nphe_from_nsb = new float[ct_Number];
+
+  float timefirst = 0.;
+  float timelast = 0.;
+
+  TArrayC *fadcValues;     //@< the analog Fadc High gain signal for pixels
+  TArrayC *fadcValuesLow;  //@< the analog Fadc Low gain signal for pixels
+
+  // allocate space for FADC info
+  fadcValues    =  new TArrayC(FADC_slices_written);
+  fadcValuesLow =  new TArrayC(FADC_slices_written);
+
+
+  // allocate space for PMTs numbers of pixels
+  float *pheinpix = new float [ct_NPixels];
+
+  for (int calevent = 0; calevent < 1; calevent++) 
+    { 
+      //
+      // Clear Trigger and Fadc
+      //
+      for(int ict = 0; ict < ct_Number; ict++)
+	{
+	  Trigger_CT[ict]->Reset() ; 
+	  Trigger_CT[ict]->ClearFirst();
+	  Trigger_CT[ict]->ClearZero();
+	  Fadc_CT[ict]->Reset() ; 
+
+	  ntotphe[ict] = ntotphot[ict] = 0;
+	  nphe_from_nsb[ict] = 0.;
+
+	}
+
+      ntshow++;
+
+      if((ntshow+1)%100 == 1)
+	log(SIGNATURE, "Event %d\n", ntshow);
+	
+      // Produce the photoelectrons
+      float phot_per_pix = 100.;
+      for(int ict = 0; ict < ct_Number; ict++)
+	produce_calib_phes( (MGeomCam*)(camgeom.UncheckedAt(ict)),
+			    phot_per_pix,
+			    Trigger_CT[ict],
+			    Fadc_CT[ict],
+			    &(ntotphe[ict]),
+			    pheinpix,
+			    &(ntotphot[ict]),
+			    ict);
+
+      // NSB simulation
+
+      if(lons && lons_outer)
+	{
+	  //  Fill trigger and fadc response in the trigger class from the NSB database
+	  
+	  for (int ict = 0; ict < ct_Number; ict++)
+	    {
+	      for (UInt_t ui = 0;
+		  ui < ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
+		  ui++)
+		{
+		  nphe_from_nsb[ict] += nsb_phepns[ict][ui];
+
+		  if (nsb_phepns[ict][ui] > 0.0)
+		    {
+		      if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD() > 
+			 (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD())
+			lons_return = lons_outer->GetResponse(nsb_phepns[ict][ui],0.01,
+							      & nsb_trigresp[0],
+							      & nsb_fadcresp[0]);
+		      else
+			lons_return = lons->GetResponse(nsb_phepns[ict][ui],0.01,
+							& nsb_trigresp[0],
+							& nsb_fadcresp[0]);
+
+		      if (lons_return == 0)
+			{
+			  cout << "Exiting.\n";
+			  exit(1);
+			}
+
+		      Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
+		    }
+		}
+	    }
+	  
+	}// end if(simulateNSB) ...
+
+
+      //++++++++++++++++++++++++++++++++++++++++++++++++++
+      // at this point we have a camera full of
+      // ph.e.s
+      //--------------------------------------------------
+	  
+      //   now the noise of the electronic 
+      //   (preamps, optical transmission,..)  is introduced. 
+      //   
+      
+      for (int ict = 0; ict < ct_Number; ict++) 
+	{
+	  if (addElecNoise)
+	    Fadc_CT[ict]->ElecNoise();
+
+	  //  now a shift in the fadc signal due to the pedestals is 
+	  //  introduced
+	  //  This is done inside the class MFadc by the method Pedestals
+
+	  Fadc_CT[ict]->Pedestals();
+
+	  // Set the trigger time. The 3 ns are such that the pulse appears well 
+	  // centered if the FADC range is the usual 50 ns.
+	  // If you want to shift the pulse position, do not change the value here,
+	  // use the option trigger_delay in the input card instead.
+
+	  Fadc_CT[ict]->TriggeredFadc(3.);
+
+	  //
+	  // Fill the event header information
+	  //
+	  EvtHeader[ict]->FillHeader((UInt_t) calevent, 0);
+
+	  McEvt[ict]->Fill( 0, 0, 0., -1.0, -1.0, -1.0, 0., 0., 0., 0., 0.,
+			    0., 0., 0., timefirst, timelast, 0., 0., 0.,
+			    0., 0., 0., 0., 0, 0, 0,
+			    (UInt_t) ntotphot[ict], 
+			    (UInt_t) ntotphe[ict],
+			    (UInt_t) nphe_from_nsb[ict]+ ntotphe[ict],
+			    0., 0., 0.);
+
+	  //
+	  // Fill the FADC information
+	  //
+	  for(UInt_t ui=0;
+	      ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
+	      ui++)
+	    {
+	      for (int jslice = 0; jslice < FADC_slices_written; jslice++)
+		{
+		  fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(ui, jslice),jslice);
+		  fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(ui,jslice),jslice);
+		}
+	      EvtData[ict]->AddPixel(ui,fadcValues,0);
+	      EvtData[ict]->AddPixel(ui,fadcValuesLow,kTRUE);
+	    }
+	}
+
+      EvtTree->Fill(); 
+
+      //    clear all
+      for(int ict=0;ict<ct_Number;ict++)
+	{
+	  EvtHeader[ict]->Clear() ;
+	  EvtData[ict]->ResetPixels(0,0);
+	  McEvt[ict]->Clear() ; 
+	}
+    }
+
+  return(0);
+}
+
+
+//------------------------------------------------------------
+// 
+// Function produce_calib_phes, A. Moralejo Oct 2004
+//    
+//------------------------------------------------------------
+
+int produce_calib_phes( MGeomCam *camgeom,      // The camera layout
+			float     phot_per_pix, // Average number of photons per inner pixel
+			MTrigger *trigger,
+			MFadc    *fadc,
+			int      *itotnphe,     // total number of produced photoelectrons		
+			float    *nphe,         // number of photoelectrons in each pixel
+			int      *nphot,        // total number of photons in all pixels
+			int       ict           // Telescope that is being analised to get the right QE.
+			)
+{
+
+  int   ipixnum;
+  float cx, cy, wl, time, phi, phi_deg, qe;
+  float **qept;
+  float radius_mm, focal_dist_mm;
+
+  // reset variables
+
+  for ( uint i = 0; i < camgeom->GetNumPixels(); i++ )
+    nphe[i] = 0.0;
+
+  *itotnphe = 0;
+  *nphot = 0;
+
+  TRandom random;
+  random.SetSeed((Int_t)(get_seeds(0)*get_seeds(1)));
+
+  //
+  // Create photons and "map" them into the pixels
+  //
+  
+  // Maximum radius of camera:
+  radius_mm = camgeom->GetMaxRadius();
+
+  // Distance from center of mirror dish to camera plane:
+  focal_dist_mm = camgeom->GetCameraDist()*1000;
+
+  // Cosine of the angle between telescope axis and line from center of mirror
+  // dish to the edge of the camera:
+  float cos_epsilon_max = cos(atan2(radius_mm, focal_dist_mm));
+
+
+  // Calculate total number of photons to be produced.
+  int total_photons = (int) (phot_per_pix * 3.14159265 * radius_mm * radius_mm / 
+			     (*camgeom)[0].GetA());
+
+  // loop over the photons
+  
+  for (int iph = 0; iph < total_photons; iph++) 
+    {
+
+      time = 0.;
+
+      // Obtain photon coordinates on the camera. We assume a point source of light placed
+      // in the center of the mirror dish.
+
+      // polar angle 
+      float psi = RandomNumber * 2 * 3.14159265;
+      // angle between the telescope axis and the photon trajectory.
+      float epsilon = acos(1.-RandomNumber*(1.-cos_epsilon_max));
+      float tanepsilon = tan(epsilon);
+
+      cx = focal_dist_mm*tanepsilon*cos(psi); // mm
+      cy = focal_dist_mm*tanepsilon*sin(psi); // mm
+
+      // Angle between photon trajectory and camera plane:
+      phi = 3.14159265/2.-epsilon; // rad
+
+      // wavelength
+	  
+      wl = 500.; // nm
+
+
+      if( (wl > WAVEBANDBOUND6) || (wl < WAVEBANDBOUND1) || 
+	  (sqrt(cx*cx + cy*cy) > radius_mm ) )
+	continue;
+
+      //
+      // Pixelization
+      //
+      ipixnum = bpoint_is_in_pix(cx, cy, camgeom);
+
+      // -1 = the photon is in none of the pixels
+      //  0 = the phton is in the central pixel, which is not used
+      if (ipixnum==-1 || ipixnum==0)
+        continue;
+
+      // increase number of photons within pixels
+      *nphot += 1;
+
+      //
+      // QE simulation
+      //    
+      // set pointer to the QE table of the relevant pixel
+      
+      qept = (float **)QE[ict][ipixnum];
+
+      // check if wl is inside table; outside the table, QE is assumed to be zero
+    
+      if((wl < qept[0][0]) || (wl > qept[0][pointsQE[ict]-1]))
+	continue;
+	    
+
+      // find data point in the QE table (-> k)
+      int k = 1; // start at 1 because the condition was already tested for 0
+      while (k < pointsQE[ict]-1 && qept[0][k] < wl)
+	k++;
+
+      // calculate the qe value between 0. and 1.
+
+      qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
+
+
+      //
+      // Apply incient angular correction due to Light Guides, plexiglas,
+      // 1st dynode collection efficiency, double crossings... etc. 
+      // This information is contained in the file Data/LightCollection.dat,
+      // and has been read into the array WC (which stands for "Winston Cones")
+      //
+      phi_deg = phi*180./3.14159265;
+
+      // find data point in the WC table (-> k)
+
+      if ( camgeom->GetPixRatio(ipixnum) < 1.) // => Pixel is outer pixel
+	{
+	  k = 0; 
+	  while (k < pointsWC-1 && WC_outer[0][k] < phi_deg)
+	    k++;
+	  // correct the qe with WC data.
+	  qe = qe*lin_interpol(WC_outer[0][k-1], WC_outer[1][k-1], 
+			       WC_outer[0][k],   WC_outer[1][k], phi_deg);
+	}
+
+      else    // => Pixel is inner pixel
+	{
+	  k = 0; 
+	  while (k < pointsWC-1 && WC[0][k] < phi_deg)
+	    k++;
+	  // correct the qe with WC data.
+	  qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], phi_deg);
+	}
+
+
+      // if random > quantum efficiency, reject it
+
+      if ( (RandomNumber) > qe ) {
+	continue;
+      }
+	  
+      //
+      // The photon has produced a photo electron
+      //
+
+      // increment the number of photoelectrons in the relevant pixel
+	  
+      nphe[ipixnum] += 1.;
+
+      // store the new photoelectron
+
+      fadc->Fill(ipixnum,(time) , 
+		 trigger->FillShow(ipixnum, time),
+		 !((*camgeom)[ipixnum].GetD()>(*camgeom)[0].GetD()));
+    
+      *itotnphe += 1;
+    }
+
+  return(0);
+
+}
+
+
 // @endcode
 
@@ -4426,4 +4808,79 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.72  2004/10/12 13:39:34  moralejo
+//
+// Lots of changes intended to make it possible to select the FADC sampling
+// frequency from the input card, through the command fadc_GHz. The most
+// important ones are the following:
+//
+//  - Removed FADC_SLICES de Mars/mmc/Mdefine.h
+//    Already defined in MFadcDefine.h!
+//
+//  - Replaced fixed numbers in array dimensions of starresponse.cxx
+//
+//  - Added in MFadc.cxx and MFadc.hxx (Float_t) casts to initializations
+//    of single phe response array
+//
+//  - IMPORTANT: Fixed MFadc::GetResponse  -> The returned single phe response
+//    had only RESPONSE_SLICES (which is actually for the trigger branch),
+//    whereas it should have RESPONSE_SLICES_MFADC. Fixed the same confusion
+//    in other two points of the code (filling of the FADC for the signal),
+//    in Fill() and FillOuter().
+//
+//  - RESPONSE_SLICES_MFADC is now eliminated, since this quantity is now
+//    decided by MFadc depending on other parameters.
+//
+//  - Fixed problem in starresponse.cxx due to which the histograms fadcresp
+//    and fadcbase in the root output were actually identical.
+//
+//  - Changed starresponse.cxx such that above 1 phe/ns/pixel the precision
+//    of the database is forced to be 0.1, and below 1, to 0.01 phe/ns/pixel
+//
+//  - In Camera/creadparam.cxx: Now unkown tokens cause the program to stop,
+//    instead of being simply skipped as it was until now.
+//
+//  - Fixed error in MStarLight::StoreHisto. TH1::SetBinContent uses bin numbers
+//    from 1 to nbins (in old code 0 to nbins-1 was assumed).
+//
+//  - In MLons.cxx: use memcpy to copy pieces of the database into the FADC and
+//    trigger branches, to (hopefully) speed up execution. For this I had to add
+//    2 getter functions in MStarLight.hxx
+//
+//  - Everywhere: changed the shape parameter for trigger and FADC response to
+//    be an Int. Changed version of NSB database from 1002 to 1003.
+//
+//  - In MTrigger.cxx, changed all initializations of SlicesFirst and
+//    SlicesSecond to 0 (instead of -50 as it was before). This controls the
+//    time of trigger. If no trigger happened (like when making pedestal files)
+//    the time was negative and the array index used to retrieve the noise from
+//    the database was negative, resulting in "discontinuities" in the noise
+//    (half-photoelectrons for instance).
+//
+//  - In MStarLight changed fBinsTrig  and fBinsFadc from Float_t to Int_t
+//
+//  - Replace WIDTH_FADC_TIMESLICE by FADC_SLICES_PER_NSEC (which is the
+//    inverse of the former)
+//
+//  - Replaced SLICES_PER_NSEC by TRIG_SLICES_PER_NSEC
+//
+//  - TRIGBINS eliminated. It depends on other two values
+//        TRIGBINS = TIMERANGE*TRIG_SLICES_PER_NSEC
+//
+//  - FADCBINS eliminated. It depends on other two values
+//        FADCBINS = TIMERANGE*FADC_SLICES_PER_NSEC
+//
+//  - MTriggerDefine.h  Changed RESPONSE_SLICES to RESPONSE_SLICES_TRIG
+//
+//  - Added to the MLons constructor an argument regarding the FADC sampling
+//    frequency
+//
+//  - MFadc: now the number of response slices for the FADC simulation is
+//    decided by the program according to the other parameters.
+//
+//  - MStarLight: Added arguments (fadc_slices_per_ns, response_slices_fadc)
+//    to constructor.
+//
+//  - creadparam.cxx, camera.cxx   Changed default value of digital_noise to 0.
+//
 // Revision 1.71  2004/09/17 09:20:52  moralejo
 //
