Ignore:
Timestamp:
01/14/03 13:40:17 (22 years ago)
Author:
blanch
Message:
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.
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1694 r1706  
    2121//
    2222// $RCSfile: camera.cxx,v $
    23 // $Revision: 1.50 $
     23// $Revision: 1.51 $
    2424// $Author: blanch $
    25 // $Date: 2003-01-07 16:33:31 $
     25// $Date: 2003-01-14 13:40:17 $
    2626//
    2727////////////////////////////////////////////////////////////////////////
    … …  
    443443  int addElecNoise;           //@< Will we add ElecNoise?
    444444
     445  float riseDiskThres;
     446  float secureDiskThres;
     447
    445448  int simulateNSB;            //@< Will we simulate NSB?
    446449  int nphe2NSB;               //@< From how many phe will we simulate NSB?
    … …  
    472475
    473476  int anaPixels;
    474    
     477
     478  int nstoredevents = 0;
     479  int flagstoring = 0;
     480
    475481  float fCorrection;          //@< Factor to apply to pixel values (def. 1.)
    476482
    … …  
    609615  Individual_Thres_Pixel = get_indi_thres_pixel();
    610616
     617  get_secure_threhold(&riseDiskThres,&secureDiskThres);
     618
    611619  get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm);
    612620
    … …  
    811819  Trigger.CheckThreshold(&qThreshold[0]);
    812820
     821  // Set flag in pixel 0 (not used for trigger) that indicates if secure pixel
     822  // is active: secureDiskThres*10000+riseDiskThres
     823
     824  if(riseDiskThres>0.0)
     825    qThreshold[0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres;
     826  else
     827    qThreshold[0]=0.0;
    813828
    814829  //  Initialise McTrig information class if we want to save trigger informtion
    … …  
    963978
    964979  RunHeader->SetMagicNumber(kMagicNumber);
    965   RunHeader->SetFormatVersion(3);
     980  RunHeader->SetFormatVersion(4);
    966981  RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
    967982  RunHeader->SetRunType(256);
    … …  
    12881303        if (reflector_file_version<6){
    12891304          fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
    1290           thetashw = mcevth.get_theta();
    1291           phishw = mcevth.get_phi();
    12921305        }
    12931306        else{
    12941307          fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile );
    1295           thetashw = mcevth_2.get_theta();
    1296           phishw = mcevth_2.get_phi();
    12971308        }
    12981309
    … …  
    13411352        m1 = sin(thetashw)*sin(phishw);
    13421353        n1 = cos(thetashw);
    1343        
    1344         // read the deviation of the telescope with respect to the shower
    1345        
    1346         if (reflector_file_version<6)
     1354
     1355        if (reflector_file_version<6){
     1356
     1357          mcevth.get_core(&coreX, &coreY);     
     1358         
     1359          // read the deviation of the telescope with respect to the shower
    13471360          mcevth.get_deviations ( &thetaCT, &phiCT );
    1348         else
    1349           mcevth_2.get_deviations ( &thetaCT, &phiCT );
    1350 
    1351         if ( (thetaCT == 0.) && (phiCT == 0.) ) {
     1361
     1362          if ( (thetaCT == 0.) && (phiCT == 0.) ) {
    13521363         
    1353           // CT was looking to the source (both lines are parallel)
    1354           // therefore, we calculate the impact parameter as the distance
    1355           // between the CT axis and the core position
    1356 
    1357           if (reflector_file_version<6)
    1358             mcevth.get_core(&coreX, &coreY);
    1359           else
    1360             mcevth_2.get_core(&coreX, &coreY);
    1361 
    1362           impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
     1364            // CT was looking to the source (both lines are parallel)
     1365            // therefore, we calculate the impact parameter as the distance
     1366            // between the CT axis and the core position
     1367
     1368            impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
     1369            thetaCT += thetashw;
     1370            phiCT   += phishw;
     1371           
     1372          } else {
     1373           
     1374            // the shower comes off-axis
     1375           
     1376            // obtain with this the final direction of the CT
     1377           
     1378            thetaCT += thetashw;
     1379            phiCT   += phishw;
     1380           
     1381            // calculate vector for telescope
     1382           
     1383            l2 = sin(thetaCT)*cos(phiCT);
     1384            m2 = sin(thetaCT)*sin(phiCT);
     1385            n2 = cos(thetaCT);
     1386           
     1387            num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
     1388            den = (SQR(l1*m2 - l2*m1) +
     1389                   SQR(m1*n2 - m2*n1) +
     1390                   SQR(n1*l2 - n2*l1));
     1391            den = sqrt(den);
     1392           
     1393            impactD = fabs(num)/den;
     1394           
     1395          }
     1396        }
     1397        else{
     1398          mcevth_2.get_core(&coreX, &coreY);   
     1399          thetaCT=mcevth_2.get_theta_CT();
     1400          phiCT=mcevth_2.get_phi_CT();
     1401          if ( (thetaCT == thetashw) && (phiCT == phishw) ) {
    13631402         
    1364         } else {
    1365          
    1366           // the shower comes off-axis
    1367          
    1368           // obtain with this the final direction of the CT
    1369          
    1370           thetaCT += thetashw;
    1371           phiCT   += phishw;
    1372          
    1373           // calculate vector for telescope
    1374          
    1375           l2 = sin(thetaCT)*cos(phiCT);
    1376           m2 = sin(thetaCT)*sin(phiCT);
    1377           n2 = cos(thetaCT);
    1378          
    1379           num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
    1380           den = (SQR(l1*m2 - l2*m1) +
    1381                  SQR(m1*n2 - m2*n1) +
    1382                  SQR(n1*l2 - n2*l1));
    1383           den = sqrt(den);
    1384          
    1385           impactD = fabs(num)/den;
    1386          
     1403            // CT was looking to the source (both lines are parallel)
     1404            // therefore, we calculate the impact parameter as the distance
     1405            // between the CT axis and the core position
     1406
     1407            impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
     1408           
     1409          } else {
     1410           
     1411            // the shower comes off-axis
     1412           
     1413            // calculate vector for telescope
     1414           
     1415            l2 = sin(thetaCT)*cos(phiCT);
     1416            m2 = sin(thetaCT)*sin(phiCT);
     1417            n2 = cos(thetaCT);
     1418           
     1419            num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
     1420            den = (SQR(l1*m2 - l2*m1) +
     1421                   SQR(m1*n2 - m2*n1) +
     1422                   SQR(n1*l2 - n2*l1));
     1423            den = sqrt(den);
     1424           
     1425            impactD = fabs(num)/den;
     1426           
     1427          }
    13871428        }
    13881429
    … …  
    14651506                             rho);
    14661507           
    1467  
    14681508          }
    14691509
    … …  
    15511591          // Set to zero the flag to know if some conditon has triggered
    15521592          btrigger=0;
     1593          flagstoring = 0;
    15531594          //  Loop over trigger threshold
    15541595          int iconcount;
    15551596          for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
    1556             for (i=0;i<CAMERA_PIXELS;i++)
     1597            for (i=0;i<CAMERA_PIXELS;i++){
    15571598              fpixelthres[i]=
    15581599                ((Float_t)(fthrescount)>=qThreshold[i])?
    15591600                (Float_t)(fthrescount):qThreshold[i];
     1601
     1602              // Rise the discrimnator threshold to avoid huge rates
     1603
     1604              if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
     1605                for(int ii=0;ii<CAMERA_PIXELS;ii++)
     1606                  if( nsb_phepns_rotated[ii]>riseDiskThres)
     1607                    fpixelthres[ii]=secureDiskThres;
     1608            }
    15601609            Trigger.SetThreshold(fpixelthres);
    15611610
    … …  
    16651714            //
    16661715
     1716            if(!flagstoring)
     1717              nstoredevents++;
     1718            flagstoring = 1;
     1719
    16671720            if (Write_McEvt) {
    16681721              Float_t ftime, ltime;
    … …  
    16811734                             mcevth.get_coreY(),
    16821735                             impactD,
    1683                              mcevth_2.get_theta_CT(),
    1684                              mcevth_2.get_phi_CT(),
     1736                             phiCT,
     1737                             thetaCT,
    16851738                             ftime,
    16861739                             ltime,
    … …  
    17011754              else{
    17021755                Float_t Nmax, t0, tmax, a, b, c, chi2;
    1703                 mcevth.get_times(&ftime, &ltime);
     1756                mcevth_2.get_times(&ftime, &ltime);
    17041757                chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
    17051758                McEvt->Fill( mcevth_2.get_evt_number(),
    … …  
    17151768                             mcevth_2.get_coreY(),
    17161769                             impactD,
     1770                             mcevth_2.get_phi_CT(),
    17171771                             mcevth_2.get_theta_CT(),
    1718                              mcevth_2.get_phi_CT(),
    17191772                             ftime,
    17201773                             ltime,
    … …  
    17611814          for (int i=0;i<CAMERA_PIXELS;i++)
    17621815            fpixelthres[i]=qThreshold[i];
     1816
     1817          // Rise the discrimnator threshold to avoid huge rates
     1818          if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
     1819            for(int ii=0;ii<CAMERA_PIXELS;ii++){
     1820              if( nsb_phepns_rotated[ii]>riseDiskThres)
     1821                fpixelthres[ii]=secureDiskThres;
     1822            }
     1823         
    17631824          Trigger.SetThreshold(fpixelthres);
    17641825                                 
    … …  
    18021863            //
    18031864            fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
     1865
     1866            nstoredevents++;
    18041867
    18051868            if (Write_McTrig){
    … …  
    18511914                             mcevth.get_coreY(),
    18521915                             impactD,
    1853                              mcevth_2.get_theta_CT(),
    1854                              mcevth_2.get_phi_CT(),
     1916                             phiCT,
     1917                             thetaCT,
    18551918                             ftime,
    18561919                             ltime,
    … …  
    18711934              else{
    18721935                Float_t Nmax, t0, tmax, a, b, c, chi2;
    1873                 mcevth.get_times(&ftime, &ltime);
     1936                mcevth_2.get_times(&ftime, &ltime);
    18741937                chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
    18751938                McEvt->Fill( mcevth_2.get_evt_number(),
    … …  
    18851948                             mcevth_2.get_coreY(),
    18861949                             impactD,
     1950                             mcevth_2.get_phi_CT(),
    18871951                             mcevth_2.get_theta_CT(),
    1888                              mcevth_2.get_phi_CT(),
    18891952                             ftime,
    18901953                             ltime,
    … …  
    20222085
    20232086  get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
    2024   get_teles_axis(&telestheta,&telesphi);
     2087  if (reflector_file_version<6){
     2088    get_teles_axis(&telestheta,&telesphi);
     2089    mcevth.get_theta_range(&shthetamin, &shthetamax);
     2090    mcevth.get_phi_range(&shphimin,&shphimax);
     2091    mcevth.get_theta_range(&shthetamin, &shthetamax);
     2092    mcevth.get_phi_range(&shphimin,&shphimax);
     2093    corsika=UInt_t(mcevth.get_VersionPGM()*1000);
     2094    for (int i=0; i< 10;i++)
     2095      heights[i]=mcevth.get_HeightLev (i);
     2096    rnum=mcevth.get_RunNumber();
     2097  }
     2098  else{
     2099    telestheta=mcevth_2.get_theta_CT();
     2100    telesphi=mcevth_2.get_phi_CT();
     2101    mcevth_2.get_theta_range(&shthetamin, &shthetamax);
     2102    mcevth_2.get_phi_range(&shphimin,&shphimax);
     2103    corsika=UInt_t(mcevth_2.get_VersionPGM()*1000);
     2104    for (int i=0; i< 10;i++)
     2105      heights[i]=mcevth_2.get_HeightLev (i);
     2106    rnum=mcevth_2.get_RunNumber();
     2107  }
     2108
    20252109  get_source_off(&sofftheta,&soffphi);
    2026   mcevth.get_theta_range(&shthetamin, &shthetamax);
    2027   mcevth.get_phi_range(&shphimin,&shphimax);
    2028   corsika=UInt_t(mcevth.get_VersionPGM()*1000);
    20292110  if(!Trigger_Loop) icontrigger=0;
    20302111  time (&ltime);
    20312112  ftime = ((Float_t)ltime)/1000;
    2032   for (int i=0; i< 10;i++)
    2033     heights[i]=mcevth.get_HeightLev (i);
    2034   rnum=mcevth.get_RunNumber();
    2035 
    2036   McRunHeader->Fill(rnum,
     2113
     2114  if (reflector_file_version<6)
     2115    McRunHeader->Fill(rnum,
    20372116                    (UInt_t) 0,
    20382117                    mcevth.get_DateRun(),
    … …  
    20472126                    CAMERA_PIXELS,
    20482127                    (UInt_t)ntshow,
    2049                     (UInt_t)ntrigger,
     2128                    (UInt_t)nstoredevents,
    20502129                    0,
    20512130                    sfRaH,
    … …  
    20702149                    heights,
    20712150                    corsika,
    2072                     (UInt_t)(REFL_VERSION_A*100),
     2151                    (UInt_t)(reflector_file_version*100),
    20732152                    (UInt_t)(VERSION*100),
    20742153                    0);
     2154  else
     2155    McRunHeader->Fill(rnum,
     2156                    (UInt_t) 0,
     2157                    mcevth_2.get_DateRun(),
     2158                    ftime,
     2159                    icontrigger,
     2160                    !Write_All_Images,
     2161                    Write_McEvt,
     2162                    Write_McTrig,
     2163                    Write_McFADC,
     2164                    Write_RawEvt,
     2165                    addElecNoise,
     2166                    CAMERA_PIXELS,
     2167                    (UInt_t)ntshow,
     2168                    (UInt_t)nstoredevents,
     2169                    0,
     2170                    sfRaH,
     2171                    sfRaM,
     2172                    sfRaS,
     2173                    sfDeD,
     2174                    sfDeM,
     2175                    sfDeS,
     2176                    meanNSB,
     2177                    telestheta,
     2178                    telesphi,
     2179                    sofftheta,
     2180                    soffphi,
     2181                    shthetamax,
     2182                    shthetamin,
     2183                    shphimax,
     2184                    shphimin,
     2185                    mcevth_2.get_CWaveLower(),
     2186                    mcevth_2.get_CWaveUpper(),
     2187                    mcevth_2.get_slope(),
     2188                    1,
     2189                    heights,
     2190                    corsika,
     2191                    (UInt_t)(reflector_file_version*100),
     2192                    (UInt_t)(VERSION*100),
     2193                    0);
     2194  // Fill some missing values for MRawRunHeader
     2195
     2196  RunHeader->SetRunNumber(rnum);
     2197  RunHeader->SetNumEvents(nstoredevents);
     2198
    20752199
    20762200  //  Store qe for each PMT in output file
    … …  
    40284152//
    40294153// $Log: not supported by cvs2svn $
     4154// Revision 1.50  2003/01/07 16:33:31  blanch
     4155// Star Field Rotation has been implemented by Raul Orduna. Now there is a
     4156// rotation for each shower. It is done by a non enter pixel rotation assuming
     4157// a circular symetry of the camera. It is not exact but it is accurate enough and
     4158// much faster.
     4159//
    40304160// Revision 1.49  2002/12/13 10:04:07  blanch
    40314161// *** empty log message ***
    … …  
    41794309//
    41804310// $Log: not supported by cvs2svn $
     4311// Revision 1.50  2003/01/07 16:33:31  blanch
     4312// Star Field Rotation has been implemented by Raul Orduna. Now there is a
     4313// rotation for each shower. It is done by a non enter pixel rotation assuming
     4314// a circular symetry of the camera. It is not exact but it is accurate enough and
     4315// much faster.
     4316//
    41814317// Revision 1.49  2002/12/13 10:04:07  blanch
    41824318// *** empty log message ***
Note: See TracChangeset for help on using the changeset viewer.