Changeset 1517 for trunk/MagicSoft

Timestamp:

09/09/02 17:00:49 (22 years ago)

Author:

blanch

Message:

Statement has been included to avoid writting to disk MParContainer and MArray.
It has also been added the effect of the WC, the actual values must be added,
once they are measured.

File:

• trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (modified) (19 diffs)

trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx

@@ -21 +21 @@
 //
 // $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 +304 @@
 //@: table of QE
 static float *QElambda;
+
+//@: table of lightguide = WC;
+static float **WC;
+
+//@: number of datapoints for the WC curve
+static int pointsWC;
+
 //!@}
 
@@ -458 +465 @@
   float zenfactor=1.0;  //  correction factor calculated from the extinction
 
-  float qTailCut;             //@< Tail Cut value
   int anaPixels;
     
@@ -499 +505 @@
   MGeomCamMagic camgeom; // structure holding the camera definition
 
+  ct_NPixels=camgeom.GetNumPixels();
+  read_QE();
+  read_WC();
 
   //!@' @#### Definition of variables for |getopt()|.
@@ -618 +627 @@
   // get different parameters of the simulation
 
-  qTailCut = get_tail_cut();
   addElecNoise = add_elec_noise(&FADC_noise, &Trigger_noise);
   simulateNSB = get_nsb( &meanNSB, &nphe2NSB );  
@@ -699 +707 @@
   
   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 +775 @@
     
   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 +1192 @@
   // allocate space for PMTs numbers of pixels
   
-  fnpix = new float [ ct_NPixels ];
+  fnpix = new float [ camgeom.GetNumPixels() ];
 
      
@@ -1490 +1502 @@
                   
                   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 +1529 @@
                       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 +1627 @@
             ++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 +1661 @@
               //   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 +1746 @@
         }
         
-        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 +2278 @@
 
 //!-----------------------------------------------------------
+// @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 +3104 @@
   static int k, ipixnum;
   static float cx, cy, wl, qe, t;
+  static float cu, cv, cw;
   static MCCphoton photon;
   static float **qept;
@@ -2944 +3149 @@
     fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
           
+
     // increase number of photons
           
@@ -3050 +3256 @@
     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 +3607 @@
 //
 // $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 +3732 @@
 //
 // $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.