Changeset 5268 for trunk/MagicSoft/Simulation/Detector/Camera
- Timestamp:
- 10/13/04 15:28:59 (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
r5247 r5268 21 21 // 22 22 // $RCSfile: camera.cxx,v $ 23 // $Revision: 1.7 2$23 // $Revision: 1.73 $ 24 24 // $Author: moralejo $ 25 // $Date: 2004-10-1 2 13:39:34$25 // $Date: 2004-10-13 14:28:59 $ 26 26 // 27 27 //////////////////////////////////////////////////////////////////////// … … 238 238 static float Spot_y=0.0; 239 239 static float Spotsigma=0.0; 240 241 static int CalibrationRun = 0; 240 242 241 243 … … 715 717 // get filenames 716 718 717 for(int ict=0;ict<ct_Number;ict++) 718 { 719 strcpy( inname_CT[ict], get_input_filename(ict) ); 720 if (strlen(inname_CT[ict]) == 0) 721 { 722 printf("\nError: missing input file name for CT id = %d. Exiting.\n\n", ict); 723 exit(1); 724 } 725 } 719 CalibrationRun = 1; 720 if (! CalibrationRun) 721 for(int ict=0;ict<ct_Number;ict++) 722 { 723 strcpy( inname_CT[ict], get_input_filename(ict) ); 724 if (strlen(inname_CT[ict]) == 0) 725 { 726 printf("\nError: missing input file name for CT id = %d. Exiting.\n\n", ict); 727 exit(1); 728 } 729 } 726 730 727 731 strcpy( starfieldname, get_starfield_filename() ); … … 1548 1552 zenfactor = pow(10., -0.4 * ext[j] ); 1549 1553 1550 for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){ 1551 nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j]; 1552 nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui]; 1553 } 1554 for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++) 1555 { 1556 nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + 1557 zenfactor * nsbrate_phepns[ui][j]; 1558 nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui]; 1559 } 1554 1560 } 1555 1561 } … … 1675 1681 } 1676 1682 1683 1684 if (CalibrationRun) 1685 { 1686 DoCalibration(Fadc_CT, Trigger_CT, camgeom, nsb_trigresp, nsb_fadcresp, 1687 &lons, &lons_outer, nsb_phepns, addElecNoise, 1688 &EvtTree, EvtHeader, McEvt, EvtData); 1689 1690 HeaderTree.Fill() ; 1691 1692 outfile_temp.Write(); 1693 outfile_temp.Close(); 1694 1695 cout << endl << "Calibration run finished!" << endl << endl; 1696 1697 return (0); 1698 } 1699 1677 1700 // 1678 1701 // Read the reflector file with the Cherenkov data … … 1681 1704 // select input file 1682 1705 1683 if ( Data_From_STDIN ) { 1684 1706 if ( Data_From_STDIN ) 1685 1707 inputfile[0] = stdin; 1686 1708 1687 } 1688 else{ 1689 for(int ict=0;ict<ct_Number;ict++){ 1690 1691 log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] ); 1692 inputfile[ict] = fopen( inname_CT[ict], "r" ); 1693 1694 if ( inputfile[ict] == NULL ) 1695 error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] ); 1696 } 1697 1698 } 1709 else 1710 { 1711 for(int ict = 0; ict < ct_Number; ict++) 1712 { 1713 log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] ); 1714 inputfile[ict] = fopen( inname_CT[ict], "r" ); 1715 1716 if ( inputfile[ict] == NULL ) 1717 error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] ); 1718 } 1719 } 1699 1720 1700 1721 // get signature, and check it 1701 1722 1702 for(int ict=0;ict<ct_Number;ict++) {1703 if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE){1704 exit(1);1705 } 1706 }1723 for(int ict=0;ict<ct_Number;ict++) 1724 { 1725 if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE) 1726 exit(1); 1727 } 1707 1728 1708 1729 // open data file … … 2463 2484 // Fill the header of this event 2464 2485 // 2465 2486 2466 2487 EvtHeader[ict] 2467 2488 ->FillHeader ( (UInt_t) (ntshow + nshow) , 20 ) ; … … 2594 2615 } 2595 2616 } 2596 // We do nnot count photons out of the camera.2617 // We do not count photons out of the camera. 2597 2618 2598 2619 … … 2985 3006 << SIGNATURE << '\n' << '\n' 2986 3007 << "Processor of the reflector output\n" 2987 << "J C Gonzalez, Jun 1998\n" 3008 << "J. C. Gonzalez, Jun 1998\n" 3009 << "O. Blanch, A. Moralejo, 2004\n" 2988 3010 << "##################################################\n\n" 2989 3011 << flush ; … … 4417 4439 4418 4440 4441 4442 //------------------------------------------------------------- 4443 // 4444 // Function DoCalibration. A. Moralejo October 2004 4445 // 4446 //------------------------------------------------------------- 4447 4448 int DoCalibration(MFadc **Fadc_CT, MTrigger **Trigger_CT, 4449 TObjArray camgeom, 4450 float *nsb_trigresp, float *nsb_fadcresp, 4451 MLons *lons, MLons *lons_outer, 4452 float **nsb_phepns, int addElecNoise, 4453 TTree *EvtTree, MRawEvtHeader **EvtHeader, 4454 MMcEvt **McEvt, MRawEvtData** EvtData) 4455 { 4456 4457 int ntshow = 0; 4458 int lons_return; 4459 4460 int *ntotphe = new int[ct_Number]; 4461 int *ntotphot = new int[ct_Number]; 4462 float *nphe_from_nsb = new float[ct_Number]; 4463 4464 float timefirst = 0.; 4465 float timelast = 0.; 4466 4467 TArrayC *fadcValues; //@< the analog Fadc High gain signal for pixels 4468 TArrayC *fadcValuesLow; //@< the analog Fadc Low gain signal for pixels 4469 4470 // allocate space for FADC info 4471 fadcValues = new TArrayC(FADC_slices_written); 4472 fadcValuesLow = new TArrayC(FADC_slices_written); 4473 4474 4475 // allocate space for PMTs numbers of pixels 4476 float *pheinpix = new float [ct_NPixels]; 4477 4478 for (int calevent = 0; calevent < 1; calevent++) 4479 { 4480 // 4481 // Clear Trigger and Fadc 4482 // 4483 for(int ict = 0; ict < ct_Number; ict++) 4484 { 4485 Trigger_CT[ict]->Reset() ; 4486 Trigger_CT[ict]->ClearFirst(); 4487 Trigger_CT[ict]->ClearZero(); 4488 Fadc_CT[ict]->Reset() ; 4489 4490 ntotphe[ict] = ntotphot[ict] = 0; 4491 nphe_from_nsb[ict] = 0.; 4492 4493 } 4494 4495 ntshow++; 4496 4497 if((ntshow+1)%100 == 1) 4498 log(SIGNATURE, "Event %d\n", ntshow); 4499 4500 // Produce the photoelectrons 4501 float phot_per_pix = 100.; 4502 for(int ict = 0; ict < ct_Number; ict++) 4503 produce_calib_phes( (MGeomCam*)(camgeom.UncheckedAt(ict)), 4504 phot_per_pix, 4505 Trigger_CT[ict], 4506 Fadc_CT[ict], 4507 &(ntotphe[ict]), 4508 pheinpix, 4509 &(ntotphot[ict]), 4510 ict); 4511 4512 // NSB simulation 4513 4514 if(lons && lons_outer) 4515 { 4516 // Fill trigger and fadc response in the trigger class from the NSB database 4517 4518 for (int ict = 0; ict < ct_Number; ict++) 4519 { 4520 for (UInt_t ui = 0; 4521 ui < ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); 4522 ui++) 4523 { 4524 nphe_from_nsb[ict] += nsb_phepns[ict][ui]; 4525 4526 if (nsb_phepns[ict][ui] > 0.0) 4527 { 4528 if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD() > 4529 (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD()) 4530 lons_return = lons_outer->GetResponse(nsb_phepns[ict][ui],0.01, 4531 & nsb_trigresp[0], 4532 & nsb_fadcresp[0]); 4533 else 4534 lons_return = lons->GetResponse(nsb_phepns[ict][ui],0.01, 4535 & nsb_trigresp[0], 4536 & nsb_fadcresp[0]); 4537 4538 if (lons_return == 0) 4539 { 4540 cout << "Exiting.\n"; 4541 exit(1); 4542 } 4543 4544 Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp); 4545 } 4546 } 4547 } 4548 4549 }// end if(simulateNSB) ... 4550 4551 4552 //++++++++++++++++++++++++++++++++++++++++++++++++++ 4553 // at this point we have a camera full of 4554 // ph.e.s 4555 //-------------------------------------------------- 4556 4557 // now the noise of the electronic 4558 // (preamps, optical transmission,..) is introduced. 4559 // 4560 4561 for (int ict = 0; ict < ct_Number; ict++) 4562 { 4563 if (addElecNoise) 4564 Fadc_CT[ict]->ElecNoise(); 4565 4566 // now a shift in the fadc signal due to the pedestals is 4567 // introduced 4568 // This is done inside the class MFadc by the method Pedestals 4569 4570 Fadc_CT[ict]->Pedestals(); 4571 4572 // Set the trigger time. The 3 ns are such that the pulse appears well 4573 // centered if the FADC range is the usual 50 ns. 4574 // If you want to shift the pulse position, do not change the value here, 4575 // use the option trigger_delay in the input card instead. 4576 4577 Fadc_CT[ict]->TriggeredFadc(3.); 4578 4579 // 4580 // Fill the event header information 4581 // 4582 EvtHeader[ict]->FillHeader((UInt_t) calevent, 0); 4583 4584 McEvt[ict]->Fill( 0, 0, 0., -1.0, -1.0, -1.0, 0., 0., 0., 0., 0., 4585 0., 0., 0., timefirst, timelast, 0., 0., 0., 4586 0., 0., 0., 0., 0, 0, 0, 4587 (UInt_t) ntotphot[ict], 4588 (UInt_t) ntotphe[ict], 4589 (UInt_t) nphe_from_nsb[ict]+ ntotphe[ict], 4590 0., 0., 0.); 4591 4592 // 4593 // Fill the FADC information 4594 // 4595 for(UInt_t ui=0; 4596 ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); 4597 ui++) 4598 { 4599 for (int jslice = 0; jslice < FADC_slices_written; jslice++) 4600 { 4601 fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(ui, jslice),jslice); 4602 fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(ui,jslice),jslice); 4603 } 4604 EvtData[ict]->AddPixel(ui,fadcValues,0); 4605 EvtData[ict]->AddPixel(ui,fadcValuesLow,kTRUE); 4606 } 4607 } 4608 4609 EvtTree->Fill(); 4610 4611 // clear all 4612 for(int ict=0;ict<ct_Number;ict++) 4613 { 4614 EvtHeader[ict]->Clear() ; 4615 EvtData[ict]->ResetPixels(0,0); 4616 McEvt[ict]->Clear() ; 4617 } 4618 } 4619 4620 return(0); 4621 } 4622 4623 4624 //------------------------------------------------------------ 4625 // 4626 // Function produce_calib_phes, A. Moralejo Oct 2004 4627 // 4628 //------------------------------------------------------------ 4629 4630 int produce_calib_phes( MGeomCam *camgeom, // The camera layout 4631 float phot_per_pix, // Average number of photons per inner pixel 4632 MTrigger *trigger, 4633 MFadc *fadc, 4634 int *itotnphe, // total number of produced photoelectrons 4635 float *nphe, // number of photoelectrons in each pixel 4636 int *nphot, // total number of photons in all pixels 4637 int ict // Telescope that is being analised to get the right QE. 4638 ) 4639 { 4640 4641 int ipixnum; 4642 float cx, cy, wl, time, phi, phi_deg, qe; 4643 float **qept; 4644 float radius_mm, focal_dist_mm; 4645 4646 // reset variables 4647 4648 for ( uint i = 0; i < camgeom->GetNumPixels(); i++ ) 4649 nphe[i] = 0.0; 4650 4651 *itotnphe = 0; 4652 *nphot = 0; 4653 4654 TRandom random; 4655 random.SetSeed((Int_t)(get_seeds(0)*get_seeds(1))); 4656 4657 // 4658 // Create photons and "map" them into the pixels 4659 // 4660 4661 // Maximum radius of camera: 4662 radius_mm = camgeom->GetMaxRadius(); 4663 4664 // Distance from center of mirror dish to camera plane: 4665 focal_dist_mm = camgeom->GetCameraDist()*1000; 4666 4667 // Cosine of the angle between telescope axis and line from center of mirror 4668 // dish to the edge of the camera: 4669 float cos_epsilon_max = cos(atan2(radius_mm, focal_dist_mm)); 4670 4671 4672 // Calculate total number of photons to be produced. 4673 int total_photons = (int) (phot_per_pix * 3.14159265 * radius_mm * radius_mm / 4674 (*camgeom)[0].GetA()); 4675 4676 // loop over the photons 4677 4678 for (int iph = 0; iph < total_photons; iph++) 4679 { 4680 4681 time = 0.; 4682 4683 // Obtain photon coordinates on the camera. We assume a point source of light placed 4684 // in the center of the mirror dish. 4685 4686 // polar angle 4687 float psi = RandomNumber * 2 * 3.14159265; 4688 // angle between the telescope axis and the photon trajectory. 4689 float epsilon = acos(1.-RandomNumber*(1.-cos_epsilon_max)); 4690 float tanepsilon = tan(epsilon); 4691 4692 cx = focal_dist_mm*tanepsilon*cos(psi); // mm 4693 cy = focal_dist_mm*tanepsilon*sin(psi); // mm 4694 4695 // Angle between photon trajectory and camera plane: 4696 phi = 3.14159265/2.-epsilon; // rad 4697 4698 // wavelength 4699 4700 wl = 500.; // nm 4701 4702 4703 if( (wl > WAVEBANDBOUND6) || (wl < WAVEBANDBOUND1) || 4704 (sqrt(cx*cx + cy*cy) > radius_mm ) ) 4705 continue; 4706 4707 // 4708 // Pixelization 4709 // 4710 ipixnum = bpoint_is_in_pix(cx, cy, camgeom); 4711 4712 // -1 = the photon is in none of the pixels 4713 // 0 = the phton is in the central pixel, which is not used 4714 if (ipixnum==-1 || ipixnum==0) 4715 continue; 4716 4717 // increase number of photons within pixels 4718 *nphot += 1; 4719 4720 // 4721 // QE simulation 4722 // 4723 // set pointer to the QE table of the relevant pixel 4724 4725 qept = (float **)QE[ict][ipixnum]; 4726 4727 // check if wl is inside table; outside the table, QE is assumed to be zero 4728 4729 if((wl < qept[0][0]) || (wl > qept[0][pointsQE[ict]-1])) 4730 continue; 4731 4732 4733 // find data point in the QE table (-> k) 4734 int k = 1; // start at 1 because the condition was already tested for 0 4735 while (k < pointsQE[ict]-1 && qept[0][k] < wl) 4736 k++; 4737 4738 // calculate the qe value between 0. and 1. 4739 4740 qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0; 4741 4742 4743 // 4744 // Apply incient angular correction due to Light Guides, plexiglas, 4745 // 1st dynode collection efficiency, double crossings... etc. 4746 // This information is contained in the file Data/LightCollection.dat, 4747 // and has been read into the array WC (which stands for "Winston Cones") 4748 // 4749 phi_deg = phi*180./3.14159265; 4750 4751 // find data point in the WC table (-> k) 4752 4753 if ( camgeom->GetPixRatio(ipixnum) < 1.) // => Pixel is outer pixel 4754 { 4755 k = 0; 4756 while (k < pointsWC-1 && WC_outer[0][k] < phi_deg) 4757 k++; 4758 // correct the qe with WC data. 4759 qe = qe*lin_interpol(WC_outer[0][k-1], WC_outer[1][k-1], 4760 WC_outer[0][k], WC_outer[1][k], phi_deg); 4761 } 4762 4763 else // => Pixel is inner pixel 4764 { 4765 k = 0; 4766 while (k < pointsWC-1 && WC[0][k] < phi_deg) 4767 k++; 4768 // correct the qe with WC data. 4769 qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], phi_deg); 4770 } 4771 4772 4773 // if random > quantum efficiency, reject it 4774 4775 if ( (RandomNumber) > qe ) { 4776 continue; 4777 } 4778 4779 // 4780 // The photon has produced a photo electron 4781 // 4782 4783 // increment the number of photoelectrons in the relevant pixel 4784 4785 nphe[ipixnum] += 1.; 4786 4787 // store the new photoelectron 4788 4789 fadc->Fill(ipixnum,(time) , 4790 trigger->FillShow(ipixnum, time), 4791 !((*camgeom)[ipixnum].GetD()>(*camgeom)[0].GetD())); 4792 4793 *itotnphe += 1; 4794 } 4795 4796 return(0); 4797 4798 } 4799 4800 4419 4801 // @endcode 4420 4802 … … 4426 4808 // 4427 4809 // $Log: not supported by cvs2svn $ 4810 // Revision 1.72 2004/10/12 13:39:34 moralejo 4811 // 4812 // Lots of changes intended to make it possible to select the FADC sampling 4813 // frequency from the input card, through the command fadc_GHz. The most 4814 // important ones are the following: 4815 // 4816 // - Removed FADC_SLICES de Mars/mmc/Mdefine.h 4817 // Already defined in MFadcDefine.h! 4818 // 4819 // - Replaced fixed numbers in array dimensions of starresponse.cxx 4820 // 4821 // - Added in MFadc.cxx and MFadc.hxx (Float_t) casts to initializations 4822 // of single phe response array 4823 // 4824 // - IMPORTANT: Fixed MFadc::GetResponse -> The returned single phe response 4825 // had only RESPONSE_SLICES (which is actually for the trigger branch), 4826 // whereas it should have RESPONSE_SLICES_MFADC. Fixed the same confusion 4827 // in other two points of the code (filling of the FADC for the signal), 4828 // in Fill() and FillOuter(). 4829 // 4830 // - RESPONSE_SLICES_MFADC is now eliminated, since this quantity is now 4831 // decided by MFadc depending on other parameters. 4832 // 4833 // - Fixed problem in starresponse.cxx due to which the histograms fadcresp 4834 // and fadcbase in the root output were actually identical. 4835 // 4836 // - Changed starresponse.cxx such that above 1 phe/ns/pixel the precision 4837 // of the database is forced to be 0.1, and below 1, to 0.01 phe/ns/pixel 4838 // 4839 // - In Camera/creadparam.cxx: Now unkown tokens cause the program to stop, 4840 // instead of being simply skipped as it was until now. 4841 // 4842 // - Fixed error in MStarLight::StoreHisto. TH1::SetBinContent uses bin numbers 4843 // from 1 to nbins (in old code 0 to nbins-1 was assumed). 4844 // 4845 // - In MLons.cxx: use memcpy to copy pieces of the database into the FADC and 4846 // trigger branches, to (hopefully) speed up execution. For this I had to add 4847 // 2 getter functions in MStarLight.hxx 4848 // 4849 // - Everywhere: changed the shape parameter for trigger and FADC response to 4850 // be an Int. Changed version of NSB database from 1002 to 1003. 4851 // 4852 // - In MTrigger.cxx, changed all initializations of SlicesFirst and 4853 // SlicesSecond to 0 (instead of -50 as it was before). This controls the 4854 // time of trigger. If no trigger happened (like when making pedestal files) 4855 // the time was negative and the array index used to retrieve the noise from 4856 // the database was negative, resulting in "discontinuities" in the noise 4857 // (half-photoelectrons for instance). 4858 // 4859 // - In MStarLight changed fBinsTrig and fBinsFadc from Float_t to Int_t 4860 // 4861 // - Replace WIDTH_FADC_TIMESLICE by FADC_SLICES_PER_NSEC (which is the 4862 // inverse of the former) 4863 // 4864 // - Replaced SLICES_PER_NSEC by TRIG_SLICES_PER_NSEC 4865 // 4866 // - TRIGBINS eliminated. It depends on other two values 4867 // TRIGBINS = TIMERANGE*TRIG_SLICES_PER_NSEC 4868 // 4869 // - FADCBINS eliminated. It depends on other two values 4870 // FADCBINS = TIMERANGE*FADC_SLICES_PER_NSEC 4871 // 4872 // - MTriggerDefine.h Changed RESPONSE_SLICES to RESPONSE_SLICES_TRIG 4873 // 4874 // - Added to the MLons constructor an argument regarding the FADC sampling 4875 // frequency 4876 // 4877 // - MFadc: now the number of response slices for the FADC simulation is 4878 // decided by the program according to the other parameters. 4879 // 4880 // - MStarLight: Added arguments (fadc_slices_per_ns, response_slices_fadc) 4881 // to constructor. 4882 // 4883 // - creadparam.cxx, camera.cxx Changed default value of digital_noise to 0. 4884 // 4428 4885 // Revision 1.71 2004/09/17 09:20:52 moralejo 4429 4886 //
Note:
See TracChangeset
for help on using the changeset viewer.