Ignore:
Timestamp:
02/18/00 17:40:35 (25 years ago)
Author:
petry
Message:
This version includes drastic changes compared to camera.cxx 1.4.
It is not yet finished and not immediately useful because the
trigger simulation is not yet re-implemented. I had to take it
out together with some other stuff in order to tidy the whole
program up. This is not meant as an insult to anyone. I needed
to do this in order to be able to work on it.

This version has been put in the repository in order to be
able to share the further development with others.

If you need something working, wait or take an earlier one.
See file README.
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r344 r365  
    2121//
    2222// $RCSfile: camera.cxx,v $
    23 // $Revision: 1.4 $
     23// $Revision: 1.5 $
    2424// $Author: petry $
    25 // $Date: 2000-01-25 08:36:23 $
     25// $Date: 2000-02-18 17:40:35 $
    2626//
    2727////////////////////////////////////////////////////////////////////////
    … …  
    5454
    5555#include "MDiag.h"
    56 
    57 #include "MTrigger.hxx"
    5856
    5957#include "MRawEvt.h"
    … …  
    8684#undef  __DEBUG__
    8785
    88 // flag for NNT in CT1 camera (default: ON )
    89 #undef  __CT1_NO_NEIGHBOURS__
    90 #define __CT1_NO_NEIGHBOURS__
    91 
    92 // flag for calculation of NSB (default: ON )
    93 #undef  __NSB__
    94 #define __NSB__
    95 
    96 // flag for calculation of QE for pixels (default: ON )
    97 #undef  __QE__
    98 #define __QE__
    99 
    100 
    101 // flag for implementation of DETAIL_TRIGGER (default: ON )
    102 //
    103 //      This is the new implementation of Trigger studies
    104 //      It relies on a better simulation of the time stucture
    105 //      of the PhotoMultiplier. For more details see the
    106 //      documentation of the --> class MTrigger <--
    107 #define __DETAIL_TRIGGER__
    108 #undef  __DETAIL_TRIGGER__
    109 
    110 // flag for implementation of TRIGGER (default: ON )
    111 #undef  __TRIGGER__
    112 #define __TRIGGER__
    113 
    114 // flag for implementation of Tail-Cut (default: ON )
    115 #undef  __TAILCUT__
    116 #define __TAILCUT__
    117 
    118 // flag for calculation of islands stat. (default: ON )
    119 #define __ISLANDS__
    120 #undef  __ISLANDS__
    121 
    122 // flag for calculation of image parameters (default: ON )
    123 #undef  __MOMENTS__
    124 #define __MOMENTS__
    125 
    126 // flag for ROOT  (default: ON )
    127 #undef  __ROOT__
    128 #define __ROOT__
    129 
    13086//!@}
    13187
    … …  
    259215// and data stored on them with information about the pixels
    260216
    261 //@: table for IJ sys.
    262 static float pixels[PIX_ARRAY_SIDE][PIX_ARRAY_SIDE][4];   
    263217
    264218//@: coordinates x,y for each pixel
    … …  
    276230//@: contents of the pixels (ph.e.) after cleanning
    277231static float *fnpixclean;
    278 
    279 //@: contents of the sum of all ph.e. in one timeslice of 1 ns
    280 static float *fnslicesum ;
    281232
    282233 
    … …  
    364315  @"*/
    365316
    366 //!@{
    367 // table of random numbers
    368 
    369 // (unused)
    370 //static double RandomNumbers[500]; 
    371 //!@}
    372 
    373 /*!@"
    374 
    375   The following is a variable to count the number of Cphotons
    376   in the different steps of the simulation.
    377   The definition is as follows:
    378   @[
    379   \mbox{CountCphotons}[ \mbox{FILTER} ] \equiv
    380   \mbox{\it Number of photons after the filter} \mbox{FILTER}
    381   @]
    382   The filters are defined and can be found in the file |camera.h|.
    383 
    384   @"*/
    385 
    386 //!@{
    387 // vector to count photons at any given step of the simulation
    388 
    389 static int CountCphotons[10]; 
    390 //!@}
    391 
    392317/*!@"
    393318
    … …  
    430355
    431356  char inname[256];           //@< input file name
    432   char outname[256];          //@< output file name
     357  char starfieldname[256];    //@< starfield input file name
    433358  char datname[256];          //@< data (ASCII) output file name
    434359  char diagname[256];         //@< diagnistic output file (ROOT format)
    … …  
    437362  char parname[256];          //@< parameters file name
    438363
    439   char sign[20];              //@< initialize sign
    440 
    441364  char flag[SIZE_OF_FLAGS + 1];  //@< flags in the .rfl file
    442365
    443   ifstream inputfile;         //@< stream for the input file
    444   ofstream outputfile;        //@< stream for the output file
     366  FILE *inputfile;            //@< stream for the input file
    445367  ofstream datafile;          //@< stream for the data file
    446368
    447369  MCEventHeader mcevth;       //@< Event Header class (MC)
    448   MCCphoton cphoton;          //@< Cherenkov Photon class (MC)
     370
     371  Photoelectron *photoe = NULL; //@< array of the photoelectrons of one event
     372  int inumphe;                //@< number of photoelectrons in an event
     373
     374  float arrtmin_ns;           //@ arrival time of the first photoelectron
     375  float arrtmax_ns;           //@ arrival time of the last photoelectron
    449376
    450377  float thetaCT, phiCT;       //@< parameters of a given shower
    … …  
    458385  int nshow=0;                //@< partial number of shower in a given run
    459386  int ntshow=0;               //@< total number of showers
    460   int ncph=0;                 //@< partial number of shower in a given run
    461   int ntcph=0;                //@< total number of showers
    462 
    463   int i, ii, j, k;            //@< simple counters
    464 
    465   float t_ini;                //@< time of the first Cphoton read in
    466   float t;                    //@< time for a single photon
    467   int   t_chan ;              //@< the bin (channel) in time of a single photon
    468  
    469   int   startchan ;           //@< the first bin with entries in the time slices
    470 
    471   float cx, cy, ci, cj;       //@< coordinates in the XY and IJ systems
    472   int ici, icj, iici, iicj;   //@< coordinates in the IJ (integers)
    473 
    474   int nPMT;                   //@< number of pixel
    475 
    476   float wl, last_wl;          //@< wavelength of the photon
    477   float qe;                   //@< quantum efficiency
    478   float **qeptr;
     387  int ncph=0;                 //@< partial number of photons in a given run
     388  int ntcph=0;                //@< total number of photons
     389
     390  int i, j, k;                //@< simple counters
    479391
    480392  int simulateNSB;            //@< Will we simulate NSB?
    481   float meanNSB;              //@< NSB mean value (per pixel)
     393  float meanNSB;              //@< diffuse NSB mean value (phe per ns per central pixel)
     394  float diffnsb_phepns[iMAXNUMPIX];  //@< diffuse NSB values for each pixel derived
     395                                     //@< from meanNSB
     396  float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS];    //@< non-diffuse nsb
     397                                                     //@< photoelectron rates
     398  float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
     399    EXTWAVEBAND1,
     400    EXTWAVEBAND2,
     401    EXTWAVEBAND3,
     402    EXTWAVEBAND4,
     403    EXTWAVEBAND5
     404  };                         
     405  float baseline_mv[iMAXNUMPIX]; //@< The baseline (mV) caused by the NSB; to be subtracted
     406                                 //@< in order to simulate the preamps' AC coupling
     407
    482408  float qThreshold;           //@< Threshold value
    483409  float qTailCut;             //@< Tail Cut value
    … …  
    488414  float fCorrection;          //@< Factor to apply to pixel values (def. 1.)
    489415
    490   float q0;                   //@< trigger threshold ( intermediate variable )
    491   float maxcharge;            //@< maximum charge in pixels
    492   int noverq0, novq0;         //@< number of pixels above threshold
    493   int ngrpq0, mxgrp;          //@< number of pixels in a group
    494 
    495416  int trigger;                //@< trigger flag
    496417  int itrigger;               //@< index of pixel fired
    497418  int ntrigger = 0;           //@< number of triggers in the whole file
    498   unsigned char triggerBits;  //@< byte for trigger condition check (MAGIC)
    499   int bit;                    //@< intermediate variable
    500419
    501420  float plateScale_cm2deg;    //@< plate scale (deg/cm)
    … …  
    506425  int still_in_loop = FALSE;
    507426
    508   char    Signature[20];
    509 
    510427  float *image_data;
    511428  int nvar, hidt;
    … …  
    513430  struct camera cam; // structure holding the camera definition
    514431
    515 #ifdef __DETAIL_TRIGGER__
    516 
    517   MTrigger  Trigger ;         //@< A instance of the Class MTrigger
    518 
    519 #endif // __DETAIL_TRIGGER__
    520432
    521433  //!@' @#### Definition of variables for |getopt()|.
    … …  
    592504
    593505  strcpy( inname, get_input_filename() );
    594   strcpy( outname, get_output_filename() );
     506  strcpy( starfieldname, get_starfield_filename() );
    595507  strcpy( datname, get_data_filename() );
    596508  strcpy( diagname, get_diag_filename() );
    … …  
    615527      "Filenames",
    616528      "In", inname,
    617       "Out", outname,
     529      "Stars", starfieldname,
     530      "CT", ct_filename,
    618531      "Data", datname,
    619532      "Diag", diagname,
    620       "ROOT",  rootname,
    621       "CT", ct_filename);
     533      "ROOT",  rootname
     534      );
    622535
    623536 
    … …  
    669582  read_ct_file();
    670583
    671   // read pixels data
     584  // read camera setup
    672585
    673586  read_pixels(&cam);
    674587
     588  // allocate memory for the photoelectrons
     589
     590  photoe = new Photoelectron[iMAXNUMPHE];
     591
    675592  // initialise ROOT
    676593
    … …  
    682599
    683600  hfile = new TFile( diagname,"RECREATE", "MAGIC Telescope MC diagnostic data");
    684  
    685601 
    686602  // Create the ROOT Tree for the diagnostic data
    … …  
    698614  tree->Branch("event", "MDiagEventobject", &event, bsize, split);
    699615
    700 #ifdef __ROOT__
     616  // Prepare the raw data output
    701617
    702618  MRawEvt *Evt   = new MRawEvt() ;
    … …  
    706622  //
    707623
    708   TFile outfile ( rootname , "RECREATE" ) ;
    709 
    710   //
     624  TFile outfile ( rootname , "RECREATE" );
     625
    711626  //      create a Tree for the Event data stream
    712   //
    713627
    714628  TTree EvtTree("EvtTree","Events of Run");
    … …  
    722636                 &McEvt, bsize, split);
    723637
    724 
    725   float  flli = 0. ;
    726   unsigned short ulli = 0 ;
    727 
    728 #endif // __ROOT__
    729638
    730639  // for safety and for dimensioning image_data: count the elements in the
    … …  
    787696  anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
    788697
    789   // open input file if we DO read data from a file
    790 
    791   if (! Data_From_STDIN) { 
    792     log( SIGNATURE, "Openning input \"rfl\" file %s\n", inname );
    793     inputfile.open( inname );
    794     if ( inputfile.bad() )
     698  // prepare the NSB simulation
     699
     700  if( simulateNSB ){
     701
     702    //
     703    // Calculate the non-diffuse NSB photoelectron rates
     704    //
     705   
     706    k = produce_nsbrates( starfieldname,
     707                          &cam,
     708                          photoe, // only a dummy here
     709                          nsbrate_phepns );
     710    if (k != 0){
     711      cout << "Error when reading starfield... \nExiting.\n";
     712      exit(1);
     713    }
     714 
     715//     for(i=0; i<cam.inumpixels; i++){
     716//       cout << i;
     717//       for(j=0; j<iNUMWAVEBANDS; j++){
     718//         cout << " " << j << " " << nsbrate_phepns[i][j];
     719//       }
     720//       cout << "\n";
     721//     }
     722
     723    // calculate diffuse rate correcting for the pixel size
     724
     725    for(i=0; i<cam.inumpixels; i++){
     726      diffnsb_phepns[i] = meanNSB *
     727        cam.dpixsizefactor[i] * cam.dpixsizefactor[i];
     728    }
     729
     730  }
     731
     732  //
     733  // Read the reflector file with the Cherenkov data
     734  //                   
     735
     736  // select input file
     737
     738  if ( Data_From_STDIN ) {
     739
     740    inputfile = stdin;
     741
     742  }
     743  else{
     744
     745    log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
     746    inputfile = fopen( inname, "r" );
     747    if ( inputfile == NULL )
    795748      error( SIGNATURE, "Cannot open input file: %s\n", inname );
     749
    796750  }
    797751 
    798752  // get signature, and check it
    799   // NOTE: this part repeats further down in the code;
    800   // if you change something here you probably want to change it
    801   // there as well
    802 
    803   strcpy(Signature, REFL_SIGNATURE);
    804 
    805   strcpy(sign, Signature);
    806 
    807   if ( Data_From_STDIN )
    808     cin.read( (char *)sign, strlen(Signature));
    809   else
    810     inputfile.read( (char *)sign, strlen(Signature));
    811 
    812   if (strcmp(sign, Signature) != 0) {
    813     cerr << "ERROR: Signature of .rfl file is not correct\n";
    814     cerr << '"' << sign << '"' << '\n';
    815     cerr << "should be: " << Signature << '\n';
     753
     754  if(check_reflector_file( inputfile )==FALSE){
    816755    exit(1);
    817756  }
    818757
    819   if ( Data_From_STDIN )
    820     cin.read( (char *)sign, 1);
    821   else
    822     inputfile.read( (char *)sign, 1);
    823 
    824   // open output file
    825  
    826   log( SIGNATURE, "Openning output \"phe\" file %s\n", outname );
    827   outputfile.open( outname );
    828  
    829   if ( outputfile.bad() )
    830     error( SIGNATURE, "Cannot open output file: %s\n", outname );
    831  
    832758  // open data file
    833759 
    834   log( SIGNATURE, "Openning data \"dat\" file %s\n", datname );
     760  log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
    835761  datafile.open( datname );
    836762 
    837   if ( outputfile.bad() )
    838     error( SIGNATURE, "Cannot open output file: %s\n", outname );
    839  
    840   // write signature
    841 
    842   outputfile.write( SIGNATURE, sizeof(SIGNATURE) );
    843 
     763  if ( datafile.bad() )
     764    error( SIGNATURE, "Cannot open data file: %s\n", datname );
     765 
    844766  // initializes flag
    845767 
    … …  
    851773  fnpixclean = new float [ ct_NPixels ];
    852774
    853 #ifdef __ROOT__
    854  
    855   fnslicesum = new float [ (2 * SLICES) ] ;
    856  
    857   float slices   [CAMERA_PIXELS][ (2 * SLICES) ] ;
    858   float slices2  [CAMERA_PIXELS][ SLICES ] ;
    859 
    860   float trans    [ SLICES ] ;
    861 #endif // __ROOT__
    862 
    863  
    864775  moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
     776
     777  // initialize baseline
     778
     779  for(i=0; i<cam.inumpixels; i++){
     780    baseline_mv[i] = 0.;
     781  }
    865782       
    866783  //!@' @#### Main loop.
    867784  //@'
    868785
    869   //begin my version
    870                                            
    871786  // get flag
    872787   
    873   if ( Data_From_STDIN )
    874     cin.read( flag, SIZE_OF_FLAGS );
    875   else
    876     inputfile.read ( flag, SIZE_OF_FLAGS );
     788  fread( flag, SIZE_OF_FLAGS, 1, inputfile );
    877789
    878790  // loop over the file
    … …  
    881793
    882794  while (
    883          ((! Data_From_STDIN) && (! inputfile.eof()))
     795         ((! Data_From_STDIN) && ( !feof(inputfile) ))
    884796         ||
    885797         (Data_From_STDIN && still_in_loop)
    … …  
    891803    }
    892804    else { // found start of run
     805
    893806      nshow=0;
    894       // read next flag
    895 
    896       if ( Data_From_STDIN )
    897         cin.read( flag, SIZE_OF_FLAGS );
    898       else
    899         inputfile.read ( flag, SIZE_OF_FLAGS );
     807
     808      fread( flag, SIZE_OF_FLAGS, 1, inputfile );
    900809
    901810      while( isA( flag, FLAG_START_OF_EVENT   )){ // while there is a next event
    902         /*!@'
    903          
    904           For the case  |FLAG\_START\_OF\_EVENT|,
    905           we read each Cherenkov photon, and follow these steps:
    906          
    907           @enumerate
    908          
    909           @- Transform XY-coordinates to IJ-coordinates.
    910          
    911           @- With this, we obtain the pixel where the photon hits.
    912          
    913           @- Use the wavelength $\lambda$ and the table of QE, and
    914           calculate the estimated (third order interpolated) quantum
    915           efficiency for that photon. The photon can be rejected.
    916          
    917           @- If accepted, then add to the pixel.
    918          
    919           @endenumerate
    920          
    921           In principle, we should stop here, and use another program to
    922           'smear' the image, to add the Night Sky Background, and to
    923           simulate the trigger logic, but we will make this program
    924           quick and dirty, and include all here.
    925          
    926           If we are reading PHE files, we jump to the point where the
    927           pixelization process already has finished.
    928          
    929         */
    930811       
    931812        ++nshow;
    … …  
    934815        // get MCEventHeader
    935816       
    936         if ( Data_From_STDIN )
    937           cin.read( (char*)&mcevth, mcevth.mysize() );
    938         else
    939           mcevth.read( inputfile );
     817        fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
    940818       
    941819        // calculate core distance and impact parameter
    … …  
    1023901        }
    1024902
    1025         // clear camera
    1026        
    1027         for ( i=0; i<ct_NPixels; ++i ){
    1028  
    1029           fnpix[i] = 0.0;
    1030 #ifdef __ROOT__
    1031           for ( ii=0; ii<(2 * SLICES); ii++ ) {
    1032             slices [i][ii] = 0 ;
    1033           }
    1034 #endif // __ROOT__
     903        // read the photons and produce the photoelectrons
     904
     905        k = produce_phes( inputfile,
     906                          &cam,
     907                          WAVEBANDBOUND1,
     908                          WAVEBANDBOUND6,
     909                          photoe,   // will be changed by the function!
     910                          &inumphe, // important for later: the size of photoe[]
     911                          fnpix,    // will be changed by the function!
     912                          &ncph,    // will be changed by the function!
     913                          &arrtmin_ns, // will be changed by the function!
     914                          &arrtmax_ns // will be changed by the function!
     915                          );
     916
     917        if( k != 0 ){ // non-zero returnvalue means error
     918          cout << "Exiting.\n";
     919          exit(1);
    1035920        }
    1036 
    1037         ntcph +=ncph;
    1038         ncph = 0;
    1039 
    1040 #ifdef __DETAIL_TRIGGER__
    1041         //
    1042         //   clear Trigger
    1043         //
    1044      
    1045         Trigger.Reset() ;
    1046 #endif // __DETAIL_TRIGGER__
    1047        
    1048         //- - - - - - - - - - - - - - - - - - - - - - - - -
    1049         // read photons and "map" them into the pixels
    1050         //--------------------------------------------------     
    1051        
    1052         // initialize CPhoton
    1053        
    1054         cphoton.fill(0., 0., 0., 0., 0., 0., 0., 0.);
    1055        
    1056         // read the photons data
    1057        
    1058         if ( Data_From_STDIN )
    1059           cin.read( flag, SIZE_OF_FLAGS );
    1060         else
    1061           inputfile.read ( flag, SIZE_OF_FLAGS );
    1062          
    1063         // loop over the photons
    1064 
    1065         t_ini = -99999;
    1066        
    1067         while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
    1068          
    1069           memcpy( (char*)&cphoton, flag, SIZE_OF_FLAGS );
    1070 
    1071           if ( Data_From_STDIN )
    1072             cin.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
    1073           else
    1074             inputfile.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
    1075          
    1076           // increase number of photons
    1077          
    1078           ncph++;
    1079          
    1080           t = cphoton.get_t() ;
    1081          
    1082           if(t_ini == -99999){ // this is the first photon we read from this event
    1083             t_ini = t; // memorize time
    1084           }
    1085          
    1086           // The photons don't come in chronological order!
    1087           // Put the first photon at the center of the array by adding the constant SLICES
    1088          
    1089           t_chan = (int) ((t - t_ini )/ WIDTH_TIMESLICE ) + SLICES ;     
    1090          
    1091           if (t_chan > (2 * SLICES)){
    1092             log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n",
    1093                 t_chan, (2 * SLICES), (2 * SLICES));
    1094             t_chan = (2 * SLICES);
    1095           }
    1096           else if(t_chan < 0){
    1097             log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n",
    1098                 t_chan, 0, 0);
    1099             t_chan = 0;
    1100           }
    1101            
    1102           // Pixelization
    1103            
    1104           cx = cphoton.get_x();
    1105           cy = cphoton.get_y();
    1106          
    1107           // get wavelength
    1108          
    1109           last_wl = wl;
    1110           wl = cphoton.get_wl();
    1111          
    1112           // check if photon has valid wavelength and is inside outermost camera radius
    1113          
    1114           if( (wl > 800.0) || (wl < 290.0) ||
    1115               (sqrt(cx*cx + cy*cy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)) ){
    1116            
    1117             // read next CPhoton
    1118             if ( Data_From_STDIN )
    1119               cin.read( flag, SIZE_OF_FLAGS );
    1120             else
    1121               inputfile.read ( flag, SIZE_OF_FLAGS );
    1122            
    1123             // go to beginning of loop, the photon is lost
    1124             continue;
    1125            
    1126           }
    1127 
    1128           // cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
    1129          
    1130           nPMT = -1;
    1131 
    1132           for(i=0; i<ct_NPixels; i++){
    1133             if( bpoint_is_in_pix( cx, cy, i, &cam) ){
    1134               nPMT = i;
    1135               break;
    1136             }
    1137           }
    1138            
    1139           if(nPMT==-1){// the photon is in none of the pixels
    1140 
    1141             // read next CPhoton
    1142             if ( Data_From_STDIN )
    1143               cin.read( flag, SIZE_OF_FLAGS );
    1144             else
    1145               inputfile.read ( flag, SIZE_OF_FLAGS );
    1146            
    1147             // go to beginning of loop, the photon is lost
    1148             continue;
    1149           }
    1150          
    1151        
    1152          
    1153 #ifdef __QE__
    1154          
    1155           //!@' @#### QE simulation.
    1156           //@'
    1157          
    1158           //+++
    1159           // QE simulation
    1160           //---
    1161          
    1162           // find data point to be used in Lagrange interpolation (-> k)
    1163          
    1164           qeptr = (float **)QE[nPMT];
    1165          
    1166           FindLagrange(qeptr,k,wl);
    1167          
    1168           // if random > quantum efficiency, reject it
    1169          
    1170           qe = Lagrange(qeptr,k,wl) / 100.0;
    1171          
    1172           // fprintf(stdout, "%f\n", qe);
    1173          
    1174           if ( RandomNumber > qe ) {
    1175            
    1176             // read next CPhoton
    1177             if ( Data_From_STDIN )
    1178               cin.read( flag, SIZE_OF_FLAGS );
    1179             else
    1180               inputfile.read ( flag, SIZE_OF_FLAGS );
    1181            
    1182             // go to beginning of loop
    1183             continue;
    1184            
    1185           }
    1186          
    1187 #endif // __QE__
    1188          
    1189           //+++
    1190           // Cphoton is accepted
    1191           //---
    1192          
    1193           // increase the number of Cphs. in the PMT, i.e.,
    1194           // increase in one unit the counter of the photons
    1195           // stored in the pixel nPMT
    1196          
    1197           fnpix[nPMT] += 1.0;
    1198          
    1199 #ifdef __ROOT__
    1200           fnslicesum[t_chan]  += 1.0 ;
    1201           slices[nPMT][t_chan] += 1.0 ;
    1202 #endif // __ROOT__       
    1203          
    1204 #ifdef __DETAIL_TRIGGER__
    1205           //
    1206           //  fill the Trigger class with this phe
    1207           //
    1208           //
    1209          
    1210           Trigger.Fill( nPMT, ( t - t_ini  ) ) ;
    1211 #endif // __DETAIL_TRIGGER__
    1212          
    1213           // read next CPhoton
    1214          
    1215           if ( Data_From_STDIN )
    1216             cin.read( flag, SIZE_OF_FLAGS );
    1217           else
    1218             inputfile.read ( flag, SIZE_OF_FLAGS );
    1219          
    1220         }  // end while, i.e. found end of event
    1221        
     921         
    1222922        log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
    1223923            ncph, ntcph);
    1224      
    1225         // show number of photons
    1226        
    1227         //cout << ncph << " photons read . . . " << endl << flush;
    1228      
     924
     925        ntcph += ncph;
     926
    1229927        // skip it ?
    1230928       
    … …  
    1245943        }
    1246944       
    1247         /*!@'
    1248          
    1249           After reading all the Cherenkov photons for a given event,
    1250           we have in the table of number of photons for each pixel
    1251           only the 'raw' amount of Cherenkov photons @$n_p@$. Now, we
    1252           should take this number as the mean value of the
    1253           distribution of photons in that pixel @$p@$, following a
    1254           Poisson distribution.
    1255          
    1256           @[ n_p \equiv \mu_p @]
    1257          
    1258           and with this number the amount of light coming from the
    1259           shower is calculated @$\hat{n}_p@$.
    1260          
    1261           Then, we calculate the amount of Night Sky Background we
    1262           must introduce in that pixel @$p@$. We calculate this using
    1263           again a Poisson distribution with mean @$\mu_\mathrm{NSB}@$
    1264           (defined in the |camera.h| file). The value of
    1265           @$\mu_\mathrm{NSB}@$ is obtained from measurements. With this
    1266           value, the amount of photons @$\hat{n}_\mathrm{NSB}@$ coming
    1267           from the Night Sky Background is calculated.
    1268          
    1269           Finally, the amount of photons for that pixels is:
    1270           @[ \hat{n}_p^\mathrm{final} = \hat{n}_p + \hat{n}_\mathrm{NSB} @]
    1271          
    1272         */
    1273        
    1274         // after reading all the photons, our camera is filled
     945        // energy cut
    1275946       
    1276947        if ( Select_Energy ) {
    … …  
    1281952          }
    1282953        }
    1283        
    1284 #ifdef __NSB__
    1285        
    1286         //!@' @#### NSB (Night Sky Background) simulation.
    1287         //@'
    1288        
    1289         //+++
     954
    1290955        // NSB simulation
    1291         //---
    1292        
    1293         // add NSB "noise"
    1294         // TO DO: make meanNSB an array and read the contents from a file!
    1295        
    1296         if ( simulateNSB )
    1297           for ( i=0; i<ct_NPixels; ++i )
    1298             fnpix[i] += (float)ignpoi( meanNSB );
    1299        
    1300 #endif // __NSB__
    1301        
    1302         // if we should apply any kind of correction, do it here.
    1303        
    1304         for ( i=0; i<ct_NPixels; ++i )
    1305           fnpix[i] *= fCorrection;
    1306        
    1307 #ifdef __DETAIL_TRIGGER__
    1308         //   Trigger.Print() ;
    1309         cout << Trigger.Diskriminate() << endl << endl ;
    1310 #endif // __DETAIL_TRIGGER__
    1311        
    1312 #ifdef __ROOT__
    1313        
    1314         //
    1315         //  Fill the header of this event
    1316         //
    1317        
    1318         Evt->FillHeader ( (UShort_t) (ntshow + nshow) ,  20 ) ;
    1319        
    1320         //  now put out the data of interest
    1321         //
    1322         //  1.  -> look for the first slice with signal
    1323         //
    1324        
    1325         for ( i=0; i<(2 * SLICES); i++ )
    1326           if ( fnslicesum[i] > 0. )
    1327             break ;
    1328        
    1329         startchan = i ;
    1330        
    1331         //
    1332         //     copy the slices out of the big array
    1333         //   
    1334         //     put the first slice with signal to slice 4
    1335         //
    1336        
    1337         for (i=0; i<ct_NPixels; i++ )
    1338           for ( ii=(startchan-3); ii < (startchan+12); ii++ ) 
    1339             slices2 [i][ii-startchan+3] = slices [i][ii] ;
    1340        
    1341        
    1342         //
    1343         //  if a pixes has a signal put it to the MRawEvt
    1344         //
    1345        
    1346         for (i=0; i<ct_NPixels; i++ ) {
    1347           if ( fnpix[i] > 0 ) {
    1348            
    1349             for ( ii=0; ii < 15; ii++ ) {
    1350               trans [ii] = slices2[i][ii] ;
    1351             }
    1352            
    1353             Evt->FillPixel ( (UShort_t) i , trans ) ;
    1354            
     956
     957        if(simulateNSB){
     958
     959          k = produce_nsb_phes( &arrtmin_ns, // will be changed by the function!
     960                                &arrtmax_ns, // will be changed by the function!
     961                                thetaCT,
     962                                &cam,
     963                                nsbrate_phepns,
     964                                diffnsb_phepns,
     965                                ext,
     966                                fnpix,   // will be changed by the function!
     967                                photoe,  // will be changed by the function!
     968                                &inumphe, // important for later: the size of photoe[]
     969                                baseline_mv // will be generated by the function
     970                                );
     971
     972          if( k != 0 ){ // non-zero returnvalue means error
     973            cout << "Exiting.\n";
     974            exit(1);
    1355975          }
    1356         }
    1357        
    1358         //
    1359         //   
    1360         //
    1361        
    1362         McEvt->Fill( (UShort_t) mcevth.get_primary() ,
    1363                      mcevth.get_energy(),
    1364                      mcevth.get_theta(),
    1365                      mcevth.get_phi(),
    1366                      mcevth.get_core(),
    1367                      mcevth.get_coreX(),
    1368                      mcevth.get_coreY(),
    1369                      flli,
    1370                      ulli, ulli, ulli, ulli, ulli ) ;
    1371        
    1372         //
    1373         //    write it out to the file outfile
    1374         //
    1375        
    1376         EvtTree.Fill() ;
    1377        
    1378         //    clear all
    1379        
    1380         Evt->Clear() ;
    1381         McEvt->Clear() ;
    1382        
    1383 #endif // __ROOT__
    1384        
    1385         //++++++++++++++++++++++++++++++++++++++++++++++++++
    1386         // at this point we have a camera full of
    1387         // ph.e.s
    1388         // we should first apply the trigger condition,
    1389         // and if there's trigger, then clean the image,
    1390         // calculate the islands statistics and the
    1391         // other parameters of the image (Hillas' parameters
    1392         // and so on).
    1393         //--------------------------------------------------
    1394        
    1395 #ifdef __TRIGGER__
    1396        
    1397         /*!@'
    1398          
    1399           @#### Trigger logic simulation.
    1400          
    1401           In the following block we look at the pixel contents, looking
    1402           for pixels fulfilling the trigger condition. This condition,
    1403           in this current version of the program, is the following:
    1404          
    1405           @itemize
    1406          
    1407           @- |CT1|: Two neighbour pixels with charge above the threshold
    1408           @$q_0@$. For the old CT1 data, however, the trigger condition
    1409           was 'any two pixels with charge above the threshold @$q_0@$'.
    1410          
    1411           @- |MAGIC|: A 'closed-packet' of four neighbour pixels, each
    1412           of them with charge above the threshold @$q_0@$.
    1413          
    1414           @enditemize
    1415          
    1416           In the following figure you can find a sort of description
    1417           about the meanning of 'closed-packet'.
    1418          
    1419           @F
    1420          
    1421           \begin{figure}[htbp]
    1422           \begin{center}
    1423           \includegraphics{closepck.eps}
    1424           \caption{Meanning of the expression ``{\it close-packet}''}
    1425           \label{fig:closepacket}
    1426           \end{center}
    1427           \end{figure}
    1428          
    1429           @F
    1430          
    1431         */
    1432        
    1433         //++
    1434         // TRIGGER Condition
    1435         //--
    1436        
    1437         //@ If the input parameter "threshold" is 0 we find the maximum
    1438         //@ trigger threshold this event can pass
    1439        
    1440         for(k=0; ( qThreshold == 0 ? (k <= iMAX_THRESHOLD_PHE) : (k == 1) ); k++){
    1441          
    1442           // is there trigger?
    1443          
    1444           noverq0 = 0;
    1445           q0 = ( qThreshold == 0. ? (float) k : qThreshold );
    1446           trigger = FALSE;
    1447           mxgrp = 0;
    1448           maxcharge = 0.0;
    1449          
    1450           // Warning! NOT all the camera is able to give trigger
    1451           // only up to 'degTrigger' degrees
    1452          
    1453           for ( i=0 ; (i<ct_NCentralPixels) && (trigger==FALSE) ; ++i ) {
    1454            
    1455             // calculate absolute maximum
    1456            
    1457             maxcharge = MAX(fnpix[i],maxcharge);
    1458            
    1459             // is this pixel above threshold ?
    1460            
    1461             if ( fnpix[i] <= q0 )   
    1462               continue;
    1463            
    1464             // it is: increment the number of pixels above threshold
    1465            
    1466             ++noverq0;
    1467            
    1468             // if the trigger already fired, just count the pixels
    1469             // above threshold
    1470            
    1471             if ( trigger == TRUE )   
    1472               continue;
    1473            
    1474             // is this pixel inside the trigger zone in the camera ?
    1475            
    1476             if ( (sqrt(SQR(pixary[i][0]) +
    1477                        SQR(pixary[i][1]))*plateScale_cm2deg) > degTriggerZone)
    1478               continue;
    1479            
    1480             // 'ngrpq0' is the number of neighbours of pixel i with q > q0
    1481            
    1482             ngrpq0 = 0;
    1483            
    1484             // look at each pixel in the neighborhood, and count
    1485             // those above threshold q0
    1486            
    1487             for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1 ; ++j )
    1488               if ( fnpix[pixneig[i][j]] > q0 )
    1489                 ++ngrpq0;
    1490            
    1491             // check whether we have trigger
    1492            
    1493             if ( ct_Type == 0 ) {
    1494              
    1495               //++ >>>>> CT1 <<<<<
    1496              
    1497 #ifdef __CT1_NO_NEIGHBOURS__
    1498              
    1499               if ( noverq0 > 1 )
    1500                 trigger = TRUE;
    1501              
    1502 #else
    1503              
    1504               if ( ngrpq0 > 0 )
    1505                 trigger = TRUE;
    1506              
    1507 #endif
    1508              
    1509               //-- >>>>> CT1 <<<<<
    1510              
    1511             } else {
    1512              
    1513               //++ >>>>> MAGIC <<<<<
    1514              
    1515               // (at least 4 packed with at least q0 phes)
    1516              
    1517               // there are 3 cases
    1518               // 1. zero, one or two neighbours have enough charge: no trigger
    1519               // 2. five or six neighbours with enough charge: trigger? sure!!
    1520               // 3. three or four neighbours with q > q0 : we must look
    1521               //    for 'closeness'.
    1522              
    1523               switch ( ngrpq0 ) {
    1524                
    1525               case 0:
    1526               case 1:
    1527               case 2:
    1528                
    1529                 trigger = FALSE;
    1530                 break;
    1531                
    1532               case 3:
    1533               case 4:
    1534                
    1535                 // if reaches this line, it means three or four neighbours
    1536                 // around the central pixel
    1537                
    1538                 triggerBits = 1;
    1539                
    1540                 for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) {
    1541                
    1542                   if ( fnpix[pixneig[i][j]] > q0 ) {
    1543                    
    1544                     if ( pixary[pixneig[i][j]][0] > pixary[i][0] ) {
    1545                      
    1546                       if ( nint(pixary[pixneig[i][j]][1]*10.0) >
    1547                            nint(pixary[i][1]*10.0) )
    1548                         bit = 2;
    1549                       else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
    1550                                 nint(pixary[i][1]*10.0) )
    1551                         bit = 6;
    1552                       else
    1553                         bit = 1;
    1554                      
    1555                     } else {
    1556                      
    1557                       if ( nint(pixary[pixneig[i][j]][1]*10.0) >
    1558                            nint(pixary[i][1]*10.0) )
    1559                         bit = 3;
    1560                       else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
    1561                                 nint(pixary[i][1]*10.0) )
    1562                         bit = 5;
    1563                       else
    1564                         bit = 4;
    1565                      
    1566                     }
    1567                    
    1568                     triggerBits |= (1<<bit);
    1569                    
    1570                   }
    1571                  
    1572                 }
    1573                
    1574                 if ( ngrpq0 == 3 ) {  // 4-fold trigger
    1575                  
    1576                   switch ( triggerBits ) {
    1577                    
    1578                   case 0x0f:          // 0 000111 1
    1579                   case 0x1d:          // 0 001110 1
    1580                   case 0x39:          // 0 011100 1
    1581                   case 0x71:          // 0 111000 1
    1582                   case 0x63:          // 0 110001 1
    1583                   case 0x47:          // 0 100011 1
    1584                    
    1585                     trigger = TRUE;
    1586                     break;
    1587                    
    1588                   default:
    1589                    
    1590                     trigger = FALSE;
    1591                    
    1592                   }
    1593                  
    1594                 } else {              // 4-fold trigger
    1595                  
    1596                   switch ( triggerBits ) {
    1597                    
    1598                   case 0x1f:          // 0 001111 1
    1599                   case 0x3d:          // 0 011110 1
    1600                   case 0x79:          // 0 111100 1
    1601                   case 0x73:          // 0 111001 1
    1602                   case 0x67:          // 0 110011 1
    1603                   case 0x4f:          // 0 100111 1
    1604                    
    1605                     trigger = TRUE;
    1606                     break;
    1607                    
    1608                   default:
    1609                    
    1610                     trigger = FALSE;
    1611                    
    1612                   }
    1613                  
    1614                 }
    1615                
    1616                 mxgrp = MAX(ngrpq0,mxgrp);
    1617                
    1618                 break;
    1619                
    1620               case 5:
    1621               case 6:
    1622                
    1623                 trigger = TRUE;
    1624                 break;
    1625                
    1626               default:
    1627                
    1628                 trigger = FALSE;
    1629                 error( SIGNATURE, "Number of neighbours > 6 !!! Exiting.\n\n");
    1630                 break;
    1631                
    1632               } // switch (ngrpq0)
    1633              
    1634             } // ct_Type
    1635            
    1636           } // for each pixel i
    1637          
    1638           if ( trigger == FALSE ) {
    1639             break;
    1640           } // end if
    1641          
    1642         } //end  for each threshold
    1643         maxtrigthr_phe = (float) k-1;  // i.e. maxtrigthr_phe < 0. means that
    1644         // the event doesn't even pass threshold 0.
    1645         //  maxtrigthr_phe >= 0 means, the event passes some threshold
    1646         //  or (in case the input parameter "threshold" was > 0), the event
    1647         //  passes the threshold given by the input parameter.
    1648         if ( maxtrigthr_phe >= 0. ) {
    1649           trigger = TRUE;
    1650         }
    1651        
    1652        
    1653         novq0 = noverq0;
     976
     977        }// end if(simulateNSB) ...
     978
     979
     980//      cout << arrtmin_ns << " " << arrtmax_ns << "\n";
     981//      for(i=0; i<cam.inumpixels; i++){
     982//        cout << i << " " << baseline_mv[i] <<"\n";
     983//      }
     984 
     985
     986        cout << "Total number of phes: " << inumphe << "\n";
     987
     988        // TRIGGER HERE
     989
     990        trigger = FALSE;
     991
     992        if( TRUE ) trigger = TRUE; // put your trigger function here
     993
     994//      for( i=0; i<inumphe; i++){
     995//        cout << "phe " << photoe[i].ipixnum << " " << photoe[i].iarrtime_ns << "\n";
     996//      }
    1654997       
    1655998        if ( trigger == TRUE ) {
    … …  
    16601003          memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels );
    16611004         
    1662 #ifdef __TAILCUT__
    1663          
    1664           //!@' @#### Tail-cut condition.
    1665           //@'
    1666          
    1667           //++
    1668           // tail-cut
    1669           //--
    1670          
    1671           // Tail-Cut = 0   : No Tail-Cut
    1672           // Tail-Cut > 0   : Make Tail-Cut
    1673           // Tail-Cut < 0   : Make Tail-Cut with t_0 = Sqrt[ maximum ]
    1674          
    1675           if (qTailCut > 0.0) {
    1676            
    1677             for ( i=0; i<ct_NPixels; ++i )
    1678               if ( fnpix[i] < qTailCut )
    1679                 fnpixclean[i] = 0.0;
    1680            
    1681           } else if (qTailCut < 0.0) {
    1682            
    1683             maxcharge = sqrt(maxcharge);
    1684             for ( i=0; i<ct_NPixels; ++i )
    1685               if ( fnpix[i] < maxcharge )
    1686                 fnpixclean[i] = 0.0;
    1687            
    1688           }
    1689          
    1690 #endif // __TAILCUT__
    16911005         
    16921006#ifdef __ISLANDS__
    … …  
    17051019#endif // __ISLANDS__
    17061020         
    1707 #ifdef __MOMENTS__
    17081021         
    17091022          //!@' @#### Calculation of parameters of the image.
    … …  
    17971110          }
    17981111         
    1799          
    1800 #endif // __MOMENTS__
    1801          
    1802          
    18031112          // revert the fnpixclean matrix into fnpix
    18041113          // (now we do this, but maybe in a future we want to
    … …  
    18891198         
    18901199        } // trigger == FALSE
    1891        
    1892 #endif // __TRIGGER__
    1893        
     1200               
    18941201        //!@' @#### Save data.
    18951202        //@'
    … …  
    19041211        // save the image to the file
    19051212        //--
    1906        
    1907         // write MCEventHeader to output file
    1908        
    1909         outputfile.write( (char *)&mcevth, mcevth.mysize() );
    1910        
    1911 #ifdef __TRIGGER__
    1912        
    1913         // save the image
    1914        
    1915         if ( (trigger == TRUE) || (Write_All_Images == TRUE) )
    1916           outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) );
    1917        
    1918 #else
    1919        
    1920         // save the image
    1921        
    1922         outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) );
    1923        
    1924 #endif // __TRIGGER__
    1925        
    1926         if ( Data_From_STDIN )
    1927           cin.read( flag, SIZE_OF_FLAGS );
    1928         else
    1929           inputfile.read ( flag, SIZE_OF_FLAGS );
    1930        
     1213
     1214
     1215
     1216        // look for the next event
     1217       
     1218        fread( flag, SIZE_OF_FLAGS, 1, inputfile );
     1219
    19311220      } // end while there is a next event
    19321221     
    … …  
    19381227        log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
    19391228       
    1940         if ( Data_From_STDIN )
    1941           cin.read( flag, SIZE_OF_FLAGS );
    1942         else
    1943           inputfile.read ( flag, SIZE_OF_FLAGS );
     1229        fread( flag, SIZE_OF_FLAGS, 1, inputfile );
    19441230       
    19451231        if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
    … …  
    19471233          still_in_loop  = FALSE;
    19481234         
    1949           if ((! Data_From_STDIN) && (! inputfile.eof())){
     1235          if ((! Data_From_STDIN) && ( !feof(inputfile) )){
    19501236           
    19511237            // we have concatenated input files.
    19521238            // get signature of the next part and check it.
    1953             // NOTE: this part repeats further up in the code;
    1954             // if you change something here you probably want to change it
    1955             // there as well
    1956            
    1957             strcpy(Signature, REFL_SIGNATURE);
    1958            
    1959             strcpy(sign, Signature);
    1960            
    1961             inputfile.read( (char *)sign, strlen(Signature));
    1962            
    1963             if (strcmp(sign, Signature) != 0) {
    1964               cerr << "ERROR: Signature of .rfl file is not correct\n";
    1965               cerr << '"' << sign << '"' << '\n';
    1966               cerr << "should be: " << Signature << '\n';
     1239
     1240            if(check_reflector_file( inputfile )==FALSE){
    19671241              exit(1);
    19681242            }
    19691243           
    1970             if ( Data_From_STDIN )
    1971               cin.read( (char *)sign, 1);
    1972             else
    1973               inputfile.read( (char *)sign, 1);
    1974            
    19751244          }
    19761245         
    19771246        } // end if found end of file
    19781247      } // end if found end of run
    1979       if ( Data_From_STDIN )
    1980         cin.read( flag, SIZE_OF_FLAGS );
    1981       else
    1982         inputfile.read ( flag, SIZE_OF_FLAGS );
     1248
     1249      fread( flag, SIZE_OF_FLAGS, 1, inputfile );
     1250
    19831251    } // end if else found start of run
    19841252  } // end big while loop
    19851253 
    1986 //!@' @#### End of program.
    1987 //@'
    1988  
    1989   //end my version
    1990 
    1991 #ifdef __ROOT__
    19921254  //++
    19931255  // put the Event to the root file
    … …  
    19981260  outfile.Close() ;
    19991261
    2000 #endif // __ROOT__
    20011262 
    20021263  // close input file
    20031264 
    2004   ntcph += ncph;
    20051265  log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
    20061266       ntshow, ntcph );
    … …  
    20121272  log( SIGNATURE, "Closing files\n" );
    20131273 
    2014   inputfile.close();
    2015   outputfile.close();
     1274  if( ! Data_From_STDIN ){
     1275    fclose( inputfile );
     1276  }
    20161277  datafile.close();
    20171278
    … …  
    20191280
    20201281  hfile->Close();
    2021 
    2022 #ifdef __DETAIL_TRIGGER__
    2023   // Output of Trigger statistics
    2024   //
    2025 
    2026   Trigger.PrintStat() ;
    2027 #endif // __DETAIL_TRIGGER__
    20281282
    20291283  // program finished
    … …  
    21381392
    21391393  //  Display the name of the function that called error
    2140   fprintf(stderr, "ERROR in %s: ", funct);
     1394  fprintf(stdout, "ERROR in %s: ", funct);
    21411395
    21421396  // Display the remainder of the message
    21431397  va_start(args, fmt);
    2144   vfprintf(stderr, fmt, args);
     1398  vfprintf(stdout, fmt, args);
    21451399  va_end(args);
    21461400
    … …  
    24251679  ifstream qefile;
    24261680  char line[LINE_MAX_LENGTH];
    2427   int n, i, j, k;
     1681  int n, i, j, icount;
    24281682  float qe;
    24291683
    … …  
    24381692
    24391693  // initialize pixel numbers
    2440 
    2441   for ( i=0; i<PIX_ARRAY_SIDE; ++i )
    2442     for ( j=0; j<PIX_ARRAY_SIDE; ++j )
    2443       pixels[i][j][PIXNUM] = -1;
    24441694
    24451695  pixary = new float* [2*ct_NCentralPixels];
    … …  
    24621712  igen_pixel_coordinates(pcam);
    24631713
    2464   // transfer coordinates to the working arrays for
    2465   // the central pixels
    2466 
    2467   for(k=0; k<ct_NCentralPixels; k++){
    2468 
    2469     i = (int) pcam->di[k];
    2470     j = (int) pcam->dj[k];
    2471 
    2472     pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k;
    2473     pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k];
    2474     pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k];
    2475    
    2476     pixary[k][0] = pcam->dxc[k];
    2477     pixary[k][1] = pcam->dyc[k];
    2478   }
    24791714
    24801715  // calculate tables of neighbours
    … …  
    25191754  // try to open the file
    25201755
    2521   log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE);
     1756  log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
    25221757 
    25231758  qefile.open( QE_FILE );
    … …  
    25331768
    25341769  i=-1;
     1770  icount = 0;
    25351771
    25361772  while ( ! qefile.eof() ) {         
    … …  
    25771813    // get the values (num-pixel, num-datapoint, QE-value)
    25781814   
    2579     sscanf(line, "%d %d %f", &i, &j, &qe);
     1815    if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )
     1816      break;
    25801817
    25811818    if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
    … …  
    25851822    }
    25861823
     1824    if ( i > ct_NPixels) break;
     1825
     1826    icount++;
     1827
    25871828  }
    25881829
     1830  if(icount/pointsQE < ct_NPixels){
     1831    error( "read_pixels", "The quantum efficiency file is faulty\n  (found only %d pixels instead of %d).\n",
     1832           icount/pointsQE, ct_NPixels );
     1833  }
     1834
    25891835  // close file
    25901836
    25911837  qefile.close();
     1838
     1839  // test QE
     1840
     1841  for(icount=0; icount< ct_NPixels; icount++){
     1842    for(i=0; i<pointsQE; i++){
     1843      if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
     1844          QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
     1845        error( "read_pixels", "The quantum efficiency file is faulty\n  pixel %d, point %d  is % f, %f\n",
     1846               icount, i, QE[icount][0][i], QE[icount][1][i] );
     1847      }
     1848    }
     1849  }
    25921850
    25931851  // end
    … …  
    28752133  } /* end if */
    28762134
    2877   /* generate the ij coordinates */
    2878 
    2879   for(i=0; i<pcam->inumpixels; i++){
    2880     pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize;
    2881     pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60;
    2882 
    2883     //    fprintf(stdout, "%d %f %f %f %f %f\n",
    2884     //      i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i],
    2885     //      pcam->dpixsizefactor[i]);
    2886 
    2887   }
    2888 
    28892135  return(pcam->inumpixels);
    28902136
    … …  
    29682214          );
    29692215}
     2216
     2217//------------------------------------------------------------
     2218// @name check_reflector_file                         
     2219//                                     
     2220// @desc check if a given reflector file has the right signature
     2221// @desc return TRUE or FALSE
     2222//
     2223// @date Mon Feb 14 16:44:21 CET 2000
     2224// @function @code
     2225//------------------------------------------------------------
     2226
     2227int check_reflector_file(FILE *infile){
     2228
     2229  char Signature[20];           // auxiliary variable
     2230  char sign[20];                // auxiliary variable
     2231
     2232  strcpy(Signature, REFL_SIGNATURE);
     2233 
     2234  strcpy(sign, Signature);
     2235 
     2236  fread( (char *)sign, strlen(Signature), 1, infile);
     2237
     2238  if (strcmp(sign, Signature) != 0) {
     2239    cout << "ERROR: Signature of .rfl file is not correct\n";
     2240    cout << '"' << sign << '"' << '\n';
     2241    cout << "should be: " << Signature << '\n';
     2242    return(FALSE);
     2243  }
     2244
     2245  fread( (char *)sign, 1, 1, infile);
     2246
     2247  return(TRUE);
     2248
     2249}
     2250
     2251//------------------------------------------------------------
     2252// @name lin_interpol                         
     2253//                                     
     2254// @desc interpolate linearly between two points returning the
     2255// @desc y-value of the result
     2256//
     2257// @date Thu Feb 17 11:31:32 CET 2000
     2258// @function @code
     2259//------------------------------------------------------------
     2260
     2261float lin_interpol( float x1, float y1, float x2, float y2, float x){
     2262
     2263  if( (x2 - x1)<=0. ){ // avoid division by zero, return average
     2264    cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
     2265    return((y1+y2)/2.);
     2266  }
     2267  else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
     2268    return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
     2269  }
     2270}
     2271
     2272
     2273//------------------------------------------------------------
     2274// @name Photoelectron                         
     2275//                                     
     2276// @desc constructor for class Photoelectron
     2277//
     2278// @date Mon Feb 15 16:44:21 CET 2000
     2279// @function @code
     2280//------------------------------------------------------------
     2281
     2282Photoelectron::Photoelectron(void){
     2283  iarrtime_ns = NOTIME;
     2284  ipixnum = -1;
     2285}
     2286
     2287//------------------------------------------------------------
     2288// @name produce_phes                         
     2289//                                     
     2290// @desc read the photons of an event, pixelize them and simulate QE
     2291// @desc return various statistics and the array of Photoelectrons
     2292//
     2293// @date Mon Feb 14 16:44:21 CET 2000
     2294// @function @code
     2295//------------------------------------------------------------
     2296
     2297int produce_phes( FILE *sp, // the input file
     2298                  struct camera *cam, // the camera layout
     2299                  float minwl_nm, // the minimum accepted wavelength
     2300                  float maxwl_nm, // the maximum accepted wavelength
     2301                  class Photoelectron phe[iMAXNUMPHE], // the generated phes
     2302                  int *itotnphe, // total number of produced photoelectrons
     2303                  float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
     2304                  int *incph,    // total number of cph read
     2305                  float *tmin_ns,   // minimum arrival time of all phes
     2306                  float *tmax_ns    // maximum arrival time of all phes
     2307                  ){
     2308
     2309  static int i, k, ipixnum;
     2310  static float cx, cy, wl, qe, t;
     2311  static MCCphoton photon;
     2312  static float **qept;
     2313  static char flag[SIZE_OF_FLAGS + 1];
     2314  static float radius;
     2315
     2316  // reset variables
     2317
     2318  for ( i=0; i<cam->inumpixels; ++i ){
     2319   
     2320    nphe[i] = 0.0;
     2321   
     2322  }
     2323
     2324  *itotnphe = 0;
     2325  *incph = 0;
     2326  *tmin_ns = NOTIME; // very big
     2327  *tmax_ns = -NOTIME; // very small
     2328
     2329  radius = cam->dxc[cam->inumpixels-1]
     2330    + 1.5*cam->dpixdiameter_cm*cam->dpixsizefactor[cam->inumpixels-1];
     2331
     2332  //- - - - - - - - - - - - - - - - - - - - - - - - -
     2333  // read photons and "map" them into the pixels
     2334  //--------------------------------------------------     
     2335 
     2336  // initialize CPhoton
     2337 
     2338  photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
     2339 
     2340  // read the photons data
     2341 
     2342  fread ( flag, SIZE_OF_FLAGS, 1, sp );
     2343 
     2344  // loop over the photons
     2345 
     2346  while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
     2347         
     2348    memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
     2349   
     2350    fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
     2351         
     2352    // increase number of photons
     2353         
     2354    (*incph)++;
     2355                 
     2356    //
     2357    // Pixelization
     2358    //
     2359           
     2360    cx = photon.get_x();
     2361    cy = photon.get_y();
     2362         
     2363    // get wavelength
     2364         
     2365    wl = photon.get_wl();
     2366
     2367    // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n";
     2368         
     2369    // check if photon has valid wavelength and is inside outermost camera radius
     2370         
     2371    if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){
     2372           
     2373      // cout << " lost first check\n";
     2374
     2375      // read next CPhoton
     2376
     2377      fread ( flag, SIZE_OF_FLAGS, 1, sp );
     2378           
     2379      // go to beginning of loop, the photon is lost
     2380      continue;
     2381           
     2382    }
     2383
     2384    ipixnum = -1;
     2385
     2386    for(i=0; i<cam->inumpixels; i++){
     2387      if( bpoint_is_in_pix( cx, cy, i, cam) ){
     2388        ipixnum = i;
     2389        break;
     2390      }
     2391    }
     2392           
     2393    if(ipixnum==-1){// the photon is in none of the pixels
     2394
     2395      // cout << " lost pixlization\n";
     2396
     2397      // read next CPhoton
     2398
     2399      fread ( flag, SIZE_OF_FLAGS, 1, sp );
     2400     
     2401      // go to beginning of loop, the photon is lost
     2402      continue;
     2403    }
     2404                 
     2405    //+++
     2406    // QE simulation
     2407    //---
     2408   
     2409    // set pointer to the QE table of the relevant pixel
     2410         
     2411    qept = (float **)QE[ipixnum];
     2412         
     2413    // check if wl is inside table; outside the table, QE is assumed to be zero
     2414
     2415    if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
     2416
     2417      // cout << " lost\n";
     2418     
     2419      // read next Photon
     2420
     2421      fread ( flag, SIZE_OF_FLAGS, 1, sp );
     2422           
     2423      // go to beginning of loop
     2424      continue;
     2425           
     2426    }
     2427
     2428    // find data point in the QE table (-> k)
     2429
     2430    k = 1; // start at 1 because the condition was already tested for 0
     2431    while (k < pointsQE-1 && qept[0][k] < wl){
     2432      k++;
     2433    }
     2434
     2435    // calculate the qe value between 0. and 1.
     2436
     2437    qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
     2438   
     2439    // if random > quantum efficiency, reject it
     2440               
     2441    if ( RandomNumber > qe ) {
     2442
     2443      // cout << " lost\n";
     2444           
     2445      // read next Photon
     2446
     2447      fread ( flag, SIZE_OF_FLAGS, 1, sp );
     2448           
     2449      // go to beginning of loop
     2450      continue;
     2451           
     2452    }
     2453         
     2454    //+++
     2455    // The photon has produced a photo electron
     2456    //---
     2457
     2458    // cout << " accepted\n";
     2459         
     2460    // increment the number of photoelectrons in the relevant pixel
     2461         
     2462    nphe[ipixnum] += 1.0;
     2463
     2464    t = photon.get_t() ;
     2465
     2466    //    cout << " t " << t;
     2467         
     2468    // find minimum and maximum arrival time
     2469
     2470    if(t < *tmin_ns){
     2471      *tmin_ns = t; // memorize time
     2472    }
     2473    if(t > *tmax_ns){
     2474      *tmax_ns = t; // memorize time
     2475    }
     2476
     2477    // store the new photoelectron
     2478
     2479    if(*itotnphe >= iMAXNUMPHE){
     2480      cout << "Error: Memory overflow. Event produces more than maximum\n";
     2481      cout << "       allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
     2482      return(1);
     2483    }
     2484
     2485    phe[*itotnphe].iarrtime_ns = (int)t;
     2486    phe[*itotnphe].ipixnum = ipixnum;
     2487         
     2488    *itotnphe += 1;
     2489
     2490    // read next Photon
     2491         
     2492    fread( flag, SIZE_OF_FLAGS, 1, sp );
     2493         
     2494  }  // end while, i.e. found end of event
     2495         
     2496  return(0);
     2497
     2498}
     2499
     2500
     2501//------------------------------------------------------------
     2502// @name produce_nsbrates
     2503//                                     
     2504// @desc read the starfield file, call produce_phes on it in,
     2505// @desc different wavebands, calculate the nsbrates
     2506//
     2507// @date Mon Feb 14 16:44:21 CET 2000
     2508// @function @code
     2509//------------------------------------------------------------
     2510
     2511int produce_nsbrates( char *iname, // the starfield input file name
     2512                      struct camera *cam, // camera layout
     2513                      class Photoelectron phe[iMAXNUMPHE], // array for photoelectrons
     2514                      float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
     2515                                               // the NSB rates in phe/ns for each pixel
     2516                      ){
     2517
     2518  int i, j, k, ii; // counters
     2519                     
     2520  static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
     2521                                            WAVEBANDBOUND2,
     2522                                            WAVEBANDBOUND3,
     2523                                            WAVEBANDBOUND4,
     2524                                            WAVEBANDBOUND5,
     2525                                            WAVEBANDBOUND6 };
     2526
     2527  FILE *infile;    // the input file
     2528  fpos_t fileposition; // marker on the input file
     2529  static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
     2530  static MCEventHeader evth; // the event header
     2531  static float nphe[iMAXNUMPIX];    // the number of photoelectrons in each pixel
     2532  int itnphe;   // total number of produced photoelectrons
     2533  int itotnphe;   // total number of produced photoelectrons after averaging
     2534  int incph;      // total number of cph read
     2535  float tmin_ns;  // minimum arrival time of all phes
     2536  float tmax_ns;  // maximum arrival time of all phes
     2537  float integtime_ns; // integration time from the starfield generator
     2538
     2539  // open input file
     2540 
     2541  log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
     2542
     2543  infile = fopen( iname, "r" );
     2544  if ( infile == NULL )
     2545    error( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
     2546 
     2547  // get signature, and check it
     2548 
     2549  if(check_reflector_file(infile)==FALSE){
     2550    exit(1);
     2551  }
     2552
     2553  // initialize flag
     2554 
     2555  strcpy( cflag, "                                        \0" );
     2556
     2557  // get flag
     2558   
     2559  fread( cflag, SIZE_OF_FLAGS, 1, infile );
     2560
     2561  if ( ! feof(infile)){
     2562
     2563    // reading .rfl file
     2564
     2565    if(!isA( cflag, FLAG_START_OF_RUN )){
     2566      error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
     2567    }
     2568    else { // found start of run
     2569     
     2570      fread( cflag, SIZE_OF_FLAGS, 1, infile );
     2571     
     2572      if( isA( cflag, FLAG_START_OF_EVENT   )){ // there is a event
     2573       
     2574        // get MCEventHeader
     2575       
     2576        fread( (char*)&evth, evth.mysize(), 1, infile );
     2577       
     2578        integtime_ns = evth.get_energy();
     2579
     2580        // memorize where we are in the file
     2581
     2582        if (fgetpos( infile, &fileposition ) != 0){
     2583          error( SIGNATURE, "Cannot position in file ...\n");
     2584        }
     2585
     2586        // loop over the wavebands
     2587
     2588        for(i=0; i<iNUMWAVEBANDS; i++){
     2589
     2590          // initialize the rate array
     2591
     2592          for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
     2593            rate_phepns[j][i] = 0.;
     2594          }
     2595
     2596          itotnphe = 0;
     2597
     2598          // read the photons and produce the photoelectrons
     2599          // - in order to average over the QE simulation, call the
     2600          // production function iNUMNSBPRODCALLS times
     2601
     2602          for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
     2603
     2604            // position the file pointer to the beginning of the photons
     2605           
     2606            fsetpos( infile, &fileposition);
     2607
     2608            // produce photoelectrons
     2609
     2610            k = produce_phes( infile,
     2611                              cam,
     2612                              wl_nm[i],
     2613                              wl_nm[i+1],
     2614                              phe, // this is a dummy here
     2615                              &itnphe,
     2616                              nphe, // we want this!
     2617                              &incph,
     2618                              &tmin_ns,
     2619                              &tmax_ns );
     2620
     2621            if( k != 0 ){ // non-zero returnvalue means error
     2622              cout << "Exiting.\n";
     2623              exit(1);
     2624            }
     2625         
     2626            for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
     2627              rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
     2628            }
     2629
     2630            itotnphe += itnphe;
     2631
     2632          } // end for(ii=0 ...
     2633
     2634          fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
     2635                  wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
     2636
     2637        } // end for(i=0 ...
     2638
     2639      }
     2640      else{
     2641        cout << "Starfield file contains no event.\nExiting.\n";
     2642        exit(1);
     2643      } // end if( isA ... event
     2644    } // end if ( isA ... run
     2645  }
     2646  else{
     2647    cout << "Starfield file contains no run.\nExiting.\n";
     2648    exit(1);
     2649  }
     2650
     2651  fclose( infile );
     2652
     2653  return(0);
     2654
     2655}
     2656
     2657
     2658//------------------------------------------------------------
     2659// @name produce_nsb_phes
     2660//                                     
     2661// @desc produce the photoelectrons from the nsbrates
     2662//
     2663// @date Thu Feb 17 17:10:40 CET 2000
     2664// @function @code
     2665//------------------------------------------------------------
     2666
     2667int produce_nsb_phes( float *atmin_ns,
     2668                      float *atmax_ns,
     2669                      float theta_rad,
     2670                      struct camera *cam,
     2671                      float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
     2672                      float difnsbr_phepns[iMAXNUMPIX],
     2673                      float extinction[iNUMWAVEBANDS],
     2674                      float fnpx[iMAXNUMPIX],
     2675                      Photoelectron photo[iMAXNUMPHE],
     2676                      int *inphe,
     2677                      float base_mv[iMAXNUMPIX]){
     2678
     2679  float simtime_ns; // NSB simulation time
     2680  int i, j, k, ii;
     2681  float zenfactor;  //  correction factor calculated from the extinction
     2682  int inumnsbphe;   //  number of photoelectrons caused by NSB
     2683
     2684  ii = *inphe; // avoid dereferencing
     2685               
     2686  // check if the arrival times are set; if not generate them
     2687
     2688  if(*atmin_ns == NOTIME){
     2689
     2690    *atmin_ns = 0.;
     2691    *atmax_ns = simtime_ns = SLICES*WIDTH_TIMESLICE;
     2692
     2693  }
     2694  else{ // extend the event time window by the given offsets
     2695
     2696    *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
     2697    *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
     2698
     2699    simtime_ns = *atmax_ns - *atmin_ns;
     2700
     2701    // make sure the simulated time is long enough for the FADC simulation
     2702
     2703    if(simtime_ns< SLICES*WIDTH_TIMESLICE){
     2704      *atmax_ns = *atmin_ns + SLICES*WIDTH_TIMESLICE;
     2705      simtime_ns = SLICES*WIDTH_TIMESLICE;
     2706    }
     2707     
     2708  }
     2709
     2710  // initialize baselines
     2711
     2712  for(i=0; i<cam->inumpixels; i++){
     2713    base_mv[i] = 0.;
     2714  }
     2715
     2716  // calculate baselines and generate phes
     2717
     2718  for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
     2719   
     2720    // calculate the effect of the atmospheric extinction
     2721   
     2722    zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
     2723   
     2724    for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
     2725     
     2726      inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS) 
     2727                          * simtime_ns );
     2728
     2729      base_mv[j] += inumnsbphe;
     2730
     2731      // randomize
     2732     
     2733      inumnsbphe =  ignpoi( inumnsbphe );
     2734     
     2735      // create the photoelectrons
     2736     
     2737      for(k=0; k < inumnsbphe; k++){
     2738
     2739        if(ii >= iMAXNUMPHE){
     2740          cout << "Error: Memory overflow. NSB simulation produces more than maximum\n";
     2741          cout << "       allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
     2742          return(1);
     2743        }
     2744
     2745        photo[ii].iarrtime_ns = (int)(RandomNumber * simtime_ns + *atmin_ns );
     2746        photo[ii].ipixnum = j;
     2747       
     2748        // cout << "Created phe " << photo[ii].iarrtime_ns << " "
     2749        //     << photo[ii].ipixnum << "\n";
     2750       
     2751        ii++; // increment total number of photoelectons
     2752       
     2753        fnpx[j] += 1.; // increment number of photoelectrons in this pixel
     2754
     2755      }
     2756     
     2757    } // end for(j=0 ...
     2758  } // end for(i=0 ...
     2759
     2760  // finish baseline calculation
     2761
     2762  for(i=0; i<cam->inumpixels; i++){
     2763    base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
     2764  }
     2765
     2766  *inphe = ii; // update the pointer
     2767
     2768  return(0);
     2769}
     2770
     2771
    29702772// @endcode
    29712773
    … …  
    29772779//
    29782780// $Log: not supported by cvs2svn $
     2781// Revision 1.4  2000/01/25 08:36:23  petry
     2782// The pixelization in previous versions was buggy.
     2783// This is the first version with a correct pixelization.
     2784//
    29792785// Revision 1.3  2000/01/20 18:22:17  petry
    29802786// Found little bug which makes camera crash if it finds a photon
    … …  
    29962802// The "rootification" was done by Dirk Petry and Harald Kornmayer.
    29972803//
    2998 // In the following you can see the README file of that version:
    2999 //
    3000 // ==================================================
    3001 //
    3002 // Fri Oct 22  1999   D.P.
    3003 //
    3004 // The MAGIC Monte Carlo System
    3005 //
    3006 // Camera Simulation Programme
    3007 // ---------------------------
    3008 //
    3009 // 1) Description
    3010 //
    3011 // This version is the result of the fusion of H.K.'s
    3012 // root_camera which is described below (section 2)
    3013 // and another version by D.P. which had a few additional
    3014 // useful features.
    3015 //
    3016 // The version compiles under Linux with ROOT 2.22 installed
    3017 // (variable ROOTSYS has to be set).
    3018 //
    3019 // Compile as before simply using "make" in the root_camera
    3020 // directory.
    3021 //
    3022 // All features of H.K.'s root_camera were retained.
    3023 //
    3024 // Additional features of this version are:
    3025 //
    3026 //   a) HBOOK is no longer used and all references are removed.
    3027 //
    3028 //   b) Instead of HBOOK, the user is given now the possibility of
    3029 //      having Diagnostic data in ROOT format as a complement
    3030 //      to the ROOT Raw data.
    3031 //
    3032 //      This data is written to the file which is determined by
    3033 //      the new input parameter "diag_file" in the camera parameter
    3034 //      file.
    3035 //
    3036 //      All source code file belonging to this part have filenames
    3037 //      starting with "MDiag".
    3038 //
    3039 //      The user can read the output file using the following commands
    3040 //      in an interactive ROOT session:
    3041 //
    3042 //              root [0] .L MDiag.so
    3043 //      root [1] new TFile("diag.root");
    3044 //      root [2] new TTreeViewer("T");
    3045 //     
    3046 //      This brings up a viewer from which all variables of the
    3047 //      TTree can be accessed and histogrammed. This example
    3048 //      assumes that you have named the file "diag.root", that
    3049 //      you are using ROOT version 2.22 or later and that you have
    3050 //      the shared object library "MDiag.so" which is produced
    3051 //      by the Makefile along with the executable "camera".
    3052 //       
    3053 //  !   The contents of the so-called diag file is not yet fixed.
    3054 //  !   At the moment it is what J.C.G. used to put into the HBOOK
    3055 //  !   ntuple. In future versions the moments calculation can be
    3056 //  !   removed and the parameter list be modified correspondingly.
    3057 //
    3058 //   c) Now concatenated reflector files can be read. This is useful
    3059 //      if you have run the reflector with different parameters but
    3060 //      you want to continue the analysis with all reflector data
    3061 //      going into ONE ROOT outputfile.
    3062 //
    3063 //      The previous camera version contained a bug which made reading
    3064 //      of two or more concatenated reflector files impossible.
    3065 //
    3066 //   d) The reflector output format was changed. It is now version
    3067 //      0.4 .
    3068 //      The change solely consists in a shortening of the flag
    3069 //      definition in the file
    3070 //
    3071 //            include-MC/MCCphoton.hxx 
    3072 //
    3073 //  !   IF YOU WANT TO READ REFLECTOR FORMAT 0.3, you can easily
    3074 //  !   do so by recompiling camera with the previous version of
    3075 //  !   include-MC/MCCphoton.hxx.
    3076 //
    3077 //      The change was necessary for saving space and better
    3078 //      debugging. From now on, this format can be frozen.
    3079 //
    3080 //  !   For producing reflector output in the new format, you
    3081 //  !   of course have to recompile your reflector with the
    3082 //  !   new include-MC/MCCphoton.hxx .
    3083 //
    3084 //   e) A first version of the pixelization with the larger
    3085 //      outer pixels is implemented. THIS IS NOT YET FULLY
    3086 //      TESTED, but first rough tests show that it works
    3087 //      at least to a good approximation.
    3088 //
    3089 //      The present version implements the camera outline
    3090 //      with 18 "gap-pixels" and 595 pixels in total as
    3091 //      shown in
    3092 //
    3093 //         http://sarastro.ifae.es/internal/home/hardware/camera/numbering.ps
    3094 //
    3095 //      This change involved
    3096 //
    3097 //      (i) The file pixels.dat is no longer needed. Instead
    3098 //           the coordinates are generated by the program itself
    3099 //           (takes maybe 1 second). In the file
    3100 //
    3101 //              pixel-coords.txt
    3102 //
    3103 //        in the same directory as this README, you find a list
    3104 //           of the coordinates generated by this new routine. It
    3105 //           has the format
    3106 //
    3107 //               number   i   j   x  y  size-factor
    3108 //
    3109 //           where i and j are J.C.G.'s so called biaxis hexagonal
    3110 //           coordinates (for internal use) and x and y are the
    3111 //           coordinates of the pixel centers in the standard camera
    3112 //           coordinate system in units of centimeters. The value
    3113 //           of "size-factor" determines the linear size of the pixel
    3114 //           relative to the central pixels.
    3115 //
    3116 //         (ii) The magic.def file has two additional parameters
    3117 //           which give the number of central pixels and the
    3118 //           number of gap pixels
    3119 //
    3120 //         (iii) In camera.h and camera.cxx several changes were
    3121 //           necessary, among them the introduction of several
    3122 //           new functions
    3123 //
    3124 //      The newly suggested outline with asymmetric Winston cones
    3125 //      will be implemented in a later version.
    3126 //
    3127 //   f) phe files can no longer be read since this contradicts
    3128 //      our philosophy that the analysis should be done with other
    3129 //      programs like e.g. EVITA and not with "camera" itself.
    3130 //      This possibility was removed.
    3131 //
    3132 //   g) ROOT is no longer invoked with an interactive interface.
    3133 //      In this way, camera can better be run as a batch program and
    3134 //      it uses less memory.
    3135 //
    3136 //   h) small changes concerning the variable "t_chan" were necessary in
    3137 //      order to avoid segmentation faults: The variable is used as an
    3138 //      index and it went sometimes outside the limits when camera
    3139 //      was reading proton data. This is because the reflector files
    3140 //      don't contain the photons in a chronological order and also
    3141 //      the timespread can be considerably longer that the foreseen
    3142 //      digitisation timespan. Please see the source code of camera.cxx
    3143 //      round about line 1090.
    3144 //
    3145 //   j) several unused variables were removed, a few warning messages
    3146 //      occur when you compile camera.cxx but these can be ignored at
    3147 //      the moment.
    3148 //
    3149 // In general the program is of course not finished. It still needs
    3150 // debugging, proper trigger simulation, simulation of the asymmetric
    3151 // version of the outer pixels, proper NSB simulation, adaption of
    3152 // the diag "ntuple" contents to our need and others small improvements.
    3153 //
    3154 // In the directory rfl-files there is now a file in reflector format 0.4
    3155 // containing a single event produced by the starfiled adder. It has
    3156 // a duration of 30 ns and represents the region around the Crab Nebula.
    3157 // Using the enclosed input parameter file, camera should process this
    3158 // file without problems.
    3159 //
    3160 // 2) The README for the previous version of root_camera
    3161 //
    3162 // README for a preliminary version of the
    3163 // root_camera program.
    3164 //
    3165 // root_camera is based on the program "camera"of Jose Carlos
    3166 // Gonzalez. It was changed in the way that only the pixelisation
    3167 // and the distibution of the phe to the FADCs works in a
    3168 // first version.
    3169 //
    3170 // Using the #undef command most possibilities of the orignal
    3171 // program are switched of.
    3172 //
    3173 // The new parts are signed by
    3174 //
    3175 // - ROOT or __ROOT__
    3176 //   nearly all  important codelines for ROOT output are enclosed
    3177 //   in structures like
    3178 //   #ifdef __ROOT__
    3179 //   
    3180 //     code
    3181 //
    3182 //   #endif __ROOT__
    3183 //
    3184 //   In same case the new lines are signed by a comment with the word
    3185 //   ROOT in it.
    3186 //
    3187 //   For timing of the pulse some variable names are changed.
    3188 //   (t0, t1, t  -->  t_ini, t_fin, t_1st, t_chan,...)
    3189 //   Look also for this changes.
    3190 //
    3191 //   For the new root-file is also a change in readparm-files
    3192 //
    3193 //
    3194 // - __DETAIL_TRIGGER__
    3195 //
    3196 //   This is for the implementation of the current work on trigger
    3197 //   studies. Because the class MTrigger is not well documented it
    3198 //   isn´t a part of this tar file. Only a dummy File exists.
    3199 //
    3200 //
    3201 //
    3202 // With all files in the archive, the root_camera program should run.
    3203 //
    3204 // A reflector file is in the directory rfl-files
    3205 //
    3206 // ==================================================
    3207 //
    3208 // From now on, use CVS for development!!!!
    3209 //
    3210 //
    3211 //
    32122804// Revision 1.3  1999/10/22 15:01:28  petry
    32132805// version sent to H.K. and N.M. on Fri Oct 22 1999
Note: See TracChangeset for help on using the changeset viewer.