Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 1515)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 1517)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.41 $
+// $Revision: 1.42 $
 // $Author: blanch $ 
-// $Date: 2002-09-04 09:57:42 $
+// $Date: 2002-09-09 16:00:49 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -304,4 +304,11 @@
 //@: table of QE
 static float *QElambda;
+
+//@: table of lightguide = WC;
+static float **WC;
+
+//@: number of datapoints for the WC curve
+static int pointsWC;
+
 //!@}
 
@@ -458,5 +465,4 @@
   float zenfactor=1.0;  //  correction factor calculated from the extinction
 
-  float qTailCut;             //@< Tail Cut value
   int anaPixels;
     
@@ -499,4 +505,7 @@
   MGeomCamMagic camgeom; // structure holding the camera definition
 
+  ct_NPixels=camgeom.GetNumPixels();
+  read_QE();
+  read_WC();
 
   //!@' @#### Definition of variables for |getopt()|.
@@ -618,5 +627,4 @@
   // get different parameters of the simulation
 
-  qTailCut = get_tail_cut();
   addElecNoise = add_elec_noise(&FADC_noise, &Trigger_noise);
   simulateNSB = get_nsb( &meanNSB, &nphe2NSB );  
@@ -699,7 +707,6 @@
   
   log(SIGNATURE,
-      "%s:\n\t%20s: %f\n\t%20s: %f %s\n",
+      "%s:\n\t%20s: %f\n",
       "Parameters",
-      "t0 (Tail-cut)", qTailCut,
       "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB));
   
@@ -768,5 +775,10 @@
     
   anaPixels = get_ana_pixels();
-  anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
+  anaPixels = (anaPixels == -1) ? camgeom.GetNumPixels() : anaPixels;
+
+  // Switched off writing TObject
+
+  MArray::Class()->IgnoreTObjectStreamer();
+  MParContainer::Class()->IgnoreTObjectStreamer();
 
   // initialise ROOT
@@ -1180,5 +1192,5 @@
   // allocate space for PMTs numbers of pixels
   
-  fnpix = new float [ ct_NPixels ];
+  fnpix = new float [ camgeom.GetNumPixels() ];
 
      
@@ -1490,11 +1502,12 @@
 		  
 		  Lev0=1;
+		  Int_t NumImages = Lev1;
 		  if(Lev1==0 && (Write_All_Images || btrigger)){
 		    btrigger= 1;
-		    Lev1=1;
+		    NumImages=1;
 		    Lev0=0;
 		  }
 
-		  for (Int_t ii=0;ii<Lev1;ii++){
+		  for (Int_t ii=0;ii<NumImages;ii++){
 		      if (Write_McTrig){
 			  McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
@@ -1516,9 +1529,11 @@
 		      EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
 		      //   fill pixel information
-		      for(int i=0;i<ct_NPixels;i++){
-			for (j=0;j<FADC_SLICES;j++){
-			  fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
+		      if (Lev1){
+			for(UInt_t i=0;i<camgeom.GetNumPixels();i++){
+			  for (j=0;j<FADC_SLICES;j++){
+			    fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
+			  }
+			  EvtData[iconcount]->AddPixel(i,fadcValues,0);
 			}
-			EvtData[iconcount]->AddPixel(i,fadcValues,0);
 		      }
 		    }
@@ -1612,11 +1627,12 @@
 	    ++ntrigger;
 	  }
+	  Int_t NumImages = Lev1;
 	  Lev0=1;
 	  if (Lev1==0 && Write_All_Images){ 
-	    Lev1=1;
+	    NumImages=1;
 	    Lev0=0;
 	  }
 
-	  for(Int_t ii=0;ii<Lev1;ii++){
+	  for(Int_t ii=0;ii<NumImages;ii++){
 	    //  Loop over different level one triggers
 
@@ -1645,9 +1661,11 @@
 	      //   fill pixel information
 	      
-	      for(int i=0;i<ct_NPixels;i++){
-		for (j=0;j<FADC_SLICES;j++){
-		  fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
+	      if (Lev1){
+		for(UInt_t i=0;i<camgeom.GetNumPixels();i++){
+		  for (j=0;j<FADC_SLICES;j++){
+		    fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
+		  }
+		  EvtData[0]->AddPixel(i,fadcValues,0);
 		}
-		EvtData[0]->AddPixel(i,fadcValues,0);
 	      }
 	    }	    
@@ -1728,5 +1746,5 @@
 	}
 	
-	for (int i=0; i<ct_NPixels; ++i) {
+	for (int i=0; i<camgeom.GetNumPixels(); ++i) {
 	  printf("%d (%d): ", i, npixneig[i]);
 	  for (j=0; j<npixneig[i]; ++i) 
@@ -2260,4 +2278,190 @@
 
 //!-----------------------------------------------------------
+// @name read_QE  
+//                          
+// @desc read QE data
+//
+// @date thu  5 17:59:57 CEST 2002
+//------------------------------------------------------------
+// @function
+
+//!@{
+void 
+read_QE(void){
+  ifstream qefile;
+  char line[LINE_MAX_LENGTH];
+  int i, j, icount;
+  float qe;
+
+  //------------------------------------------------------------
+  // second, pixels' QE
+
+  // try to open the file
+
+  log("read_QE", "Opening the file \"%s\" . . .\n", QE_FILE);
+  
+  qefile.open( QE_FILE );
+  
+  // if it is wrong or does not exist, exit
+  
+  if ( qefile.bad() )
+    error( "read_QE", "Cannot open \"%s\". Exiting.\n", QE_FILE );
+  
+  // read file
+
+  log("read_QE", "Reading data . . .\n");
+
+  i=-1;
+  icount = 0;
+
+  while ( ! qefile.eof() ) {          
+
+    // get line from the file
+
+    qefile.getline(line, LINE_MAX_LENGTH);
+
+    // skip if comment
+
+    if ( *line == '#' )
+      continue;
+
+    // if it is the first valid value, it is the number of QE data points
+
+    if ( i < 0 ) {
+
+      // get the number of datapoints 
+
+      sscanf(line, "%d", &pointsQE);
+      
+      // allocate memory for the table of QEs
+      
+      QE = new float ** [ct_NPixels];
+
+      for ( i=0; i<ct_NPixels; ++i ) {
+        QE[i] = new float * [2];
+        QE[i][0] = new float[pointsQE];
+        QE[i][1] = new float[pointsQE];
+      }
+      
+      QElambda = new float [pointsQE];
+
+      for ( i=0; i<pointsQE; ++i ) {
+        qefile.getline(line, LINE_MAX_LENGTH);
+        sscanf(line, "%f", &QElambda[i]);
+      }
+
+      i=0;
+
+      continue;
+    }
+
+    // get the values (num-pixel, num-datapoint, QE-value)
+    
+    if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 ) 
+      break;
+
+    if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
+         ((j-1) < pointsQE)   && ((j-1) > -1) ) {
+      QE[i-1][0][j-1] = QElambda[j-1];
+      QE[i-1][1][j-1] = qe;
+    }
+
+    if ( i > ct_NPixels) break;
+
+    icount++;
+
+  }
+
+  if(icount/pointsQE < ct_NPixels){
+    error( "read_QE", "The quantum efficiency file is faulty\n  (found only %d pixels instead of %d).\n",
+	   icount/pointsQE, ct_NPixels );
+  }
+
+  // close file
+
+  qefile.close();
+
+  // test QE
+
+  for(icount=0; icount< ct_NPixels; icount++){
+    for(i=0; i<pointsQE; i++){
+      if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
+	  QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
+	error( "read_QE", "The quantum efficiency file is faulty\n  pixel %d, point %d  is % f, %f\n",
+	       icount, i, QE[icount][0][i], QE[icount][1][i] );
+      }
+    }
+  }
+
+  // end
+
+  log("read_QE", "Done.\n");
+}
+//!@}
+
+//!-----------------------------------------------------------
+// @name read_WC  
+//                          
+// @desc read WC data
+//
+// @date thu  5 17:59:57 CEST 2002
+//------------------------------------------------------------
+// @function
+
+//!@{
+void 
+read_WC(void){
+  ifstream wcfile;
+  char line[LINE_MAX_LENGTH];
+  int i;
+
+  //------------------------------------------------------------
+  // Read Light Guides data
+
+  // try to open the file
+
+  log("read_QE", "Opening the file \"%s\" . . .\n", WC_FILE);
+  
+  wcfile.open( WC_FILE );
+  
+  // if it is wrong or does not exist, exit
+  
+  if ( wcfile.bad() )
+    error( "read_WC", "Cannot open \"%s\". Exiting.\n", QE_FILE );
+  
+  // read file
+
+  log("read_WC", "Reading data . . .\n");
+
+  // get line from the file
+
+  wcfile.getline(line, LINE_MAX_LENGTH);
+  
+  // get the number of datapoints 
+    
+  sscanf(line, "%d", &pointsWC);
+    
+  // allocate memory for the table of QEs
+    
+  WC = new float * [2];
+  WC[0] = new float[pointsWC];
+  WC[1] = new float[pointsWC];
+
+  for ( i=0; i<pointsWC; ++i ) {
+    wcfile.getline(line, LINE_MAX_LENGTH);
+    sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
+  }
+    
+  // close file
+
+  wcfile.close();
+
+  // read
+
+  log("read_WC", "Done.\n");
+}
+//!@}
+
+//!-----------------------------------------------------------
 // @name read_pixels  
 //                          
@@ -2900,4 +3104,5 @@
   static int k, ipixnum;
   static float cx, cy, wl, qe, t;
+  static float cu, cv, cw;
   static MCCphoton photon;
   static float **qept;
@@ -2944,4 +3149,5 @@
     fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
 	  
+
     // increase number of photons
 	  
@@ -3050,5 +3256,28 @@
     qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
 
-    // if random > quantum efficiency, reject it
+    // Apply incient angular correction due to Light Guides
+
+    cu=photon.get_u();
+    cv=photon.get_v();
+    cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*pi/180.0;
+    cout<<" Hola "<<cu<<" "<<cv<<" "<<cw<<endl;
+
+    // find data point in the QE table (-> k)
+
+    k = 0; 
+    while (k < pointsWC-1 && WC[0][k] < cw){
+      k++;
+    }
+
+    cout<<" Hola 2 "<<k<<endl;
+
+     // correct the qe with WC data.
+
+    cout<<" Hola 3 "<<qe<<" ";
+
+    qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
+
+    cout<<qe<<endl;
+   // if random > quantum efficiency, reject it
 
     if ( (RandomNumber) > qe ) {
@@ -3378,4 +3607,7 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.41  2002/09/04 09:57:42  blanch
+// Modifications done to use MGeomCam from MARS.
+//
 // Revision 1.40  2002/07/16 16:15:22  blanch
 // A first implementation of the Star field rotation has been done.
@@ -3500,4 +3732,7 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.41  2002/09/04 09:57:42  blanch
+// Modifications done to use MGeomCam from MARS.
+//
 // Revision 1.40  2002/07/16 16:15:22  blanch
 // A first implementation of the Star field rotation has been done.