Changeset 1512


Ignore:
Timestamp:
09/04/02 10:57:42 (22 years ago)
Author:
blanch
Message:
Modifications done to use MGeomCam from MARS.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx

    r1419 r1512  
    2121//
    2222// $RCSfile: camera.cxx,v $
    23 // $Revision: 1.40 $
     23// $Revision: 1.41 $
    2424// $Author: blanch $
    25 // $Date: 2002-07-16 16:15:22 $
     25// $Date: 2002-09-04 09:57:42 $
    2626//
    2727////////////////////////////////////////////////////////////////////////
     
    7272#include "MMcTrigHeader.hxx"
    7373#include "MMcFadcHeader.hxx"
     74#include "MGeomCamMagic.h"
     75#include "MGeomPix.h"
    7476
    7577/*!@"
     
    495497  UInt_t corsika = 5200 ;
    496498
    497 
    498   struct camera cam; // structure holding the camera definition
     499  MGeomCamMagic camgeom; // structure holding the camera definition
    499500
    500501
     
    759760
    760761  read_ct_file();
    761 
    762   // read camera setup
    763 
    764   read_pixels(&cam);
    765762
    766763  Int_t Lev0, Lev1;
     
    11121109   
    11131110    k = produce_nsbrates( starfieldname,
    1114                           &cam,
     1111                          &camgeom,
    11151112                          nsbrate_phepns);
    11161113
     
    11221119    // calculate diffuse rate correcting for the pixel size
    11231120
    1124     for(i=0; i<cam.inumpixels; i++){
    1125       diffnsb_phepns[i] = meanNSB *
    1126         cam.dpixsizefactor[i] * cam.dpixsizefactor[i];
     1121    const double size_inner = camgeom[0].GetR();
     1122     
     1123    for(UInt_t ui=0; ui<camgeom.GetNumPixels(); ui++){
     1124      const double actual_size = camgeom[ui].GetR();
     1125      diffnsb_phepns[ui] = meanNSB *
     1126        actual_size*actual_size/(size_inner*size_inner);
    11271127    }
    11281128
     
    11341134        zenfactor = pow(10., -0.4 * ext[j] );
    11351135       
    1136         for(i=0; i<cam.inumpixels;i++){
    1137             nsb_phepns[i]+=diffnsb_phepns[i]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[i][j];
     1136        for(UInt_t ui=0; ui<camgeom.GetNumPixels();ui++){
     1137            nsb_phepns[ui]+=diffnsb_phepns[ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
    11381138        }
    11391139    }
     
    13271327
    13281328        k = produce_phes( inputfile,
    1329                           &cam,
     1329                          &camgeom,
    13301330                          WAVEBANDBOUND1,
    13311331                          WAVEBANDBOUND6,
     
    13491349
    13501350          //  Fill trigger and fadc response in the trigger class from the database
    1351           for(i=0;i<cam.inumpixels;i++)
    1352             if(nsb_phepns[i]>0.0){
    1353               k=lons.GetResponse(nsb_phepns[i],0.1,
     1351          for(UInt_t ui=0;ui<camgeom.GetNumPixels();ui++)
     1352            if(nsb_phepns[ui]>0.0){
     1353              k=lons.GetResponse(nsb_phepns[ui],0.1,
    13541354                                 & nsb_trigresp[0],& nsb_fadcresp[0]);
    13551355              if(k==0){
     
    13571357                exit(1);
    13581358              }
    1359               Trigger.AddNSB(i,nsb_trigresp);
    1360               fadc.AddSignal(i,nsb_fadcresp);
     1359              Trigger.AddNSB(ui,nsb_trigresp);
     1360              fadc.AddSignal(ui,nsb_fadcresp);
    13611361            }
    13621362        }// end if(simulateNSB && nphe2NSB<=inumphe) ...
     
    13871387        cout << "Total number of phes: " << inumphe<<" + ";
    13881388        inumphensb=0;
    1389         for (i=0;i<cam.inumpixels;i++){
    1390           inumphensb+=nsb_phepns[i]*TOTAL_TRIGGER_TIME;
     1389        for (UInt_t ui=0;ui<camgeom.GetNumPixels();ui++){
     1390          inumphensb+=nsb_phepns[ui]*TOTAL_TRIGGER_TIME;
    13911391        }
    13921392        cout<<inumphensb<<endl;
     
    27512751
    27522752/******** bpoint_is_in_pix() ***************************************/
    2753 /*
    2754 int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
    2755   // return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system)
    2756   // the pixel is assumed to be a "closed set"
    2757 
    2758   double a, b; // a = length of one of the edges of one pixel, b = half the width of one pixel
    2759   double c, xx, yy; // auxiliary variable
    2760 
    2761   b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
    2762   a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
    2763   c = 1./sqrt(3.);
    2764 
    2765   if((ipixnum < 0)||(ipixnum >= pcam->inumpixels)){
    2766     fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
    2767     fprintf(stderr, "Exiting.\n");
    2768     exit(203);
    2769   }
    2770   xx = dx - pcam->dxc[ipixnum];
    2771   yy = dy - pcam->dyc[ipixnum];
    2772 
    2773   if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) ||
    2774      ((0. <  xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a)))   ){
    2775     return(TRUE); // point is inside
    2776   }
    2777   else{
    2778     return(FALSE); // point is outside
    2779   }
    2780 }
    2781 */
    27822753
    27832754#define sqrt13 0.577350269 // = 1./sqrt(3.)
    27842755#define sqrt3  1.732050807 // = sqrt(3.)
    27852756
    2786 int bpoint_is_in_pix(double dx, double dy, struct camera *pcam)
     2757int bpoint_is_in_pix(double dx, double dy, MGeomCamMagic *pgeomcam)
    27872758{
    27882759    /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
     
    27942765     */
    27952766
    2796     const int    num  = pcam->inumpixels;
    2797     const double diam = pcam->dpixdiameter_cm;
    2798 
    2799     const double diam2 = diam/2;
    2800     const double diam3 = diam/sqrt3;
    2801 
    2802     for (int i=0; i<num; i++)
     2767    const int    numN  = pgeomcam->GetNumPixels();
     2768
     2769    for (int i=0; i<numN; i++)
    28032770    {
    2804         const double size = pcam->dpixsizefactor[i];
    2805 
    2806         const double b = diam2 * size;
    2807         const double a = diam3 * size;
    2808 
    2809         const double xx = dx - pcam->dxc[i];
    2810         const double yy = dy - pcam->dyc[i];
    2811 
    2812         if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
    2813            ((0. <  xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a)))   ){
    2814 
    2815             return i; // inside i
     2771      MGeomPix &pixel = (*pgeomcam)[i];
     2772      const double b = pixel.GetR()/2;
     2773      const double a = pixel.GetR()/sqrt3;
     2774     
     2775      const double xx = dx - pixel.GetX();
     2776      const double yy = dy - pixel.GetY();
     2777     
     2778      if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
     2779         ((0. <  xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a)))   ){
     2780
     2781        return i; // inside i
    28162782        }
    28172783
    28182784        //   return -1; // outside
    28192785    }
     2786
    28202787    return -1; // outside
    28212788}
     
    29172884
    29182885int produce_phes( FILE *sp, // the input file
    2919                   struct camera *cam, // the camera layout
     2886                  MGeomCamMagic *camgeom, // the camera layout
    29202887                  float minwl_nm, // the minimum accepted wavelength
    29212888                  float maxwl_nm, // the maximum accepted wavelength
     
    29302897                   ){
    29312898
    2932   static int i, k, ipixnum;
     2899  static uint i;
     2900  static int k, ipixnum;
    29332901  static float cx, cy, wl, qe, t;
    29342902  static MCCphoton photon;
    29352903  static float **qept;
    29362904  static char flag[SIZE_OF_FLAGS + 1];
    2937   static float radius;
     2905  static float radius_mm;
    29382906  static float SFR  , C3 , C2 , C1;
    29392907  const double pi = 3.14159;
     
    29412909  // reset variables
    29422910
    2943   for ( i=0; i<cam->inumpixels; ++i ){
     2911  for ( i=0; i<camgeom->GetNumPixels(); ++i ){
    29442912   
    29452913    nphe[i] = 0.0;
     
    29492917  *incph = 0;
    29502918
    2951   radius = cam->dxc[cam->inumpixels-1]
    2952     + 1.5*cam->dpixdiameter_cm*cam->dpixsizefactor[cam->inumpixels-1];
     2919  radius_mm = camgeom->GetMaxRadius();
     2920
    29532921
    29542922  //- - - - - - - - - - - - - - - - - - - - - - - - -
     
    30433011    // check if photon has valid wavelength and is inside outermost camera radius
    30443012         
    3045     if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){
     3013    if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy)*10 > radius_mm ) ){
    30463014      continue;
    30473015           
    30483016    }
    30493017
    3050     ipixnum = bpoint_is_in_pix(cx, cy, cam);
     3018    ipixnum = bpoint_is_in_pix(cx*10, cy*10, camgeom);
    30513019
    30523020    // -1 = the photon is in none of the pixels
     
    31213089
    31223090int produce_nsbrates( char *iname, // the starfield input file name
    3123                       struct camera *cam, // camera layout
     3091                      MGeomCamMagic *camgeom, // camera layout
    31243092                      float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
    31253093                                               // the NSB rates in phe/ns for each pixel
    31263094                      ){
    31273095
    3128   int i, j, k, ii; // counters
     3096  uint i, j; // counters
     3097  int k, ii; // counters
    31293098
    31303099  MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
     
    32123181          // initialize the rate array
    32133182
    3214           for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
     3183          for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
    32153184            rate_phepns[j][i] = 0.;
    32163185          }
     
    32313200
    32323201            k = produce_phes( infile,
    3233                               cam,
     3202                              camgeom,
    32343203                              wl_nm[i],
    32353204                              wl_nm[i+1],
     
    32483217            }
    32493218         
    3250             for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
     3219            for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels
    32513220              rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
    32523221            }
     
    34093378//
    34103379// $Log: not supported by cvs2svn $
     3380// Revision 1.40  2002/07/16 16:15:22  blanch
     3381// A first implementation of the Star field rotation has been done.
     3382//
    34113383// Revision 1.39  2002/04/27 10:48:39  blanch
    34123384// Some unused varibles have been removed.
     
    35283500//
    35293501// $Log: not supported by cvs2svn $
     3502// Revision 1.40  2002/07/16 16:15:22  blanch
     3503// A first implementation of the Star field rotation has been done.
     3504//
    35303505// Revision 1.39  2002/04/27 10:48:39  blanch
    35313506// Some unused varibles have been removed.
Note: See TracChangeset for help on using the changeset viewer.