Changeset 1694 for trunk/MagicSoft/Simulation/Detector

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:

• 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.49 $
+// $Revision: 1.50 $
 // $Author: blanch $ 
-// $Date: 2002-12-13 10:04:07 $
+// $Date: 2003-01-07 16:33:31 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -380 +380 @@
 // switch on  starfield rotation
 static int Starfield_rotate = TRUE;
-static float zenith = 0.0;//zenith angle for starfield rotation
-static float azimutal = 90.0;//azimuth angle for starfield rotation
+
 //=-----------------------------------------------------------
 // @subsection Main program.
@@ -451 +450 @@
   float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS];    //@< non-diffuse nsb 
                                                      //@< photoelectron rates 
-  float nsb_phepns[iMAXNUMPIX];      //@< NSB values for each pixel
+  float nsb_phepns[iMAXNUMPIX];                 //@< NSB values for each pixel
+  float nsb_phepns_rotated[iMAXNUMPIX];         //@< NSB values for each pixel after rotation
 
   Float_t nsb_trigresp[TRIGGER_TIME_SLICES];    //@< array to write the trigger
@@ -509 +509 @@
   UInt_t corsika = 5200 ;
 
+  // Star Field Rotation variables
+  Float_t zenith = 0.0;
+  Float_t azimutal = 90.0;
+  Float_t rho  , C3 , C2 , C1; 
+
   MGeomCamMagic camgeom; // structure holding the camera definition
 
@@ -525 +530 @@
   for(i=0;i<iMAXNUMPIX;i++){
     nsb_phepns[i]=0;
+    nsb_phepns_rotated[i]=0;
     for(j=0;j<iNUMWAVEBANDS;j++)
       nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
@@ -640 +646 @@
   Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
  
-  //parameters for starfield rotation
+  //Parameters for starfield rotation
+
     get_teles_axis(&zenith, &azimutal); 
     Starfield_rotate = get_starfield_rotate();
@@ -1163 +1170 @@
         for(UInt_t ui=0; ui<camgeom.GetNumPixels();ui++){
             nsb_phepns[ui]+=diffnsb_phepns[ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
+            nsb_phepns_rotated[ui]=nsb_phepns[ui];
         }
     }
@@ -1278 +1286 @@
         
         // get MCEventHeader
-        if (reflector_file_version<6)
+        if (reflector_file_version<6){
           fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
-        else
+          thetashw = mcevth.get_theta();
+          phishw = mcevth.get_phi();
+        }
+        else{
           fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile );
-          
+          thetashw = mcevth_2.get_theta();
+          phishw = mcevth_2.get_phi();
+        }
 
         // calculate impact parameter (shortest distance betwee the original
@@ -1412 +1425 @@
                           &ncph,    // will be changed by the function!
                           &arrtmin_ns, // will be changed by the function!
