Changeset 1078 for trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
- Timestamp:
- 11/14/01 17:38:23 (23 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
r1001 r1078 21 21 // 22 22 // $RCSfile: camera.cxx,v $ 23 // $Revision: 1.2 8$23 // $Revision: 1.29 $ 24 24 // $Author: blanch $ 25 // $Date: 2001-1 0-26 16:31:45$25 // $Date: 2001-11-14 17:38:23 $ 26 26 // 27 27 //////////////////////////////////////////////////////////////////////// … … 69 69 #include "MMcEvt.hxx" 70 70 #include "MMcTrig.hxx" 71 #include "MMcRunHeader.hxx" 71 72 #include "MMcTrigHeader.hxx" 72 73 #include "MMcFadcHeader.hxx" … … 431 432 432 433 int simulateNSB; //@< Will we simulate NSB? 433 int nphe2NSB; //@< From how many ohe will we simulate NSB?434 int nphe2NSB; //@< From how many phe will we simulate NSB? 434 435 float meanNSB; //@< diffuse NSB mean value (phe per ns per central pixel) 435 436 float diffnsb_phepns[iMAXNUMPIX]; //@< diffuse NSB values for each pixel … … 455 456 EXTWAVEBAND5 456 457 }; 458 float zenfactor=1.0; // correction factor calculated from the extinction 457 459 458 460 float qTailCut; //@< Tail Cut value … … 473 475 float fpixelthres[TRIGGER_PIXELS]; //@< Threshold values 474 476 475 TArrayC *fadcValues; //@< the analog Fadc sig anl for pixels477 TArrayC *fadcValues; //@< the analog Fadc signal for pixels 476 478 477 479 float plateScale_cm2deg; //@< plate scale (deg/cm) … … 479 481 480 482 int still_in_loop = FALSE; 483 484 //@< Variables to fill the McRunHeader 485 Int_t sfRaH = 5; 486 Int_t sfRaM = 34; 487 Int_t sfRaS = 32; 488 Int_t sfDeD = 22; 489 Int_t sfDeM = 00; 490 Int_t sfDeS = 52; 491 Float_t shthetamax = 0.0; 492 Float_t shthetamin = 0.0; 493 Float_t shphimax = 0.0; 494 Float_t shphimin = 0.0; 495 Float_t telestheta = 0.0; 496 Float_t telesphi = 0.0; 497 Float_t sofftheta = 0.0; 498 Float_t soffphi = 0.0; 499 UInt_t corsika = 5200 ; 500 481 501 482 502 struct camera cam; // structure holding the camera definition … … 743 763 // Initialise McTrig information class if we want to save trigger informtion 744 764 745 MMcTrig **McTrig ;746 MMcTrigHeader **HeaderTrig ;747 MMcFadcHeader **HeaderFadc ;765 MMcTrig **McTrig = NULL; 766 MMcTrigHeader **HeaderTrig = NULL; 767 MMcFadcHeader **HeaderFadc = NULL; 748 768 749 769 if (Write_McTrig){ … … 796 816 // Header Tree 797 817 798 MRawRunHeader **RunHeader; 799 800 RunHeader = new MRawRunHeader * [1]; 801 RunHeader[0] = new MRawRunHeader(); 818 MRawRunHeader *RunHeader= new MRawRunHeader(); 819 MMcRunHeader *McRunHeader = new MMcRunHeader(); 802 820 803 821 // Header branch 804 822 805 MRawEvtHeader **EvtHeader ;823 MRawEvtHeader **EvtHeader = NULL; 806 824 807 825 if (Write_RawEvt) { … … 815 833 // Data branch 816 834 817 MRawEvtData **EvtData ; // Data branch835 MRawEvtData **EvtData = NULL; // Data branch 818 836 819 837 if (Write_RawEvt) { … … 822 840 for (i=0;i<icontrigger;i++) { 823 841 EvtData[i] = new MRawEvtData(); 842 EvtData[i]->Init(RunHeader); // We need thr RunHeader to read 843 // number of pixels 824 844 } 825 845 } … … 833 853 TFile outfile_temp ( rootname , "RECREATE" ); 834 854 835 Int_t bsize=128000; Int_t split=1;836 837 855 // create a Tree for the Header Event 838 856 TTree HeaderTree("RunHeaders","Headers of Run"); … … 843 861 844 862 HeaderTree.Branch("MRawRunHeader","MRawRunHeader", 845 &RunHeader[0], bsize, split); 863 &RunHeader); 864 865 HeaderTree.Branch("MMcRunHeader","MMcRunHeader", 866 &McRunHeader); 846 867 847 868 if(!Trigger_Loop && Write_McTrig){ 848 869 849 870 HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader", 850 &HeaderTrig[0], bsize, split);871 &HeaderTrig[0]); 851 872 } 852 873 if (Trigger_Loop && Write_McTrig){ … … 858 879 strcat (branchname, "."); 859 880 HeaderTree.Branch(branchname,"MMcTrigHeader", 860 &HeaderTrig[i], bsize, split);881 &HeaderTrig[i]); 861 882 } 862 883 } … … 865 886 866 887 HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader", 867 &HeaderFadc[0], bsize, split);888 &HeaderFadc[0]); 868 889 } 869 890 if (Trigger_Loop && Write_McFADC){ … … 875 896 strcat (branchname, "."); 876 897 HeaderTree.Branch(branchname,"MMcFadcHeader", 877 &HeaderFadc[i], bsize, split);898 &HeaderFadc[i]); 878 899 } 879 900 } … … 881 902 // Fill branches for MRawRunHeader 882 903 883 RunHeader[0]->SetMagicNumber(kMagicNumber); 884 RunHeader[0]->SetFormatVersion(2); 885 RunHeader[0]->SetSoftVersion((UShort_t) (VERSION*10)); 886 RunHeader[0]->SetRunType(256); 887 RunHeader[0]->SetRunNumber(0); 888 RunHeader[0]->SetNumSamples(0,FADC_SLICES); 904 RunHeader->SetMagicNumber(kMagicNumber); 905 RunHeader->SetFormatVersion(2); 906 RunHeader->SetSoftVersion((UShort_t) (VERSION*10)); 907 RunHeader->SetRunType(256); 908 RunHeader->SetRunNumber(0); 909 RunHeader->SetNumSamples(0,FADC_SLICES); 910 RunHeader->SetNumCrates(1); 911 RunHeader->SetNumPixInCrate(CAMERA_PIXELS); 889 912 890 913 … … 956 979 957 980 // Fill the Header Tree with the current leaves of each branch 958 HeaderTree.Fill() ;981 // HeaderTree.Fill() ; 959 982 960 983 … … 965 988 966 989 EvtTree.Branch("MMcEvt","MMcEvt", 967 &McEvt , bsize, split);990 &McEvt); 968 991 } 969 992 … … 972 995 if (Write_RawEvt){ 973 996 EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", 974 &EvtHeader[0], bsize, split);997 &EvtHeader[0]); 975 998 EvtTree.Branch("MRawEvtData","MRawEvtData", 976 &EvtData[0] , bsize, split);999 &EvtData[0]); 977 1000 } 978 1001 if (Write_McTrig){ 979 1002 EvtTree.Branch("MMcTrig","MMcTrig", 980 &McTrig[0] , bsize, split);1003 &McTrig[0]); 981 1004 } 982 1005 } … … 990 1013 strcat (branchname, "."); 991 1014 EvtTree.Branch(branchname,"MMcTrig", 992 &McTrig[i] , bsize, split);1015 &McTrig[i]); 993 1016 } 994 1017 } … … 1003 1026 strcat (branchname, "."); 1004 1027 EvtTree.Branch(branchname,"MRawEvtHeader", 1005 &EvtHeader[i] , bsize, split);1028 &EvtHeader[i]); 1006 1029 } 1007 1030 for(char branchname[15],i=0;i<icontrigger;i++){ … … 1012 1035 strcat (branchname, "."); 1013 1036 EvtTree.Branch(branchname,"MRawEvtData", 1014 &EvtData[i], bsize, split); 1015 } 1016 } 1017 1018 unsigned short ulli = 0 ; 1037 &EvtData[i]); 1038 } 1039 } 1019 1040 1020 1041 TApplication theAppTrigger("App", &argc, argv); 1021 1042 1022 if (gROOT->IsBatch()) {1023 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);1024 // return 1;1025 }1026 1027 1043 if(FADC_Scan){ 1028 TApplication theAppFadc("App", &argc, argv); 1029 1030 if (gROOT->IsBatch()) { 1031 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); 1032 // return 1; 1033 } 1044 if (gROOT->IsBatch()) { 1045 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); 1046 // return 1; 1047 } 1048 } 1049 if(FADC_Scan){ 1050 //TApplication theAppFadc("App", &argc, argv); 1051 1052 if (gROOT->IsBatch()) { 1053 fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]); 1054 // return 1; 1055 } 1034 1056 } 1035 1057 … … 1060 1082 k = produce_nsbrates( starfieldname, 1061 1083 &cam, 1062 nsbrate_phepns 1084 nsbrate_phepns); 1063 1085 1064 1086 if (k != 0){ … … 1075 1097 1076 1098 // calculate nsb rate including diffuse and starlight 1077 for(i=0; i<cam.inumpixels;i++){ 1078 nsb_phepns[i]= diffnsb_phepns[i]; 1079 for(j=0;j<iNUMWAVEBANDS;j++) 1080 nsb_phepns[i]=nsb_phepns[i]+ nsbrate_phepns[i][j]; 1099 // we also include the extinction effect 1100 for(j=0;j<iNUMWAVEBANDS;j++){ 1101 // calculate the effect of the atmospheric extinction 1102 1103 zenfactor = pow(10., -0.4 * ext[i] ); 1104 1105 for(i=0; i<cam.inumpixels;i++){ 1106 nsb_phepns[i]+=diffnsb_phepns[i]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[i][j]; 1107 } 1081 1108 } 1082 1109 } … … 1357 1384 fadc.ElecNoise() ; 1358 1385 1386 // We should simulate the AC coupling behaviour: 1387 // For the FADC it is only done for the NSB while producinr 1388 // the StarResponse database. 1389 // For the trigger is done in the Trigger.Diskriminate(), which 1390 // is called later (it should be seperated to speed up the program 1391 1359 1392 // now a shift in the fadc signal due to the pedestlas is 1360 1393 // intorduced … … 1406 1439 // Start the First Level Trigger simulation 1407 1440 // 1408 Lev1=Trigger.FirstLevel(); 1441 if(Lev0!=0) 1442 Lev1=Trigger.FirstLevel(); 1409 1443 if (Write_McTrig) 1410 1444 McTrig[iconcount]->SetFirstLevel (Lev1); … … 1700 1734 } // end if else found start of run 1701 1735 } // end big while loop 1736 1737 //<@ Finally we should fill th McRunHeader 1738 1739 1740 get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS); 1741 get_teles_axis(&telestheta,&telesphi); 1742 get_source_off(&sofftheta,&soffphi); 1743 mcevth.get_theta_range(&shthetamin, &shthetamax); 1744 mcevth.get_phi_range(&shphimin,&shphimax); 1745 corsika=get_corsika_ver(); 1746 if(!Trigger_Loop) icontrigger=0; 1747 McRunHeader->Fill(icontrigger, 1748 !Write_All_Images, 1749 Write_McEvt, 1750 Write_McTrig, 1751 Write_McFADC, 1752 CAMERA_PIXELS, 1753 (UInt_t)ntshow, 1754 (UInt_t)ntrigger, 1755 sfRaH, 1756 sfRaM, 1757 sfRaS, 1758 sfDeD, 1759 sfDeM, 1760 sfDeS, 1761 meanNSB, 1762 telestheta, 1763 telesphi, 1764 sofftheta, 1765 soffphi, 1766 shthetamax, 1767 shthetamin, 1768 shphimax, 1769 shphimin, 1770 corsika, 1771 (UInt_t)(REFL_VERSION*100), 1772 (UInt_t)(VERSION*100)); 1773 1774 // Fill the Header Tree with the current leaves of each branch 1775 HeaderTree.Fill() ; 1702 1776 1703 1777 //++ … … 2621 2695 2622 2696 /******** bpoint_is_in_pix() ***************************************/ 2623 2697 /* 2624 2698 int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){ 2625 / * return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */2626 / * the pixel is assumed to be a "closed set" */2627 2628 double a, b; / * a = length of one of the edges of one pixel, b = half the width of one pixel */2629 double c, xx, yy; / * auxiliary variable */2699 // return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) 2700 // the pixel is assumed to be a "closed set" 2701 2702 double a, b; // a = length of one of the edges of one pixel, b = half the width of one pixel 2703 double c, xx, yy; // auxiliary variable 2630 2704 2631 2705 b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum]; … … 2643 2717 if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) || 2644 2718 ((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){ 2645 return(TRUE); / * point is inside */2719 return(TRUE); // point is inside 2646 2720 } 2647 2721 else{ 2648 return(FALSE); /* point is outside */ 2649 } 2722 return(FALSE); // point is outside 2723 } 2724 } 2725 */ 2726 2727 #define sqrt13 0.577350269 // = 1./sqrt(3.) 2728 #define sqrt3 1.732050807 // = sqrt(3.) 2729 2730 int bpoint_is_in_pix(double dx, double dy, struct camera *pcam) 2731 { 2732 /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */ 2733 /* the pixel is assumed to be a "closed set" */ 2734 2735 /* 2736 a = length of one of the edges of one pixel, 2737 b = half the width of one pixel 2738 */ 2739 2740 const int num = pcam->inumpixels; 2741 const double diam = pcam->dpixdiameter_cm; 2742 2743 const double diam2 = diam/2; 2744 const double diam3 = diam/sqrt3; 2745 2746 for (int i=0; i<num; i++) 2747 { 2748 const double size = pcam->dpixsizefactor[i]; 2749 2750 const double b = diam2 * size; 2751 const double a = diam3 * size; 2752 2753 const double xx = dx - pcam->dxc[i]; 2754 const double yy = dy - pcam->dyc[i]; 2755 2756 if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) || 2757 ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){ 2758 2759 return i; // inside i 2760 } 2761 2762 // return -1; // outside 2763 } 2764 return -1; // outside 2650 2765 } 2651 2766 … … 2789 2904 // read the photons data 2790 2905 2791 fread ( flag, SIZE_OF_FLAGS, 1, sp );2792 2906 2793 2907 // loop over the photons 2794 2908 2795 while ( !isA( flag, FLAG_END_OF_EVENT ) ) { 2909 while (1) { 2910 2911 fread ( flag, SIZE_OF_FLAGS, 1, sp ); 2912 2913 if (isA( flag, FLAG_END_OF_EVENT )) 2914 break; 2796 2915 2797 2916 memcpy( (char*)&photon, flag, SIZE_OF_FLAGS ); … … 2808 2927 2809 2928 if (t-*tmin_ns>TOTAL_TRIGGER_TIME) { 2810 2811 // read next Photon 2812 2813 fread( flag, SIZE_OF_FLAGS, 1, sp ); 2814 2815 // go to beginning of loop, the photon is lost 2816 continue; 2929 continue; 2817 2930 } 2818 2931 … … 2833 2946 2834 2947 if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){ 2835 2836 // cout << " lost first check\n";2837 2838 // read next CPhoton2839 2840 fread ( flag, SIZE_OF_FLAGS, 1, sp );2841 2842 // go to beginning of loop, the photon is lost2843 2948 continue; 2844 2949 2845 2950 } 2846 2951 2847 ipixnum = -1; 2848 2849 for(i=0; i<cam->inumpixels; i++){ 2850 if( bpoint_is_in_pix( cx, cy, i, cam) ){ 2851 ipixnum = i; 2852 break; 2853 } 2854 } 2855 2856 if(ipixnum==-1){// the photon is in none of the pixels 2857 2858 // cout << " lost pixlization\n"; 2859 2860 // read next CPhoton 2861 2862 fread ( flag, SIZE_OF_FLAGS, 1, sp ); 2863 2864 // go to beginning of loop, the photon is lost 2865 continue; 2866 } 2867 2868 if(ipixnum==0) {// the phton is in the central pixel, which is not used for trigger 2869 // read next CPhoton 2870 2871 fread ( flag, SIZE_OF_FLAGS, 1, sp ); 2872 2873 // go to beginning of loop, the photon is lost 2874 continue; 2875 } 2952 ipixnum = bpoint_is_in_pix(cx, cy, cam); 2953 2954 // -1 = the photon is in none of the pixels 2955 // 0 = the phton is in the central pixel, which is not used for trigger 2956 if (ipixnum==-1 || ipixnum==0) { 2957 continue; 2958 } 2959 2876 2960 //+++ 2877 2961 // QE simulation … … 2885 2969 2886 2970 if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){ 2887 2888 // cout << " lost\n";2889 2890 // read next Photon2891 2892 fread ( flag, SIZE_OF_FLAGS, 1, sp );2893 2894 // go to beginning of loop2895 2971 continue; 2896 2972 … … 2911 2987 2912 2988 if ( (RandomNumber) > qe ) { 2913 2914 // cout << " lost\n";2915 2916 // read next Photon2917 2918 fread ( flag, SIZE_OF_FLAGS, 1, sp );2919 2920 // go to beginning of loop2921 2989 continue; 2922 2923 2990 } 2924 2991 … … 2938 3005 2939 3006 *itotnphe += 1; 2940 2941 // read next Photon 2942 2943 fread( flag, SIZE_OF_FLAGS, 1, sp ); 2944 2945 } // end while, i.e. found end of event 3007 } 2946 3008 2947 3009 return(0); … … 3248 3310 // 3249 3311 // $Log: not supported by cvs2svn $ 3312 // Revision 1.28 2001/10/26 16:31:45 blanch 3313 // Removing several warnings. 3314 // 3250 3315 // Revision 1.27 2001/09/05 10:04:33 blanch 3251 3316 // *** empty log message *** … … 3317 3382 // 3318 3383 // $Log: not supported by cvs2svn $ 3384 // Revision 1.28 2001/10/26 16:31:45 blanch 3385 // Removing several warnings. 3386 // 3319 3387 // Revision 1.27 2001/09/05 10:04:33 blanch 3320 3388 // *** empty log message ***
Note:
See TracChangeset
for help on using the changeset viewer.