//!///////////////////////////////////////////////////////////////////// // // camera // // @file camera.cxx // @title Camera simulation // @subtitle Code for the simulation of the camera phase // @desc Code for the simulation of the camera of CT1 and MAGIC // @author J C Gonzalez // @email gonzalez@mppmu.mpg.de // @date Thu May 7 16:24:22 1998 // //---------------------------------------------------------------------- // // Created: Thu May 7 16:24:22 1998 // Author: Jose Carlos Gonzalez // Purpose: Program for reflector simulation // Notes: See files README for details // //---------------------------------------------------------------------- // // $RCSfile: timecam.cxx,v $ // $Revision: 1.1.1.1 $ // $Author: harald $ // $Date: 2000-02-08 15:13:44 $ // //////////////////////////////////////////////////////////////////////// // @tableofcontents @coverpage //=----------------------------------------------------------- //!@section Source code of |camera.cxx|. /*!@" In this section we show the (commented) code of the program for the read-out of the output files generated by the simulator of the reflector, |reflector 0.3|. @"*/ //=----------------------------------------------------------- //!@subsection Includes and Global variables definition. //!@{ // includes for ROOT // BEWARE: the order matters! #include "TROOT.h" #include "TApplication.h" #include "TFile.h" #include "TTree.h" #include "TBranch.h" #include "TCanvas.h" #include "MTrigger.hxx" #include "MFadc.hxx" #include "MRawEvt.h" #include "MMcEvt.h" #include "MMcTrig.hxx" /*!@" All the defines are located in the file |camera.h|. @"*/ #include "timecam.h" //!@} /*!@" The following set of flags are used in time of compilation. They do not affect directly the behaviour of the program at run-time (though, of course, if you disconnected the option for implementation of the Trigger logic, you will not be able to use any trigger at all. The 'default' values mean default in the sense of what you got from the server when you obtained this program. @"*/ //!@{ // flag for debugging (default: OFF ) #define __DEBUG__ #undef __DEBUG__ // flag for NNT in CT1 camera (default: ON ) #undef __CT1_NO_NEIGHBOURS__ #define __CT1_NO_NEIGHBOURS__ // flag for calculation of NSB (default: ON ) #undef __NSB__ #define __NSB__ // flag for calculation of QE for pixels (default: ON ) #undef __QE__ #define __QE__ // flag for implementation of DETAIL_TRIGGER (default: ON ) // // This is the new implementation of Trigger studies // It relies on a better simulation of the time stucture // of the PhotoMultiplier. For more details see the // documentation of the --> class MTrigger <-- #undef __DETAIL_TRIGGER__ #define __DETAIL_TRIGGER__ // flag for implementation of TRIGGER (default: ON ) #define __TRIGGER__ #undef __TRIGGER__ // flag for implementation of Tail-Cut (default: ON ) #define __TAILCUT__ #undef __TAILCUT__ // flag for calculation of islands stat. (default: ON ) #define __ISLANDS__ #undef __ISLANDS__ // flag for calculation of image parameters (default: ON ) #define __MOMENTS__ #undef __MOMENTS__ // flag for ROOT (default: ON ) #undef __ROOT__ #define __ROOT__ // flag for INTERAKTIV (default: OFF ) #undef __INTERAKTIV__ #define __INTERAKTIV__ //!@} //=----------------------------------------------------------- //!@subsection Definition of global variables. /*!@" Now we define some global variables with data about the telescope, such as "focal distance", number of pixels/mirrors, "size of the camera", and so on. @"*/ /*!@" Depending on the telescope we are using (CT1 or MAGIC), the information stored in the definition file is different. The variable |ct_Type| has the value 0 when we use CT1, and 1 when we use MAGIC. @"*/ //!@{ static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC //!@} /*!@" And this is the information about the whole telescope. @"*/ //!@{ // parameters of the CT (from the CT definition file) ////@: Focal distances [cm] //static float *ct_Focal; //@: Mean Focal distances [cm] static float ct_Focal_mean; //@: STDev. Focal distances [cm] static float ct_Focal_std; //@: Mean Point Spread function [cm] static float ct_PSpread_mean; //@: STDev. Point Spread function [cm] static float ct_PSpread_std; //@: STDev. Adjustmente deviation [cm] static float ct_Adjustment_std; //@: Radius of the Black Spot in mirror [cm] static float ct_BlackSpot_rad; //@: Radius of one mirror [cm] static float ct_RMirror; //@: Camera width [cm] static float ct_CameraWidth; //@: Pixel width [cm] static float ct_PixelWidth; //@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60) static float ct_PixelWidth_corner_2_corner; //@: ct_PixelWidth_corner_2_corner / 2 static float ct_PixelWidth_corner_2_corner_half; //@: Number of mirrors static int ct_NMirrors = 0; //@: Number of pixels static int ct_NPixels; //@: Number of pixels static int ct_NCentralPixels; //@: Number of pixels static int ct_NGapPixels; //@: ct_Apot = ct_PixelWidth / 2 static float ct_Apot; //@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth static float ct_2Apot; //@: name of the CT definition file to use static char ct_filename[256]; //@: list of showers to be skipped static int *Skip; //@: number of showers to be skipped static int nSkip=0; //@: flag: TRUE: data come from STDIN; FALSE: from file static int Data_From_STDIN = FALSE; //@: flag: TRUE: write all images to output; FALSE: only triggered showers static int Write_All_Images = FALSE; static int Write_McEvt = TRUE; static int Write_McTrig = FALSE; static int Write_RawEvt = FALSE; //@: flag: TRUE: selection on the energy static int Select_Energy = TRUE; //@: Lower edge of the selected energy range (in GeV) static float Select_Energy_le = 0.0; //@: Upper edge of the selected energy range (in GeV) static float Select_Energy_ue = 100000.0; //!@} /*!@" The following double-pointer is a 2-dimensional table with information about each pixel. The routine read_pixels will generate the information for filling it using igen_pixel_coordinates(). @"*/ //!@{ // Pointer to a tables/Arrays with information about the pixels // and data stored on them with information about the pixels //@: table for IJ sys. static float pixels[PIX_ARRAY_SIDE][PIX_ARRAY_SIDE][4]; //@: coordinates x,y for each pixel static float **pixary; //@: indexes of pixels neighbours of a given one static int **pixneig; //@: number of neighbours a pixel have static int *npixneig; //@: contents of the pixels (ph.e.) static float *fnpix; //@: contents of the pixels (ph.e.) after cleanning static float *fnpixclean; //!@} /*!@" The following double-pointer is a 2-dimensional table with the Quantum Efficiency @$QE@$ of each pixel in the camera, as a function of the wavelength @$\lambda@$. The routine |read_pixels()| will read also this information from the file |qe.dat|. @"*/ //!@{ // Pointer to a table with QE, number of datapoints, and wavelengths //@: table of QE static float ***QE; //@: number of datapoints for the QE curve static int pointsQE; //@: table of QE static float *QElambda; //!@} /*!@" The following double-pointer is a 2-dimensional table with information about each mirror in the dish. The routine |read_ct_file()| will read this information from the CT definition file. @"*/ //!@{ // Pointer to a table with the following info.: static float **ct_data; /* * TYPE=0 (CT1) * i s rho theta x y z thetan phin xn yn zn * * i : number of the mirror * s : arc length [cm] * rho : polar rho of the position of the center of the mirror [cm] * theta : polar angle of the position of the center of the mirror [cm] * x : x coordinate of the center of the mirror [cm] * y : y coordinate of the center of the mirror [cm] * z : z coordinate of the center of the mirror [cm] * thetan : polar theta angle of the direction where the mirror points to * phin : polar phi angle of the direction where the mirror points to * xn : xn coordinate of the normal vector in the center (normalized) * yn : yn coordinate of the normal vector in the center (normalized) * zn : zn coordinate of the normal vector in the center (normalized) * * TYPE=1 (MAGIC) * i f sx sy x y z thetan phin * * i : number of the mirror * f : focal distance of that mirror * sx : curvilinear coordinate of mirror's center in X[cm] * sy : curvilinear coordinate of mirror's center in X[cm] * x : x coordinate of the center of the mirror [cm] * y : y coordinate of the center of the mirror [cm] * z : z coordinate of the center of the mirror [cm] * thetan : polar theta angle of the direction where the mirror points to * phin : polar phi angle of the direction where the mirror points to * xn : xn coordinate of the normal vector in the center (normalized) * yn : yn coordinate of the normal vector in the center (normalized) * zn : zn coordinate of the normal vector in the center (normalized) */ //!@} /*!@" We define a table into where random numbers will be stored. The routines used for random number generation are provided by |RANLIB| (taken from NETLIB, |www.netlib.org|), and by the routine |double drand48(void)| (prototype defined in |stdlib.h|) through the macro |RandomNumber| defined in |camera.h|. @"*/ //!@{ // table of random numbers // (unused) //static double RandomNumbers[500]; //!@} /*!@" The following is a variable to count the number of Cphotons in the different steps of the simulation. The definition is as follows: @[ \mbox{CountCphotons}[ \mbox{FILTER} ] \equiv \mbox{\it Number of photons after the filter} \mbox{FILTER} @] The filters are defined and can be found in the file |camera.h|. @"*/ //!@{ // vector to count photons at any given step of the simulation //static int CountCphotons[10]; //!@} /*!@" The following are the set of parameters calculated for each image. The routines for their calculations are in |moments.cxx|. @"*/ //!@{ // parameters of the images //static Moments_Info *moments_ptr; //static LenWid_Info *lenwid_ptr; //static float *maxs; //static int *nmaxs; //static float length, width, dist, xdist, azw, miss, alpha, *conc; //static float phiasym, asymx, asymy; //static float charge, smax, maxtrigthr_phe; //!@} extern char FileName[]; //=----------------------------------------------------------- // @subsection Main program. //!@{ //++++++++++++++++++++++++++++++++++++++++ // MAIN PROGRAM //---------------------------------------- int main(int argc, char **argv) { //!@' @#### Definition of variables. //@' char inname[256]; //@< input file name char rootname[256] ; //@< ROOT file name char parname[256]; //@< parameters file name char sign[20]; //@< initialize sign char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file ifstream inputfile; //@< stream for the input file MCEventHeader mcevth; //@< Event Header class (MC) MCCphoton cphoton; //@< Cherenkov Photon class (MC) float thetaCT, phiCT; //@< parameters of a given shower float thetashw, phishw; //@< parameters of a given shower float coreD, coreX, coreY; //@< core position and distance float impactD; //@< impact parameter float l1, m1, n1; //@< auxiliary variables float l2, m2, n2; //@< auxiliary variables float num, den; //@< auxiliary variables int nshow=0; //@< partial number of shower in a given run int ntshow=0; //@< total number of showers int ncph_in=0; //@< number of (input) photons in one shower int ntcph_in=0; //@< total number of input photons int ncph_out=0; //@< number of (output) photons in one shower int ntcph_out=0; //@< total number of output photons int i, ii, k; //@< simple counters float t_ini; //@< time of the first Cphoton read in float t; //@< time for a single photon float t_first ; //@< time of the first cerenkov particle reaching the ground float t_last ; //@< time of the last cerenkov particle reaching the ground int t_chan ; //@< the bin (channel) in time of a single photon int startchan ; //@< the first bin with entries in the time slices float cx, cy, ci, cj; //@< coordinates in the XY and IJ systems int ici, icj, iici, iicj; //@< coordinates in the IJ (integers) int nPMT; //@< number of pixel float wl, last_wl; //@< wavelength of the photon float qe; //@< quantum efficiency float **qeptr; int simulateNSB; //@< Will we simulate NSB? float meanNSB; //@< NSB mean value (per pixel and ns) float meanPois; //@< NSB mean value for simulation float timeNSB ; //@< time of NSB photon int nIslandsCut; //@< Islands Cut value int countIslands; //@< Will we count the islands? int anaPixels; float fCorrection; //@< Factor to apply to pixel values (def. 1.) // int ntrigger = 0; //@< number of triggers in the whole file float plateScale_cm2deg; //@< plate scale (deg/cm) float degTriggerZone; //@< trigger area in the camera (radius, in deg.) float dtheta, dphi; //@< deviations of CT from shower axis int still_in_loop = FALSE; char Signature[20]; float *image_data; int nvar; // int hidt; struct camera cam; // structure holding the camera definition //!@' @#### Definition of variables for |getopt()|. //@' int ch, errflg = 0; //@< used by getopt /*!@' @#### Beginning of the program. We start with the main program. First we (could) make some presentation, and follows the reading of the parameters file (now from the |stdin|), the reading of the CT parameters file, and the creation of the output file, where the processed data will be stored. */ //++ // START //-- // make unbuffered output cout.setf ( ios::stdio ); // parse command line options (see reflector.h) parname[0] = '\0'; optarg = NULL; while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) ) switch (ch) { case 'f': strcpy(parname, optarg); break; case 'h': usage(); break; default : errflg++; } // show help if error if ( errflg>0 ) usage(); // make some sort of presentation present(); // read parameters file if ( strlen(parname) < 1 ) readparam(NULL); else readparam(parname); // read data from file or from STDIN? Data_From_STDIN = get_data_from_stdin(); // write all images, even those without trigger? Write_All_Images = get_write_all_images(); Write_McEvt = get_write_McEvt() ; Write_McTrig = get_write_McTrig() ; Write_RawEvt = get_write_RawEvt() ; // get filenames strcpy( inname, get_input_filename() ); strcpy( rootname, get_root_filename() ); strcpy( ct_filename, get_ct_filename() ); // get different parameters of the simulation simulateNSB = get_nsb( &meanNSB ); countIslands = get_islands_cut( &nIslandsCut ); // get selections on the parameters Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue); // log filenames information log(SIGNATURE, "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n", "Filenames", "In", inname, "ROOT", rootname, "CT", ct_filename); // log flags information log(SIGNATURE, "%s:\n\t%20s: %s\n\t%20s: %s\n", "Flags", "Data_From_STDIN", ONoff(Data_From_STDIN), "Write_All_Images", ONoff(Write_All_Images)); // log flags information log(SIGNATURE, "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n", "Root output", "Write_McEvt", ONoff(Write_McEvt), "Write_McTrig", ONoff(Write_McTrig), "Write_RawEvt", ONoff(Write_RawEvt)); // log parameters information log(SIGNATURE, "%s:\n\t%20s: %f %s\n\t%20s: %f %s\n", "Parameters", "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB), "i0 (Islands-cut)", nIslandsCut, ONoff(countIslands)); // log selections log(SIGNATURE, "%s:\n\t%20s: %s (%f:%f)\n", "Selections:", "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue); // set all random numbers seeds setall( get_seeds(0), get_seeds(1) ); // get list of showers to evt. skip nSkip = get_nskip_showers(); if (nSkip > 0) { Skip = new int[ nSkip ]; get_skip_showers( Skip ); log(SIGNATURE, "There are some showers to skip:\n"); for (i=0; iIsBatch()) { fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); // return 1; } #endif // __INTERAKTIV__ // for safety and for dimensioning image_data: count the elements in the // diagnostic data branch i=0; i++; // "n" i++; // "primary" i++; // "energy" i++; // "cored" i++; // "impact" i++; // "xcore" i++; // "ycore" i++; // "theta" i++; // "phi" i++; // "deviations" i++; // "dtheta" i++; // "dphi" i++; // "trigger" i++; // "ncphs" i++; // "maxpassthr_phe" i++; // "nphes" i++; // "nphes2" i++; // "length" i++; // "width" i++; // "dist" i++; // "xdist" i++; // "azw" i++; // "miss" i++; // "alpha" i++; // "conc2" i++; // "conc3" i++; // "conc4" i++; // "conc5" i++; // "conc6" i++; // "conc7" i++; // "conc8" i++; // "conc9" i++; // "conc10" i++; // "asymx" i++; // "asymy" i++; // "phiasym" nvar = i; image_data = new float[nvar]; // set plate scale (deg/cm) and trigger area (deg) plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381; if ( ! get_trigger_radius( °TriggerZone ) ) degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0); if ( ! get_correction( &fCorrection ) ) fCorrection = 1.0; // number of pixels for parameters anaPixels = get_ana_pixels(); anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels; // open input file if we DO read data from a file if (! Data_From_STDIN) { log( SIGNATURE, "Openning input \"rfl\" file %s\n", inname ); inputfile.open( inname ); if ( inputfile.bad() ) error( SIGNATURE, "Cannot open input file: %s\n", inname ); } // get signature, and check it // NOTE: this part repeats further down in the code; // if you change something here you probably want to change it // there as well strcpy(Signature, REFL_SIGNATURE); strcpy(sign, Signature); if ( Data_From_STDIN ) cin.read( (char *)sign, strlen(Signature)); else inputfile.read( (char *)sign, strlen(Signature)); if (strcmp(sign, Signature) != 0) { cerr << "ERROR: Signature of .rfl file is not correct\n"; cerr << '"' << sign << '"' << '\n'; cerr << "should be: " << Signature << '\n'; exit(1); } if ( Data_From_STDIN ) cin.read( (char *)sign, 1); else inputfile.read( (char *)sign, 1); // initializes flag strcpy( flag, " \0" ); // allocate space for PMTs numbers of pixels fnpix = new float [ ct_NPixels ]; fnpixclean = new float [ ct_NPixels ]; #ifdef __ROOT__ // construct MFADC #endif // __ROOT__ // moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 ); //!@' @#### Main loop. //@' //begin my version // get flag if ( Data_From_STDIN ) cin.read( flag, SIZE_OF_FLAGS ); else inputfile.read ( flag, SIZE_OF_FLAGS ); // loop over the file still_in_loop = TRUE; while ( ((! Data_From_STDIN) && (! inputfile.eof())) || (Data_From_STDIN && still_in_loop) ) { // reading .rfl files if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file log(SIGNATURE, "End of file . . .\n"); still_in_loop = FALSE; if ( Data_From_STDIN ) cin.read( (char *)sign, 1); else inputfile.read( (char *)sign, 1); } else if(!isA( flag, FLAG_START_OF_RUN )){ error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag ); } else { // found start of run nshow=0; // read next flag if ( Data_From_STDIN ) cin.read( flag, SIZE_OF_FLAGS ); else inputfile.read ( flag, SIZE_OF_FLAGS ); while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event /*!@' For the case |FLAG\_START\_OF\_EVENT|, we read each Cherenkov photon, and follow these steps: @enumerate @- Transform XY-coordinates to IJ-coordinates. @- With this, we obtain the pixel where the photon hits. @- Use the wavelength $\lambda$ and the table of QE, and calculate the estimated (third order interpolated) quantum efficiency for that photon. The photon can be rejected. @- If accepted, then add to the pixel. @endenumerate In principle, we should stop here, and use another program to 'smear' the image, to add the Night Sky Background, and to simulate the trigger logic, but we will make this program quick and dirty, and include all here. If we are reading PHE files, we jump to the point where the pixelization process already has finished. */ ++nshow; if ( fmod ( nshow, 1000. ) == 0. ) log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow); // get MCEventHeader if ( Data_From_STDIN ) cin.read( (char*)&mcevth, mcevth.mysize() ); else mcevth.read( inputfile ); // calculate core distance and impact parameter coreD = mcevth.get_core(&coreX, &coreY); // calculate impact parameter (shortest distance betwee the original // trajectory of the primary (assumed shower-axis) and the // direction where the telescope points to // // we use the following equation, given that the shower core position // is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1), // and the telescope position and orientation are (x2,y2,z2)=(0,0,0) // and (l2,m2,n2) // // | | // | x1-x2 y1-y2 z1-z2 | // | | // + | l1 m1 n1 | // - | | // | l2 m2 n2 | // | | // dist = ------------------------------------ ( > 0 ) // [ |l1 m1|2 |m1 n1|2 |n1 l1|2 ]1/2 // [ | | + | | + | | ] // [ |l2 m2| |m2 n2| |n2 l2| ] // // playing a little bit, we get this reduced for in our case: // // // dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) / // [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 - // 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2) // read the direction of the incoming shower thetashw = mcevth.get_theta(); phishw = mcevth.get_phi(); // calculate vector for shower l1 = sin(thetashw)*cos(phishw); m1 = sin(thetashw)*sin(phishw); n1 = cos(thetashw); // read the deviation of the telescope with respect to the shower mcevth.get_deviations ( &thetaCT, &phiCT ); if ( (thetaCT == 0.) && (phiCT == 0.) ) { // CT was looking to the source (both lines are parallel) // therefore, we calculate the impact parameter as the distance // between the CT axis and the core position impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. ); } else { // the shower comes off-axis // obtain with this the final direction of the CT thetaCT += thetashw; phiCT += phishw; // calculate vector for telescope l2 = sin(thetaCT)*cos(phiCT); m2 = sin(thetaCT)*sin(phiCT); n2 = cos(thetaCT); num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY); den = (SQR(l1*m2 - l2*m1) + SQR(m1*n2 - m2*n1) + SQR(n1*l2 - n2*l1)); den = sqrt(den); impactD = fabs(num)/den; // fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ", // thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2, // num, den); } // clear camera for ( i=0; i 600.0) || (wl < 290.0) ) break; // check if photon is inside outermost camera radius if(sqrt(cx*cx + cy*cy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)){ // read next CPhoton if ( Data_From_STDIN ) cin.read( flag, SIZE_OF_FLAGS ); else inputfile.read ( flag, SIZE_OF_FLAGS ); // go to beginning of loop, the photon is lost continue; } // cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl; ci = floor( (cx - cy*COS60/SIN60)/ ct_2Apot + 0.5); cj = floor( (cy/SIN60) / ct_2Apot + 0.5); ici = (int)(ci); icj = (int)(cj); iici = ici+PIX_ARRAY_HALF_SIDE; iicj = icj+PIX_ARRAY_HALF_SIDE; // is it inside the array? if ( (iici > 0) && (iici < PIX_ARRAY_SIDE) && (iicj > 0) && (iicj < PIX_ARRAY_SIDE) ) { // try to put into pixel // obtain the pixel number for this photon nPMT = (int) pixels[ici+PIX_ARRAY_HALF_SIDE][icj+PIX_ARRAY_HALF_SIDE][PIXNUM]; } else{ nPMT = -1; } // check if outside the central camera if ( (nPMT < 0) || (nPMT >= ct_NCentralPixels) ) { // check the outer pixels nPMT = -1; for(i=ct_NCentralPixels; i k) qeptr = (float **)QE[nPMT]; FindLagrange(qeptr,k,wl); // if random > quantum efficiency, reject it qe = Lagrange(qeptr,k,wl) / 100.0; // fprintf(stdout, "%f\n", qe); if ( RandomNumber > qe ) { // read next CPhoton if ( Data_From_STDIN ) cin.read( flag, SIZE_OF_FLAGS ); else inputfile.read ( flag, SIZE_OF_FLAGS ); // go to beginning of loop continue; } #endif // __QE__ //+++ // Cphoton is accepted //--- ncph_out++ ; // increase the number of Cphs. in the PMT, i.e., // increase in one unit the counter of the photons // stored in the pixel nPMT fnpix[nPMT] += 1.0; #ifdef __DETAIL_TRIGGER__ // // fill the Trigger class with this phe // // Trigger.Fill( nPMT, ( t - t_first ) ) ; // fadc.Fill( nPMT, ( t - t_first ), Trigger.Fill( nPMT, ( t - t_first ) ) ) ; #endif // __DETAIL_TRIGGER__ // read next CPhoton if ( Data_From_STDIN ) cin.read( flag, SIZE_OF_FLAGS ); else inputfile.read ( flag, SIZE_OF_FLAGS ); } // end while, i.e. found end of event if ( fmod ( nshow, 1000. ) == 0. ) log(SIGNATURE, "End of this event: in: %d cphs(+%d). out: %d cphs(+%d). .\n", ncph_in, ntcph_in, ncph_out, ntcph_out); // show number of photons //cout << ncph_in << " photons read . . . " << endl << flush; // skip it ? for ( i=0; i Select_Energy_ue )) { log(SIGNATURE, "select_energy: shower rejected.\n"); continue; } } #ifdef __NSB__ //!@' @#### NSB (Night Sky Background) simulation. //@' //+++ // NSB simulation //--- // add NSB "noise" // TO DO: make meanNSB an array and read the contents from a file! // if ( simulateNSB ) // for ( i=0; iSetZeroLevel( Lev0 = (Short_t) Trigger.Diskriminate() ) ; Lev1 = Lev2 = 0 ; // // Start the First Level Trigger simulation // if ( Lev0 > 0 ) { McTrig->SetFirstLevel ( Lev1 = Trigger.FirstLevel() ) ; } #endif // __DETAIL_TRIGGER__ #ifdef __ROOT__ // // Fill the header of this event // Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ; // // fill the MMcEvt with all information // McEvt->Fill( (UShort_t) mcevth.get_primary() , mcevth.get_energy(), mcevth.get_theta(), mcevth.get_phi(), mcevth.get_core(), mcevth.get_coreX(), mcevth.get_coreY(), impactD, ulli, ulli, (UShort_t) ncph_in, ulli, (UShort_t) ncph_out ) ; // // write it out to the file outfile // EvtTree.Fill() ; #endif // __ROOT__ // // if a first level trigger occurred, then // 1. do some other stuff (not implemented) // 2. start the gui tool if ( Lev1 > 0 ) { // fadc.Scan( Trigger.GetFirstLevelTime(0) ) ; #ifdef __INTERAKTIV__ Trigger.ShowSignal(McEvt) ; #endif } #ifdef __ROOT__ // clear all Evt->Clear() ; McEvt->Clear() ; McTrig->Clear() ; #endif // __ROOT__ //++++++++++++++++++++++++++++++++++++++++++++++++++ // at this point we have a camera full of // ph.e.s // we should first apply the trigger condition, // and if there's trigger, then clean the image, // calculate the islands statistics and the // other parameters of the image (Hillas' parameters // and so on). //-------------------------------------------------- #ifdef __DEBUG__ printf("\n"); for ( ici=0; ici -1 ) { if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) { printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow, (int)pixels[ici][icj][PIXNUM], pixels[ici][icj][PIXX], pixels[ici][icj][PIXY], fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj); } } } } for (i=0; i (0:CT1 Ķ 1:MAGIC) // get focal distance sscanf(line, "%s %d", token, &ct_Type); log( "read_ct_file", ": %s\n", ((ct_Type==0) ? "CT1" : "MAGIC") ); break; case focal_distance: // [cm] // get focal distance sscanf(line, "%s %f", token, &ct_Focal_mean); log( "read_ct_file", ": %f cm\n", ct_Focal_mean ); break; case focal_std: // s(focal distance) [cm] // get focal distance sscanf(line, "%s %f", token, &ct_Focal_std); log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std ); break; case point_spread: // [cm] // get point spread sscanf(line, "%s %f", token, &ct_PSpread_mean); log( "read_ct_file", ": %f cm\n", ct_PSpread_mean ); break; case point_std: // s(point spread) [cm] // get point spread sscanf(line, "%s %f", token, &ct_PSpread_std); log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std ); break; case adjustment_dev: // s(adjustment_dev) [cm] // get point spread sscanf(line, "%s %f", token, &ct_Adjustment_std); log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std ); break; case black_spot: // radius of the black spot in the center [cm] // get black spot radius sscanf(line, "%s %f", token, &ct_BlackSpot_rad); log( "read_ct_file", "Radius of the black spots: %f cm\n", ct_BlackSpot_rad); break; case r_mirror: // radius of the mirrors [cm] // get radius of mirror sscanf(line, "%s %f", token, &ct_RMirror); log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror ); break; case n_mirrors: // number of mirrors // get the name of the output_file from the line sscanf(line, "%s %d", token, &ct_NMirrors); log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors ); break; case camera_width: // camera width [cm] // get the name of the ct_file from the line sscanf(line, "%s %f", token, &ct_CameraWidth); log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth ); break; case n_pixels: // number of pixels // get the name of the output_file from the line sscanf(line, "%s %d", token, &ct_NPixels); log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels ); break; case n_centralpixels: // number of central pixels // get the name of the output_file from the line sscanf(line, "%s %d", token, &ct_NCentralPixels); log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels ); break; case n_gappixels: // number of gap pixels // get the name of the output_file from the line sscanf(line, "%s %d", token, &ct_NGapPixels); log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels ); break; case pixel_width: // pixel width [cm] // get the name of the ct_file from the line sscanf(line, "%s %f", token, &ct_PixelWidth); ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0)); ct_PixelWidth_corner_2_corner_half = ct_PixelWidth_corner_2_corner * 0.50; ct_Apot = ct_PixelWidth / 2; ct_2Apot = ct_Apot * 2.0; log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth ); break; case define_mirrors: // read table with the parameters of the mirrors log( "read_ct_file", "Table of mirrors data:\n" ); // check whether the number of mirrors was already set if ( ct_NMirrors == 0 ) error( "read_ct_file", "NMirrors was not set.\n" ); // allocate memory for paths list log( "read_ct_file", "Allocating memory for ct_data\n" ); ct_data = new float*[ct_NMirrors]; for (i=0; i> ct_data[i][j]; break; } // switch ( i ) } // end while // end log( "read_ct_file", "done.\n" ); return; } //!@} //!----------------------------------------------------------- // @name read_pixels // // @desc read pixels data // // @date Fri Mar 12 16:33:34 MET 1999 //------------------------------------------------------------ // @function //!@{ void read_pixels(struct camera *pcam) { ifstream qefile; char line[LINE_MAX_LENGTH]; int n, i, j, k; float qe; //------------------------------------------------------------ // first, pixels' coordinates pcam->inumpixels = ct_NPixels; pcam->inumcentralpixels = ct_NCentralPixels; pcam->inumgappixels = ct_NGapPixels; pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels; pcam->dpixdiameter_cm = ct_PixelWidth; // initialize pixel numbers for ( i=0; idi[k]; j = (int) pcam->dj[k]; pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k; pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k]; pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k]; pixary[k][0] = pcam->dxc[k]; pixary[k][1] = pcam->dyc[k]; } // calculate tables of neighbours #ifdef __DEBUG__ for ( n=0 ; n -1) && ((j-1) < pointsQE) && ((j-1) > -1) ) { QE[i-1][0][j-1] = QElambda[j-1]; QE[i-1][1][j-1] = qe; } } // close file qefile.close(); // end log("read_pixels", "Done.\n"); } //!@} //!----------------------------------------------------------- // @name pixels_are_neig // // @desc check whether two pixels are neighbours // // @var pix1 Number of the first pixel // @var pix2 Number of the second pixel // @return TRUE: both pixels are neighbours; FALSE: oth. // // @date Wed Sep 9 17:58:37 MET DST 1998 //------------------------------------------------------------ // @function //!@{ int pixels_are_neig(int pix1, int pix2) { if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) + SQR( pixary[pix1][1] - pixary[pix2][1] ) ) > ct_PixelWidth_corner_2_corner ) return ( FALSE ); else return ( TRUE ); } //!@} //!----------------------------------------------------------- // @name igen_pixel_coordinates // // @desc generate the pixel center coordinates // // @var *pcam structure camera containing all the // camera information // @return total number of pixels // // DP // // @date Thu Oct 14 10:41:03 CEST 1999 //------------------------------------------------------------ // @function //!@{ /******** igen_pixel_coordinates() *********************************/ int igen_pixel_coordinates(struct camera *pcam) { /* generate pixel coordinates, return value is number of pixels */ int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment; float fsegment_fract; double dtsize; double dhsize; double dpsize; double dxfirst_pix; double dyfirst_pix; double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6; double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6; double dstartx, dstarty; /* for the gap pixels and outer pixels */ int j, nrow; dpsize = pcam->dpixdiameter_cm; dtsize = dpsize * sqrt(3.) / 2.; dhsize = dpsize / 2.; /* Loop over central pixels to generate co-ordinates */ for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){ /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */ pcam->dpixsizefactor[ipixno] = 1.; in = 0; i = 0; itot_inside_ring = 0; iring_no = 0; /* Calculate the number of pixels out to and including the ring containing pixel number */ /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */ while (itot_inside_ring == 0){ iN = 3*(i*(i+1)) + 1; if (ipixno <= iN){ iring_no = i; itot_inside_ring = iN; } i++; } /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */ ipix_in_ring = 0; for (i = 0; i < iring_no; ++i){ ipix_in_ring = ipix_in_ring + 6; } /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */ /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */ /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */ isegment = 0; fsegment_fract = 0.; if (iring_no > 0) { isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */ fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ; } /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */ /* pixel width (flat to flat) dpsize. */ dxfirst_pix = dpsize*iring_no; dyfirst_pix = 0.; /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */ ddxseg1 = - dhsize*iring_no; ddyseg1 = dtsize*iring_no; ddxseg2 = -dpsize*iring_no; ddyseg2 = 0.; ddxseg3 = ddxseg1; ddyseg3 = -ddyseg1; ddxseg4 = -ddxseg1; ddyseg4 = -ddyseg1; ddxseg5 = -ddxseg2; ddyseg5 = 0.; ddxseg6 = -ddxseg1; ddyseg6 = ddyseg1; /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */ /* anti-clockwise around the ring by adding the segment to segment vectors. */ switch (isegment) { case 0: pcam->dxc[ipixno-1] = 0.; pcam->dyc[ipixno-1] = 0.; case 1: pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract; pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract; break; case 2: pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract; pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.; break; case 3: pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract; pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract; break; case 4: pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract; pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract; break; case 5: pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract; pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.; break; case 6: pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract; pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract; break; default: fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno); return(0); } /* end switch */ } /* end for */ dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize; dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize; if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */ j = pcam->inumcentralpixels; for(i=0; iinumgappixels; i=i+6){ pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize; pcam->dyc[j + i ] = dstarty; pcam->dpixsizefactor[j + i] = 1.; pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.; pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1]; pcam->dpixsizefactor[j + i + 1] = 1.; pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1]; pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1]; pcam->dpixsizefactor[j + i+ 2] = 1.; pcam->dxc[j + i + 3] = - pcam->dxc[j + i]; pcam->dyc[j + i + 3] = dstarty; pcam->dpixsizefactor[j + i+ 3] = 1.; pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2]; pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2]; pcam->dpixsizefactor[j + i+ 4] = 1.; pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1]; pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1]; pcam->dpixsizefactor[j + i + 5] = 1.; } /* end for */ } /* end if */ /* generate positions of the outer pixels */ if( pcam->inumbigpixels > 0 ){ j = pcam->inumcentralpixels + pcam->inumgappixels; for(i=0; iinumbigpixels; i++){ pcam->dpixsizefactor[j + i] = 2.; } in = 0; nrow = (int) ceil(dstartx / 2. / dpsize); while(in < pcam->inumbigpixels){ pcam->dxc[j + in] = dstartx + dpsize; pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.); pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.; pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.); pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in]; pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in]; pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in]; pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in]; pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow]; pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow]; pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in]; pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in]; for(i=1; idxc[j + in + i] = pcam->dxc[j + in] - i * dpsize; pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.); pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize; pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow]; pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i]; pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i]; pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i]; pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i]; pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow]; pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow]; pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i]; pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i]; } in = in + 6 * nrow; dstartx = dstartx + 2. * dpsize; nrow = nrow + 1; } /* end while */ } /* end if */ /* generate the ij coordinates */ for(i=0; iinumpixels; i++){ pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize; pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60; // fprintf(stdout, "%d %f %f %f %f %f\n", // i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i], // pcam->dpixsizefactor[i]); } return(pcam->inumpixels); } //!@} //!----------------------------------------------------------- // @name bpoint_is_in_pix // // @desc check if a point (x,y) in camera coordinates is inside a given pixel // // @var *pcam structure camera containing all the // camera information // @var dx, dy point coordinates in centimeters // @var ipixnum pixel number (starting at 0) // @return TRUE if the point is inside the pixel, FALSE otherwise // // DP // // @date Thu Oct 14 16:59:04 CEST 1999 //------------------------------------------------------------ // @function //!@{ /******** bpoint_is_in_pix() ***************************************/ int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){ /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */ /* the pixel is assumed to be a "closed set" */ double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */ double c, xx, yy; /* auxiliary variable */ b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum]; a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum]; c = 1. - 1./sqrt(3.); if((ipixnum < 0)||(ipixnum >= pcam->inumpixels)){ fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum); fprintf(stderr, "Exiting.\n"); exit(203); } xx = dx - pcam->dxc[ipixnum]; yy = dy - pcam->dyc[ipixnum]; if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) || ((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){ return(TRUE); /* point is inside */ } else{ return(FALSE); /* point is outside */ } } //!@} //------------------------------------------------------------ // @name dist_r_P // // @desc distance straight line r - point P // // @date Sat Jun 27 05:58:56 MET DST 1998 // @function @code //------------------------------------------------------------ // dist_r_P // // distance straight line r - point P //------------------------------------------------------------ float dist_r_P(float a, float b, float c, float l, float m, float n, float x, float y, float z) { return ( sqrt((SQR((a-x)*m-(b-y)*l) + SQR((b-y)*n-(c-z)*m) + SQR((c-z)*l-(a-x)*n))/ (SQR(l)+SQR(m)+SQR(n)) ) ); } // @endcode //=------------------------------------------------------------ //!@subsection Log of this file. //!@{ // // $Log: not supported by cvs2svn $ // Revision 1.1.1.1 1999/11/05 11:59:31 harald // This the starting point for CVS controlled further developments of the // camera program. The program was originally written by Jose Carlos. // But here you can find a "rootified" version to the program. This means // that there is no hbook stuff in it now. Also the output of the // program changed to the MagicRawDataFormat. // // The "rootification" was done by Dirk Petry and Harald Kornmayer. // // In the following you can see the README file of that version: // // ================================================== // // Fri Oct 22 1999 D.P. // // The MAGIC Monte Carlo System // // Camera Simulation Programme // --------------------------- // // 1) Description // // This version is the result of the fusion of H.K.'s // root_camera which is described below (section 2) // and another version by D.P. which had a few additional // useful features. // // The version compiles under Linux with ROOT 2.22 installed // (variable ROOTSYS has to be set). // // Compile as before simply using "make" in the root_camera // directory. // // All features of H.K.'s root_camera were retained. // // Additional features of this version are: // // a) HBOOK is no longer used and all references are removed. // // b) Instead of HBOOK, the user is given now the possibility of // having Diagnostic data in ROOT format as a complement // to the ROOT Raw data. // // This data is written to the file which is determined by // the new input parameter "diag_file" in the camera parameter // file. // // All source code file belonging to this part have filenames // starting with "MDiag". // // The user can read the output file using the following commands // in an interactive ROOT session: // // root [0] .L MDiag.so // root [1] new TFile("diag.root"); // root [2] new TTreeViewer("T"); // // This brings up a viewer from which all variables of the // TTree can be accessed and histogrammed. This example // assumes that you have named the file "diag.root", that // you are using ROOT version 2.22 or later and that you have // the shared object library "MDiag.so" which is produced // by the Makefile along with the executable "camera". // // ! The contents of the so-called diag file is not yet fixed. // ! At the moment it is what J.C.G. used to put into the HBOOK // ! ntuple. In future versions the moments calculation can be // ! removed and the parameter list be modified correspondingly. // // c) Now concatenated reflector files can be read. This is useful // if you have run the reflector with different parameters but // you want to continue the analysis with all reflector data // going into ONE ROOT outputfile. // // The previous camera version contained a bug which made reading // of two or more concatenated reflector files impossible. // // d) The reflector output format was changed. It is now version // 0.4 . // The change solely consists in a shortening of the flag // definition in the file // // include-MC/MCCphoton.hxx // // ! IF YOU WANT TO READ REFLECTOR FORMAT 0.3, you can easily // ! do so by recompiling camera with the previous version of // ! include-MC/MCCphoton.hxx. // // The change was necessary for saving space and better // debugging. From now on, this format can be frozen. // // ! For producing reflector output in the new format, you // ! of course have to recompile your reflector with the // ! new include-MC/MCCphoton.hxx . // // e) A first version of the pixelization with the larger // outer pixels is implemented. THIS IS NOT YET FULLY // TESTED, but first rough tests show that it works // at least to a good approximation. // // The present version implements the camera outline // with 18 "gap-pixels" and 595 pixels in total as // shown in // // http://sarastro.ifae.es/internal/home/hardware/camera/numbering.ps // // This change involved // // (i) The file pixels.dat is no longer needed. Instead // the coordinates are generated by the program itself // (takes maybe 1 second). In the file // // pixel-coords.txt // // in the same directory as this README, you find a list // of the coordinates generated by this new routine. It // has the format // // number i j x y size-factor // // where i and j are J.C.G.'s so called biaxis hexagonal // coordinates (for internal use) and x and y are the // coordinates of the pixel centers in the standard camera // coordinate system in units of centimeters. The value // of "size-factor" determines the linear size of the pixel // relative to the central pixels. // // (ii) The magic.def file has two additional parameters // which give the number of central pixels and the // number of gap pixels // // (iii) In camera.h and camera.cxx several changes were // necessary, among them the introduction of several // new functions // // The newly suggested outline with asymmetric Winston cones // will be implemented in a later version. // // f) phe files can no longer be read since this contradicts // our philosophy that the analysis should be done with other // programs like e.g. EVITA and not with "camera" itself. // This possibility was removed. // // g) ROOT is no longer invoked with an interactive interface. // In this way, camera can better be run as a batch program and // it uses less memory. // // h) small changes concerning the variable "t_chan" were necessary in // order to avoid segmentation faults: The variable is used as an // index and it went sometimes outside the limits when camera // was reading proton data. This is because the reflector files // don't contain the photons in a chronological order and also // the timespread can be considerably longer that the foreseen // digitisation timespan. Please see the source code of camera.cxx // round about line 1090. // // j) several unused variables were removed, a few warning messages // occur when you compile camera.cxx but these can be ignored at // the moment. // // In general the program is of course not finished. It still needs // debugging, proper trigger simulation, simulation of the asymmetric // version of the outer pixels, proper NSB simulation, adaption of // the diag "ntuple" contents to our need and others small improvements. // // In the directory rfl-files there is now a file in reflector format 0.4 // containing a single event produced by the starfiled adder. It has // a duration of 30 ns and represents the region around the Crab Nebula. // Using the enclosed input parameter file, camera should process this // file without problems. // // 2) The README for the previous version of root_camera // // README for a preliminary version of the // root_camera program. // // root_camera is based on the program "camera"of Jose Carlos // Gonzalez. It was changed in the way that only the pixelisation // and the distibution of the phe to the FADCs works in a // first version. // // Using the #undef command most possibilities of the orignal // program are switched of. // // The new parts are signed by // // - ROOT or __ROOT__ // nearly all important codelines for ROOT output are enclosed // in structures like // #ifdef __ROOT__ // // code // // #endif // __ROOT__ // // In same case the new lines are signed by a comment with the word // ROOT in it. // // For timing of the pulse some variable names are changed. // (t0, t1, t --> t_ini, t_fin, t_1st, t_chan,...) // Look also for this changes. // // For the new root-file is also a change in readparm-files // // // - __DETAIL_TRIGGER__ // // This is for the implementation of the current work on trigger // studies. Because the class MTrigger is not well documented it // isnīt a part of this tar file. Only a dummy File exists. // // // // With all files in the archive, the root_camera program should run. // // A reflector file is in the directory rfl-files // // ================================================== // // From now on, use CVS for development!!!! // // // // Revision 1.3 1999/10/22 15:01:28 petry // version sent to H.K. and N.M. on Fri Oct 22 1999 // // Revision 1.2 1999/10/22 09:44:23 petry // first synthesized version which compiles and runs without crashing; // // Revision 1.1.1.1 1999/10/21 16:35:10 petry // first synthesised version // // Revision 1.13 1999/03/15 14:59:05 gonzalez // camera-1_1 // // Revision 1.12 1999/03/02 09:56:10 gonzalez // *** empty log message *** // // //!@} //=EOF