Ignore:
Timestamp:
11/14/01 17:38:23 (23 years ago)
Author:
blanch
Message:
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
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx

    r1001 r1078  
    2121//
    2222// $RCSfile: camera.cxx,v $
    23 // $Revision: 1.28 $
     23// $Revision: 1.29 $
    2424// $Author: blanch $
    25 // $Date: 2001-10-26 16:31:45 $
     25// $Date: 2001-11-14 17:38:23 $
    2626//
    2727////////////////////////////////////////////////////////////////////////
    … …  
    6969#include "MMcEvt.hxx"
    7070#include "MMcTrig.hxx"
     71#include "MMcRunHeader.hxx"
    7172#include "MMcTrigHeader.hxx"
    7273#include "MMcFadcHeader.hxx"
    … …  
    431432
    432433  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?
    434435  float meanNSB;              //@< diffuse NSB mean value (phe per ns per central pixel)
    435436  float diffnsb_phepns[iMAXNUMPIX];  //@< diffuse NSB values for each pixel 
    … …  
    455456    EXTWAVEBAND5
    456457  };                         
     458  float zenfactor=1.0;  //  correction factor calculated from the extinction
    457459
    458460  float qTailCut;             //@< Tail Cut value
    … …  
    473475  float fpixelthres[TRIGGER_PIXELS];         //@< Threshold values
    474476
    475   TArrayC *fadcValues;  //@< the analog Fadc siganl for pixels
     477  TArrayC *fadcValues;  //@< the analog Fadc signal for pixels
    476478
    477479  float plateScale_cm2deg;    //@< plate scale (deg/cm)
    … …  
    479481
    480482  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
    481501
    482502  struct camera cam; // structure holding the camera definition
    … …  
    743763  //  Initialise McTrig information class if we want to save trigger informtion
    744764
    745   MMcTrig **McTrig;
    746   MMcTrigHeader **HeaderTrig;
    747   MMcFadcHeader **HeaderFadc;
     765  MMcTrig **McTrig = NULL;
     766  MMcTrigHeader **HeaderTrig = NULL;
     767  MMcFadcHeader **HeaderFadc = NULL;
    748768
    749769  if (Write_McTrig){
    … …  
    796816       // Header Tree
    797817
    798   MRawRunHeader **RunHeader;
    799 
    800   RunHeader = new MRawRunHeader * [1];
    801   RunHeader[0] = new MRawRunHeader();
     818  MRawRunHeader *RunHeader= new MRawRunHeader();
     819  MMcRunHeader  *McRunHeader = new MMcRunHeader();
    802820
    803821       // Header branch
    804822
    805   MRawEvtHeader **EvtHeader;
     823  MRawEvtHeader **EvtHeader = NULL;
    806824
    807825  if (Write_RawEvt) {
    … …  
    815833       // Data branch
    816834
    817   MRawEvtData **EvtData;    // Data branch
     835  MRawEvtData **EvtData = NULL;    // Data branch
    818836
    819837  if (Write_RawEvt) {
    … …  
    822840    for (i=0;i<icontrigger;i++) {
    823841      EvtData[i] = new MRawEvtData();
     842      EvtData[i]->Init(RunHeader);     //  We need thr RunHeader to read
     843                                       //  number of pixels
    824844    }
    825845  }
    … …  
    833853  TFile outfile_temp ( rootname , "RECREATE" );
    834854
    835   Int_t bsize=128000; Int_t split=1;
    836 
    837855  //      create a Tree for the Header Event
    838856  TTree HeaderTree("RunHeaders","Headers of Run");
    … …  
    843861
    844862  HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
    845                     &RunHeader[0], bsize, split);
     863                    &RunHeader);
     864
     865  HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
     866                    &McRunHeader);
    846867
    847868  if(!Trigger_Loop && Write_McTrig){
    848869   
    849870    HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader",
    850                  &HeaderTrig[0], bsize, split);   
     871                      &HeaderTrig[0]);
    851872  }
    852873  if (Trigger_Loop && Write_McTrig){
    … …  
    858879      strcat (branchname, ".");
    859880      HeaderTree.Branch(branchname,"MMcTrigHeader",
    860                      &HeaderTrig[i], bsize, split);
     881                        &HeaderTrig[i]);
    861882    }
    862883  } 
    … …  
    865886   
    866887    HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader",
    867                  &HeaderFadc[0], bsize, split);   
     888                      &HeaderFadc[0]);
    868889  }
    869890  if (Trigger_Loop && Write_McFADC){
    … …  
    875896      strcat (branchname, ".");
    876897      HeaderTree.Branch(branchname,"MMcFadcHeader",
    877                      &HeaderFadc[i], bsize, split);
     898                        &HeaderFadc[i]);
    878899    }
    879900  } 
    … …  
    881902  //  Fill branches for MRawRunHeader
    882903
    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);
    889912
    890913
    … …  
    956979
    957980  //  Fill the Header Tree with the current leaves of each branch
    958   HeaderTree.Fill() ;
     981  // HeaderTree.Fill() ;
    959982           
    960983
    … …  
    965988
    966989    EvtTree.Branch("MMcEvt","MMcEvt",
    967                    &McEvt, bsize, split); 
     990                   &McEvt);
    968991  }
    969992
    … …  
    972995    if (Write_RawEvt){
    973996      EvtTree.Branch("MRawEvtHeader","MRawEvtHeader",
    974              &EvtHeader[0], bsize, split);
     997                     &EvtHeader[0]);
    975998      EvtTree.Branch("MRawEvtData","MRawEvtData",
    976                      &EvtData[0], bsize, split);
     999                     &EvtData[0]);
    9771000    }
    9781001    if (Write_McTrig){
    9791002      EvtTree.Branch("MMcTrig","MMcTrig",
    980                      &McTrig[0], bsize, split);
     1003                     &McTrig[0]);
    9811004    }   
    9821005  }
    … …  
    9901013        strcat (branchname, ".");
    9911014        EvtTree.Branch(branchname,"MMcTrig",
    992                        &McTrig[i], bsize, split);
     1015                       &McTrig[i]);
    9931016      }
    9941017    }
    … …  
    10031026      strcat (branchname, ".");
    10041027      EvtTree.Branch(branchname,"MRawEvtHeader",
    1005                      &EvtHeader[i], bsize, split);
     1028                     &EvtHeader[i]);
    10061029    }
    10071030    for(char branchname[15],i=0;i<icontrigger;i++){
    … …  
    10121035      strcat (branchname, ".");
    10131036      EvtTree.Branch(branchname,"MRawEvtData",
    1014                      &EvtData[i], bsize, split);
    1015     }
    1016   }
    1017 
    1018   unsigned short ulli = 0 ;
     1037                     &EvtData[i]);
     1038    }
     1039  }
    10191040
    10201041  TApplication theAppTrigger("App", &argc, argv);
    10211042
    1022   if (gROOT->IsBatch()) {
    1023     fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
    1024     //    return 1;
    1025   }
    1026 
    10271043  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      }
    10341056  }
    10351057 
    … …  
    10601082    k = produce_nsbrates( starfieldname,
    10611083                          &cam,
    1062                           nsbrate_phepns );
     1084                          nsbrate_phepns);
    10631085
    10641086    if (k != 0){
    … …  
    10751097
    10761098    // 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        }
    10811108    }
    10821109  }
    … …  
    13571384        fadc.ElecNoise() ;
    13581385
     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
    13591392        //  now a shift in the fadc signal due to the pedestlas is
    13601393        //  intorduced
    … …  
    14061439                  //   Start the First Level Trigger simulation
    14071440                  //
    1408                   Lev1=Trigger.FirstLevel();
     1441                  if(Lev0!=0)
     1442                      Lev1=Trigger.FirstLevel();
    14091443                  if (Write_McTrig)
    14101444                    McTrig[iconcount]->SetFirstLevel (Lev1);
    … …  
    17001734    } // end if else found start of run
    17011735  } // 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() ;
    17021776 
    17031777  //++
    … …  
    26212695
    26222696/******** bpoint_is_in_pix() ***************************************/
    2623 
     2697/*
    26242698int 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
    26302704
    26312705  b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
    … …  
    26432717  if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) ||
    26442718     ((0. <  xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a)))   ){
    2645     return(TRUE); /* point is inside */
     2719    return(TRUE); // point is inside
    26462720  }
    26472721  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
     2730int 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
    26502765}
    26512766
    … …  
    27892904  // read the photons data
    27902905 
    2791   fread ( flag, SIZE_OF_FLAGS, 1, sp );
    27922906 
    27932907  // loop over the photons
    27942908 
    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;
    27962915         
    27972916    memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
    … …  
    28082927
    28092928    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;
    28172930    }
    28182931         
    … …  
    28332946         
    28342947    if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){
    2835            
    2836       // cout << " lost first check\n";
    2837 
    2838       // read next CPhoton
    2839 
    2840       fread ( flag, SIZE_OF_FLAGS, 1, sp );
    2841            
    2842       // go to beginning of loop, the photon is lost
    28432948      continue;
    28442949           
    28452950    }
    28462951
    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
    28762960    //+++
    28772961    // QE simulation
    … …  
    28852969   
    28862970    if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
    2887 
    2888       // cout << " lost\n";
    2889      
    2890       // read next Photon
    2891 
    2892       fread ( flag, SIZE_OF_FLAGS, 1, sp );
    2893            
    2894       // go to beginning of loop
    28952971      continue;
    28962972           
    … …  
    29112987
    29122988    if ( (RandomNumber) > qe ) {
    2913 
    2914       // cout << " lost\n";
    2915 
    2916       // read next Photon
    2917 
    2918       fread ( flag, SIZE_OF_FLAGS, 1, sp );
    2919            
    2920       // go to beginning of loop
    29212989      continue;
    2922            
    29232990    }
    29242991         
    … …  
    29383005   
    29393006    *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  }
    29463008         
    29473009  return(0);
    … …  
    32483310//
    32493311// $Log: not supported by cvs2svn $
     3312// Revision 1.28  2001/10/26 16:31:45  blanch
     3313// Removing several warnings.
     3314//
    32503315// Revision 1.27  2001/09/05 10:04:33  blanch
    32513316// *** empty log message ***
    … …  
    33173382//
    33183383// $Log: not supported by cvs2svn $
     3384// Revision 1.28  2001/10/26 16:31:45  blanch
     3385// Removing several warnings.
     3386//
    33193387// Revision 1.27  2001/09/05 10:04:33  blanch
    33203388// *** empty log message ***
Note: See TracChangeset for help on using the changeset viewer.