Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 5255) +++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (revision 5268) @@ -21,7 +21,7 @@ // // $RCSfile: camera.cxx,v $ -// $Revision: 1.72 $ +// $Revision: 1.73 $ // $Author: moralejo $ -// $Date: 2004-10-12 13:39:34 $ +// $Date: 2004-10-13 14:28:59 $ // //////////////////////////////////////////////////////////////////////// @@ -238,4 +238,6 @@ static float Spot_y=0.0; static float Spotsigma=0.0; + +static int CalibrationRun = 0; @@ -715,13 +717,15 @@ // get filenames - for(int ict=0;ictGetNumPixels();ui++){ - nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j]; - nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui]; - } + for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++) + { + nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + + zenfactor * nsbrate_phepns[ui][j]; + nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui]; + } } } @@ -1675,4 +1681,21 @@ } + + if (CalibrationRun) + { + DoCalibration(Fadc_CT, Trigger_CT, camgeom, nsb_trigresp, nsb_fadcresp, + &lons, &lons_outer, nsb_phepns, addElecNoise, + &EvtTree, EvtHeader, McEvt, EvtData); + + HeaderTree.Fill() ; + + outfile_temp.Write(); + outfile_temp.Close(); + + cout << endl << "Calibration run finished!" << endl << endl; + + return (0); + } + // // Read the reflector file with the Cherenkov data @@ -1681,28 +1704,26 @@ // select input file - if ( Data_From_STDIN ) { - + if ( Data_From_STDIN ) inputfile[0] = stdin; - } - else{ - for(int ict=0;ictFillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ; @@ -2594,5 +2615,5 @@ } } - // We don not count photons out of the camera. + // We do not count photons out of the camera. @@ -2985,5 +3006,6 @@ << SIGNATURE << '\n' << '\n' << "Processor of the reflector output\n" - << "J C Gonzalez, Jun 1998\n" + << "J. C. Gonzalez, Jun 1998\n" + << "O. Blanch, A. Moralejo, 2004\n" << "##################################################\n\n" << flush ; @@ -4417,4 +4439,364 @@ + +//------------------------------------------------------------- +// +// Function DoCalibration. A. Moralejo October 2004 +// +//------------------------------------------------------------- + +int DoCalibration(MFadc **Fadc_CT, MTrigger **Trigger_CT, + TObjArray camgeom, + float *nsb_trigresp, float *nsb_fadcresp, + MLons *lons, MLons *lons_outer, + float **nsb_phepns, int addElecNoise, + TTree *EvtTree, MRawEvtHeader **EvtHeader, + MMcEvt **McEvt, MRawEvtData** EvtData) +{ + + int ntshow = 0; + int lons_return; + + int *ntotphe = new int[ct_Number]; + int *ntotphot = new int[ct_Number]; + float *nphe_from_nsb = new float[ct_Number]; + + float timefirst = 0.; + float timelast = 0.; + + TArrayC *fadcValues; //@< the analog Fadc High gain signal for pixels + TArrayC *fadcValuesLow; //@< the analog Fadc Low gain signal for pixels + + // allocate space for FADC info + fadcValues = new TArrayC(FADC_slices_written); + fadcValuesLow = new TArrayC(FADC_slices_written); + + + // allocate space for PMTs numbers of pixels + float *pheinpix = new float [ct_NPixels]; + + for (int calevent = 0; calevent < 1; calevent++) + { + // + // Clear Trigger and Fadc + // + for(int ict = 0; ict < ct_Number; ict++) + { + Trigger_CT[ict]->Reset() ; + Trigger_CT[ict]->ClearFirst(); + Trigger_CT[ict]->ClearZero(); + Fadc_CT[ict]->Reset() ; + + ntotphe[ict] = ntotphot[ict] = 0; + nphe_from_nsb[ict] = 0.; + + } + + ntshow++; + + if((ntshow+1)%100 == 1) + log(SIGNATURE, "Event %d\n", ntshow); + + // Produce the photoelectrons + float phot_per_pix = 100.; + for(int ict = 0; ict < ct_Number; ict++) + produce_calib_phes( (MGeomCam*)(camgeom.UncheckedAt(ict)), + phot_per_pix, + Trigger_CT[ict], + Fadc_CT[ict], + &(ntotphe[ict]), + pheinpix, + &(ntotphot[ict]), + ict); + + // NSB simulation + + if(lons && lons_outer) + { + // Fill trigger and fadc response in the trigger class from the NSB database + + for (int ict = 0; ict < ct_Number; ict++) + { + for (UInt_t ui = 0; + ui < ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); + ui++) + { + nphe_from_nsb[ict] += nsb_phepns[ict][ui]; + + if (nsb_phepns[ict][ui] > 0.0) + { + if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD() > + (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD()) + lons_return = lons_outer->GetResponse(nsb_phepns[ict][ui],0.01, + & nsb_trigresp[0], + & nsb_fadcresp[0]); + else + lons_return = lons->GetResponse(nsb_phepns[ict][ui],0.01, + & nsb_trigresp[0], + & nsb_fadcresp[0]); + + if (lons_return == 0) + { + cout << "Exiting.\n"; + exit(1); + } + + Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp); + } + } + } + + }// end if(simulateNSB) ... + + + //++++++++++++++++++++++++++++++++++++++++++++++++++ + // at this point we have a camera full of + // ph.e.s + //-------------------------------------------------- + + // now the noise of the electronic + // (preamps, optical transmission,..) is introduced. + // + + for (int ict = 0; ict < ct_Number; ict++) + { + if (addElecNoise) + Fadc_CT[ict]->ElecNoise(); + + // now a shift in the fadc signal due to the pedestals is + // introduced + // This is done inside the class MFadc by the method Pedestals + + Fadc_CT[ict]->Pedestals(); + + // Set the trigger time. The 3 ns are such that the pulse appears well + // centered if the FADC range is the usual 50 ns. + // If you want to shift the pulse position, do not change the value here, + // use the option trigger_delay in the input card instead. + + Fadc_CT[ict]->TriggeredFadc(3.); + + // + // Fill the event header information + // + EvtHeader[ict]->FillHeader((UInt_t) calevent, 0); + + McEvt[ict]->Fill( 0, 0, 0., -1.0, -1.0, -1.0, 0., 0., 0., 0., 0., + 0., 0., 0., timefirst, timelast, 0., 0., 0., + 0., 0., 0., 0., 0, 0, 0, + (UInt_t) ntotphot[ict], + (UInt_t) ntotphe[ict], + (UInt_t) nphe_from_nsb[ict]+ ntotphe[ict], + 0., 0., 0.); + + // + // Fill the FADC information + // + for(UInt_t ui=0; + ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); + ui++) + { + for (int jslice = 0; jslice < FADC_slices_written; jslice++) + { + fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(ui, jslice),jslice); + fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(ui,jslice),jslice); + } + EvtData[ict]->AddPixel(ui,fadcValues,0); + EvtData[ict]->AddPixel(ui,fadcValuesLow,kTRUE); + } + } + + EvtTree->Fill(); + + // clear all + for(int ict=0;ictClear() ; + EvtData[ict]->ResetPixels(0,0); + McEvt[ict]->Clear() ; + } + } + + return(0); +} + + +//------------------------------------------------------------ +// +// Function produce_calib_phes, A. Moralejo Oct 2004 +// +//------------------------------------------------------------ + +int produce_calib_phes( MGeomCam *camgeom, // The camera layout + float phot_per_pix, // Average number of photons per inner pixel + MTrigger *trigger, + MFadc *fadc, + int *itotnphe, // total number of produced photoelectrons + float *nphe, // number of photoelectrons in each pixel + int *nphot, // total number of photons in all pixels + int ict // Telescope that is being analised to get the right QE. + ) +{ + + int ipixnum; + float cx, cy, wl, time, phi, phi_deg, qe; + float **qept; + float radius_mm, focal_dist_mm; + + // reset variables + + for ( uint i = 0; i < camgeom->GetNumPixels(); i++ ) + nphe[i] = 0.0; + + *itotnphe = 0; + *nphot = 0; + + TRandom random; + random.SetSeed((Int_t)(get_seeds(0)*get_seeds(1))); + + // + // Create photons and "map" them into the pixels + // + + // Maximum radius of camera: + radius_mm = camgeom->GetMaxRadius(); + + // Distance from center of mirror dish to camera plane: + focal_dist_mm = camgeom->GetCameraDist()*1000; + + // Cosine of the angle between telescope axis and line from center of mirror + // dish to the edge of the camera: + float cos_epsilon_max = cos(atan2(radius_mm, focal_dist_mm)); + + + // Calculate total number of photons to be produced. + int total_photons = (int) (phot_per_pix * 3.14159265 * radius_mm * radius_mm / + (*camgeom)[0].GetA()); + + // loop over the photons + + for (int iph = 0; iph < total_photons; iph++) + { + + time = 0.; + + // Obtain photon coordinates on the camera. We assume a point source of light placed + // in the center of the mirror dish. + + // polar angle + float psi = RandomNumber * 2 * 3.14159265; + // angle between the telescope axis and the photon trajectory. + float epsilon = acos(1.-RandomNumber*(1.-cos_epsilon_max)); + float tanepsilon = tan(epsilon); + + cx = focal_dist_mm*tanepsilon*cos(psi); // mm + cy = focal_dist_mm*tanepsilon*sin(psi); // mm + + // Angle between photon trajectory and camera plane: + phi = 3.14159265/2.-epsilon; // rad + + // wavelength + + wl = 500.; // nm + + + if( (wl > WAVEBANDBOUND6) || (wl < WAVEBANDBOUND1) || + (sqrt(cx*cx + cy*cy) > radius_mm ) ) + continue; + + // + // Pixelization + // + ipixnum = bpoint_is_in_pix(cx, cy, camgeom); + + // -1 = the photon is in none of the pixels + // 0 = the phton is in the central pixel, which is not used + if (ipixnum==-1 || ipixnum==0) + continue; + + // increase number of photons within pixels + *nphot += 1; + + // + // QE simulation + // + // set pointer to the QE table of the relevant pixel + + qept = (float **)QE[ict][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[ict]-1])) + continue; + + + // find data point in the QE table (-> k) + int k = 1; // start at 1 because the condition was already tested for 0 + while (k < pointsQE[ict]-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, plexiglas, + // 1st dynode collection efficiency, double crossings... etc. + // This information is contained in the file Data/LightCollection.dat, + // and has been read into the array WC (which stands for "Winston Cones") + // + phi_deg = phi*180./3.14159265; + + // find data point in the WC table (-> k) + + if ( camgeom->GetPixRatio(ipixnum) < 1.) // => Pixel is outer pixel + { + k = 0; + while (k < pointsWC-1 && WC_outer[0][k] < phi_deg) + k++; + // correct the qe with WC data. + qe = qe*lin_interpol(WC_outer[0][k-1], WC_outer[1][k-1], + WC_outer[0][k], WC_outer[1][k], phi_deg); + } + + else // => Pixel is inner pixel + { + k = 0; + while (k < pointsWC-1 && WC[0][k] < phi_deg) + 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], phi_deg); + } + + + // if random > quantum efficiency, reject it + + if ( (RandomNumber) > qe ) { + continue; + } + + // + // The photon has produced a photo electron + // + + // increment the number of photoelectrons in the relevant pixel + + nphe[ipixnum] += 1.; + + // store the new photoelectron + + fadc->Fill(ipixnum,(time) , + trigger->FillShow(ipixnum, time), + !((*camgeom)[ipixnum].GetD()>(*camgeom)[0].GetD())); + + *itotnphe += 1; + } + + return(0); + +} + + // @endcode @@ -4426,4 +4808,79 @@ // // $Log: not supported by cvs2svn $ +// Revision 1.72 2004/10/12 13:39:34 moralejo +// +// Lots of changes intended to make it possible to select the FADC sampling +// frequency from the input card, through the command fadc_GHz. The most +// important ones are the following: +// +// - Removed FADC_SLICES de Mars/mmc/Mdefine.h +// Already defined in MFadcDefine.h! +// +// - Replaced fixed numbers in array dimensions of starresponse.cxx +// +// - Added in MFadc.cxx and MFadc.hxx (Float_t) casts to initializations +// of single phe response array +// +// - IMPORTANT: Fixed MFadc::GetResponse -> The returned single phe response +// had only RESPONSE_SLICES (which is actually for the trigger branch), +// whereas it should have RESPONSE_SLICES_MFADC. Fixed the same confusion +// in other two points of the code (filling of the FADC for the signal), +// in Fill() and FillOuter(). +// +// - RESPONSE_SLICES_MFADC is now eliminated, since this quantity is now +// decided by MFadc depending on other parameters. +// +// - Fixed problem in starresponse.cxx due to which the histograms fadcresp +// and fadcbase in the root output were actually identical. +// +// - Changed starresponse.cxx such that above 1 phe/ns/pixel the precision +// of the database is forced to be 0.1, and below 1, to 0.01 phe/ns/pixel +// +// - In Camera/creadparam.cxx: Now unkown tokens cause the program to stop, +// instead of being simply skipped as it was until now. +// +// - Fixed error in MStarLight::StoreHisto. TH1::SetBinContent uses bin numbers +// from 1 to nbins (in old code 0 to nbins-1 was assumed). +// +// - In MLons.cxx: use memcpy to copy pieces of the database into the FADC and +// trigger branches, to (hopefully) speed up execution. For this I had to add +// 2 getter functions in MStarLight.hxx +// +// - Everywhere: changed the shape parameter for trigger and FADC response to +// be an Int. Changed version of NSB database from 1002 to 1003. +// +// - In MTrigger.cxx, changed all initializations of SlicesFirst and +// SlicesSecond to 0 (instead of -50 as it was before). This controls the +// time of trigger. If no trigger happened (like when making pedestal files) +// the time was negative and the array index used to retrieve the noise from +// the database was negative, resulting in "discontinuities" in the noise +// (half-photoelectrons for instance). +// +// - In MStarLight changed fBinsTrig and fBinsFadc from Float_t to Int_t +// +// - Replace WIDTH_FADC_TIMESLICE by FADC_SLICES_PER_NSEC (which is the +// inverse of the former) +// +// - Replaced SLICES_PER_NSEC by TRIG_SLICES_PER_NSEC +// +// - TRIGBINS eliminated. It depends on other two values +// TRIGBINS = TIMERANGE*TRIG_SLICES_PER_NSEC +// +// - FADCBINS eliminated. It depends on other two values +// FADCBINS = TIMERANGE*FADC_SLICES_PER_NSEC +// +// - MTriggerDefine.h Changed RESPONSE_SLICES to RESPONSE_SLICES_TRIG +// +// - Added to the MLons constructor an argument regarding the FADC sampling +// frequency +// +// - MFadc: now the number of response slices for the FADC simulation is +// decided by the program according to the other parameters. +// +// - MStarLight: Added arguments (fadc_slices_per_ns, response_slices_fadc) +// to constructor. +// +// - creadparam.cxx, camera.cxx Changed default value of digital_noise to 0. +// // Revision 1.71 2004/09/17 09:20:52 moralejo //