-                          &arrtmax_ns, // will be changed by the function!
-                          -1 ); // photons from the shower
+                          &arrtmax_ns); // will be changed by the function!
+
 
         if( k != 0 ){ // non-zero returnvalue means error
@@ -1424 +1437 @@
         if(simulateNSB && nphe2NSB<=inumphe){
 
+          if(Starfield_rotate){
+            
+            // Introduction rho angle
+            
+            zenith = thetashw;
+            azimutal = phishw;
+            C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
+            C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
+            C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
+            rho = acos( C1/C2 );
+           
+            if ( sin(azimutal) < 0)
+              {
+                rho = 2 * 3.14159 - rho;
+              }   
+            else
+              {
+                rho = rho;
+              }
+            rho = rho*180/3.14159;
+
+            // Rottion of the NSB
+                
+            k = size_rotated(
+                             &nsb_phepns_rotated[0],
+                             nsb_phepns,
+                             rho);
+            
+ 
+          }
+
           //  Fill trigger and fadc response in the trigger class from the database
           for(UInt_t ui=0;ui<camgeom.GetNumPixels();ui++)
-            if(nsb_phepns[ui]>0.0){
+            if(nsb_phepns_rotated[ui]>0.0){
               k=lons.GetResponse(nsb_phepns[ui],0.1,
                                  & nsb_trigresp[0],& nsb_fadcresp[0]);
@@ -3400 +3444 @@
 }
 
+//------------------------------------------------------------
+// @name size_rotated
+//                                     
+// @desc It rotates the NSB
+//
+// @date Tue Jan  7 17:09:25 CET 2003
+// @function @code 
+//------------------------------------------------------------
+
+int size_rotated(
+    float *rotated,
+    float nsb[iMAXNUMPIX],
+    float rho)
+{
+    int r=0; 
+    float size_rotated;
+    float Num_Pixels;
+    float PixNum=0;
+    float rho_pixel;
+    int j=0;
+    int k=0;
+
+    for(int i=1; i<iMAXNUMPIX;i++){
+// Substituir per Int_t radius[iMAXNUMPIX]={0,1,1,1,1,1,1,2,2,2,2, ...}
+        Num_Pixels=0;
+        for (int ii=1; ii<17;ii++){
+            if (Num_Pixels >= i){
+                r=ii-1;
+                break;
+            }
+                
+            
+            
+            if (ii<12){
+                Num_Pixels+=ii*6;
+                PixNum=6*ii;
+                //      size_rotated = (360/PixNum);
+                //rho_pixel = rho/size_rotated;  
+            }   
+            
+            
+            else
+            {
+                Num_Pixels+=(ii-6)*6;
+                PixNum=(ii-6)*6;
+                // size_rotated = (360/PixNum);
+                // rho_pixel = rho/size_rotated;
+            }
+        }            
+            
+            
+            
+        size_rotated = (360/PixNum);
+        rho_pixel = rho/size_rotated;
+        
+        //Buscar j i k que no canvin de r!!!
+        
+        j=i-int(rho_pixel);
+        
+        //cout<<"j inicial "<<j<<endl;
+        int  MINPIXinner[13]= {0,1,7,19,37,61,91,127,169,217,271,331,397};
+        int  MINPIXouter[5]={397,433,475,523,578};
+            
+        if( r<12 ){
+            if (j < MINPIXinner[r])
+            {
+                j = i+6 * r - int(rho_pixel);
+            }
+            
+            k=j-1;
+            if (k < MINPIXinner[r])
+            {
+                k =  MINPIXinner[r+1]-1;
+            }
+        }
+        if(  r > 11 && r < 16){    
+            if (j < MINPIXouter[r-12])
+            {
+                j = i + (r - 6) * 6 - int(rho_pixel);
+                
+            }
+            
+            k=j-1;
+            if ( k < MINPIXouter[r-12])
+            {
+                k = MINPIXouter[r-11] - 1;
+            }
+        }
+        rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k];
+    }
+    return (int(rho_pixel));
+    
+}
+
 
 //------------------------------------------------------------
@@ -3412 +3550 @@
 
 int produce_phes( FILE *sp, // the input file
-                  MGeomCamMagic *camgeom, // the camera layout
+                  class MGeomCamMagic *camgeom, // the camera layout
                   float minwl_nm, // the minimum accepted wavelength
                   float maxwl_nm, // the maximum accepted wavelength
@@ -3421 +3559 @@
                   int *incph,    // total number of cph read
                   float *tmin_ns,   // minimum arrival time of all phes
-                  float *tmax_ns,    // maximum arrival time of all phes
-                  int star          // star=0 starfield rotation star=-1 not rotation
+                  float *tmax_ns    // maximum arrival time of all phes
                    ){
 
@@ -3433 +3570 @@
   static char flag[SIZE_OF_FLAGS + 1];
   static float radius_mm;
-  static float SFR  , C3 , C2 , C1; 
-  const double pi = 3.14159;
 
   // reset variables
@@ -3490 +3625 @@
     //
         
-   // Where the photon provide(0 from starfield, -1 from Mmcs ???
-    
-       if ( star == 0   && Starfield_rotate == TRUE )
-
-       {
-       get_teles_axis(&zenith, &azimutal); 
-         // Introduction SFR angle
-       
-        zenith = (pi/180) * zenith;
-
-        azimutal =( pi/180) * azimutal;
-
-        C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
-   
-        C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
-    
-        C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 );
-
-        SFR = acos( C1/C2 );
-
-     if ( sin(azimutal) < 0)
-            {
-               SFR = 2 * pi - SFR;
-            }   
-         else
-            {
-               SFR = SFR;
-            }
-   
-//final of introduction SFR angle//
-
-       cx = photon.get_x() * cos(SFR) - photon.get_y() * sin(SFR);
-       cy = photon.get_y() * cos(SFR) + photon.get_x() * sin(SFR); 
-
-              }
-        else
-     { 
-         cx = photon.get_x();
-         cy = photon.get_y(); 
-
-     }
- 
+    cx = photon.get_x();
+    cy = photon.get_y(); 
   
     // get wavelength
@@ -3584 +3679 @@
     cu=photon.get_u();
     cv=photon.get_v();
-    cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*pi/180.0;
+    cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*3.14159/180.0;
 
     // find data point in the QE table (-> k)
@@ -3764 +3859 @@
                               &incph,
                               &tmin_ns,
-                              &tmax_ns,
-                              0 ); // photons from starfield
+                              &tmax_ns
+                              ); // photons from starfield
 
             if( k != 0 ){ // non-zero returnvalue means error
@@ -3933 +4028 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.49  2002/12/13 10:04:07  blanch
+// *** empty log message ***
+//
 // Revision 1.48  2002/12/12 17:40:50  blanch
 // *** empty log message ***
@@ -4081 +4179 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.49  2002/12/13 10:04:07  blanch
+// *** empty log message ***
+//
 // Revision 1.48  2002/12/12 17:40:50  blanch
 // *** empty log message ***