Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 364) +++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 365) @@ -21,7 +21,7 @@ // // $RCSfile: camera.cxx,v $ -// $Revision: 1.4 $ +// $Revision: 1.5 $ // $Author: petry $ -// $Date: 2000-01-25 08:36:23 $ +// $Date: 2000-02-18 17:40:35 $ // //////////////////////////////////////////////////////////////////////// @@ -54,6 +54,4 @@ #include "MDiag.h" - -#include "MTrigger.hxx" #include "MRawEvt.h" @@ -86,46 +84,4 @@ #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 <-- -#define __DETAIL_TRIGGER__ -#undef __DETAIL_TRIGGER__ - -// flag for implementation of TRIGGER (default: ON ) -#undef __TRIGGER__ -#define __TRIGGER__ - -// flag for implementation of Tail-Cut (default: ON ) -#undef __TAILCUT__ -#define __TAILCUT__ - -// flag for calculation of islands stat. (default: ON ) -#define __ISLANDS__ -#undef __ISLANDS__ - -// flag for calculation of image parameters (default: ON ) -#undef __MOMENTS__ -#define __MOMENTS__ - -// flag for ROOT (default: ON ) -#undef __ROOT__ -#define __ROOT__ - //!@} @@ -259,6 +215,4 @@ // 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 @@ -276,7 +230,4 @@ //@: contents of the pixels (ph.e.) after cleanning static float *fnpixclean; - -//@: contents of the sum of all ph.e. in one timeslice of 1 ns -static float *fnslicesum ; @@ -364,30 +315,4 @@ @"*/ -//!@{ -// 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]; -//!@} - /*!@" @@ -430,5 +355,5 @@ char inname[256]; //@< input file name - char outname[256]; //@< output file name + char starfieldname[256]; //@< starfield input file name char datname[256]; //@< data (ASCII) output file name char diagname[256]; //@< diagnistic output file (ROOT format) @@ -437,14 +362,16 @@ 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 - ofstream outputfile; //@< stream for the output file + FILE *inputfile; //@< stream for the input file ofstream datafile; //@< stream for the data file MCEventHeader mcevth; //@< Event Header class (MC) - MCCphoton cphoton; //@< Cherenkov Photon class (MC) + + Photoelectron *photoe = NULL; //@< array of the photoelectrons of one event + int inumphe; //@< number of photoelectrons in an event + + float arrtmin_ns; //@ arrival time of the first photoelectron + float arrtmax_ns; //@ arrival time of the last photoelectron float thetaCT, phiCT; //@< parameters of a given shower @@ -458,26 +385,25 @@ int nshow=0; //@< partial number of shower in a given run int ntshow=0; //@< total number of showers - int ncph=0; //@< partial number of shower in a given run - int ntcph=0; //@< total number of showers - - int i, ii, j, k; //@< simple counters - - float t_ini; //@< time of the first Cphoton read in - float t; //@< time for a single photon - 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 ncph=0; //@< partial number of photons in a given run + int ntcph=0; //@< total number of photons + + int i, j, k; //@< simple counters int simulateNSB; //@< Will we simulate NSB? - float meanNSB; //@< NSB mean value (per pixel) + float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel) + float diffnsb_phepns[iMAXNUMPIX]; //@< diffuse NSB values for each pixel derived + //@< from meanNSB + float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS]; //@< non-diffuse nsb + //@< photoelectron rates + float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband + EXTWAVEBAND1, + EXTWAVEBAND2, + EXTWAVEBAND3, + EXTWAVEBAND4, + EXTWAVEBAND5 + }; + float baseline_mv[iMAXNUMPIX]; //@< The baseline (mV) caused by the NSB; to be subtracted + //@< in order to simulate the preamps' AC coupling + float qThreshold; //@< Threshold value float qTailCut; //@< Tail Cut value @@ -488,14 +414,7 @@ float fCorrection; //@< Factor to apply to pixel values (def. 1.) - float q0; //@< trigger threshold ( intermediate variable ) - float maxcharge; //@< maximum charge in pixels - int noverq0, novq0; //@< number of pixels above threshold - int ngrpq0, mxgrp; //@< number of pixels in a group - int trigger; //@< trigger flag int itrigger; //@< index of pixel fired int ntrigger = 0; //@< number of triggers in the whole file - unsigned char triggerBits; //@< byte for trigger condition check (MAGIC) - int bit; //@< intermediate variable float plateScale_cm2deg; //@< plate scale (deg/cm) @@ -506,6 +425,4 @@ int still_in_loop = FALSE; - char Signature[20]; - float *image_data; int nvar, hidt; @@ -513,9 +430,4 @@ struct camera cam; // structure holding the camera definition -#ifdef __DETAIL_TRIGGER__ - - MTrigger Trigger ; //@< A instance of the Class MTrigger - -#endif // __DETAIL_TRIGGER__ //!@' @#### Definition of variables for |getopt()|. @@ -592,5 +504,5 @@ strcpy( inname, get_input_filename() ); - strcpy( outname, get_output_filename() ); + strcpy( starfieldname, get_starfield_filename() ); strcpy( datname, get_data_filename() ); strcpy( diagname, get_diag_filename() ); @@ -615,9 +527,10 @@ "Filenames", "In", inname, - "Out", outname, + "Stars", starfieldname, + "CT", ct_filename, "Data", datname, "Diag", diagname, - "ROOT", rootname, - "CT", ct_filename); + "ROOT", rootname + ); @@ -669,8 +582,12 @@ read_ct_file(); - // read pixels data + // read camera setup read_pixels(&cam); + // allocate memory for the photoelectrons + + photoe = new Photoelectron[iMAXNUMPHE]; + // initialise ROOT @@ -682,5 +599,4 @@ hfile = new TFile( diagname,"RECREATE", "MAGIC Telescope MC diagnostic data"); - // Create the ROOT Tree for the diagnostic data @@ -698,5 +614,5 @@ tree->Branch("event", "MDiagEventobject", &event, bsize, split); -#ifdef __ROOT__ + // Prepare the raw data output MRawEvt *Evt = new MRawEvt() ; @@ -706,9 +622,7 @@ // - TFile outfile ( rootname , "RECREATE" ) ; - - // + TFile outfile ( rootname , "RECREATE" ); + // create a Tree for the Event data stream - // TTree EvtTree("EvtTree","Events of Run"); @@ -722,9 +636,4 @@ &McEvt, bsize, split); - - float flli = 0. ; - unsigned short ulli = 0 ; - -#endif // __ROOT__ // for safety and for dimensioning image_data: count the elements in the @@ -787,59 +696,72 @@ 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() ) + // prepare the NSB simulation + + if( simulateNSB ){ + + // + // Calculate the non-diffuse NSB photoelectron rates + // + + k = produce_nsbrates( starfieldname, + &cam, + photoe, // only a dummy here + nsbrate_phepns ); + if (k != 0){ + cout << "Error when reading starfield... \nExiting.\n"; + exit(1); + } + +// for(i=0; i (2 * SLICES)){ - log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n", - t_chan, (2 * SLICES), (2 * SLICES)); - t_chan = (2 * SLICES); - } - else if(t_chan < 0){ - log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n", - t_chan, 0, 0); - t_chan = 0; - } - - // Pixelization - - cx = cphoton.get_x(); - cy = cphoton.get_y(); - - // get wavelength - - last_wl = wl; - wl = cphoton.get_wl(); - - // check if photon has valid wavelength and is inside outermost camera radius - - if( (wl > 800.0) || (wl < 290.0) || - (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; - - nPMT = -1; - - for(i=0; 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 - //--- - - // 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 __ROOT__ - fnslicesum[t_chan] += 1.0 ; - slices[nPMT][t_chan] += 1.0 ; -#endif // __ROOT__ - -#ifdef __DETAIL_TRIGGER__ - // - // fill the Trigger class with this phe - // - // - - Trigger.Fill( nPMT, ( t - t_ini ) ) ; -#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 - + log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n", ncph, ntcph); - - // show number of photons - - //cout << ncph << " photons read . . . " << endl << flush; - + + ntcph += ncph; + // skip it ? @@ -1245,32 +943,5 @@ } - /*!@' - - After reading all the Cherenkov photons for a given event, - we have in the table of number of photons for each pixel - only the 'raw' amount of Cherenkov photons @$n_p@$. Now, we - should take this number as the mean value of the - distribution of photons in that pixel @$p@$, following a - Poisson distribution. - - @[ n_p \equiv \mu_p @] - - and with this number the amount of light coming from the - shower is calculated @$\hat{n}_p@$. - - Then, we calculate the amount of Night Sky Background we - must introduce in that pixel @$p@$. We calculate this using - again a Poisson distribution with mean @$\mu_\mathrm{NSB}@$ - (defined in the |camera.h| file). The value of - @$\mu_\mathrm{NSB}@$ is obtained from measurements. With this - value, the amount of photons @$\hat{n}_\mathrm{NSB}@$ coming - from the Night Sky Background is calculated. - - Finally, the amount of photons for that pixels is: - @[ \hat{n}_p^\mathrm{final} = \hat{n}_p + \hat{n}_\mathrm{NSB} @] - - */ - - // after reading all the photons, our camera is filled + // energy cut if ( Select_Energy ) { @@ -1281,375 +952,47 @@ } } - -#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; iFillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ; - - // now put out the data of interest - // - // 1. -> look for the first slice with signal - // - - for ( i=0; i<(2 * SLICES); i++ ) - if ( fnslicesum[i] > 0. ) - break ; - - startchan = i ; - - // - // copy the slices out of the big array - // - // put the first slice with signal to slice 4 - // - - for (i=0; i 0 ) { - - for ( ii=0; ii < 15; ii++ ) { - trans [ii] = slices2[i][ii] ; - } - - Evt->FillPixel ( (UShort_t) i , trans ) ; - + + if(simulateNSB){ + + k = produce_nsb_phes( &arrtmin_ns, // will be changed by the function! + &arrtmax_ns, // will be changed by the function! + thetaCT, + &cam, + nsbrate_phepns, + diffnsb_phepns, + ext, + fnpix, // will be changed by the function! + photoe, // will be changed by the function! + &inumphe, // important for later: the size of photoe[] + baseline_mv // will be generated by the function + ); + + if( k != 0 ){ // non-zero returnvalue means error + cout << "Exiting.\n"; + exit(1); } - } - - // - // - // - - 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(), - flli, - ulli, ulli, ulli, ulli, ulli ) ; - - // - // write it out to the file outfile - // - - EvtTree.Fill() ; - - // clear all - - Evt->Clear() ; - McEvt->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 __TRIGGER__ - - /*!@' - - @#### Trigger logic simulation. - - In the following block we look at the pixel contents, looking - for pixels fulfilling the trigger condition. This condition, - in this current version of the program, is the following: - - @itemize - - @- |CT1|: Two neighbour pixels with charge above the threshold - @$q_0@$. For the old CT1 data, however, the trigger condition - was 'any two pixels with charge above the threshold @$q_0@$'. - - @- |MAGIC|: A 'closed-packet' of four neighbour pixels, each - of them with charge above the threshold @$q_0@$. - - @enditemize - - In the following figure you can find a sort of description - about the meanning of 'closed-packet'. - - @F - - \begin{figure}[htbp] - \begin{center} - \includegraphics{closepck.eps} - \caption{Meanning of the expression ``{\it close-packet}''} - \label{fig:closepacket} - \end{center} - \end{figure} - - @F - - */ - - //++ - // TRIGGER Condition - //-- - - //@ If the input parameter "threshold" is 0 we find the maximum - //@ trigger threshold this event can pass - - for(k=0; ( qThreshold == 0 ? (k <= iMAX_THRESHOLD_PHE) : (k == 1) ); k++){ - - // is there trigger? - - noverq0 = 0; - q0 = ( qThreshold == 0. ? (float) k : qThreshold ); - trigger = FALSE; - mxgrp = 0; - maxcharge = 0.0; - - // Warning! NOT all the camera is able to give trigger - // only up to 'degTrigger' degrees - - for ( i=0 ; (i degTriggerZone) - continue; - - // 'ngrpq0' is the number of neighbours of pixel i with q > q0 - - ngrpq0 = 0; - - // look at each pixel in the neighborhood, and count - // those above threshold q0 - - for ( j=0 ; j-1 ; ++j ) - if ( fnpix[pixneig[i][j]] > q0 ) - ++ngrpq0; - - // check whether we have trigger - - if ( ct_Type == 0 ) { - - //++ >>>>> CT1 <<<<< - -#ifdef __CT1_NO_NEIGHBOURS__ - - if ( noverq0 > 1 ) - trigger = TRUE; - -#else - - if ( ngrpq0 > 0 ) - trigger = TRUE; - -#endif - - //-- >>>>> CT1 <<<<< - - } else { - - //++ >>>>> MAGIC <<<<< - - // (at least 4 packed with at least q0 phes) - - // there are 3 cases - // 1. zero, one or two neighbours have enough charge: no trigger - // 2. five or six neighbours with enough charge: trigger? sure!! - // 3. three or four neighbours with q > q0 : we must look - // for 'closeness'. - - switch ( ngrpq0 ) { - - case 0: - case 1: - case 2: - - trigger = FALSE; - break; - - case 3: - case 4: - - // if reaches this line, it means three or four neighbours - // around the central pixel - - triggerBits = 1; - - for ( j=0 ; j-1; ++j ) { - - if ( fnpix[pixneig[i][j]] > q0 ) { - - if ( pixary[pixneig[i][j]][0] > pixary[i][0] ) { - - if ( nint(pixary[pixneig[i][j]][1]*10.0) > - nint(pixary[i][1]*10.0) ) - bit = 2; - else if ( nint(pixary[pixneig[i][j]][1]*10.0) < - nint(pixary[i][1]*10.0) ) - bit = 6; - else - bit = 1; - - } else { - - if ( nint(pixary[pixneig[i][j]][1]*10.0) > - nint(pixary[i][1]*10.0) ) - bit = 3; - else if ( nint(pixary[pixneig[i][j]][1]*10.0) < - nint(pixary[i][1]*10.0) ) - bit = 5; - else - bit = 4; - - } - - triggerBits |= (1< 6 !!! Exiting.\n\n"); - break; - - } // switch (ngrpq0) - - } // ct_Type - - } // for each pixel i - - if ( trigger == FALSE ) { - break; - } // end if - - } //end for each threshold - maxtrigthr_phe = (float) k-1; // i.e. maxtrigthr_phe < 0. means that - // the event doesn't even pass threshold 0. - // maxtrigthr_phe >= 0 means, the event passes some threshold - // or (in case the input parameter "threshold" was > 0), the event - // passes the threshold given by the input parameter. - if ( maxtrigthr_phe >= 0. ) { - trigger = TRUE; - } - - - novq0 = noverq0; + + }// end if(simulateNSB) ... + + +// cout << arrtmin_ns << " " << arrtmax_ns << "\n"; +// for(i=0; i 0 : Make Tail-Cut - // Tail-Cut < 0 : Make Tail-Cut with t_0 = Sqrt[ maximum ] - - if (qTailCut > 0.0) { - - for ( i=0; iClose(); - -#ifdef __DETAIL_TRIGGER__ - // Output of Trigger statistics - // - - Trigger.PrintStat() ; -#endif // __DETAIL_TRIGGER__ // program finished @@ -2138,9 +1392,9 @@ // Display the name of the function that called error - fprintf(stderr, "ERROR in %s: ", funct); + fprintf(stdout, "ERROR in %s: ", funct); // Display the remainder of the message va_start(args, fmt); - vfprintf(stderr, fmt, args); + vfprintf(stdout, fmt, args); va_end(args); @@ -2425,5 +1679,5 @@ ifstream qefile; char line[LINE_MAX_LENGTH]; - int n, i, j, k; + int n, i, j, icount; float qe; @@ -2438,8 +1692,4 @@ // 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 @@ -2519,5 +1754,5 @@ // try to open the file - log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE); + log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE); qefile.open( QE_FILE ); @@ -2533,4 +1768,5 @@ i=-1; + icount = 0; while ( ! qefile.eof() ) { @@ -2577,5 +1813,6 @@ // get the values (num-pixel, num-datapoint, QE-value) - sscanf(line, "%d %d %f", &i, &j, &qe); + if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 ) + break; if ( ((i-1) < ct_NPixels) && ((i-1) > -1) && @@ -2585,9 +1822,30 @@ } + if ( i > ct_NPixels) break; + + icount++; + } + if(icount/pointsQE < ct_NPixels){ + error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n", + icount/pointsQE, ct_NPixels ); + } + // close file qefile.close(); + + // test QE + + for(icount=0; icount< ct_NPixels; icount++){ + for(i=0; i 1000. || + QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){ + error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n", + icount, i, QE[icount][0][i], QE[icount][1][i] ); + } + } + } // end @@ -2875,16 +2133,4 @@ } /* 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); @@ -2968,4 +2214,560 @@ ); } + +//------------------------------------------------------------ +// @name check_reflector_file +// +// @desc check if a given reflector file has the right signature +// @desc return TRUE or FALSE +// +// @date Mon Feb 14 16:44:21 CET 2000 +// @function @code +//------------------------------------------------------------ + +int check_reflector_file(FILE *infile){ + + char Signature[20]; // auxiliary variable + char sign[20]; // auxiliary variable + + strcpy(Signature, REFL_SIGNATURE); + + strcpy(sign, Signature); + + fread( (char *)sign, strlen(Signature), 1, infile); + + if (strcmp(sign, Signature) != 0) { + cout << "ERROR: Signature of .rfl file is not correct\n"; + cout << '"' << sign << '"' << '\n'; + cout << "should be: " << Signature << '\n'; + return(FALSE); + } + + fread( (char *)sign, 1, 1, infile); + + return(TRUE); + +} + +//------------------------------------------------------------ +// @name lin_interpol +// +// @desc interpolate linearly between two points returning the +// @desc y-value of the result +// +// @date Thu Feb 17 11:31:32 CET 2000 +// @function @code +//------------------------------------------------------------ + +float lin_interpol( float x1, float y1, float x2, float y2, float x){ + + if( (x2 - x1)<=0. ){ // avoid division by zero, return average + cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n"; + return((y1+y2)/2.); + } + else{ // note: the check whether x1 < x < x2 is omitted for speed reasons + return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) ); + } +} + + +//------------------------------------------------------------ +// @name Photoelectron +// +// @desc constructor for class Photoelectron +// +// @date Mon Feb 15 16:44:21 CET 2000 +// @function @code +//------------------------------------------------------------ + +Photoelectron::Photoelectron(void){ + iarrtime_ns = NOTIME; + ipixnum = -1; +} + +//------------------------------------------------------------ +// @name produce_phes +// +// @desc read the photons of an event, pixelize them and simulate QE +// @desc return various statistics and the array of Photoelectrons +// +// @date Mon Feb 14 16:44:21 CET 2000 +// @function @code +//------------------------------------------------------------ + +int produce_phes( FILE *sp, // the input file + struct camera *cam, // the camera layout + float minwl_nm, // the minimum accepted wavelength + float maxwl_nm, // the maximum accepted wavelength + class Photoelectron phe[iMAXNUMPHE], // the generated phes + int *itotnphe, // total number of produced photoelectrons + float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel + int *incph, // total number of cph read + float *tmin_ns, // minimum arrival time of all phes + float *tmax_ns // maximum arrival time of all phes + ){ + + static int i, k, ipixnum; + static float cx, cy, wl, qe, t; + static MCCphoton photon; + static float **qept; + static char flag[SIZE_OF_FLAGS + 1]; + static float radius; + + // reset variables + + for ( i=0; iinumpixels; ++i ){ + + nphe[i] = 0.0; + + } + + *itotnphe = 0; + *incph = 0; + *tmin_ns = NOTIME; // very big + *tmax_ns = -NOTIME; // very small + + radius = cam->dxc[cam->inumpixels-1] + + 1.5*cam->dpixdiameter_cm*cam->dpixsizefactor[cam->inumpixels-1]; + + //- - - - - - - - - - - - - - - - - - - - - - - - - + // read photons and "map" them into the pixels + //-------------------------------------------------- + + // initialize CPhoton + + photon.fill(0., 0., 0., 0., 0., 0., 0., 0.); + + // read the photons data + + fread ( flag, SIZE_OF_FLAGS, 1, sp ); + + // loop over the photons + + while ( !isA( flag, FLAG_END_OF_EVENT ) ) { + + memcpy( (char*)&photon, flag, SIZE_OF_FLAGS ); + + fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp ); + + // increase number of photons + + (*incph)++; + + // + // Pixelization + // + + cx = photon.get_x(); + cy = photon.get_y(); + + // get wavelength + + wl = photon.get_wl(); + + // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n"; + + // check if photon has valid wavelength and is inside outermost camera radius + + if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){ + + // cout << " lost first check\n"; + + // read next CPhoton + + fread ( flag, SIZE_OF_FLAGS, 1, sp ); + + // go to beginning of loop, the photon is lost + continue; + + } + + ipixnum = -1; + + for(i=0; iinumpixels; i++){ + if( bpoint_is_in_pix( cx, cy, i, cam) ){ + ipixnum = i; + break; + } + } + + if(ipixnum==-1){// the photon is in none of the pixels + + // cout << " lost pixlization\n"; + + // read next CPhoton + + fread ( flag, SIZE_OF_FLAGS, 1, sp ); + + // go to beginning of loop, the photon is lost + continue; + } + + //+++ + // QE simulation + //--- + + // set pointer to the QE table of the relevant pixel + + qept = (float **)QE[ipixnum]; + + // check if wl is inside table; outside the table, QE is assumed to be zero + + if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){ + + // cout << " lost\n"; + + // read next Photon + + fread ( flag, SIZE_OF_FLAGS, 1, sp ); + + // go to beginning of loop + continue; + + } + + // find data point in the QE table (-> k) + + k = 1; // start at 1 because the condition was already tested for 0 + while (k < pointsQE-1 && qept[0][k] < wl){ + k++; + } + + // calculate the qe value between 0. and 1. + + qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0; + + // if random > quantum efficiency, reject it + + if ( RandomNumber > qe ) { + + // cout << " lost\n"; + + // read next Photon + + fread ( flag, SIZE_OF_FLAGS, 1, sp ); + + // go to beginning of loop + continue; + + } + + //+++ + // The photon has produced a photo electron + //--- + + // cout << " accepted\n"; + + // increment the number of photoelectrons in the relevant pixel + + nphe[ipixnum] += 1.0; + + t = photon.get_t() ; + + // cout << " t " << t; + + // find minimum and maximum arrival time + + if(t < *tmin_ns){ + *tmin_ns = t; // memorize time + } + if(t > *tmax_ns){ + *tmax_ns = t; // memorize time + } + + // store the new photoelectron + + if(*itotnphe >= iMAXNUMPHE){ + cout << "Error: Memory overflow. Event produces more than maximum\n"; + cout << " allowed number of photoelectrons (" << iMAXNUMPHE << ").\n"; + return(1); + } + + phe[*itotnphe].iarrtime_ns = (int)t; + phe[*itotnphe].ipixnum = ipixnum; + + *itotnphe += 1; + + // read next Photon + + fread( flag, SIZE_OF_FLAGS, 1, sp ); + + } // end while, i.e. found end of event + + return(0); + +} + + +//------------------------------------------------------------ +// @name produce_nsbrates +// +// @desc read the starfield file, call produce_phes on it in, +// @desc different wavebands, calculate the nsbrates +// +// @date Mon Feb 14 16:44:21 CET 2000 +// @function @code +//------------------------------------------------------------ + +int produce_nsbrates( char *iname, // the starfield input file name + struct camera *cam, // camera layout + class Photoelectron phe[iMAXNUMPHE], // array for photoelectrons + float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function: + // the NSB rates in phe/ns for each pixel + ){ + + int i, j, k, ii; // counters + + static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1, + WAVEBANDBOUND2, + WAVEBANDBOUND3, + WAVEBANDBOUND4, + WAVEBANDBOUND5, + WAVEBANDBOUND6 }; + + FILE *infile; // the input file + fpos_t fileposition; // marker on the input file + static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable + static MCEventHeader evth; // the event header + static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel + int itnphe; // total number of produced photoelectrons + int itotnphe; // total number of produced photoelectrons after averaging + int incph; // total number of cph read + float tmin_ns; // minimum arrival time of all phes + float tmax_ns; // maximum arrival time of all phes + float integtime_ns; // integration time from the starfield generator + + // open input file + + log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname ); + + infile = fopen( iname, "r" ); + if ( infile == NULL ) + error( SIGNATURE, "Cannot open starfield input file: %s\n", iname ); + + // get signature, and check it + + if(check_reflector_file(infile)==FALSE){ + exit(1); + } + + // initialize flag + + strcpy( cflag, " \0" ); + + // get flag + + fread( cflag, SIZE_OF_FLAGS, 1, infile ); + + if ( ! feof(infile)){ + + // reading .rfl file + + if(!isA( cflag, FLAG_START_OF_RUN )){ + error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag ); + } + else { // found start of run + + fread( cflag, SIZE_OF_FLAGS, 1, infile ); + + if( isA( cflag, FLAG_START_OF_EVENT )){ // there is a event + + // get MCEventHeader + + fread( (char*)&evth, evth.mysize(), 1, infile ); + + integtime_ns = evth.get_energy(); + + // memorize where we are in the file + + if (fgetpos( infile, &fileposition ) != 0){ + error( SIGNATURE, "Cannot position in file ...\n"); + } + + // loop over the wavebands + + for(i=0; iinumpixels; j++){ // loop over pixels + rate_phepns[j][i] = 0.; + } + + itotnphe = 0; + + // read the photons and produce the photoelectrons + // - in order to average over the QE simulation, call the + // production function iNUMNSBPRODCALLS times + + for(ii=0; iiinumpixels; j++){ // loop over pixels + rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS; + } + + itotnphe += itnphe; + + } // end for(ii=0 ... + + fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n", + wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns); + + } // end for(i=0 ... + + } + else{ + cout << "Starfield file contains no event.\nExiting.\n"; + exit(1); + } // end if( isA ... event + } // end if ( isA ... run + } + else{ + cout << "Starfield file contains no run.\nExiting.\n"; + exit(1); + } + + fclose( infile ); + + return(0); + +} + + +//------------------------------------------------------------ +// @name produce_nsb_phes +// +// @desc produce the photoelectrons from the nsbrates +// +// @date Thu Feb 17 17:10:40 CET 2000 +// @function @code +//------------------------------------------------------------ + +int produce_nsb_phes( float *atmin_ns, + float *atmax_ns, + float theta_rad, + struct camera *cam, + float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS], + float difnsbr_phepns[iMAXNUMPIX], + float extinction[iNUMWAVEBANDS], + float fnpx[iMAXNUMPIX], + Photoelectron photo[iMAXNUMPHE], + int *inphe, + float base_mv[iMAXNUMPIX]){ + + float simtime_ns; // NSB simulation time + int i, j, k, ii; + float zenfactor; // correction factor calculated from the extinction + int inumnsbphe; // number of photoelectrons caused by NSB + + ii = *inphe; // avoid dereferencing + + // check if the arrival times are set; if not generate them + + if(*atmin_ns == NOTIME){ + + *atmin_ns = 0.; + *atmax_ns = simtime_ns = SLICES*WIDTH_TIMESLICE; + + } + else{ // extend the event time window by the given offsets + + *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS; + *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS; + + simtime_ns = *atmax_ns - *atmin_ns; + + // make sure the simulated time is long enough for the FADC simulation + + if(simtime_ns< SLICES*WIDTH_TIMESLICE){ + *atmax_ns = *atmin_ns + SLICES*WIDTH_TIMESLICE; + simtime_ns = SLICES*WIDTH_TIMESLICE; + } + + } + + // initialize baselines + + for(i=0; iinumpixels; i++){ + base_mv[i] = 0.; + } + + // calculate baselines and generate phes + + for(i=0; iinumpixels; j++){ // loop over the pixels + + inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS) + * simtime_ns ); + + base_mv[j] += inumnsbphe; + + // randomize + + inumnsbphe = ignpoi( inumnsbphe ); + + // create the photoelectrons + + for(k=0; k < inumnsbphe; k++){ + + if(ii >= iMAXNUMPHE){ + cout << "Error: Memory overflow. NSB simulation produces more than maximum\n"; + cout << " allowed number of photoelectrons (" << iMAXNUMPHE << ").\n"; + return(1); + } + + photo[ii].iarrtime_ns = (int)(RandomNumber * simtime_ns + *atmin_ns ); + photo[ii].ipixnum = j; + + // cout << "Created phe " << photo[ii].iarrtime_ns << " " + // << photo[ii].ipixnum << "\n"; + + ii++; // increment total number of photoelectons + + fnpx[j] += 1.; // increment number of photoelectrons in this pixel + + } + + } // end for(j=0 ... + } // end for(i=0 ... + + // finish baseline calculation + + for(i=0; iinumpixels; i++){ + base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns; + } + + *inphe = ii; // update the pointer + + return(0); +} + + // @endcode @@ -2977,4 +2779,8 @@ // // $Log: not supported by cvs2svn $ +// Revision 1.4 2000/01/25 08:36:23 petry +// The pixelization in previous versions was buggy. +// This is the first version with a correct pixelization. +// // Revision 1.3 2000/01/20 18:22:17 petry // Found little bug which makes camera crash if it finds a photon @@ -2996,218 +2802,4 @@ // 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