Changeset 5268 for trunk/MagicSoft/Simulation/Detector

Timestamp:

10/13/04 15:28:59 (20 years ago)

Author:

moralejo

Message:

*** empty log message ***

File:

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

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

@@ -21 +21 @@
 //
 // $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 +238 @@
 static  float Spot_y=0.0;
 static  float Spotsigma=0.0;
+
+static  int CalibrationRun = 0;
 
   
@@ -715 +717 @@
   // 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 +1552 @@
         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 +1681 @@
   }
 
+
+  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 +1704 @@
   // 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 +2484 @@
                 //  Fill the header of this event 
                 //
-                
+
                 EvtHeader[ict]
                   ->FillHeader ( (UInt_t) (ntshow + nshow) ,  20 ) ; 
@@ -2594 +2615 @@
               }
             }
-            //   We don not count photons out of the camera.    
+            //   We do not count photons out of the camera.     
               
             
@@ -2985 +3006 @@
        <<  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 +4439 @@
 
 
+
+//-------------------------------------------------------------
+//
+// 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 +4808 @@
 //
 // $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
 //