Changeset 383


Ignore:
Timestamp:
03/24/00 18:10:46 (25 years ago)
Author:
blanch
Message:
A first FADC simulation and a trigger simulation are already implemented.
The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
A loop over trigger threshold and some corretcions in the time range where it looks for a trigger will be implemented as soon as possible.
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r379 r383  
    2121//
    2222// $RCSfile: camera.cxx,v $
    23 // $Revision: 1.6 $
     23// $Revision: 1.7 $
    2424// $Author: blanch $
    25 // $Date: 2000-03-20 18:35:11 $
     25// $Date: 2000-03-24 18:10:46 $
    2626//
    2727////////////////////////////////////////////////////////////////////////
    … …  
    4848
    4949#include "TROOT.h"
     50
     51#include "TApplication.h"
    5052
    5153#include "TFile.h"
    … …  
    5456#include "TCanvas.h"
    5557
    56 #include "MDiag.h"
     58#include "MTrigger.hxx"
     59#include "MFadc.hxx"
    5760
    5861#include "MRawEvt.h"
    … …  
    6063#include "MMcTrig.hxx"
    6164
    62 #include "MTrigger.hxx"
    63 
    6465/*!@"
    6566
    … …  
    6970
    7071#include "camera.h"
     72
    7173//!@}
    7274
    … …  
    8890#undef  __DEBUG__
    8991
     92
    9093//!@}
    9194
    … …  
    196199static int Write_All_Data = FALSE;
    197200
     201static int Write_McEvt  = TRUE;
     202static int Write_McTrig = FALSE;
     203static int Write_RawEvt = FALSE;
     204
    198205//@: flag: TRUE: selection on the energy
    199206static int Select_Energy = TRUE;
    … …  
    204211//@: Upper edge of the selected energy range (in GeV)
    205212static float Select_Energy_ue = 100000.0;
     213
     214//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
     215static int FADC_Scan = FALSE;
     216
     217//@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't
     218static int Trigger_Scan = FALSE;
    206219
    207220//!@}
    … …  
    319332  @"*/
    320333
    321 /*!@"
    322 
    323   The following are the set of parameters calculated for each image.
    324   The routines for their calculations are in |moments.cxx|.
    325 
    326   @"*/
    327 
    328 //!@{
    329 // parameters of the images
    330 
    331 static Moments_Info *moments_ptr;
    332 static LenWid_Info *lenwid_ptr;
    333 
    334 static float *maxs;
    335 static int *nmaxs;
    336 static float length, width, dist, xdist, azw, miss, alpha, *conc;
    337 static float phiasym, asymx, asymy;
    338 static float charge, smax, maxtrigthr_phe;
    339334
    340335//!@}
    … …  
    360355  char inname[256];           //@< input file name
    361356  char starfieldname[256];    //@< starfield input file name
     357
    362358  char datname[256];          //@< data (ASCII) output file name
    363359  char diagname[256];         //@< diagnistic output file (ROOT format)
     360
    364361  char rootname[256] ;        //@< ROOT file name
    365362
    366363  char parname[256];          //@< parameters file name
     364
     365  char sign[20];              //@< initialize sign
    367366
    368367  char flag[SIZE_OF_FLAGS + 1];  //@< flags in the .rfl file
    … …  
    372371
    373372  MCEventHeader mcevth;       //@< Event Header class (MC)
     373  MCCphoton cphoton;          //@< Cherenkov Photon class (MC)
    374374
    375375  Photoelectron *photoe = NULL; //@< array of the photoelectrons of one event
    … …  
    421421  int itrigger;               //@< index of pixel fired
    422422  int ntrigger = 0;           //@< number of triggers in the whole file
    423   //  MTrigger  Trigger ;         //@< A instance of the Class MTrigger
    424 
    425   //    MMcTrig *McTrig   = new MMcTrig() ;
    426423
    427424  float plateScale_cm2deg;    //@< plate scale (deg/cm)
    … …  
    508505  Write_All_Data = get_write_all_data();
    509506
     507  Write_McEvt  = get_write_McEvt()  ;
     508  Write_McTrig = get_write_McTrig() ;
     509  Write_RawEvt = get_write_RawEvt() ;
     510
     511  FADC_Scan = get_FADC_Scan();
     512  Trigger_Scan = get_Trigger_Scan();
     513
    510514  // get filenames
    511515
    … …  
    550554      "Write_All_Images",  ONoff(Write_All_Images),
    551555      "Write_All_Data",    ONoff(Write_All_Data));
     556
     557  // log flags information
     558
     559  log(SIGNATURE,
     560      "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
     561      "Root output",
     562      "Write_McEvt",   ONoff(Write_McEvt), 
     563      "Write_McTrig",  ONoff(Write_McTrig), 
     564      "Write_RawEvt",  ONoff(Write_RawEvt));
    552565     
    553566  // log parameters information
    … …  
    597610  photoe = new Photoelectron[iMAXNUMPHE];
    598611
     612  Int_t Lev0, Lev1, Lev2 ;
     613
    599614  // initialise ROOT
    600615
    601616  TROOT simple("simple", "MAGIC Telescope Monte Carlo");
    602617
    603   // prepare ROOT tree for the diagnostic data
    604  
    605   TFile *hfile;
    606 
    607   hfile = new TFile( diagname,"RECREATE", "MAGIC Telescope MC diagnostic data");
    608  
    609   // Create the ROOT Tree for the diagnostic data
    610  
    611   TTree *tree = new TTree("T","MAGIC Telescope MC diagnostic data");
    612   tree->SetAutoSave(100000000);
    613  
    614   Int_t split = 1;
    615   Int_t bsize = 64000;
    616   MDiagEventobject  *event = 0;
    617  
    618   // Create one branch. If splitlevel is set, event is a superbranch
    619   // creating a sub branch for each data member of the Eventobject event.
    620  
    621   tree->Branch("event", "MDiagEventobject", &event, bsize, split);
     618  // initialise instance of Trigger and FADC classes
     619
     620  MTrigger  Trigger ;         //@< A instance of the Class MTrigger
     621
     622  MMcTrig *McTrig   = new MMcTrig() ;
     623
     624  MFadc fadc ;                //@< A instance of the Class MFadc
     625
    622626
    623627  // Prepare the raw data output
    … …  
    635639  TTree EvtTree("EvtTree","Events of Run");
    636640
    637   bsize=128000; split=1;
     641  Int_t bsize=128000; Int_t split=1;
    638642
    639643  EvtTree.Branch("MRawEvt","MRawEvt",
    … …  
    643647                 &McEvt, bsize, split);
    644648
    645 
     649  EvtTree.Branch("MMcTrig","MMcTrig",
     650                   &McTrig, bsize, split);
     651
     652  unsigned short ulli = 0 ;
     653
     654  TApplication theAppTrigger("App", &argc, argv);
     655
     656  if (gROOT->IsBatch()) {
     657    fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
     658    //    return 1;
     659  }
     660
     661  if(FADC_Scan){
     662  TApplication theAppFadc("App", &argc, argv);
     663
     664  if (gROOT->IsBatch()) {
     665    fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
     666    //    return 1;
     667  }
     668  }
     669 
    646670  // for safety and for dimensioning image_data: count the elements in the
    647671  // diagnostic data branch
    … …  
    720744    }
    721745 
    722 //     for(i=0; i<cam.inumpixels; i++){
    723 //       cout << i;
    724 //       for(j=0; j<iNUMWAVEBANDS; j++){
    725 //         cout << " " << j << " " << nsbrate_phepns[i][j];
    726 //       }
    727 //       cout << "\n";
    728 //     }
    729 
    730746    // calculate diffuse rate correcting for the pixel size
    731747
    … …  
    780796  fnpixclean = new float [ ct_NPixels ];
    781797
    782   moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
    783 
    784798  // initialize baseline
    785799
    … …  
    787801    baseline_mv[i] = 0.;
    788802  }
    789 
    790   //  Instance of the Trigger class
    791 
    792   MTrigger  Trigger ;         //@< A instance of the Class MTrigger
    793803       
    794804  //!@' @#### Main loop.
    … …  
    820830
    821831      while( isA( flag, FLAG_START_OF_EVENT   )){ // while there is a next event
    822        
     832
     833        //
     834        // Clear Trigger and Fadc
     835        //
     836        Trigger.Reset() ;
     837        fadc.Reset() ;
     838
    823839        ++nshow;
    824840        log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
    … …  
    9991015
    10001016        //   
    1001         //   remove the prvious values of the analog signal for
    1002         //   each pixel and put the new ones from photoe array
     1017        //   Put values of the analog signal for
     1018        //   each pixel from photoe array
    10031019        //
    1004         Trigger.Reset();
    10051020       
    10061021        for(i=0;i<inumphe;i++){
    1007           Trigger.FillShow(photoe[i].ipixnum,float((photoe[i].iarrtime_ns-arrtmin_ns+TOTAL_TRIGGER_TIME/10.0)));
    1008         }
     1022          Trigger.FillShow(photoe[i].ipixnum,float((photoe[i].iarrtime_ns-arrtmin_ns)));
     1023          fadc.Fill( photoe[i].ipixnum,(photoe[i].iarrtime_ns-arrtmin_ns) , Trigger.FillShow(photoe[i].ipixnum,float((photoe[i].iarrtime_ns-arrtmin_ns))));
     1024}
    10091025
    10101026        //
    … …  
    10141030        //   
    10151031        Trigger.ElecNoise() ;
    1016        
     1032        fadc.ElecNoise() ;
     1033       
     1034
     1035        Trigger.Diskriminate() ;
     1036
    10171037        //
    10181038        //   look if in all the signals in the trigger signal branch
    … …  
    10241044        //
    10251045
    1026         Trigger.Diskriminate();
    1027 
    1028         //McTrig->SetZeroLevel(trigger = (Short_t) Trigger.ZeroLevel()) ;
    1029         trigger = (Short_t) Trigger.ZeroLevel();
     1046        McTrig->SetZeroLevel( Lev0 = (Short_t) Trigger.ZeroLevel() ) ;
     1047
     1048        Lev1 = Lev2 = 0 ;
    10301049       
    10311050        //
    10321051        //   Start the First Level Trigger simulation
    10331052        //
    1034         if ( trigger > 0 ) {
    1035           //      McTrig->SetFirstLevel ( trigger = Trigger.FirstLevel() )  ;
    1036           trigger = Trigger.FirstLevel();
    1037         }
    1038         if ( trigger>0 ) {
    1039           cout << "TRIG"<<"  "<<thetashw <<"  "<< trigger  << endl;
     1053
     1054        if ( Lev0 > 0 ) {
     1055          McTrig->SetFirstLevel ( Lev1 = Trigger.FirstLevel() )  ;
    10401056        }
    10411057
    … …  
    10441060//      }
    10451061       
    1046         if ( trigger>0 ) {
    1047          
     1062        //if ( trigger>0 ) {
     1063        if (Lev1>0){
     1064       
    10481065          itrigger = i;
    10491066          ++ntrigger;
    10501067         
    10511068          memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels );
    1052          
    1053          
    1054 #ifdef __ISLANDS__
    1055          
    1056           //!@' @#### Islands algorithm.
    1057           //@'
    1058          
    1059           //++
    1060           // islands counting, and cleanning
    1061           //--
    1062          
    1063           if ( countIslands )
    1064             do_islands( ct_NPixels, fnpixclean, pixneig, npixneig,
    1065                         countIslands, nIslandsCut);
    1066          
    1067 #endif // __ISLANDS__
    1068          
    1069          
    1070           //!@' @#### Calculation of parameters of the image.
    1071           //@'
    1072          
    1073           //++
    1074           // moments calculation
    1075           //--
    1076          
    1077           // calculate moments and other things
    1078          
    1079           moments_ptr = moments( anaPixels, fnpixclean, pixary,
    1080                                  plateScale_cm2deg, 0 );
    1081          
    1082           charge = moments_ptr->charge ;
    1083           smax   = moments_ptr->smax   ;
    1084           maxs   = moments_ptr->maxs   ;
    1085           nmaxs  = moments_ptr->nmaxs  ;
    1086           length = moments_ptr->length ;
    1087           width  = moments_ptr->width  ;
    1088           dist   = moments_ptr->dist   ;
    1089           xdist  = moments_ptr->xdist  ;
    1090           azw    = moments_ptr->azw    ;
    1091           miss   = moments_ptr->miss   ;
    1092           alpha  = moments_ptr->alpha  ;
    1093           conc   = moments_ptr->conc   ;
    1094           asymx  = moments_ptr->asymx  ;
    1095           asymx  = moments_ptr->asymx  ;
    1096           phiasym= moments_ptr->phi;
    1097          
    1098           lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary,
    1099                                plateScale_cm2deg,
    1100                                ct_PixelWidth_corner_2_corner_half);
    1101          
    1102          
    1103           // fill the diagnostic Tree
    1104          
    1105           event = new MDiagEventobject();
    1106          
    1107           i=0;
    1108           image_data[i] =       event->n       = hidt/10; i++;
    1109           image_data[i] =       event->primary = mcevth.get_primary(); i++;
    1110           image_data[i] =       event->energy  = mcevth.get_energy(); i++;
    1111           image_data[i] =       event->cored   = coreD;  i++;
    1112           image_data[i] =       event->impact  = impactD; i++;
    1113           image_data[i] =       event->xcore   = coreX; i++;
    1114           image_data[i] =       event->ycore   = coreY; i++;
    1115           image_data[i] =       event->theta   = mcevth.get_theta(); i++;
    1116           image_data[i] =       event->phi     = mcevth.get_phi(); i++;
    1117           image_data[i] =       event->deviations = mcevth.get_deviations (&dtheta, &dphi); i++;
    1118           image_data[i] =       event->dtheta  = dtheta; i++;
    1119           image_data[i] =       event->dphi    = dphi; i++;
    1120           image_data[i] =       event->trigger = trigger; i++;
    1121           image_data[i] =       event->ncphs   = ncph; i++;
    1122           image_data[i] =       event->maxpassthr_phe   = maxtrigthr_phe; i++;
    1123           image_data[i] =       event->nphes   = charge; i++;
    1124           image_data[i] =       event->nphes2  = smax; i++;
    1125           image_data[i] =       event->length  = length; i++;
    1126           image_data[i] =       event->width   = width; i++;
    1127           image_data[i] =       event->dist    = dist; i++;
    1128           image_data[i] =       event->xdist   = xdist; i++;
    1129           image_data[i] =       event->azw     = azw; i++;
    1130           image_data[i] =       event->miss    = miss; i++;
    1131           image_data[i] =       event->alpha   = alpha; i++;
    1132           image_data[i] =       event->conc2   = conc[0]; i++;
    1133           image_data[i] =       event->conc3   = conc[1]; i++;
    1134           image_data[i] =       event->conc4   = conc[2]; i++;
    1135           image_data[i] =       event->conc5   = conc[3]; i++;
    1136           image_data[i] =       event->conc6   = conc[4]; i++;
    1137           image_data[i] =       event->conc7   = conc[5]; i++;
    1138           image_data[i] =       event->conc8   = conc[6]; i++;
    1139           image_data[i] =       event->conc9   = conc[7]; i++;
    1140           image_data[i] =       event->conc10  = conc[8]; i++;
    1141           image_data[i] =       event->asymx   = asymx; i++;
    1142           image_data[i] =       event->asymy   = asymy; i++;
    1143           image_data[i] =       event->phiasym = phiasym; i++;
    1144          
    1145           // there should be "nvar" variables
    1146          
    1147           if ( i != nvar )
    1148             error( SIGNATURE, "Wrong entry number for diagnostic data.\n" );
    1149          
    1150           tree->Fill();
    1151           delete event;
    1152          
    1153           // put information in the data file,
    1154          
    1155           datafile << ntrigger;
    1156           for(i=0;i<nvar;i++) {
    1157             datafile << ' ' << image_data[i];
     1069        }
     1070        //
     1071        //  Fill the header of this event
     1072        //
     1073       
     1074        Evt->FillHeader ( (UShort_t) (ntshow + nshow) ,  20 ) ;
     1075       
     1076        //
     1077        //   fill the MMcEvt with all information 
     1078        //
     1079       
     1080        McEvt->Fill( (UShort_t) mcevth.get_primary() ,
     1081                     mcevth.get_energy(),
     1082                     mcevth.get_theta(),
     1083                     mcevth.get_phi(),
     1084                     mcevth.get_core(),
     1085                     mcevth.get_coreX(),
     1086                     mcevth.get_coreY(),
     1087                     impactD,
     1088                     ulli, ulli,
     1089                     (UShort_t) ncph,
     1090                     ulli,
     1091                     (UShort_t) ncph) ;
     1092
     1093        //   We don not count phtons out of the camera.
     1094
     1095        //
     1096        //    write it out to the file outfile
     1097        //
     1098       
     1099        EvtTree.Fill() ;
     1100       
     1101
     1102        //
     1103        //    if a first level trigger occurred, then
     1104        //       1. do some other stuff (not implemented)
     1105        //       2. start the gui tool
     1106
     1107        if(FADC_Scan){
     1108          if ( Lev0 > 0 ) {
     1109            fadc.ShowSignal( McEvt, (Float_t) 60. ) ;
    11581110          }
    1159          
    1160           // revert the fnpixclean matrix into fnpix
    1161           // (now we do this, but maybe in a future we want to
    1162           // use both fnpix and fnpixclean for different things
    1163          
    1164           memcpy( fnpix, fnpixclean, sizeof(float) * ct_NPixels );
    1165          
    1166           // put this information in the data file,
    1167          
    1168           if ( Write_All_Data ) {
    1169             datafile << ' ' << -9999;
    1170             for ( i=0; i<ct_NPixels; ++i )
    1171               datafile << ' ' << fnpix[i];
     1111        }
     1112
     1113        if(Trigger_Scan){
     1114          if ( Lev0 > 0 ) {
     1115            Trigger.ShowSignal(McEvt) ;
    11721116          }
    1173          
    1174           datafile << endl; 
    1175          
    1176           mcevth.set_trigger( TRUE );
    1177          
    1178           log(SIGNATURE, "TRIGGER\n");
    1179          
    1180         } else { // ( trigger == FALSE )
    1181          
    1182           event = new MDiagEventobject();
    1183          
    1184           i=0;
    1185           image_data[i] =       event->n       = hidt/10; i++;
    1186           image_data[i] =       event->primary = mcevth.get_primary(); i++;
    1187           image_data[i] =       event->energy  = mcevth.get_energy(); i++;
    1188           image_data[i] =       event->cored   = coreD = mcevth.get_core(&coreX, &coreY); i++;
    1189           image_data[i] =       event->impact  = coreD; i++;
    1190           image_data[i] =       event->xcore   = coreX; i++;
    1191           image_data[i] =       event->ycore   = coreY; i++;
    1192           image_data[i] =       event->theta   = mcevth.get_theta(); i++;
    1193           image_data[i] =       event->phi     = mcevth.get_phi(); i++;
    1194           image_data[i] =       event->deviations = mcevth.get_deviations(&dtheta, &dphi); i++;
    1195           image_data[i] =       event->dtheta = dtheta; i++;
    1196           image_data[i] =       event->dphi = dphi; i++;
    1197           image_data[i] =       event->trigger = trigger; i++;
    1198           image_data[i] =       event->ncphs   = ncph; i++;
    1199           image_data[i] =       event->maxpassthr_phe   = maxtrigthr_phe; i++;
    1200           image_data[i] =       -1.; i++;
    1201           image_data[i] =       -1.; i++;
    1202           image_data[i] =       -1.; i++;
    1203           image_data[i] =       -1.; i++;
    1204           image_data[i] =       -1.; i++;
    1205           image_data[i] =       -1.; i++;
    1206           image_data[i] =       -1.; i++;
    1207           image_data[i] =       -1.; i++;
    1208           image_data[i] =       -1.; i++;
    1209           image_data[i] =       -1.; i++;
    1210           image_data[i] =       -1.; i++;
    1211           image_data[i] =       -1.; i++;
    1212           image_data[i] =       -1.; i++;
    1213           image_data[i] =       -1.; i++;
    1214           image_data[i] =       -1.; i++;
    1215           image_data[i] =       -1.; i++;
    1216           image_data[i] =       -1.; i++;
    1217           image_data[i] =       -1.; i++;
    1218           image_data[i] =       -1.; i++;
    1219           image_data[i] =       -1.; i++;
    1220           image_data[i] =       -1.; i++;
    1221          
    1222           // there should be "nvar" variables
    1223          
    1224           if ( i != nvar )
    1225             error( SIGNATURE, "Wrong entry length for Ntuple.\n" );
    1226          
    1227           tree->Fill();
    1228           delete event;
    1229          
    1230           // put this information in the data file,
    1231          
    1232           if ( Write_All_Data ) {
     1117        }
     1118
     1119        //    clear all
     1120        Evt->Clear() ;
     1121        McEvt->Clear() ;
     1122        McTrig->Clear() ;
     1123     
     1124       
     1125        //++++++++++++++++++++++++++++++++++++++++++++++++++
     1126        // at this point we have a camera full of
     1127        // ph.e.s
     1128        // we should first apply the trigger condition,
     1129        // and if there's trigger, then clean the image,
     1130        // calculate the islands statistics and the
     1131        // other parameters of the image (Hillas' parameters
     1132        // and so on).
     1133        //--------------------------------------------------
     1134       
     1135#ifdef __DEBUG__ 
     1136        printf("\n");
     1137       
     1138        for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
     1139         
     1140          for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
    12331141           
    1234             datafile << ntrigger;
    1235             for ( i=0; i<nvar; ++i )
    1236               datafile << ' ' << image_data[i];
    1237            
    1238             datafile << -9999;
    1239             for ( i=0; i<ct_NPixels; ++i )
    1240               datafile << ' ' << fnpix[i];
    1241            
    1242             datafile << endl; 
     1142            if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
     1143             
     1144              if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
     1145               
     1146                printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
     1147                        (int)pixels[ici][icj][PIXNUM],
     1148                        pixels[ici][icj][PIXX],
     1149                        pixels[ici][icj][PIXY],
     1150                        fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
     1151               
     1152              }
     1153            } 
    12431154          }
    1244          
    1245           mcevth.set_trigger( FALSE );
    1246          
    1247         } // trigger == FALSE
    1248                
     1155        }
     1156       
     1157        for (i=0; i<ct_NPixels; ++i) {
     1158          printf("%d (%d): ", i, npixneig[i]);
     1159          for (j=0; j<npixneig[i]; ++i)
     1160            printf(" %d", pixneig[i][j]);
     1161          printf("\n");
     1162        }
     1163       
     1164#endif // __DEBUG__
     1165                         
     1166       
    12491167        //!@' @#### Save data.
    12501168        //@'
    … …  
    13241242  }
    13251243  datafile.close();
    1326 
    1327   hfile->Write();
    1328 
    1329   hfile->Close();
    13301244
    13311245  // program finished
    … …  
    28272741//
    28282742// $Log: not supported by cvs2svn $
     2743// Revision 1.6  2000/03/20 18:35:11  blanch
     2744// The trigger is already implemented but it does not save the trigger information in any file as it is implemented in timecam. In the next days there will be a version which also creates the files with the trigger information. It is going to be a mixing of the current camera and timecam programs.
     2745//
    28292746// Revision 1.5  2000/02/18 17:40:35  petry
    28302747// This version includes drastic changes compared to camera.cxx 1.4.
Note: See TracChangeset for help on using the changeset viewer.