Changeset 1694


Ignore:
Timestamp:
01/07/03 16:33:31 (22 years ago)
Author:
blanch
Message:
Star Field Rotation has been implemented by Raul Orduna. Now there is a
rotation for each shower. It is done by a non enter pixel rotation assuming
a circular symetry of the camera. It is not exact but it is accurate enough and
much faster.
File:
1 edited

Legend:

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

    r1685 r1694  
    2121//
    2222// $RCSfile: camera.cxx,v $
    23 // $Revision: 1.49 $
     23// $Revision: 1.50 $
    2424// $Author: blanch $
    25 // $Date: 2002-12-13 10:04:07 $
     25// $Date: 2003-01-07 16:33:31 $
    2626//
    2727////////////////////////////////////////////////////////////////////////
     
    380380// switch on  starfield rotation
    381381static int Starfield_rotate = TRUE;
    382 static float zenith = 0.0;//zenith angle for starfield rotation
    383 static float azimutal = 90.0;//azimuth angle for starfield rotation
     382
    384383//=-----------------------------------------------------------
    385384// @subsection Main program.
     
    451450  float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS];    //@< non-diffuse nsb
    452451                                                     //@< photoelectron rates
    453   float nsb_phepns[iMAXNUMPIX];      //@< NSB values for each pixel
     452  float nsb_phepns[iMAXNUMPIX];                 //@< NSB values for each pixel
     453  float nsb_phepns_rotated[iMAXNUMPIX];         //@< NSB values for each pixel after rotation
    454454
    455455  Float_t nsb_trigresp[TRIGGER_TIME_SLICES];    //@< array to write the trigger
     
    509509  UInt_t corsika = 5200 ;
    510510
     511  // Star Field Rotation variables
     512  Float_t zenith = 0.0;
     513  Float_t azimutal = 90.0;
     514  Float_t rho  , C3 , C2 , C1;
     515
    511516  MGeomCamMagic camgeom; // structure holding the camera definition
    512517
     
    525530  for(i=0;i<iMAXNUMPIX;i++){
    526531    nsb_phepns[i]=0;
     532    nsb_phepns_rotated[i]=0;
    527533    for(j=0;j<iNUMWAVEBANDS;j++)
    528534      nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
     
    640646  Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
    641647 
    642   //parameters for starfield rotation
     648  //Parameters for starfield rotation
     649
    643650    get_teles_axis(&zenith, &azimutal);
    644651    Starfield_rotate = get_starfield_rotate();
     
    11631170        for(UInt_t ui=0; ui<camgeom.GetNumPixels();ui++){
    11641171            nsb_phepns[ui]+=diffnsb_phepns[ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
     1172            nsb_phepns_rotated[ui]=nsb_phepns[ui];
    11651173        }
    11661174    }
     
    12781286       
    12791287        // get MCEventHeader
    1280         if (reflector_file_version<6)
     1288        if (reflector_file_version<6){
    12811289          fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
    1282         else
     1290          thetashw = mcevth.get_theta();
     1291          phishw = mcevth.get_phi();
     1292        }
     1293        else{
    12831294          fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile );
    1284          
     1295          thetashw = mcevth_2.get_theta();
     1296          phishw = mcevth_2.get_phi();
     1297        }
    12851298
    12861299        // calculate impact parameter (shortest distance betwee the original
     
    14121425                          &ncph,    // will be changed by the function!
    14131426                          &arrtmin_ns, // will be changed by the function!
    1414                           &arrtmax_ns, // will be changed by the function!
    1415                           -1 ); // photons from the shower
     1427                          &arrtmax_ns); // will be changed by the function!
     1428
    14161429
    14171430        if( k != 0 ){ // non-zero returnvalue means error
     
    14241437        if(simulateNSB && nphe2NSB<=inumphe){
    14251438
     1439          if(Starfield_rotate){
     1440           
     1441            // Introduction rho angle
     1442           
     1443            zenith = thetashw;
     1444            azimutal = phishw;
     1445            C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
     1446            C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
     1447            C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
     1448            rho = acos( C1/C2 );
     1449           
     1450            if ( sin(azimutal) < 0)
     1451              {
     1452                rho = 2 * 3.14159 - rho;
     1453              }   
     1454            else
     1455              {
     1456                rho = rho;
     1457              }
     1458            rho = rho*180/3.14159;
     1459
     1460            // Rottion of the NSB
     1461               
     1462            k = size_rotated(
     1463                             &nsb_phepns_rotated[0],
     1464                             nsb_phepns,
     1465                             rho);
     1466           
     1467 
     1468          }
     1469
    14261470          //  Fill trigger and fadc response in the trigger class from the database
    14271471          for(UInt_t ui=0;ui<camgeom.GetNumPixels();ui++)
    1428             if(nsb_phepns[ui]>0.0){
     1472            if(nsb_phepns_rotated[ui]>0.0){
    14291473              k=lons.GetResponse(nsb_phepns[ui],0.1,
    14301474                                 & nsb_trigresp[0],& nsb_fadcresp[0]);
     
    34003444}
    34013445
     3446//------------------------------------------------------------
     3447// @name size_rotated
     3448//                                     
     3449// @desc It rotates the NSB
     3450//
     3451// @date Tue Jan  7 17:09:25 CET 2003
     3452// @function @code
     3453//------------------------------------------------------------
     3454
     3455int size_rotated(
     3456    float *rotated,
     3457    float nsb[iMAXNUMPIX],
     3458    float rho)
     3459{
     3460    int r=0;
     3461    float size_rotated;
     3462    float Num_Pixels;
     3463    float PixNum=0;
     3464    float rho_pixel;
     3465    int j=0;
     3466    int k=0;
     3467
     3468    for(int i=1; i<iMAXNUMPIX;i++){
     3469// Substituir per Int_t radius[iMAXNUMPIX]={0,1,1,1,1,1,1,2,2,2,2, ...}
     3470        Num_Pixels=0;
     3471        for (int ii=1; ii<17;ii++){
     3472            if (Num_Pixels >= i){
     3473                r=ii-1;
     3474                break;
     3475            }
     3476               
     3477           
     3478           
     3479            if (ii<12){
     3480                Num_Pixels+=ii*6;
     3481                PixNum=6*ii;
     3482                //      size_rotated = (360/PixNum);
     3483                //rho_pixel = rho/size_rotated; 
     3484            }   
     3485           
     3486           
     3487            else
     3488            {
     3489                Num_Pixels+=(ii-6)*6;
     3490                PixNum=(ii-6)*6;
     3491                // size_rotated = (360/PixNum);
     3492                // rho_pixel = rho/size_rotated;
     3493            }
     3494        }           
     3495           
     3496           
     3497           
     3498        size_rotated = (360/PixNum);
     3499        rho_pixel = rho/size_rotated;
     3500       
     3501        //Buscar j i k que no canvin de r!!!
     3502       
     3503        j=i-int(rho_pixel);
     3504       
     3505        //cout<<"j inicial "<<j<<endl;
     3506        int  MINPIXinner[13]= {0,1,7,19,37,61,91,127,169,217,271,331,397};
     3507        int  MINPIXouter[5]={397,433,475,523,578};
     3508           
     3509        if( r<12 ){
     3510            if (j < MINPIXinner[r])
     3511            {
     3512                j = i+6 * r - int(rho_pixel);
     3513            }
     3514           
     3515            k=j-1;
     3516            if (k < MINPIXinner[r])
     3517            {
     3518                k =  MINPIXinner[r+1]-1;
     3519            }
     3520        }
     3521        if(  r > 11 && r < 16){   
     3522            if (j < MINPIXouter[r-12])
     3523            {
     3524                j = i + (r - 6) * 6 - int(rho_pixel);
     3525               
     3526            }
     3527           
     3528            k=j-1;
     3529            if ( k < MINPIXouter[r-12])
     3530            {
     3531                k = MINPIXouter[r-11] - 1;
     3532            }
     3533        }
     3534        rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k];
     3535    }
     3536    return (int(rho_pixel));
     3537   
     3538}
     3539
    34023540
    34033541//------------------------------------------------------------
     
    34123550
    34133551int produce_phes( FILE *sp, // the input file
    3414                   MGeomCamMagic *camgeom, // the camera layout
     3552                  class MGeomCamMagic *camgeom, // the camera layout
    34153553                  float minwl_nm, // the minimum accepted wavelength
    34163554                  float maxwl_nm, // the maximum accepted wavelength
     
    34213559                  int *incph,    // total number of cph read
    34223560                  float *tmin_ns,   // minimum arrival time of all phes
    3423                   float *tmax_ns,    // maximum arrival time of all phes
    3424                   int star          // star=0 starfield rotation star=-1 not rotation
     3561                  float *tmax_ns    // maximum arrival time of all phes
    34253562                   ){
    34263563
     
    34333570  static char flag[SIZE_OF_FLAGS + 1];
    34343571  static float radius_mm;
    3435   static float SFR  , C3 , C2 , C1;
    3436   const double pi = 3.14159;
    34373572
    34383573  // reset variables
     
    34903625    //
    34913626       
    3492    // Where the photon provide(0 from starfield, -1 from Mmcs ???
    3493    
    3494        if ( star == 0   && Starfield_rotate == TRUE )
    3495 
    3496        {
    3497        get_teles_axis(&zenith, &azimutal);
    3498          // Introduction SFR angle
    3499        
    3500         zenith = (pi/180) * zenith;
    3501 
    3502         azimutal =( pi/180) * azimutal;
    3503 
    3504         C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
    3505    
    3506         C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
    3507    
    3508         C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
    3509 
    3510         SFR = acos( C1/C2 );
    3511 
    3512      if ( sin(azimutal) < 0)
    3513             {
    3514                SFR = 2 * pi - SFR;
    3515             }   
    3516          else
    3517             {
    3518                SFR = SFR;
    3519             }
    3520    
    3521 //final of introduction SFR angle//
    3522 
    3523        cx = photon.get_x() * cos(SFR) - photon.get_y() * sin(SFR);
    3524        cy = photon.get_y() * cos(SFR) + photon.get_x() * sin(SFR);
    3525 
    3526               }
    3527         else
    3528      {
    3529          cx = photon.get_x();
    3530          cy = photon.get_y();
    3531 
    3532      }
    3533  
     3627    cx = photon.get_x();
     3628    cy = photon.get_y();
    35343629 
    35353630    // get wavelength
     
    35843679    cu=photon.get_u();
    35853680    cv=photon.get_v();
    3586     cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*pi/180.0;
     3681    cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*3.14159/180.0;
    35873682
    35883683    // find data point in the QE table (-> k)
     
    37643859                              &incph,
    37653860                              &tmin_ns,
    3766                               &tmax_ns,
    3767                               0 ); // photons from starfield
     3861                              &tmax_ns
     3862                              ); // photons from starfield
    37683863
    37693864            if( k != 0 ){ // non-zero returnvalue means error
     
    39334028//
    39344029// $Log: not supported by cvs2svn $
     4030// Revision 1.49  2002/12/13 10:04:07  blanch
     4031// *** empty log message ***
     4032//
    39354033// Revision 1.48  2002/12/12 17:40:50  blanch
    39364034// *** empty log message ***
     
    40814179//
    40824180// $Log: not supported by cvs2svn $
     4181// Revision 1.49  2002/12/13 10:04:07  blanch
     4182// *** empty log message ***
     4183//
    40834184// Revision 1.48  2002/12/12 17:40:50  blanch
    40844185// *** empty log message ***
Note: See TracChangeset for help on using the changeset viewer.