//!///////////////////////////////////////////////////////////////////// // // 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: camera.cxx,v $ // $Revision: 1.57 $ // $Author: blanch $ // $Date: 2003-07-17 18:02:46 $ // //////////////////////////////////////////////////////////////////////// // @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 "TRandom.h" #include "TApplication.h" #include "TFile.h" #include "TTree.h" #include "TBranch.h" #include "TCanvas.h" #include "TArrayC.h" #include "MTrigger.hxx" #include "MFadc.hxx" #include "MLons.hxx" #include "MRawRunHeader.h" #include "MRawEvtData.h" #include "MRawEvtHeader.h" #include "MRawCrateData.h" #include "MMcEvt.hxx" #include "MMcTrig.hxx" #include "MMcRunHeader.hxx" #include "MMcConfigRunHeader.h" #include "MMcCorsikaRunHeader.h" #include "MMcTrigHeader.hxx" #include "MMcFadcHeader.hxx" #include "MGeomCamMagic.h" #include "MGeomPix.h" /*!@" All the defines are located in the file |camera.h|. @"*/ #include "camera.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__ //!@} //=----------------------------------------------------------- //!@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; //@: 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; //@: 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; //@: flag: TRUE: write all data to output; FALSE: only triggered showers static int Write_All_Data = FALSE; static int Write_McEvt = TRUE; static int Write_McTrig = TRUE; static int Write_McFADC = TRUE; 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; //@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't static int FADC_Scan = FALSE; //@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't static int Trigger_Scan = FALSE; //@: flag: TRUE: loop trigger analysis over several thresholds, multiplicities and topologies; FALSE: a single trigger configuration static int Trigger_Loop = FALSE; //@: flag: TRUE: Different threshold for each pixel ; FALSE: same threshold for all pixels static int Individual_Thres_Pixel = FALSE; //@: Properties of the trigger static float Trigger_gate_length = 6.0; static float Trigger_response_ampl = 1.0; static float Trigger_response_fwhm = 2.0; static float Trigger_overlaping_time= 0.25; static float Trigger_noise= 0.3; //@: Properties of the FADC static float FADC_response_ampl = MFADC_RESPONSE_AMPLITUDE; static float FADC_response_fwhm = MFADC_RESPONSE_FWHM; static float FADC_resp_ampl_out = MFADC_RESPONSE_AMPLITUDE; static float FADC_resp_fwhm_out = MFADC_RESPONSE_FWHM; static float FADC_noise = 2.0; //@: Trigger conditions for a single trigger mode static float qThreshold[CAMERA_PIXELS]; static int Trigger_multiplicity = 4; static int Trigger_topology = 2; //@: Upper and lower edges of the trigger loop static float Trigger_loop_lthres = 2.0; static float Trigger_loop_uthres = 10.0; static float Trigger_loop_sthres = 1.0; static int Trigger_loop_lmult = 2; static int Trigger_loop_umult = 10; static int Trigger_loop_ltop = 0; static int Trigger_loop_utop = 2; //!@} /*!@" 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 //@: 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; //!@} /*!@" 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; //@: table of lightguide = WC; static float **WC; //@: number of datapoints for the WC curve static int pointsWC; //!@} /*!@" 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|. @"*/ //!@} extern char FileName[]; // switch on starfield rotation static int Starfield_rotate = TRUE; //=----------------------------------------------------------- // @subsection Main program. //!@{ //++++++++++++++++++++++++++++++++++++++++ // MAIN PROGRAM //---------------------------------------- int main(int argc, char **argv) { //!@' @#### Definition of variables. //@' char inname[256]; //@< input file name char starfieldname[256]; //@< starfield input file name char qe_filename[256]; //@< qe input file name char datname[256]; //@< data (ASCII) output file name char rootname[256] ; //@< ROOT file name char rootname_loop[256] ; //@< ROOT file name char parname[256]; //@< parameters file name char nsbpathname[256]; //@< directory with the dataabse for th NSB char nsbpath_outer[256]; //@< directory with the dataabse for outer pixels char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file char flag_new[SIZE_OF_FLAGS + 1]; //@< New flag for the run header in the .rfl file int reflector_file_version; //@< vrsion of he reflector file FILE *inputfile; //@< stream for the input file ofstream datafile; //@< stream for the data file MCRunHeader mcrunh; //@< Run Header class (MC) MCEventHeader mcevth; //@< Event Header class (MC) MCEventHeader_2 mcevth_2; //@< Event Header class (MC) for reflector > 0.6 MCCphoton cphoton; //@< Cherenkov Photon class (MC) int inumphe; //@< number of photoelectrons in an event from showers float inumphensb; //@< number of photoelectrons in an event from nsb 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 float thetashw, phishw; //@< parameters of a given shower float coreX, coreY; //@< core position 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=0; //@< partial number of photons in a given run int ntcph=0; //@< total number of photons int i, j, k; //@< simple counters int addElecNoise; //@< Will we add ElecNoise? float riseDiskThres; float secureDiskThres; int simulateNSB; //@< Will we simulate NSB? int nphe2NSB; //@< From how many phe will we simulate NSB? 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 nsb_phepns[iMAXNUMPIX]; //@< NSB values for each pixel float nsb_phepns_rotated[iMAXNUMPIX]; //@< NSB values for each pixel after rotation Float_t nsb_trigresp[TRIGGER_TIME_SLICES]; //@< array to write the trigger //@< response from the database Float_t nsb_fadcresp[(Int_t) SLICES_MFADC]; //@< array to write the fadc //@< response from the database Byte_t trigger_map[((Int_t)(CAMERA_PIXELS/8))+1]; //@< Pixels on when the //@< camera triggers Float_t fadc_elecnoise[CAMERA_PIXELS]; //@< Electronic niose for each pixel Float_t fadc_pedestals[CAMERA_PIXELS]; //@< array for fadc pedestals values float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband EXTWAVEBAND1, EXTWAVEBAND2, EXTWAVEBAND3, EXTWAVEBAND4, EXTWAVEBAND5 }; float zenfactor=1.0; // correction factor calculated from the extinction int anaPixels; int nstoredevents = 0; int flagstoring = 0; int ntrigger = 0; //@< number of triggers in the whole file int btrigger = 0; //@< trigger flag int ithrescount; //@< counter for loop over threshold trigger float fthrescount; //@< value for loop over threshold trigger int imulticount; //@< counter for loop over multiplicity trigger int itopocount; //@< counter for loop over topology trigger int isorttopo[3]; //@< sorting the topologies int icontrigger; //@< number of trigger conditions to be analised float fpixelthres[CAMERA_PIXELS]; //@< Threshold values TArrayC *fadcValues; //@< the analog Fadc High gain signal for pixels TArrayC *fadcValuesLow; //@< the analog Fadc Low gain signal for pixels int still_in_loop = FALSE; //@< Variables to fill the McRunHeader Int_t sfRaH = 5; Int_t sfRaM = 34; Int_t sfRaS = 32; Int_t sfDeD = 22; Int_t sfDeM = 00; Int_t sfDeS = 52; Float_t shthetamax = 0.0; Float_t shthetamin = 0.0; Float_t shphimax = 0.0; Float_t shphimin = 0.0; Float_t telestheta = 0.0; Float_t telesphi = 0.0; Float_t sofftheta = 0.0; Float_t soffphi = 0.0; UInt_t corsika = 5200 ; // Star Field Rotation variables Float_t zenith = 0.0; Float_t azimutal = 90.0; Float_t rho , C3 , C2 , C1; MGeomCamMagic camgeom; // structure holding the camera definition ct_NPixels=camgeom.GetNumPixels(); //!@' @#### Definition of variables for |getopt()|. //@' int ch, errflg = 0; //@< used by getopt //!@' @#### Initialisation of some variables //@' for(i=0;i0 ) 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(); // read WC and QE files strcpy( qe_filename, get_qe_filename()); read_QE(qe_filename); read_WC(); // write all images, even those without trigger? Write_All_Images = get_write_all_events(); // write all data (i.e., ph.e.s in pixels) Write_All_Data = get_write_all_data(); Write_McEvt = get_write_McEvt() ; Write_McTrig = get_write_McTrig() ; Write_McFADC = get_write_McFadc() ; Write_RawEvt = get_write_RawEvt() ; FADC_Scan = get_FADC_Scan(); Trigger_Scan = get_Trigger_Scan(); Individual_Thres_Pixel = get_indi_thres_pixel(); get_secure_threhold(&riseDiskThres,&secureDiskThres); get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm, &FADC_resp_ampl_out, &FADC_resp_fwhm_out); get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm); Trigger_Loop = get_Trigger_Loop(&Trigger_loop_lthres, &Trigger_loop_uthres, &Trigger_loop_sthres, &Trigger_loop_lmult, &Trigger_loop_umult, &Trigger_loop_ltop, &Trigger_loop_utop); icontrigger =((int)((Trigger_loop_uthres-Trigger_loop_lthres) /Trigger_loop_sthres)+1)* (Trigger_loop_umult-Trigger_loop_lmult+1)* (Trigger_loop_utop-Trigger_loop_ltop+1); if (!Trigger_Loop){ get_Trigger_Single (qThreshold, &Trigger_multiplicity, &Trigger_topology); icontrigger=1; } else get_threshold(qThreshold); // get filenames strcpy( inname, get_input_filename() ); strcpy( starfieldname, get_starfield_filename() ); strcpy( datname, get_data_filename() ); strcpy( rootname, get_root_filename() ); strcpy( rootname_loop, get_loop_filename() ); strcpy( ct_filename, get_ct_filename() ); strcpy( nsbpathname, get_nsb_directory() ); strcpy( nsbpath_outer, get_nsb_directory_outer() ); // get different parameters of the simulation addElecNoise = add_elec_noise(&FADC_noise, &Trigger_noise); simulateNSB = get_nsb( &meanNSB, &nphe2NSB ); // get selections on the parameters Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue); //Parameters for starfield rotation Starfield_rotate = get_starfield_rotate(); log(SIGNATURE, "%s:\n\t%20s: %s\n", "Starfield Rotate", "Rotate Starfield", ONoff(Starfield_rotate) ); // log filenames information log(SIGNATURE, "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n", "Filenames", "In", inname, "Stars", starfieldname, "NSB database (inner pixels)", nsbpathname, "NSB database (outer pixels)", nsbpath_outer, "CT", ct_filename, "Data", datname, "ROOT", rootname ); // log Trigger information if (Trigger_Loop) { log(SIGNATURE, "%s:\n\t%20s: from %5.2f to %5.2f by %5.2f step\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s\n", "Trigger Loop mode", "Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres, "Multiplicity",Trigger_loop_lmult,Trigger_loop_umult, "Topology",Trigger_loop_ltop,Trigger_loop_utop, rootname_loop); } else if (Individual_Thres_Pixel == FALSE){ log(SIGNATURE, "%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n", "Single Trigger mode", "Threshold",qThreshold[0], "Multiplicity",Trigger_multiplicity, "Topology",Trigger_topology); } else{ log(SIGNATURE, "%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n", "Single Trigger mode", "Threshold","Different for each pixel", "Multiplicity",Trigger_multiplicity, "Topology",Trigger_topology); } // log flags information log(SIGNATURE, "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n", "Flags", "Data_From_STDIN", ONoff(Data_From_STDIN), "Write_All_Events", ONoff(Write_All_Images), "Write_All_Data", ONoff(Write_All_Data)); // log flags information log(SIGNATURE, "%s:\n\t%20s: %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_McFADC", ONoff(Write_McFADC), "Write_RawEvt", ONoff(Write_RawEvt)); // log parameters information log(SIGNATURE, "%s:\n\t%20s: %f\n", "Parameters", "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB)); // log selections log(SIGNATURE, "%s:\n\t%20s: %s (%f:%f)\n", "Selections:", "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue); // Definition and initialization of array to save trigger statistics int ***ntriggerloop; ntriggerloop= new int ** [(int)((Trigger_loop_uthres-Trigger_loop_lthres) /Trigger_loop_sthres)]; for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;fthrescount+=Trigger_loop_sthres, ithrescount++){ ntriggerloop[ithrescount]= new int * [Trigger_loop_umult-Trigger_loop_lmult+1]; for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){ ntriggerloop[ithrescount][imulticount]= new int [Trigger_loop_utop-Trigger_loop_ltop+1]; for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){ ntriggerloop[ithrescount][imulticount][itopocount]=0; } } } // We should be careful that topologies are sort from // the less to the more restrictive one. if (Trigger_loop_utop==Trigger_loop_ltop) for(int is=0; is<3;is++) isorttopo[is]=is; else { isorttopo[0]=1; isorttopo[1]=0; isorttopo[2]=2; } // 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; iIgnoreTObjectStreamer(); MParContainer::Class()->IgnoreTObjectStreamer(); // initialise ROOT TROOT simple("simple", "MAGIC Telescope Monte Carlo"); // initialise instance of Trigger and FADC classes MTrigger Trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm); //@< A instance of the Class MTrigger // Generate databse for trigger electronic noise Trigger.SetElecNoise(Trigger_noise) ; // Set Right Discriminator threshold, taking into account trigger pixels Trigger.CheckThreshold(&qThreshold[0]); // Set flag in pixel 0 (not used for trigger) that indicates if secure pixel // is active: secureDiskThres*10000+riseDiskThres if(riseDiskThres>0.0) qThreshold[0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres; else qThreshold[0]=0.0; // Initialise McTrig information class if we want to save trigger informtion MMcTrig **McTrig = NULL; MMcTrigHeader **HeaderTrig = NULL; MMcFadcHeader **HeaderFadc = NULL; if (Write_McTrig){ McTrig = new MMcTrig * [icontrigger]; for (i=0;iInit(RunHeader); // We need thr RunHeader to read // number of pixels } } MMcEvt *McEvt = new MMcEvt (); // // initalize a temporal ROOT file // TFile outfile_temp ( rootname , "RECREATE" ); // create a Tree for the Header Event TTree HeaderTree("RunHeaders","Headers of Run"); // define branches of Header Tree char help[4]; HeaderTree.Branch("MRawRunHeader","MRawRunHeader", &RunHeader); HeaderTree.Branch("MMcRunHeader","MMcRunHeader", &McRunHeader); HeaderTree.Branch("MMcConfigRunHeader","MMcConfigRunHeader", &McConfigRunHeader); HeaderTree.Branch("MMcCorsikaRunHeader","MMcCorsikaRunHeader", &McCorsikaRunHeader); if(!Trigger_Loop && Write_McTrig){ HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader", &HeaderTrig[0]); } if (Trigger_Loop && Write_McTrig){ for(char branchname[20],i=0;iSetMagicNumber(kMagicNumber); RunHeader->SetFormatVersion(4); RunHeader->SetSoftVersion((UShort_t) (VERSION*10)); RunHeader->SetRunType(256); RunHeader->SetRunNumber(0); RunHeader->SetNumSamples(0,FADC_SLICES); RunHeader->SetNumCrates(1); RunHeader->SetNumPixInCrate(CAMERA_PIXELS); // Fill branches for MMcTrigHeader if(!Trigger_Loop && Write_McTrig){ HeaderTrig[0]->SetTopology((Short_t) Trigger_topology); HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity); HeaderTrig[0]->SetThreshold(qThreshold); HeaderTrig[0]->SetAmplitud(Trigger_response_ampl); HeaderTrig[0]->SetFwhm(Trigger_response_fwhm); HeaderTrig[0]->SetOverlap(Trigger_overlaping_time); HeaderTrig[0]->SetGate(Trigger_gate_length); HeaderTrig[0]->SetElecNoise(Trigger_noise); } if(Trigger_Loop && Write_McTrig){ int iconcount; for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++,fthrescount+=Trigger_loop_sthres){ for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){ for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){ HeaderTrig[iconcount]->SetTopology((Short_t) isorttopo[itopocount+Trigger_loop_ltop]); HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult); for(int i=0;i=qThreshold[i])? (Float_t)(fthrescount):qThreshold[i]; } HeaderTrig[iconcount]->SetThreshold( fpixelthres); HeaderTrig[iconcount]->SetAmplitud(Trigger_response_ampl); HeaderTrig[iconcount]->SetFwhm(Trigger_response_fwhm); HeaderTrig[iconcount]->SetOverlap(Trigger_overlaping_time); HeaderTrig[iconcount]->SetGate(Trigger_gate_length); HeaderTrig[iconcount]->SetElecNoise(Trigger_noise); iconcount++; } } } } // Fill branches for MMcFadcHeader for(i=0;iSetShape(0.0); HeaderFadc[0]->SetAmplitud(FADC_response_ampl); HeaderFadc[0]->SetFwhm(FADC_response_fwhm); HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],anaPixels); HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0],anaPixels); } if(Trigger_Loop && Write_McFADC){ int iconcount; for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){ for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){ for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){ HeaderFadc[iconcount]->SetShape(0.0); HeaderFadc[iconcount]->SetAmplitud(Trigger_response_ampl); HeaderFadc[iconcount]->SetFwhm(Trigger_response_fwhm); HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0],anaPixels); HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],anaPixels); iconcount++; } } } } // Fill the Header Tree with the current leaves of each branch // HeaderTree.Fill() ; // create a Tree for the Event data stream TTree EvtTree("Events","Normal Triggered Events"); if (Write_McEvt){ EvtTree.Branch("MMcEvt","MMcEvt", &McEvt); } if(!Trigger_Loop){ if (Write_RawEvt){ EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", &EvtHeader[0]); EvtTree.Branch("MRawEvtData","MRawEvtData", &EvtData[0]); } if (Write_McTrig){ EvtTree.Branch("MMcTrig","MMcTrig", &McTrig[0]); } } else{ if (Write_McTrig){ for(char branchname[10],i=0;iIsBatch()) { fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); // return 1; } } if(FADC_Scan){ //TApplication theAppFadc("App", &argc, argv); if (gROOT->IsBatch()) { fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); // return 1; } } // prepare the NSB simulation // Instance of the Mlons class MLons lons(Trigger_response_ampl, Trigger_response_fwhm, FADC_response_ampl,FADC_response_fwhm); lons.SetPath(nsbpathname); // Instance of the Mlons class MLons lons_outer(Trigger_response_ampl, Trigger_response_fwhm, FADC_resp_ampl_out,FADC_resp_fwhm_out); lons_outer.SetPath(nsbpath_outer); if( simulateNSB ){ // // Calculate the non-diffuse NSB photoelectron rates // k = produce_nsbrates( starfieldname, &camgeom, nsbrate_phepns); if (k != 0){ cout << "Error when reading starfield... \nExiting.\n"; exit(1); } // calculate diffuse rate correcting for the pixel size const double size_inner = camgeom[0].GetR(); for(UInt_t ui=0; ui 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 if (reflector_file_version<6){ thetashw = mcevth.get_theta(); phishw = mcevth.get_phi(); } else{ thetashw = mcevth_2.get_theta(); phishw = mcevth_2.get_phi(); } // calculate vector for shower l1 = sin(thetashw)*cos(phishw); m1 = sin(thetashw)*sin(phishw); n1 = cos(thetashw); if (reflector_file_version<6){ mcevth.get_core(&coreX, &coreY); // 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. ); thetaCT += thetashw; phiCT += phishw; } 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; } } else{ mcevth_2.get_core(&coreX, &coreY); thetaCT=mcevth_2.get_theta_CT(); phiCT=mcevth_2.get_phi_CT(); if ( (thetaCT == thetashw) && (phiCT == phishw) ) { // 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 // 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; } } // energy cut if ( Select_Energy ) { if (reflector_file_version<6) if (( mcevth.get_energy() < Select_Energy_le ) || ( mcevth.get_energy() > Select_Energy_ue )) { log(SIGNATURE, "select_energy: shower rejected.\n"); continue; } else if (( mcevth_2.get_energy() < Select_Energy_le ) || ( mcevth_2.get_energy() > Select_Energy_ue )) { log(SIGNATURE, "select_energy: shower rejected.\n"); continue; } } // Read first and last time and put inumphe to 0 if (reflector_file_version<6) mcevth.get_times(&arrtmin_ns,&arrtmax_ns); else mcevth_2.get_times(&arrtmin_ns,&arrtmax_ns); inumphe=0; // read the photons and produce the photoelectrons k = produce_phes( inputfile, &camgeom, WAVEBANDBOUND1, WAVEBANDBOUND6, &Trigger, // will be changed by the function! &fadc, // will be changed by the function! &inumphe, // important for later: the size of photoe[] fnpix, // will be changed by the function! &ncph, // will be changed by the function! &arrtmin_ns, // will be changed by the function! &arrtmax_ns); // will be changed by the function! if( k != 0 ){ // non-zero returnvalue means error cout << "Exiting.\n"; exit(1); } // NSB simulation if(simulateNSB && nphe2NSB<=inumphe){ if(Starfield_rotate){ // Introduction rho angle zenith = thetashw; azimutal = phishw; C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal); C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal)); C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + C3 * C3 ); rho = acos( C1/C2 ); if ( sin(azimutal) < 0) { rho = 2 * 3.14159 - rho; } else { rho = rho; } rho = rho*180/3.14159; // Rottion of the NSB k = size_rotated( &nsb_phepns_rotated[0], nsb_phepns, rho); } // Fill trigger and fadc response in the trigger class from the database for(UInt_t ui=0;ui0.0){ if(camgeom[ui].GetR()>camgeom[0].GetR()){ k=lons_outer.GetResponse(nsb_phepns_rotated[ui],0.1, & nsb_trigresp[0],& nsb_fadcresp[0]); } else{ k=lons.GetResponse(nsb_phepns_rotated[ui],0.1, & nsb_trigresp[0],& nsb_fadcresp[0]); } if(k==0){ cout << "Exiting.\n"; exit(1); } Trigger.AddNSB(ui,nsb_trigresp); fadc.AddSignal(ui,nsb_fadcresp); } }// end if(simulateNSB && nphe2NSB<=inumphe) ... inumphensb=0; for (UInt_t ui=0;ui=qThreshold[i])? (Float_t)(fthrescount):qThreshold[i]; // Rise the discrimnator threshold to avoid huge rates if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe) for(int ii=0;iiriseDiskThres) fpixelthres[ii]=secureDiskThres; } Trigger.SetThreshold(fpixelthres); Trigger.Diskriminate(); // // look if in all the signals in the trigger signal branch // is a possible Trigger. Therefore we habe to diskriminate all // the simulated analog signals (Method Diskriminate in class // MTrigger). We look simultanously for the moments at which // there are more than TRIGGER_MULTI pixels above the // CHANNEL_THRESHOLD. // // Set trigger flags to zero Lev1=0; // loop over multiplicity of trigger configuration for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){ Trigger.SetMultiplicity(imulticount); Trigger.ClearZero(); Lev0=(Short_t) Trigger.ZeroLevel(); if (Lev0>0 || Write_All_Images || btrigger){ // loop over topologies for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){ Lev1=0; if(itopocount==0 && imulticount>7) continue; //COBB if(itopocount==2 && imulticount<3) continue; // It only makes to look for a different topology // if there are 3 or more N pixels. if(imulticount<3) Trigger.SetTopology(1); else { // We should be careful that topologies are sort from // the less to the more restrictive one. Trigger.SetTopology(isorttopo[itopocount]); } Trigger.ClearFirst(); // // Start the First Level Trigger simulation // if(Lev0!=0) Lev1=Trigger.FirstLevel(); if(Lev1>0) { btrigger= 1; ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++; } Lev0=1; Int_t NumImages = Lev1; if(Lev1==0 && (Write_All_Images || btrigger)){ btrigger= 1; NumImages=1; Lev0=0; } for (Int_t ii=0;iiSetFirstLevel ((ii+1)*Lev0); McTrig[iconcount]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1); Trigger.GetMapDiskriminator(trigger_map); McTrig[iconcount]->SetMapPixels(trigger_map,ii); } // // fill inside the class fadc the member output // fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii)); if( Write_RawEvt ){ // // Fill the header of this event // EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0); // fill pixel information if (Lev1){ for(UInt_t i=0;iAddAt(fadc.GetFadcSignal(i,j),j); fadcValuesLow->AddAt(fadc.GetFadcLowGainSignal(i,j),j); } EvtData[iconcount]->AddPixel(i,fadcValues,0); EvtData[iconcount]->AddPixel(i,fadcValuesLow,kTRUE); } } } } // // Increase counter of analised trigger conditions // iconcount++; } } else{ break; } } if (!btrigger) break; } if (btrigger){ // // fill the MMcEvt with all information // if(!flagstoring) nstoredevents++; flagstoring = 1; if (Write_McEvt) { Float_t ftime, ltime; if (reflector_file_version<6){ mcevth.get_times(&ftime, <ime); McEvt->Fill( 0, (UShort_t) mcevth.get_primary() , mcevth.get_energy(), -1.0, -1.0, -1.0, mcevth.get_theta(), mcevth.get_phi(), mcevth.get_core(), mcevth.get_coreX(), mcevth.get_coreY(), impactD, phiCT, thetaCT, ftime, ltime, 0, 0, 0, 0, 0, 0, 0, (UInt_t)mcevth.get_CORSIKA(), (UInt_t)mcevth.get_AtmAbs(), (UInt_t)(mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot()), (UInt_t) ncph, (UInt_t) inumphe, (UInt_t) inumphensb+inumphe, -1.0, -1.0, -1.0); } else{ Float_t Nmax, t0, tmax, a, b, c, chi2; mcevth_2.get_times(&ftime, <ime); chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c); McEvt->Fill((UInt_t) mcevth_2.get_evt_number(), (UShort_t) mcevth_2.get_primary() , mcevth_2.get_energy(), mcevth_2.get_thick0(), mcevth_2.get_first_target(), mcevth_2.get_z_first_int(), mcevth_2.get_theta(), mcevth_2.get_phi(), mcevth_2.get_core(), mcevth_2.get_coreX(), mcevth_2.get_coreY(), impactD, mcevth_2.get_phi_CT(), mcevth_2.get_theta_CT(), ftime, ltime, Nmax, t0, tmax, a, b, c, chi2, (UInt_t)mcevth_2.get_CORSIKA(), (UInt_t)mcevth_2.get_AtmAbs(), (UInt_t)(mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot()), (UInt_t) ncph, (UInt_t) inumphe, (UInt_t) inumphensb+inumphe, mcevth_2.get_ElecFraction(), mcevth_2.get_MuonFraction(), mcevth_2.get_OtherFraction()); } } // Fill the Tree with the current leaves of each branch i=EvtTree.Fill() ; // Clear the branches if(Write_McTrig){ for(i=0;iClear() ; } } if( Write_RawEvt ){ for(i=0;iClear() ; EvtData[i]->DeletePixels(); } } if (Write_McEvt) McEvt->Clear() ; } } // We study a single trigger condition else { // Setting trigger conditions Trigger.SetMultiplicity(Trigger_multiplicity); Trigger.SetTopology(Trigger_topology); for (int i=0;i0.0 && simulateNSB && nphe2NSB<=inumphe) for(int ii=0;iiriseDiskThres) fpixelthres[ii]=secureDiskThres; } Trigger.SetThreshold(fpixelthres); Trigger.Diskriminate() ; // // look if in all the signals in the trigger signal branch // is a possible Trigger. Therefore we habe to diskriminate all // the simulated analog signals (Method Diskriminate in class // MTrigger). We look simultanously for the moments at which // there are more than TRIGGER_MULTI pixels above the // CHANNEL_THRESHOLD. // Lev0 = (Short_t) Trigger.ZeroLevel() ; Lev1 = 0 ; // // Start the First Level Trigger simulation // if ( Lev0 > 0 || Write_All_Images) { Lev1= Trigger.FirstLevel(); } if (Lev1>0){ ++ntrigger; } Int_t NumImages = Lev1; Lev0=1; if (Lev1==0 && Write_All_Images){ NumImages=1; Lev0=0; } for(Int_t ii=0;iiSetFirstLevel ((ii+1)*Lev0); McTrig[0]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1); Trigger.GetMapDiskriminator(trigger_map); McTrig[0]->SetMapPixels(trigger_map,ii); } // Fill Evt information if (Write_RawEvt){ // // Fill the header of this event // EvtHeader[0]->FillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ; // fill pixel information if (Lev1){ for(UInt_t i=0;iAddAt(fadc.GetFadcSignal(i,j),j); fadcValuesLow->AddAt(fadc.GetFadcLowGainSignal(i,j),j); } EvtData[0]->AddPixel(i,fadcValues,0); EvtData[0]->AddPixel(i,fadcValuesLow,kTRUE); } } } // // fill the MMcEvt with all information // if (Write_McEvt){ Float_t ftime, ltime; if (reflector_file_version<6){ mcevth.get_times(&ftime, <ime); McEvt->Fill( 0, (UShort_t) mcevth.get_primary() , mcevth.get_energy(), -1.0, -1.0, -1.0, mcevth.get_theta(), mcevth.get_phi(), mcevth.get_core(), mcevth.get_coreX(), mcevth.get_coreY(), impactD, phiCT, thetaCT, ftime, ltime, 0, 0, 0, 0, 0, 0, 0, (UInt_t)mcevth.get_CORSIKA(), (UInt_t)mcevth.get_AtmAbs(), (UInt_t)(mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot()), (UInt_t) ncph, (UInt_t) inumphe, (UInt_t) inumphensb+inumphe, -1.0, -1.0, -1.0); } else{ Float_t Nmax, t0, tmax, a, b, c, chi2; mcevth_2.get_times(&ftime, <ime); chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c); McEvt->Fill( (UInt_t) mcevth_2.get_evt_number(), (UShort_t) mcevth_2.get_primary() , mcevth_2.get_energy(), mcevth_2.get_thick0(), mcevth_2.get_first_target(), mcevth_2.get_z_first_int(), mcevth_2.get_theta(), mcevth_2.get_phi(), mcevth_2.get_core(), mcevth_2.get_coreX(), mcevth_2.get_coreY(), impactD, mcevth_2.get_phi_CT(), mcevth_2.get_theta_CT(), ftime, ltime, Nmax, t0, tmax, a, b, c, chi2, (UInt_t)mcevth_2.get_CORSIKA(), (UInt_t)mcevth_2.get_AtmAbs(), (UInt_t) (mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot()), (UInt_t) ncph, (UInt_t) inumphe, (UInt_t) inumphensb+inumphe, mcevth_2.get_ElecFraction(), mcevth_2.get_MuonFraction(), mcevth_2.get_OtherFraction()); } } // We don not count photons out of the camera. // // write it out to the file outfile // EvtTree.Fill() ; // // if a first level trigger occurred, then // 1. do some other stuff (not implemented) // 2. start the gui tool if(FADC_Scan){ if ( Lev0 > 0 ) { fadc.ShowSignal( McEvt, (Float_t) 60. ) ; } } if(Trigger_Scan){ if ( Lev0 > 0 ) { Trigger.ShowSignal(McEvt) ; } } // clear all if (Write_RawEvt) EvtHeader[0]->Clear() ; if (Write_RawEvt) EvtData[0]->DeletePixels(); if (Write_McEvt) McEvt->Clear() ; } if (Write_McTrig) McTrig[0]->Clear() ; } #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 (int i=0; iFill(rnum, (UInt_t) 0, mcevth.get_DateRun(), ftime, icontrigger, !Write_All_Images, Write_McEvt, Write_McTrig, Write_McFADC, Write_RawEvt, addElecNoise, CAMERA_PIXELS, (UInt_t)ntshow, (UInt_t)nstoredevents, 0, sfRaH, sfRaM, sfRaS, sfDeD, sfDeM, sfDeS, meanNSB, telestheta, telesphi, sofftheta, soffphi, shthetamax, shthetamin, shphimax, shphimin, mcevth.get_CWaveLower(), mcevth.get_CWaveUpper(), mcevth.get_slope(), 1, heights, corsika, (UInt_t)(reflector_file_version*100), (UInt_t)(VERSION*100), 0); else McRunHeader->Fill(rnum, (UInt_t) 0, mcevth_2.get_DateRun(), ftime, icontrigger, !Write_All_Images, Write_McEvt, Write_McTrig, Write_McFADC, Write_RawEvt, addElecNoise, CAMERA_PIXELS, (UInt_t)ntshow, (UInt_t)nstoredevents, 0, sfRaH, sfRaM, sfRaS, sfDeD, sfDeM, sfDeS, meanNSB, telestheta, telesphi, sofftheta, soffphi, shthetamax, shthetamin, shphimax, shphimin, mcevth_2.get_CWaveLower(), mcevth_2.get_CWaveUpper(), mcevth_2.get_slope(), 1, heights, corsika, (UInt_t)(reflector_file_version*100), (UInt_t)(VERSION*100), 0); // Fill some missing values for MRawRunHeader RunHeader->SetRunNumber(rnum); RunHeader->SetNumEvents(nstoredevents); // Fill MMcCorsikaRunHeader Float_t constantC[50]; mcrunh.get_constantC(&constantC[0]); Float_t constantCKA[40]; mcrunh.get_constantCKA(&constantCKA[0]); Float_t constantCETA[5]; mcrunh.get_constantCETA(&constantCETA[0]); Float_t constantCSTRBA[11]; mcrunh.get_constantCSTRBA(&constantCSTRBA[0]); Float_t constantAATM[5]; mcrunh.get_constantAATM(&constantAATM[0]); Float_t constantBATM[5]; mcrunh.get_constantBATM(&constantBATM[0]); Float_t constantCATM[5]; mcrunh.get_constantCATM(&constantCATM[0]); Float_t constantNFL[4]; mcrunh.get_constantNFL(&constantNFL[0]); if(reflector_file_version>5) McCorsikaRunHeader->Fill(rnum, mcrunh.get_date(), corsika, 1, heights, mcevth_2.get_slope(), mcrunh.get_ELow(), mcrunh.get_EUpp(), mcrunh.get_EGS4(), mcrunh.get_NKG(), mcrunh.get_Ecutoffh(), mcrunh.get_Ecutoffm(), mcrunh.get_Ecutoffe(), mcrunh.get_Ecutoffg(), constantC, constantCKA, constantCETA, constantCSTRBA, constantAATM, constantBATM, constantCATM, constantNFL, mcrunh.get_wobble(), mcrunh.get_atmophere() ); // Store qe for each PMT in output file TArrayF qe_pmt; TArrayF wav_pmt; McConfigRunHeader->InitSizePMTs(ct_NPixels); for(int i=0; iAddPMT(i); MGeomPMT &pmt = McConfigRunHeader->GetPMT(i); qe_pmt.Set(pointsQE,QE[i][1]); wav_pmt.Set(pointsQE,QElambda); pmt.SetArraySize(pointsQE); pmt.SetPMTContent(i,wav_pmt,qe_pmt); } // Store Light Guides factors in the output file TArrayF theta_lg; TArrayF factor_lg; theta_lg.Set(pointsWC,WC[0]); factor_lg.Set(pointsWC,WC[1]); McConfigRunHeader->SetLightGuides(theta_lg,factor_lg); // Fill the Header Tree with the current leaves of each branch HeaderTree.Fill() ; //++ // put the Event to the root file //-- outfile_temp.Write() ; outfile_temp.Close() ; // close input file log( SIGNATURE, "%d event(s), with a total of %d C.photons\n", ntshow, ntcph ); datafile< (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)); 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_QE // // @desc read QE data // // @date thu 5 17:59:57 CEST 2002 //------------------------------------------------------------ // @function //!@{ void read_QE(char fname[256]){ ifstream qefile; char line[LINE_MAX_LENGTH]; int i, j, icount; float qe; //------------------------------------------------------------ // second, pixels' QE // try to open the file log("read_QE", "Opening the file \"%s\" . . .\n", fname); qefile.open( fname ); // if it is wrong or does not exist, exit if ( qefile.bad() ) error( "read_QE", "Cannot open \"%s\". Exiting.\n", fname ); // read file log("read_QE", "Reading data . . .\n"); i=-1; icount = 0; while ( ! qefile.eof() ) { // get line from the file qefile.getline(line, LINE_MAX_LENGTH); // skip if comment if ( *line == '#' ) continue; // if it is the first valid value, it is the number of QE data points if ( i < 0 ) { // get the number of datapoints sscanf(line, "%d", &pointsQE); // allocate memory for the table of QEs QE = new float ** [ct_NPixels]; for ( i=0; i -1) && ((j-1) < pointsQE) && ((j-1) > -1) ) { QE[i-1][0][j-1] = QElambda[j-1]; QE[i-1][1][j-1] = qe; } if ( i > ct_NPixels) break; icount++; } if(icount/pointsQE < ct_NPixels){ error( "read_QE", "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_QE", "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 log("read_QE", "Done.\n"); } //!@} //!----------------------------------------------------------- // @name read_WC // // @desc read WC data // // @date thu 5 17:59:57 CEST 2002 //------------------------------------------------------------ // @function //!@{ void read_WC(void){ ifstream wcfile; char line[LINE_MAX_LENGTH]; int i; //------------------------------------------------------------ // Read Light Guides data // try to open the file log("read_WC", "Opening the file \"%s\" . . .\n", WC_FILE); wcfile.open( WC_FILE ); // if it is wrong or does not exist, exit if ( wcfile.bad() ) error( "read_WC", "Cannot open \"%s\". Exiting.\n", WC_FILE ); // read file log("read_WC", "Reading data . . .\n"); // get line from the file wcfile.getline(line, LINE_MAX_LENGTH); // get the number of datapoints sscanf(line, "%d", &pointsWC); // allocate memory for the table of QEs WC = new float * [2]; WC[0] = new float[pointsWC]; WC[1] = new float[pointsWC]; for ( i=0; iinumpixels = 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 pixary = new float* [2*ct_NCentralPixels]; for ( i=0; i<2*ct_NCentralPixels; ++i ) pixary[i] = new float[2]; pixneig = new int* [ct_NCentralPixels]; for ( i=0; i -1) && ((j-1) < pointsQE) && ((j-1) > -1) ) { QE[i-1][0][j-1] = QElambda[j-1]; QE[i-1][1][j-1] = qe; } 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 log("read_pixels", "Done.\n"); } //!@} //!----------------------------------------------------------- // @name read_ascii // // @desc read ascii configuration files used by the reflector // // @date tue dec 10 17:14:10 CET 2002 //------------------------------------------------------------ // @function //!@{ void read_ascii(FILE *sp, MMcConfigRunHeader *config) { Float_t radius = -1.0; Float_t focal = -1.0; Float_t stdfocal = -1.0; Float_t point = -1.0; Float_t stdpoint = -1.0; Float_t adjust = -1.0; Float_t spot = -1.0; Float_t camwidth = -1.0; Int_t imir; Float_t f,sx,sy,x,y,z,thetan,phin,xn,yn,zn; Float_t dx,dy; Int_t nummir, numref; Float_t wav,ref; Char_t token[40]; Char_t line[511]; Char_t flag; while(1){ if((flag=fgetc(sp))==EOF) break; fgets(&line[1],500,sp); line[0]=flag; if (line[0]== '#' && line[2]=='n' && line[3]=='u' && line[4]=='m' && line[5]== 'b' && line[12]=='d' && line[13]=='a' && line[14]=='t'){ fgets(line,500,sp); sscanf(line,"%i",&numref); for(int i=0; iInitSizeMirror(nummir); continue; } if (strstr(line, "r_mirror") == line){ sscanf(line,"%s %f",token,&radius); continue; } if (strstr(line, "define_mirrors") == line){ for(int i=0;iAddMirror(i); MGeomMirror &mirror = config->GetMirror(i); mirror.SetMirrorContent(imir,f,sx,sy,x,y,z,thetan,phin,xn,yn,zn); } fgets(line,500,sp); while(line[2]!='w' || line[3]!='i' || line[4]!='t' || line[5]!='h'){ fgets(line,500,sp); } fgets(line,500,sp); for(int i=0;iGetMirror(i); mirror.SetMirrorDeviations(dx,dy); } continue; } } config->SetMagicDef(radius, focal, stdfocal, point, stdpoint, adjust, spot, camwidth); } //!----------------------------------------------------------- // @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; double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5; 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] = 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.; /* 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 */ 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() ***************************************/ #define sqrt13 0.577350269 // = 1./sqrt(3.) #define sqrt3 1.732050807 // = sqrt(3.) int bpoint_is_in_pix(double dx, double dy, MGeomCamMagic *pgeomcam) { /* 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" */ /* a = length of one of the edges of one pixel, b = half the width of one pixel */ const int numN = pgeomcam->GetNumPixels(); for (int i=0; i=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 size_rotated // // @desc It rotates the NSB // // @date Tue Jan 7 17:09:25 CET 2003 // @function @code //------------------------------------------------------------ int size_rotated( float *rotated, float nsb[iMAXNUMPIX], float rho) { int r=0; float size_rotated; float Num_Pixels; float PixNum=0; float rho_pixel; int j=0; int k=0; for(int i=1; i= i){ r=ii-1; break; } if (ii<12){ Num_Pixels+=ii*6; PixNum=6*ii; // size_rotated = (360/PixNum); //rho_pixel = rho/size_rotated; } else { Num_Pixels+=(ii-6)*6; PixNum=(ii-6)*6; // size_rotated = (360/PixNum); // rho_pixel = rho/size_rotated; } } size_rotated = (360/PixNum); rho_pixel = rho/size_rotated; //Buscar j i k que no canvin de r!!! j=i-int(rho_pixel); //cout<<"j inicial "< 11 && r < 16){ if (j < MINPIXouter[r-12]) { j = i + (r - 6) * 6 - int(rho_pixel); } k=j-1; if ( k < MINPIXouter[r-12]) { k = MINPIXouter[r-11] - 1; } } rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k]; } return (int(rho_pixel)); } //------------------------------------------------------------ // @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 class MGeomCamMagic *camgeom, // the camera layout float minwl_nm, // the minimum accepted wavelength float maxwl_nm, // the maximum accepted wavelength class MTrigger *trigger, // the generated phes class MFadc *fadc, 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 uint i; static int k, ipixnum; static float cx, cy, wl, qe, t; static float cu, cv, cw; static MCCphoton photon; static float **qept; static char flag[SIZE_OF_FLAGS + 1]; static float radius_mm; // reset variables for ( i=0; iGetNumPixels(); ++i ){ nphe[i] = 0.0; } *incph = 0; radius_mm = camgeom->GetMaxRadius(); //- - - - - - - - - - - - - - - - - - - - - - - - - // read photons and "map" them into the pixels //-------------------------------------------------- // initialize CPhoton photon.fill(0., 0., 0., 0., 0., 0., 0., 0.); // read the photons data // loop over the photons while (1) { fread ( flag, SIZE_OF_FLAGS, 1, sp ); if (isA( flag, FLAG_END_OF_EVENT )) break; 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)++; // Chceck if photon is inside trigger time range t = photon.get_t() ; if (t-*tmin_ns>TOTAL_TRIGGER_TIME) { continue; } // // 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)*10 > radius_mm ) ){ continue; } ipixnum = bpoint_is_in_pix(cx*10, cy*10, camgeom); // -1 = the photon is in none of the pixels // 0 = the phton is in the central pixel, which is not used for trigger if (ipixnum==-1 || ipixnum==0) { 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])){ 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; // Apply incient angular correction due to Light Guides cu=photon.get_u(); cv=photon.get_v(); cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*3.14159/180.0; // find data point in the QE table (-> k) k = 0; while (k < pointsWC-1 && WC[0][k] < cw){ k++; } // correct the qe with WC data. qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw); // if random > quantum efficiency, reject it if ( (RandomNumber) > qe ) { continue; } //+++ // The photon has produced a photo electron //--- // cout << " accepted\n"; // increment the number of photoelectrons in the relevant pixel nphe[ipixnum] += 1.0; // store the new photoelectron fadc->Fill(ipixnum,(t-*tmin_ns) , trigger->FillShow(ipixnum,t-*tmin_ns), !((*camgeom)[ipixnum].GetR()>(*camgeom)[0].GetR())); *itotnphe += 1; } 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 MGeomCamMagic *camgeom, // camera layout float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function: // the NSB rates in phe/ns for each pixel ){ uint i, j; // counters int k, ii; // counters MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm); MFadc flashadc; 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 MCEventHeader evth_2; // the event header static float nphe[iMAXNUMPIX]; // the number of photoelectrons in each pixel int reflector_file_version; 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 char flag_new[4]; // open input file log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname ); infile = fopen( iname, "r" ); // check if the satrfield input file exists if ( infile == NULL ) { log( SIGNATURE, "Cannot open starfield input file: %s\n", iname ); log( SIGNATURE, "There is not NSB from the Stars\n"); return (0); } // get signature, and check it if((reflector_file_version=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 if (reflector_file_version > 5){ fread( flag_new, 4, 1, infile ); if(!isA( flag_new, FLAG_START_OF_HEADER)){ // We break the main loop cout<<"Warning: Expected start of run header flag, but found:"<GetNumPixels(); 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; iiGetNumPixels(); 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], MTrigger *trigger, MFadc *fadc, 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 float t; TRandom random; // Random numbers generator random.SetSeed ((UInt_t) (RandomNumber*1000)); ii = *inphe; // avoid dereferencing // check if the arrival times are set; if not generate them if(*atmin_ns *atmax_ns){ *atmin_ns = 0.; *atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME; } 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 and not too long if(simtime_ns< TOTAL_TRIGGER_TIME){ *atmin_ns = *atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2; *atmax_ns = *atmin_ns + TOTAL_TRIGGER_TIME; simtime_ns = TOTAL_TRIGGER_TIME; } if(simtime_ns> TOTAL_TRIGGER_TIME){ *atmax_ns =*atmin_ns + TOTAL_TRIGGER_TIME; simtime_ns = TOTAL_TRIGGER_TIME; } } // 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 if (inumnsbphe>0.0){ inumnsbphe = random.Poisson (inumnsbphe ); } // create the photoelectrons for(k=0; k < inumnsbphe; k++){ t=(RandomNumber * simtime_ns); (*fadc).Fill(j,t ,(*trigger).FillNSB(j,t)); 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 //=------------------------------------------------------------ //!@subsection Log of this file. //!@{ // // $Log: not supported by cvs2svn $ // Revision 1.55 2003/02/12 12:22:10 blanch // *** empty log message *** // // Revision 1.54 2003/02/12 11:55:01 blanch // *** empty log message *** // // Revision 1.53 2003/01/23 18:35:21 blanch // *** empty log message *** // // Revision 1.52 2003/01/20 17:19:20 blanch // It fills the MMcCorsikaRun. // // Revision 1.51 2003/01/14 13:40:17 blanch // MRawRunHeader::fNumEvents has been filled with the number of events in // this file. // Problems in fImpact computation have been solved. // Option to set a dc value to rise the discriminator threshold has been added. // // Revision 1.50 2003/01/07 16:33:31 blanch // Star Field Rotation has been implemented by Raul Orduna. Now there is a // rotation for each shower. It is done by a non enter pixel rotation assuming // a circular symetry of the camera. It is not exact but it is accurate enough and // much faster. // // Revision 1.49 2002/12/13 10:04:07 blanch // *** empty log message *** // // Revision 1.48 2002/12/12 17:40:50 blanch // *** empty log message *** // // Revision 1.47 2002/12/10 17:19:31 blanch // *** empty log message *** // // Revision 1.46 2002/11/08 17:51:00 blanch // *** empty log message *** // // Revision 1.45 2002/10/29 17:15:27 blanch // Reading several reflector versions. // // Revision 1.44 2002/10/18 16:53:03 blanch // Modification to read several reflector files. // // Revision 1.43 2002/09/13 10:53:39 blanch // Minor change to remove some undisired comments. // // Revision 1.42 2002/09/09 16:00:49 blanch // Statement has been included to avoid writting to disk MParContainer and MArray. // It has also been added the effect of the WC, the actual values must be added, // once they are measured. // // Revision 1.41 2002/09/04 09:57:42 blanch // Modifications done to use MGeomCam from MARS. // // Revision 1.40 2002/07/16 16:15:22 blanch // A first implementation of the Star field rotation has been done. // // Revision 1.39 2002/04/27 10:48:39 blanch // Some unused varibles have been removed. // // Revision 1.38 2002/03/18 18:44:29 blanch // Small modificatin to set the electronic Noise in the MMcTrigHeader class. // // Revision 1.37 2002/03/18 16:42:20 blanch // The data member fProductionSite of the MMcRunHeader has been set to 0, // which means that the prodution site is unknown. // // Revision 1.35 2002/03/15 15:15:52 blanch // Several modifications to simulate the actual trigger zone. // // Revision 1.34 2002/03/13 18:13:56 blanch // Some changes to fill correctly the new format of MMcRunHeader. // // Revision 1.33 2002/03/04 17:21:48 blanch // Small and not important changes. // // Revision 1.32 2002/02/28 15:04:52 blanch // A small back has been solved. Before, while not using the option // writte_all_images, not all triggered showers were stored. Now it is solved. // For that it is importatn taht the less restrictive trigger option is // checked first. // A new facility has been introduced and now one can choose the step size in // trigger loop mode for the discriminator threshold. // The close-compact topology for less than 3 pixels does not make sense. Before // the program was ignoring that, now it switch to simple neighbour condition. // // Revision 1.31 2002/01/18 17:41:02 blanch // The option of adding noise to all pixels or to not adding the noise // has been added. // We removed the pixels larger than 577. When there were more than one // trigger in one shower, the pixel number was increasing. Now it is // flagged by the variable MMcTrig::fFirstLvlTrig. // // Revision 1.30 2001/11/27 09:49:54 blanch // Fixing bug which was treating wrongly the extension of star photons. // // Revision 1.29 2001/11/14 17:38:23 blanch // Sveral changes have been done: // - bpoint_in_in_pixel has been dodified to speed up the program // - Some minor changes have been done to remove warnings // - buffer size and split version of the Branches have been removed // - Some modifications were needed to adat the program to the new // MRawEvtData::DeletePixels // // Revision 1.28 2001/10/26 16:31:45 blanch // Removing several warnings. // // Revision 1.27 2001/09/05 10:04:33 blanch // *** empty log message *** // // Revision 1.26 2001/07/19 09:29:53 blanch // Different threshold for each pixel can be used. // // Revision 1.25 2001/05/08 08:07:54 blanch // New numbering for branches from different trigger conditions has been // implemented. Now, they are calles: ClassName;1., ClasNema;2., ... // The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to // the RunHeaders tree. Also the MRawRunHeader is thera with some of its // information already filled. // // Revision 1.24 2001/04/06 16:48:09 magicsol // New camera version able to read the new format of the reflector output: // reflector 0.4 // // Revision 1.23 2001/03/28 16:13:41 blanch // While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has // been solved. // // Revision 1.22 2001/03/20 18:52:43 blanch // *** empty log message *** // // Revision 1.21 2001/03/19 19:53:03 blanch // Some print outs have been removed. // // Revision 1.20 2001/03/19 19:30:06 magicsol // Minor changes have been done to improve the FADC pedestals treatment. // // Revision 1.19 2001/03/05 11:14:41 magicsol // I changed the position of readinf a parameter. It is a minnor change. // // Revision 1.18 2001/03/05 10:36:52 blanch // A branch with information about the FADC simulation (MMcFadcHeader) is writen // in the McHeaders tree of the aoutput root file. // The NSB contribution is only applied if the the number of phe form the shower // are above a given number. // // Revision 1.17 2001/02/23 11:05:57 magicsol // Small changes due to slightly different output format and the introduction of // pedesals for teh FADC. // // Revision 1.16 2001/01/18 18:44:40 magicsol // *** empty log message *** // // Revision 1.15 2001/01/17 09:32:27 harald // The changes are neccessary to have the same name for trees in MC and in // data. So now there should be no differences in MC data and real data from // FADC system. // // Revision 1.14 2001/01/15 12:33:34 magicsol // Some modifications have been done to use the new (Dec'2000) Raw data format. // There are also some minnors modifications to adapt the improvements in the // MTrigger class (overlaping time and trigger cells). // // Revision 1.13 2000/10/25 08:14:23 magicsol // The routine that produce poisson random numbers to decide how many phe // form NSB are emmited in each pixel has been replaced. Now a ROOT routine // is used. // // Revision 1.12 2000/09/22 17:40:18 harald // Added a lot of changes done by oscar. // // Revision 1.11 2000/09/21 11:47:33 harald // Oscar found some smaller errors in the calculation of the pixel shape and // corrected it. // // $Log: not supported by cvs2svn $ // Revision 1.55 2003/02/12 12:22:10 blanch // *** empty log message *** // // Revision 1.54 2003/02/12 11:55:01 blanch // *** empty log message *** // // Revision 1.53 2003/01/23 18:35:21 blanch // *** empty log message *** // // Revision 1.52 2003/01/20 17:19:20 blanch // It fills the MMcCorsikaRun. // // Revision 1.51 2003/01/14 13:40:17 blanch // MRawRunHeader::fNumEvents has been filled with the number of events in // this file. // Problems in fImpact computation have been solved. // Option to set a dc value to rise the discriminator threshold has been added. // // Revision 1.50 2003/01/07 16:33:31 blanch // Star Field Rotation has been implemented by Raul Orduna. Now there is a // rotation for each shower. It is done by a non enter pixel rotation assuming // a circular symetry of the camera. It is not exact but it is accurate enough and // much faster. // // Revision 1.49 2002/12/13 10:04:07 blanch // *** empty log message *** // // Revision 1.48 2002/12/12 17:40:50 blanch // *** empty log message *** // // Revision 1.47 2002/12/10 17:19:31 blanch // *** empty log message *** // // Revision 1.46 2002/11/08 17:51:00 blanch // *** empty log message *** // // Revision 1.45 2002/10/29 17:15:27 blanch // Reading several reflector versions. // // Revision 1.44 2002/10/18 16:53:03 blanch // Modification to read several reflector files. // // Revision 1.43 2002/09/13 10:53:39 blanch // Minor change to remove some undisired comments. // // Revision 1.42 2002/09/09 16:00:49 blanch // Statement has been included to avoid writting to disk MParContainer and MArray. // It has also been added the effect of the WC, the actual values must be added, // once they are measured. // // Revision 1.41 2002/09/04 09:57:42 blanch // Modifications done to use MGeomCam from MARS. // // Revision 1.40 2002/07/16 16:15:22 blanch // A first implementation of the Star field rotation has been done. // // Revision 1.39 2002/04/27 10:48:39 blanch // Some unused varibles have been removed. // // Revision 1.38 2002/03/18 18:44:29 blanch // Small modificatin to set the electronic Noise in the MMcTrigHeader class. // // Revision 1.37 2002/03/18 16:42:20 blanch // The data member fProductionSite of the MMcRunHeader has been set to 0, // which means that the prodution site is unknown. // // Revision 1.35 2002/03/15 15:15:52 blanch // Several modifications to simulate the actual trigger zone. // // Revision 1.34 2002/03/13 18:13:56 blanch // Some changes to fill correctly the new format of MMcRunHeader. // // Revision 1.33 2002/03/04 17:21:48 blanch // Small and not important changes. // // Revision 1.32 2002/02/28 15:04:52 blanch // A small back has been solved. Before, while not using the option // writte_all_images, not all triggered showers were stored. Now it is solved. // For that it is importatn taht the less restrictive trigger option is // checked first. // A new facility has been introduced and now one can choose the step size in // trigger loop mode for the discriminator threshold. // The close-compact topology for less than 3 pixels does not make sense. Before // the program was ignoring that, now it switch to simple neighbour condition. // // Revision 1.31 2002/01/18 17:41:02 blanch // The option of adding noise to all pixels or to not adding the noise // has been added. // We removed the pixels larger than 577. When there were more than one // trigger in one shower, the pixel number was increasing. Now it is // flagged by the variable MMcTrig::fFirstLvlTrig. // // Revision 1.30 2001/11/27 09:49:54 blanch // Fixing bug which was treating wrongly the extension of star photons. // // Revision 1.29 2001/11/14 17:38:23 blanch // Sveral changes have been done: // - bpoint_in_in_pixel has been dodified to speed up the program // - Some minor changes have been done to remove warnings // - buffer size and split version of the Branches have been removed // - Some modifications were needed to adat the program to the new // MRawEvtData::DeletePixels // // Revision 1.28 2001/10/26 16:31:45 blanch // Removing several warnings. // // Revision 1.27 2001/09/05 10:04:33 blanch // *** empty log message *** // // Revision 1.26 2001/07/19 09:29:53 blanch // Different threshold for each pixel can be used. // // Revision 1.25 2001/05/08 08:07:54 blanch // New numbering for branches from different trigger conditions has been // implemented. Now, they are calles: ClassName;1., ClasNema;2., ... // The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to // the RunHeaders tree. Also the MRawRunHeader is thera with some of its // information already filled. // // Revision 1.24 2001/04/06 16:48:09 magicsol // New camera version able to read the new format of the reflector output: // reflector 0.4 // // Revision 1.23 2001/03/28 16:13:41 blanch // While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has // been solved. // // Revision 1.22 2001/03/20 18:52:43 blanch // *** empty log message *** // // Revision 1.21 2001/03/19 19:53:03 blanch // Some print outs have been removed. // // Revision 1.20 2001/03/19 19:30:06 magicsol // Minor changes have been done to improve the FADC pedestals treatment. // // Revision 1.19 2001/03/05 11:14:41 magicsol // I changed the position of readinf a parameter. It is a minnor change. // // Revision 1.18 2001/03/05 10:36:52 blanch // A branch with information about the FADC simulation (MMcFadcHeader) is writen // in the McHeaders tree of the aoutput root file. // The NSB contribution is only applied if the the number of phe form the shower // are above a given number. // // Revision 1.17 2001/02/23 11:05:57 magicsol // Small changes due to slightly different output format and the introduction of // pedesals for teh FADC. // // Revision 1.16 2001/01/18 18:44:40 magicsol // *** empty log message *** // // Revision 1.15 2001/01/17 09:32:27 harald // The changes are neccessary to have the same name for trees in MC and in // data. So now there should be no differences in MC data and real data from // FADC system. // // Revision 1.14 2001/01/15 12:33:34 magicsol // Some modifications have been done to use the new (Dec'2000) Raw data format. // There are also some minnors modifications to adapt the improvements in the // MTrigger class (overlaping time and trigger cells). // // Revision 1.13 2000/10/25 08:14:23 magicsol // The routine that produce poisson random numbers to decide how many phe // form NSB are emmited in each pixel has been replaced. Now a ROOT routine // is used. // // Revision 1.10 2000/07/04 14:10:20 MagicSol // Some changes have been done in the root output file. The RawEvt tree is only // stored in the single trigger mode. // The trigger input parameters are also given by the general input card. // The diffuse NSB and the star NSB have been decoupled. Now the contribution of // each one can be studied seperately. // // Revision 1.9 2000/06/13 13:25:24 blanch // The multiple arrays have been replaced, since they do not work // in alpha machines. Now we are using pointers and the command new // to allocate memory. // // Revision 1.8 2000/05/11 13:57:27 blanch // The option to loop over several trigger configurations has been included. // Some non-sense with triggertime range has been solved. // Montecarlo information and ADC counts are saved in a root file. // There was a problem with the maximum number of phe that the program could analyse. Now there is not limit. // // Revision 1.7 2000/03/24 18:10:46 blanch // A first FADC simulation and a trigger simulation are already implemented. // The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software. // A loop over trigger threshold and some corretcions in the time range where it looks for a trigger will be implemented as soon as possible. // // Revision 1.6 2000/03/20 18:35:11 blanch // The trigger is already implemented but it does not save the trigger information in any file as it is implemented in timecam. In the next days there will be a version which also creates the files with the trigger information. It is going to be a mixing of the current camera and timecam programs. // // Revision 1.5 2000/02/18 17:40:35 petry // 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. // // 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 // of invalid wavelength. This bug is now fixed and the range // of valid wavelengths extended to 290 - 800 nm. // This is in preparation for the NSB simulation to come. // Dirk // // Revision 1.2 1999/11/19 08:40:42 harald // Now it is possible to compile the camera programm under osf1. // // 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. // // 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