source: trunk/MagicSoft/Simulation/Corsika/Checkmc/checkmc.cxx@ 9591

Last change on this file since 9591 was 328, checked in by magiccvs, 26 years ago
This is the starting point for further development of the checkmc program. This program is based on the work of JC Gonzales, D Petry & P Morawitz. I hope it is running under linux and osf/alphas. You have to set the right values in the config.mk.* file and the Makefile. Now this program is also under CVS control.
File size: 41.4 KB
Line 
1//////////////////////////////////////////////////////////////////
2//
3// checkmc
4//
5// @file checkmc.cxx
6// @title Check of the MC generated data
7// @author J C Gonz\'alez, P. Morawitz, D. Petry
8// @email gonzalez@mppmu.mpg.de, morawitz@cern.ch, petry@ifae.es
9//
10// @maintitle
11//_______________________________________________________________
12//
13// Created: Thu May 7 16:24:22 1998
14// Author: Jose Carlos Gonzales et al.
15// Purpose: Check of the MC generated data
16// Notes: version 27-10-98
17//
18//
19
20// Below you can find a sort of, incomplete for sure, documentation
21// concerning this program. Please, do not hesiate to contact the
22// author in case of problems.
23
24// @T \newpage
25
26// @section Source code of {\tt checkmc.cxx}
27
28/* @text
29In this section we show the (commented) code of the program
30for for the check of the MC generated data, in the current version.
31@endtext */
32
33// @subsection Includes and Global variables definition
34
35/* @text
36All the defines are located in the file {\tt checkmc.h}.
37@endtext */
38
39// @code
40#include "checkmc.h"
41// @endcode
42
43// @subsubsection Definition of global variables.
44
45/* @text
46Here we define the global variables where the values from
47the parameters file are stored. See the section
48\ref{paramdesc} in page \pageref{paramdesc} for a more
49complete description of the available commands in the
50parameters file.
51@endtext */
52
53// @code
54static char **Paths_list; // list of paths to get data from
55static int Num_of_paths = 1; // number of paths
56static char Output_filename[PATH_MAX_LENGTH]; // output filename
57static char CT_filename[PATH_MAX_LENGTH]; // name of the CT def. file
58static VerboseLevel verbose = VERBOSE_DEFAULT; // verboseness
59static float fixed_Theta; // fixed target theta and phi
60static float fixed_Phi;
61static float low_Ecut = 0.0, high_Ecut = 100000.0;
62static int is_Fixed_Target = FALSE; // are we using fixed target?
63static int max_Events = 100000; // maximum number of events to read
64
65static float gammas = 1.;
66static float protons = 14.;
67static float helium = 402.;
68
69
70inline int max(int a, int b)
71{if (a>b) return a;
72 else return b;}
73
74inline int min(int a, int b)
75{if (a<b) return a;
76 else return b;}
77
78// @endcode
79
80/* @text
81Earth's radius (in cm).
82@endtext */
83
84// @code
85static const float REarth = 6348.0e5;
86// @endcode
87
88// @subsection Main program
89
90/* @text
91Let's start!
92@endtext */
93
94// @code
95//++++++++++++++++++++++++++++++++++++++++
96// MAIN PROGRAM
97//----------------------------------------
98
99void
100main(int argc, char **argv)
101{
102
103 FILE *f;
104 DIR *directory; // directory of data
105 struct dirent *de; // directory entry
106
107 char pathname[256]; // directory and file names
108 char cername[256];
109 char staname[256];
110 char outname[256];
111 char refname[256];
112 char name[256];
113
114 ifstream cerfile;
115 ofstream outputfile;
116
117 COREventHeader evth; // Event Header class (from CORSIKA)
118 CORParticle photon; // Particle (photon) class (from CORSIKA)
119 CORStatfile stat; // Statistics file class (from CORSIKA)
120
121 float thetaCT, phiCT, xiCT; // parameters of a given shower
122 float coreD, coreX, coreY; // core position and distance
123
124 float a, b, c, t, tt; // intermediate variables
125
126 float mind, d; // minimum distance trayectory-mirror center
127
128 float wl; // wavelength of the Cphoton
129 float reflec; // reflectivity for a Cphoton
130
131 int ncer; // number of the shower;
132 int np, nf; // number of path, number of cer-file
133 int i, j, k, n; // simple counters
134 int i_mirror=-1; // number of a given mirror
135
136 int nCphotons; // number of a Cphotons written on output file
137
138 float TR; // atmospheric transmitance
139 int nbeforeTR, nafterTR; // number of Cph. before and after transmission
140 int num_cer_files; // number of cer files;
141 int max_num_cer_files; // maximum number of cer files to read
142
143 float t1, t2; // first and last time
144
145 float x0, xa, xb, v; // variables for the
146 float vmax, vmaxa, vmaxb; // estimation of Xmax
147 float xmax;
148
149 int ntotfiles;
150 int kbin;
151
152 int hprimary;
153 static const int nprimaries=4;
154 char* primary_name[nprimaries] ={
155 "Photons",
156 "Protons",
157 "Helium",
158 "Heavy Nuclei"};
159 int nevents[nprimaries];
160 char htit[256];
161
162 int firstevent,lastevent;
163 int ndiag=-1;
164 int rndiag=-1;
165 char diag[256][20];
166 char rdiag[256][100];
167 char ctmp[256];
168 float par[3],vp2[4],vp3[5],vp4[6];
169 float sigpar;
170 float chi2;
171 float rnd[3][3];
172
173 // now i/p the constants to be checked
174 const float mon_HeightLev=220000.;
175 const float mon_SlopeSpec[nprimaries] = {-1.5,-2.75,-2.62,-100000.};
176 const float mon_ElowLim [nprimaries] = {1.,10.,10.,10.};
177 const float mon_EUppLim[nprimaries] = {30000.,100000.,100000.,100000.};
178 const float mon_Ecutoffh[nprimaries] = {0.3,0.3,0.3,0.3};
179 const float mon_Ecutoffm[nprimaries] = {0.3,0.3,0.3,0.3};
180 const float mon_Ecutoffe[nprimaries] = {0.02,0.02,0.02,0.02};
181 const float mon_Ecutoffg[nprimaries] = {0.02,0.02,0.02,0.02};
182 const float mon_EGS4yn[nprimaries] = {1.,1.,1.,1.};
183// the EGS4yn flag should be automatically
184 const float mon_NKGyn[nprimaries] = {0.,0.,0.,0.};
185 const float mon_GHEISHAyn[nprimaries] = {1.,1.,1.,1.};
186 const float mon_VENUSyn[nprimaries] = {1.,1.,1.,1.};
187 const float mon_CERENKOVyn[nprimaries] = {1.,1.,1.,1.};
188 const float mon_ThetaMin[nprimaries] = {5.,0.,0.,0.};
189 const float mon_ThetaMax[nprimaries] = {25.,30.,30.,30.};
190 const float mon_PhiMin[nprimaries] = {0.,0.,0.,0.};
191 const float mon_PhiMax[nprimaries] = {360.,360.,360.,360.};
192 const float mon_CWaveLower[nprimaries] = {290.,290.,290.,290.};
193 const float mon_CWaveUpper[nprimaries] = {600.,600.,600.,600.};
194
195 // @endcode
196
197 // @text
198 // Variables for the primary particle N-Tuple
199 // @endtext
200
201 // @code
202 int hid=1,istat=0,icycle=0;
203 const int nvar=44;
204 char chTags[nvar][8]={
205 "primary",
206// primary particle id
207 "energy",
208// primary energy
209 "cored",
210// Core distance in cm
211 "theta",
212 "phi",
213 "thick",
214// thickness (g/cm-2), e.g. height z converted, see below
215 "height",
216// height, z coordinate, cm, of first interaction
217 "target1",
218// number of first taget ?
219 "px",
220 "py",
221 "pz",
222// the momentum of the primary particle
223 "tfirst",
224// time of first cerenkov photon on ground
225 "tlast",
226// time of last cerenkov photon on ground
227 "xmax",
228// xmax in g/cm**2; "height" of cerenkov shower maximum
229// => xmax/dichte == penetration depth
230// depth defined as == 120km=0
231// (120km - penetration depth - observation level (2km)) == real height of interaction
232//
233 "xcore",
234// x position of core on ground in cm
235 "ycore",
236// y position of core on ground in cm
237 "Nph",
238// number of cerenkov photons
239 "RunNum",
240 "EvtNum",
241 "DateRun",
242 "Version",
243 "HeigLev",
244 "ESlope",
245 "ELowLim",
246 "EUppLim",
247//
248 "Ecutofh",
249 "Ecutofm",
250 "Ecutofe",
251 "Ecutofg",
252//
253 "EGS4yn",
254 "NKGyn",
255 "GHEISHA",
256 "VENUS",
257 "CERENKO",
258 "NEUTRIN",
259 "HORIZON",
260 "COMPUTE",
261 "ThetMin",
262 "ThetMax",
263 "PhiMin",
264 "PhiMax",
265//
266 "StepLen",
267 "CWavLow",
268 "CWavUpp"
269 };
270 // char chTags[nvar * 10];
271 float rec[nvar];
272 int record_size=4096;
273 // @endcode
274
275
276 // @text
277 // Variables for the cerenkov photon N-Tuple
278 // @endtext
279
280 // @code
281 int hid2=2,istat2=0,icycle2=0;
282 const int nvar2=11;
283 char chTags2[nvar2][8]={
284//
285// First info about primary
286 "energy",
287 "cored",
288 "theta",
289//
290// Now cerenkov photon info
291 "Wavel",
292// wavelength in nm
293 "xph",
294 "yph",
295 "uph",
296 "vph",
297 "tph",
298 "hph",
299 "EvtNum"
300 };
301 // char chTags2[nvar2 * 10];
302 float rec2[nvar2];
303 int record_size2=4096;
304 // @endcode
305
306
307 // @text
308 // We start in the main program. First we (could) make some presentation,
309 // and follows the reading of the parameters file (now from the {\tt stdin}),
310 // the reading of the CT parameters file, and the creation of the
311 // output file, where the processed data will be stored.
312 // @endtext
313 // @code
314
315 //++
316 // START
317 //--
318
319 // make unbuffered output
320
321 cout.setf ( ios::stdio );
322
323 // HBOOK core memory
324
325 HLIMIT(PAWC_SIZE);
326
327 /*
328 strcpy( chTags, "primary:" );
329 strcat( chTags, "energy:" );
330 strcat( chTags, "cored:" );
331 strcat( chTags, "theta:" );
332 strcat( chTags, "phi:" );
333 strcat( chTags, "thick:" );
334 strcat( chTags, "height:" );
335 strcat( chTags, "target1:" );
336 strcat( chTags, "px:" );
337 strcat( chTags, "py:" );
338 strcat( chTags, "pz:" );
339 strcat( chTags, "tfirst:" );
340 strcat( chTags, "tlast:" );
341 strcat( chTags, "xmax" );
342 */
343
344 // make some sort of presentation
345
346 present();
347
348 // read parameters from the stdin
349
350 readparam();
351 if ( verbose >= VERBOSE_MINIMAL )
352 log( SIGNATURE, "... reading parameters OK\n");
353 // get name of the output file
354
355 strcpy( outname, get_output_filename() );
356
357 // open HBOOK file
358
359 if ( verbose >= VERBOSE_MINIMAL )
360 log( SIGNATURE, "Openning N-Tuple file %s\n", outname );
361
362 // TFile hfile(outname,"RECREATE","Check MC generated data");
363
364
365 HROPEN(1, "CHECKMC", outname, "N", record_size, istat);
366 if (istat != 0)
367 {fprintf(stderr, "\n Output Histo couldn't be opened. Adios.\n");
368 exit(999);}
369
370 // profile (final) histograms
371
372 // TH1F *htime = new TH1F("htime", "htime", 60, 0., 60.);
373 // TH1F *helec = new TH1F("helec", "helec", NTHSTEP, 0., 790.);
374 // TH1F *hlight = new TH1F("hlight", "hlight", NTHSTEP, 0., 790.);
375
376 double fstime[18][60], fs2time[18][60];
377
378 double ftstime[18][60], fts2time[18][60];
379 double ftselec[18][NTHSTEP], fts2elec[18][60];
380 double ftslight[18][NTHSTEP], fts2light[18][60];
381
382 for (i=0; i<60; ++i)
383 for (k=0; k<18; ++k)
384 ftstime[k][i] = fts2time[k][i] = 0.0;
385
386 for (i=0; i<NTHSTEP; ++i)
387 for (k=0; k<18; ++k)
388 ftselec[k][i] = fts2elec[k][i] =
389 ftslight[k][i] = fts2light[k][i] = 0.0;
390
391 // @endcode
392 // @text
393 // After reading the parameters file (from {\tt stdin}, and the
394 // CT definition file, we start the analysis. Basically, there are
395 // three loops. The outermost loop is a loop in the {\it data directories\/}
396 // we introduced in the parameters file; the middle loop follows each
397 // {\tt cer*} file in the directory; and the innermost loop gets
398 // each Cherenkov from these files.
399 // @endtext
400 // @code
401
402 // TNtuple *ntuple = new TNtuple("ntuple","Check MC data", chTags);
403 // HBOOKN(hid, pathname, nvar, "CHECKMC", 4096, chTags);
404 // HBOOKN(hid2, pathname, nvar2, "CHECKMC", 4096, chTags2);
405
406
407
408 // First count the number of events for each primary
409 ntotfiles = 0;
410 // for each path of data files
411 for (np=0; np<get_num_of_paths(); np++) {
412 strcpy(pathname, get_path_name(np));
413 // open directory
414 if ( verbose >= VERBOSE_MINIMAL )
415 log( SIGNATURE, "Openning directory %s\n", pathname );
416 if ( (directory = opendir(pathname)) == NULL )
417 error( SIGNATURE,
418 "Cannot open directory %s\n", pathname );
419 // trace the number of events (cer* files)
420 num_cer_files = 0;
421 max_num_cer_files = get_max_events();
422 // for each directory entry (files)
423 while ( ( (de = readdir( directory )) != NULL) &&
424 (num_cer_files < max_num_cer_files) ) {
425 // skip removed files
426 if ( de->d_ino == 0)
427 continue;
428 // keep only cer* files
429 if ( strstr(de->d_name, "cer") != de->d_name )
430 continue;
431 // it is a real cer* file
432 ++num_cer_files;
433 if ( verbose >= VERBOSE_NORMAL )
434 log( SIGNATURE, "cerfile: %s\n", de->d_name );
435 // get cer* file number (number of the shower)
436 // increments the pointer in 3 to calculate the number
437 ncer = atoi(de->d_name + 3);
438 // full names of the input files cer* and sta*
439 sprintf(cername, "%s/cer%06d", pathname, ncer);
440 sprintf(staname, "%s/sta%06d", pathname, ncer);
441 // try to open cer* file
442 if ( verbose >= VERBOSE_MAXIMAL )
443 log( SIGNATURE, "Openning %s\n", cername );
444 cerfile.open( cername );
445 if ( cerfile.bad() )
446 error( SIGNATURE,
447 "Cannot open input file: %s\n", cername );
448 // try to open the statistics file
449 if ( verbose >= VERBOSE_MAXIMAL )
450 log( SIGNATURE, "Openning %s\n", staname );
451 stat.openfile ( staname );
452 // reading data from this sta* file
453 if ( verbose >= VERBOSE_MAXIMAL )
454 log( SIGNATURE, "Reading %s\n", staname );
455 stat.read();
456 // read the Event Header from the cer* file
457 evth.read( cerfile );
458 if (evth.PrimaryID == gammas)
459 ++nevents[0];
460 else if (evth.PrimaryID == protons)
461 ++nevents[1];
462 else if (evth.PrimaryID == helium)
463 ++nevents[2];
464 else
465 ++nevents[3];
466 cerfile.close();
467 stat.closefile();}}
468
469
470 cout<<"NMAXBINS"<<NMAXBINS<<endl;
471 // Now book the histogram files
472 for (i=0; i<nprimaries; ++i)
473 if (nevents[i]>0)
474// Primary histos
475 {sprintf(htit, "Log(Energy) of primary %s;log(E/GeV)", primary_name[i]);
476 HBOOK1(100+i,htit,min(nevents[i]/25,NMAXBINS),0.,4.5,0.);
477
478 sprintf(htit, "Core distance (horizontal plane) of primary %s;(m)", primary_name[i]);
479 HBOOK1(110+i,htit,min(nevents[i]/25,NMAXBINS),0.,500.,0.);
480
481 sprintf(htit, "Theta of primary %s;(degrees)", primary_name[i]);
482 HBOOK1(120+i,htit,min(nevents[i]/15,NMAXBINS),0.,70.,0.);
483
484 sprintf(htit, "Phi of primary %s;(degrees)", primary_name[i]);
485 HBOOK1(130+i,htit,min(nevents[i]/25,NMAXBINS),0.,360.,0.);
486
487 sprintf(htit, "Core position for primary %s;x (m);y (m)", primary_name[i]);
488 HBOOK2(140+i,htit,50,-500.,500.,50,-500.,500.,0.);
489
490 sprintf(htit, "Height of first interaction of primary %s;(km)", primary_name[i]);
491 HBOOK1(150+i,htit,min(nevents[i]/25,NMAXBINS),0.,100.,0.);
492
493 sprintf(htit, "Time (last-first) of C-Photons for %s;(ns)", primary_name[i]);
494 HBOOK1(160+i,htit,min(nevents[i]/25,NMAXBINS),-1.,160.,0.);
495
496 sprintf(htit, "Shower maximum XMAX for %s;(g cm^-2!)", primary_name[i]);
497 HBOOK1(170+i,htit,min(nevents[i]/10,NTHSTEP),0.,790.,0.);
498
499 sprintf(htit, "Log(Num of C-photons produced) for primary %s below 300 GeV;log(N)", primary_name[i]);
500 HBOOK1(180+i,htit,min(nevents[i]/25,NMAXBINS),-.5,6.,0.);
501
502 sprintf(htit, "Log(Num of C-photons produced) for primary %s above 300GeV;log(N)", primary_name[i]);
503 HBOOK1(190+i,htit,min(nevents[i]/25,NMAXBINS),-.5,6.,0.);
504
505// C-Photon histos
506 sprintf(htit, "C-Photon position in horiz. plane for prim. %s;x (m);y (m)", primary_name[i]);
507 HBOOK2(300+i,htit,50,-15.,15.,50,-15.,15.,0.);
508
509 sprintf(htit, "C-Photon distance to (0,0) (horizontal plane) for %s;(m)", primary_name[i]);
510 HBOOK1(310+i,htit,(min(nevents[i]/5,NMAXBINS)),0.,15.,0.);
511
512 sprintf(htit, "C-Photon spectrum for primary %s; [l] (nm)", primary_name[i]);
513 HBOOK1(320+i,htit,(min(nevents[i]/5,NMAXBINS)),200.,700.,0.);
514
515 sprintf(htit, "Production height of C-Photon for primary %s;(km)", primary_name[i]);
516 HBOOK1(330+i,htit,(min(nevents[i]/5,NMAXBINS)),0.,30.,0.);
517 }
518
519 ntotfiles = 0;
520
521 // for each path of data files
522
523 for (np=0; np<get_num_of_paths(); np++) {
524
525 strcpy(pathname, get_path_name(np));
526
527 // open directory
528
529 if ( verbose >= VERBOSE_MINIMAL )
530 log( SIGNATURE, "Openning directory %s\n", pathname );
531 if ( (directory = opendir(pathname)) == NULL )
532 error( SIGNATURE,
533 "Cannot open directory %s\n", pathname );
534
535 // trace the number of events (cer* files)
536
537 num_cer_files = 0;
538
539 // book n-tuple for this directory
540
541 // hid = np+1;
542
543 // if ( verbose >= VERBOSE_MINIMAL )
544 // log( SIGNATURE, "Create N-Tuple #%d for path %s . . .\n",
545 // hid, pathname );
546
547 // HBOOKN(hid, pathname, nvar, " ", 4096, chTags);
548
549 max_num_cer_files = get_max_events();
550
551 // for each directory entry (files)
552
553 firstevent=-100;
554 while ( ( (de = readdir( directory )) != NULL) &&
555 (num_cer_files < max_num_cer_files) ) {
556
557 // skip removed files
558
559 if ( de->d_ino == 0)
560 continue;
561
562 // keep only cer* files
563
564 if ( strstr(de->d_name, "cer") != de->d_name )
565 continue;
566
567 // it is a real cer* file
568
569 ++num_cer_files;
570
571 if ( verbose >= VERBOSE_NORMAL )
572 log( SIGNATURE, "cerfile: %s\n", de->d_name );
573
574 // get cer* file number (number of the shower)
575
576 // increments the pointer in 3 to calculate the number
577
578 ncer = atoi(de->d_name + 3);
579
580 // full names of the input files cer* and sta*
581
582 sprintf(cername, "%s/cer%06d", pathname, ncer);
583 sprintf(staname, "%s/sta%06d", pathname, ncer);
584
585 // try to open cer* file
586
587 if ( verbose >= VERBOSE_MAXIMAL )
588 log( SIGNATURE, "Openning %s\n", cername );
589
590 cerfile.open( cername );
591
592 if ( cerfile.bad() )
593 error( SIGNATURE,
594 "Cannot open input file: %s\n", cername );
595
596 // @endcode
597 // @text
598 // Each shower has associated three files: {\em Cherenkov-photons file\/}
599 // ({\tt cer*}), {\em Particles file} ({\tt dat*}) and
600 // {\em Statistics file} ({\tt sta*}). First we read the data from
601 // this last file, which is explained below, and then we read the
602 // {\em Event-Header} block from the {\tt cer*} file.
603 // @endtext
604 // @code
605
606 // try to open the statistics file
607
608 if ( verbose >= VERBOSE_MAXIMAL )
609 log( SIGNATURE, "Openning %s\n", staname );
610
611 stat.openfile ( staname );
612
613 // reading data from this sta* file
614
615 if ( verbose >= VERBOSE_MAXIMAL )
616 log( SIGNATURE, "Reading %s\n", staname );
617
618 stat.read();
619
620 // read the Event Header from the cer* file
621
622 evth.read( cerfile );
623
624 // Distinquish between gamma/proton/helium/+heavier nuclei
625 if (evth.PrimaryID == gammas)
626 hprimary=0;
627 else if (evth.PrimaryID == protons)
628 hprimary=1;
629 else if (evth.PrimaryID == helium)
630 hprimary=2;
631 else
632 hprimary=3;
633
634 // set bin for theta/energy
635
636 kbin = ((int)floor(evth.get_theta()/10.)*6 +
637 (int)floor(log10(evth.get_energy())));
638
639 // get core distance and position
640
641 coreD = evth.get_core(&coreX, &coreY);
642
643 t1 = 100000;
644 t2 = 0.0;
645
646 ++ntotfiles;
647
648 // initialize time-histogram
649
650 for (i=0; i<60; ++i)
651 fstime[kbin][i] = 0.0;
652
653 // read each and every single photon (particle) from cer* file
654
655 rec2[0] = evth.Etotal;
656 rec2[1] = coreD;
657 rec2[2] = evth.Theta;
658 rec2[10] = evth.EvtNumber;
659
660 if (firstevent == -100)
661 {firstevent = evth.EvtNumber;
662 // save random number seed of first event per run
663 for (i=0;i<3;++i)
664 for (j=0;j<3;++j)
665 rnd[i][j]=evth.RndData[i][j];}
666
667 lastevent = evth.EvtNumber;
668
669 nCphotons = 0;
670 while ( ! cerfile.eof() ) {
671
672 // read photon data from the file
673
674 photon.read( cerfile );
675
676 wl = photon.get_wl();
677
678 HF2(300+hprimary,(photon.x-coreX)/100.,(photon.y-coreY)/100.,1.);
679 HF1(310+hprimary,sqrt(pow(((photon.x-coreX)/100.),2.)+pow((photon.y-coreY)/100.,2.)),1.);
680 HF1(320+hprimary,wl,1.);
681 HF1(330+hprimary,photon.h/100000.,1.);
682
683 rec2[3] = wl;
684 rec2[4] = photon.x;
685 rec2[5] = photon.y;
686 rec2[6] = photon.u;
687 rec2[7] = photon.v;
688 rec2[8] = photon.t;
689 rec2[9] = photon.h;
690// HFN(hid2,rec2);
691
692
693 if ( wl < 1.0 )
694 break;
695
696 // fill times
697
698 t = photon.get_t() - stat.get_tfirst();
699
700 if ( (t>0.0) && (t<60.0)) {
701 i = (int)t/1.;
702 ++fstime[kbin][i];
703 }
704
705 // hdmytime->Fill(t);
706
707 // increase number of Cphotons written
708
709 ++nCphotons;
710
711 } // while still there are photons left
712
713 // calculate errors for time distribution
714
715 for (i=0;i<60;++i)
716 fs2time[kbin][i] = sqrt(fstime[kbin][i]);
717
718 // save intermediate vectors in final vector
719
720 for (i=0;i<60;++i) {
721 ftstime[kbin][i] += fstime[kbin][i];
722 fts2time[kbin][i] += fs2time[kbin][i];
723 }
724
725 for (i=0;i<NTHSTEP;++i) {
726 ftselec[kbin][i] += (stat.plong[2][i] + stat.plong[1][i]);
727 fts2elec[kbin][i] += sqrt(stat.plong[2][i] + stat.plong[1][i]);
728 ftslight[kbin][i] += stat.plong[8][i];
729 fts2light[kbin][i] += sqrt(stat.plong[8][i]);
730 }
731
732 // estimate position of Xmax
733
734 vmax = -1.0;
735 for (i=0; i<NTHSTEP; ++i) {
736 v = (stat.plong[2][i] + stat.plong[1][i]); // e- & e+
737 if (v > vmax) {
738 vmax = v;
739 j = i;
740 }
741 }
742
743 // to calculate a rough estimate of Xmax, we make a weigthed mean
744 // of the values at maximum, and the previous and next bins
745 xmax=-100.;
746 if (j > 0) {
747 x0 = j*10. + 5.;
748 xa = (j-1)*10. + 5.;
749 xb = (j+1)*10. + 5.;
750 vmaxa = (stat.plong[2][j-1] + stat.plong[1][j-1]); // e- & e+
751 vmaxb = (stat.plong[2][j+1] + stat.plong[1][j+1]); // e- & e+
752
753 xmax = (xa*vmaxa + x0*vmax + xb*vmaxb) / (vmaxa + vmax + vmaxb);}
754 // else there was no e+e- shower in this event, skip ...
755
756 // fill N-Tuple record
757
758 if (evth.Etotal>0)
759 HF1(100+hprimary,log10(evth.Etotal),1.);
760 HF1(110+hprimary,coreD/100.,1.);
761 HF1(120+hprimary,evth.Theta/3.14*180.,1.);
762 HF1(130+hprimary,evth.Phi/3.14*180.,1.);
763 HF2(140+hprimary,coreX/100.,coreY/100.,1.);
764 HF1(150+hprimary,evth.zFirstInt/100000.,1.);
765 if (stat.get_tlast()-stat.get_tfirst()>0)
766 HF1(160+hprimary,stat.get_tlast()-stat.get_tfirst(),1.);
767 else
768 HF1(160+hprimary,-5.,1.);
769
770 HF1(170+hprimary,xmax,1.);
771 if (evth.Etotal<300){
772 if (nCphotons>0)
773 HF1(180+hprimary,log10(nCphotons),1.);
774 else
775 HF1(180+hprimary,-1.,1.);
776 }
777 else{
778 if (nCphotons>0)
779 HF1(190+hprimary,log10(nCphotons),1.);
780 else
781 HF1(190+hprimary,-1.,1.);
782 }
783 rec[0] = evth.PrimaryID;
784 rec[1] = evth.Etotal;
785 rec[2] = coreD;
786 rec[3] = evth.Theta;
787 rec[4] = evth.Phi;
788 rec[5] = thick(evth.zFirstInt, evth.Theta);
789 rec[6] = evth.zFirstInt;
790 rec[7] = evth.FirstTarget;
791 rec[8] = evth.p[0];
792 rec[9] = evth.p[1];
793 rec[10] = evth.p[2];
794 rec[11] = stat.get_tfirst();
795 rec[12] = stat.get_tlast();
796 rec[13] = xmax;
797 rec[14] = coreX;
798 rec[15] = coreY;
799 rec[16] = nCphotons;
800
801 rec[17] = evth.RunNumber;
802 rec[18] = evth.EvtNumber;
803 rec[19] = evth.DateRun;
804 rec[20] = evth.VersionPGM;
805 rec[21] = evth.HeightLev[0];
806 rec[22] = evth.SlopeSpec;
807 rec[23] = evth.ELowLim;
808 rec[24] = evth.EUppLim;
809
810
811 rec[25] = evth.Ecutoffh;
812 rec[26] = evth.Ecutoffm;
813 rec[27] = evth.Ecutoffe;
814 rec[28] = evth.Ecutoffg;
815
816 rec[29] = evth.EGS4yn;
817 rec[30] = evth.NKGyn;
818 rec[31] = evth.GHEISHAyn;
819 rec[32] = evth.VENUSyn;
820 rec[33] = evth.CERENKOVyn;
821 rec[34] = evth.NEUTRINOyn;
822 rec[35] = evth.HORIZONTyn;
823 rec[36] = evth.COMPUTER;
824 rec[37] = evth.ThetaMin;
825 rec[38] = evth.ThetaMax;
826 rec[39] = evth.PhiMin;
827 rec[40] = evth.PhiMax;
828
829 rec[41] = evth.StepLength;
830 rec[42] = evth.CWaveLower;
831 rec[43] = evth.CWaveUpper;
832
833 // write record to the file
834
835// HFN(hid,rec);
836
837 // ntuple->Fill( rec );
838
839 // close files
840
841 cerfile.close();
842 stat.closefile();
843
844 // show how many photons were written
845
846 if ( verbose >= VERBOSE_MAXIMAL )
847 log( SIGNATURE, "%d C-photons written.\n", nCphotons );
848
849
850 } // while there are still directory entries left
851
852 // check parameters for each RUN...
853 // Monitor run statistics:
854
855 ++rndiag;
856 sprintf(ctmp, "%s\n", "-------------------------------------------------------------");
857 for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
858
859 ++rndiag;
860 sprintf(ctmp,"%s%f%s%f%s%f\n", "Run: ",evth.RunNumber,
861 " Date: ",evth.DateRun, " CorsikaV: ",evth.VersionPGM);
862 for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
863
864 ++rndiag;
865 sprintf(ctmp, "%s%s\n", " Primary Particle: ",
866 primary_name[hprimary]);
867 for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
868
869 ++rndiag;
870 sprintf(ctmp, "%s%d%s%d\n", " First/Last Event Number monitored: ",
871 firstevent," / ",lastevent);
872 for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
873
874 ++rndiag;
875 sprintf(ctmp, "%s %f %f %f\n", " Random number seeds: ",
876 rnd[0][0],rnd[0][1],rnd[0][2]);
877 for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
878
879 ++rndiag;
880 sprintf(ctmp, "%s %f %f %f\n", " : ",
881 rnd[1][0],rnd[1][1],rnd[1][2]);
882 for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
883
884 ++rndiag;
885 sprintf(ctmp, "%s %f %f %f\n", " : ",
886 rnd[2][0],rnd[2][1],rnd[2][2]);
887 for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
888
889 // check ranges of i/p params
890 //
891 monitor_event(diag,&ndiag,primary_name[hprimary],
892 evth.HeightLev[0],mon_HeightLev,
893 ": Observation level is ");
894 monitor_event(diag,&ndiag,primary_name[hprimary],
895 evth.SlopeSpec,mon_SlopeSpec[hprimary],
896 ": Energy-spectrum slope is ");
897 monitor_event(diag,&ndiag,primary_name[hprimary],
898 evth.ELowLim,mon_ElowLim[hprimary],
899 ": Energy-low lim is ");
900 monitor_event(diag,&ndiag,primary_name[hprimary],
901 evth.EUppLim,mon_EUppLim[hprimary],
902 ": Energy-up lim is ");
903 monitor_event(diag,&ndiag,primary_name[hprimary],
904 evth.Ecutoffh,mon_Ecutoffh[hprimary],
905 ": Energy-cutoff hadrons is ");
906 monitor_event(diag,&ndiag,primary_name[hprimary],
907 evth.Ecutoffm,mon_Ecutoffm[hprimary],
908 ": Energy-cutoff muons is ");
909
910 monitor_event(diag,&ndiag,primary_name[hprimary],
911 evth.Ecutoffe,mon_Ecutoffe[hprimary],
912 ": Energy-cutoff electrons is ");
913
914 monitor_event(diag,&ndiag,primary_name[hprimary],
915 evth.Ecutoffg,mon_Ecutoffg[hprimary],
916 ": Energy-cutoff gammas is ");
917 monitor_event(diag,&ndiag,primary_name[hprimary],
918 evth.EGS4yn,mon_EGS4yn[hprimary],
919 ": EGS4yn flag is ");
920 monitor_event(diag,&ndiag,primary_name[hprimary],
921 evth.NKGyn,mon_NKGyn[hprimary],
922 ": NKGyn flag is ");
923 monitor_event(diag,&ndiag,primary_name[hprimary],
924 evth.GHEISHAyn,mon_GHEISHAyn[hprimary],
925 ": GHEISHAyn flag is ");
926 monitor_event(diag,&ndiag,primary_name[hprimary],
927 evth.VENUSyn,mon_VENUSyn[hprimary],
928 ": VENUSyn flag is ");
929 monitor_event(diag,&ndiag,primary_name[hprimary],
930 evth.CERENKOVyn,mon_CERENKOVyn[hprimary],
931 ": CERENKOVyn flag is ");
932 monitor_event(diag,&ndiag,primary_name[hprimary],
933 evth.ThetaMin,mon_ThetaMin[hprimary],
934 ": ThetaMin is ");
935 monitor_event(diag,&ndiag,primary_name[hprimary],
936 evth.ThetaMax,mon_ThetaMax[hprimary],
937 ": ThetaMax is ");
938 monitor_event(diag,&ndiag,primary_name[hprimary],
939 evth.PhiMin,mon_PhiMin[hprimary],
940 ": PhiMin is ");
941 monitor_event(diag,&ndiag,primary_name[hprimary],
942 evth.PhiMax,mon_PhiMax[hprimary],
943 ": PhiMax is ");
944 monitor_event(diag,&ndiag,primary_name[hprimary],
945 evth.CWaveLower,mon_CWaveLower[hprimary],
946 ": CWaveLower is ");
947 monitor_event(diag,&ndiag,primary_name[hprimary],
948 evth.CWaveUpper,mon_CWaveUpper[hprimary],
949 ": CWaveUpper is ");
950
951 } // for each data directory
952
953 // store contents of temporary 1D histograms to profiles
954
955 /*
956 Int_t bin;
957 Stat_t value;
958
959 for (kbin=0; kbin<18; ++kbin) {
960
961 for (i=0; i<60; ++i) {
962 ftstime[kbin][i] /= ntotfiles;
963 fts2time[kbin][i] = sqrt(fts2time[kbin][i]/ntotfiles -
964 SQR(ftstime[kbin][i]));
965 bin=i+1, value=ftstime[kbin][i];
966 cout << bin << ' ' << value << endl << flush;
967 htime->SetBinContent( bin, value );
968 htime->SetBinError( i+1, fts2time[kbin][i] );
969 }
970
971 for (i=0; i<NTHSTEP; ++i) {
972 ftselec[kbin][i] /= ntotfiles;
973 fts2elec[kbin][i] = sqrt(fts2elec[kbin][i]/ntotfiles -
974 SQR(ftselec[kbin][i]));
975 helec->SetBinContent( i+1, ftselec[kbin][i] );
976 helec->SetBinError( i+1, fts2elec[kbin][i] );
977 }
978
979 for (i=0; i<NTHSTEP; ++i) {
980 ftslight[kbin][i] /= ntotfiles;
981 fts2light[kbin][i] = sqrt(fts2light[kbin][i]/ntotfiles -
982 SQR(ftslight[kbin][i]));
983 hlight->SetBinContent( i+1, ftslight[kbin][i] );
984 hlight->SetBinError( i+1, fts2light[kbin][i] );
985 }
986
987 TH1F h1;
988 TH1F h2;
989 TH1F h3;
990
991 htime->Copy( h1 );
992 helec->Copy( h2 );
993 hlight->Copy( h3 );
994
995 hfile.Write();
996
997 }
998 */
999
1000 // close output file
1001
1002 if ( verbose >= VERBOSE_MINIMAL )
1003 log( SIGNATURE, "Closing N-Tuple file %s\n", outname );
1004
1005 // hfile.Write();
1006 // hfile.Close();
1007
1008 HROUT(0,icycle," ");
1009
1010 fprintf(stdout, "%s\n", "-------------------------------------------------------------");
1011 fprintf(stdout, "%s\n", "---------------------- RUN STATISTICS -----------------------");
1012
1013 for (i=0;i<=rndiag;++i)
1014 {char tmp[256]; int j;
1015 {for (j=0;j<256;++j) tmp[j]=rdiag[j][i];}
1016 fprintf(stdout,"%s",tmp);}
1017
1018 fprintf(stdout, "%s\n", "-------------------------------------------------------------");
1019 fprintf(stdout, "%s\n", "----------------------- DIAGNOSTICS -------------------------");
1020 fprintf(stdout, "%s\n", "-------------------------------------------------------------");
1021
1022 if (ndiag<0)
1023 {fprintf(stdout, "%s\n", " None, everything's splendidly perfect ...");
1024 fprintf(stdout, "%s\n", " Nada, todo fue bien, hacemos fiesta ...");
1025 fprintf(stdout, "%s\n", " Keine Fehler gefunden ...");}
1026 else
1027 {for (i=0;i<=ndiag;++i)
1028 {char tmp[256]; int j;
1029 {for (j=0;j<256;++j) tmp[j]=diag[j][i];}
1030 fprintf(stdout,"%s",tmp);}}
1031
1032
1033 HREND("CHECKMC");
1034
1035 // program finished
1036
1037 if ( verbose >= VERBOSE_MINIMAL )
1038 log( SIGNATURE, "Done.\n");
1039
1040}
1041// @endcode
1042
1043// @T \newpage
1044
1045// @subsection Functions definition
1046
1047//------------------------------------------------------------
1048// @name present
1049//
1050// @desc Make some presentation
1051//
1052// @date Sat Jun 27 05:58:56 MET DST 1998
1053// @function @code
1054//------------------------------------------------------------
1055// present
1056//
1057// make some presentation
1058//------------------------------------------------------------
1059
1060void
1061present(void)
1062{
1063cout << "##################################################\n"
1064<< SIGNATURE << '\n' << '\n'
1065<< "Check of the MC data generated by MMCS\n"
1066<< "J C Gonzalez, P Morawitz, D Petry, 27 Oct 1998\n"
1067<< "##################################################\n\n";
1068}
1069// @endcode
1070
1071//------------------------------------------------------------
1072// @name log
1073//
1074// @desc function to send log information
1075//
1076// @date Sat Jun 27 05:58:56 MET DST 1998
1077// @function @code
1078//------------------------------------------------------------
1079// log
1080//
1081// function to send log information
1082//------------------------------------------------------------
1083
1084void
1085log(const char *funct, char *fmt, ...)
1086{
1087 va_list args;
1088
1089 // Display the name of the function that called error
1090 printf("[%s]: ", funct);
1091
1092 // Display the remainder of the message
1093 va_start(args, fmt);
1094 vprintf(fmt, args);
1095 va_end(args);
1096}
1097// @endcode
1098
1099//------------------------------------------------------------
1100// @name error
1101//
1102// @desc function to send an error message, and abort the program
1103//
1104// @date Sat Jun 27 05:58:56 MET DST 1998
1105// @function @code
1106//------------------------------------------------------------
1107// error
1108//
1109// function to send an error message, and abort the program
1110//------------------------------------------------------------
1111
1112void
1113error(const char *funct, char *fmt, ...)
1114{
1115 va_list args;
1116
1117 // Display the name of the function that called error
1118 fprintf(stderr, "ERROR in %s: ", funct);
1119
1120 // Display the remainder of the message
1121 va_start(args, fmt);
1122 vfprintf(stderr, fmt, args);
1123 va_end(args);
1124
1125 fprintf(stderr, "\nTerminating program.\n");
1126 exit(999);
1127}
1128// @endcode
1129
1130
1131//------------------------------------------------------------
1132// @name readparam
1133//
1134// @desc function to read the parameters of the check
1135//
1136// @date Sat Jun 27 05:58:56 MET DST 1998
1137// @function @code
1138//------------------------------------------------------------
1139// readparam
1140//
1141// function to read the parameters of the check
1142//------------------------------------------------------------
1143
1144void
1145readparam(void)
1146{
1147 char sign[] = GLUE_postp( PROGRAM, VERSION ); // initialize sign
1148 char line[LINE_MAX_LENGTH]; // line to get from the stdin
1149 char token[ITEM_MAX_LENGTH]; // a single token
1150 int i, j; // dummy counters
1151
1152 // get signature
1153 cout << "\nExpecting input from stdin (see example input file) ...\n";
1154 cin.getline(sign, strlen(SIGNATURE) + 1);
1155 // cin.getline(line, LINE_MAX_LENGTH);
1156 if (strcmp(sign, SIGNATURE) != 0) {
1157 cerr << "ERROR: Signature of parameters file is not correct\n";
1158 cerr << '"' << sign << '"' << '\n';
1159 cerr << "should be: " << SIGNATURE << '\n';
1160 exit(1);
1161 }
1162
1163 // loop till the "end" directive is reached
1164 int is_end = FALSE;
1165 while (! is_end) {
1166
1167 // get line from stdin
1168 cin.getline(line, LINE_MAX_LENGTH);
1169
1170 // look for each item at the beginning of the line
1171 for (i=0; i<=end_file; i++)
1172 if (strstr(line, PAR_ITEM_NAMES[i]) == line)
1173 break;
1174
1175 // if it is not a valid line, just ignore it
1176 if (i == end_file+1) {
1177 cerr << " Skipping unknown token in [" << line << "]\n";
1178 continue;
1179 }
1180
1181 // case block for each directive
1182 switch ( i ) {
1183
1184 case data_paths: // beginning of the list of valid paths
1185
1186 //get number of paths to read
1187 sscanf(line, "%s %d", token, &Num_of_paths);
1188
1189 if ( verbose == TRUE )
1190 cout << " " << Num_of_paths << " path(s) will be used.\n";
1191
1192 // allocate memory for paths list
1193 Paths_list = (char**)malloc(Num_of_paths * sizeof(char*));
1194 for (i=0; i<Num_of_paths; i++)
1195 Paths_list[i] = (char*)malloc(PATH_MAX_LENGTH * sizeof(char));
1196
1197 // read each single path
1198 for (i=0; i<Num_of_paths; i++) {
1199
1200 cin.getline(Paths_list[i], PATH_MAX_LENGTH);
1201
1202 // remove empty spaces at the end
1203 for (j=0; j<strlen(Paths_list[i]); j++) {
1204 if (Paths_list[i][j] == ' ') {
1205 Paths_list[i][j] = '\0';
1206 break;
1207 }
1208 }
1209
1210 // show the path
1211 if ( verbose == TRUE )
1212 cout << " " << '[' << Paths_list[i] << ']' << '\n';
1213 }
1214
1215 break;
1216
1217 case output_file: // name of the output file
1218
1219 // get the name of the output_file from the line
1220 sscanf(line, "%s %s", token, Output_filename);
1221
1222 break;
1223
1224 case verbose_level: // verbose level 0-4
1225
1226 // get the verbose level
1227 sscanf(line, "%s %d", token, &verbose);
1228
1229 break;
1230
1231 case max_events: // maximum number of events
1232
1233 // get the number of events
1234 sscanf(line, "%s %d", token, &max_Events);
1235
1236 break;
1237
1238 case end_file: // end of the parameters file
1239
1240 // show a short message
1241 is_end = TRUE;
1242
1243 break;
1244
1245 } // switch ( i )
1246
1247 } // while (! is_end)
1248
1249 // after the loop is finished, return to the main function
1250 return;
1251}
1252// @endcode
1253
1254
1255//------------------------------------------------------------
1256// @name get_num_of_paths
1257//
1258// @desc get number of defined paths
1259//
1260// @date Sat Jun 27 05:58:56 MET DST 1998
1261// @function @code
1262//--------------------------------------------------
1263// get_num_of_paths
1264//
1265// get number of defined paths
1266//--------------------------------------------------
1267int
1268get_num_of_paths(void)
1269{
1270 return (Num_of_paths);
1271}
1272// @endcode
1273
1274
1275//------------------------------------------------------------
1276// @name get_path_name
1277//
1278// @desc get path name number "i" from the list of paths
1279//
1280// @date Sat Jun 27 05:58:56 MET DST 1998
1281// @function @code
1282//--------------------------------------------------
1283// get_path_name
1284//
1285// get path name number "i" from the list of paths
1286//--------------------------------------------------
1287char *
1288get_path_name(int i)
1289{
1290 return (Paths_list[i]);
1291}
1292// @endcode
1293
1294
1295//------------------------------------------------------------
1296// @name get_output_filename
1297//
1298// @desc get name of the output file
1299//
1300// @date Sat Jun 27 05:58:56 MET DST 1998
1301// @function @code
1302//--------------------------------------------------
1303// get_output_filename
1304//
1305// get name of the output file
1306//--------------------------------------------------
1307char *
1308get_output_filename(void)
1309{
1310 return (Output_filename);
1311}
1312// @endcode
1313
1314
1315
1316
1317//------------------------------------------------------------
1318// @name get_verbose
1319//
1320// @desc get the "verbose" status
1321//
1322// @date Sat Jun 27 05:58:56 MET DST 1998
1323// @function @code
1324//--------------------------------------------------
1325// get_verbose
1326//
1327// get the "verbose" status
1328//--------------------------------------------------
1329int
1330get_verbose(void)
1331{
1332 return (verbose);
1333}
1334// @endcode
1335
1336
1337//------------------------------------------------------------
1338// @name get_max_event
1339//
1340// @desc get max events number
1341//
1342// @date Sat Jun 27 05:58:56 MET DST 1998
1343// @function @code
1344//--------------------------------------------------
1345// get_max_event
1346//
1347// get max events number
1348//--------------------------------------------------
1349int
1350get_max_events(void)
1351{
1352 return( max_Events );
1353}
1354// @endcode
1355
1356
1357//------------------------------------------------------------
1358// @name thick
1359//
1360// @desc converts height [cm] -> thickness [g cm-2]
1361//
1362// @date Sat Jun 27 05:58:56 MET DST 1998
1363// @function @code
1364//--------------------------------------------------
1365// thick
1366//
1367// converts height [cm] -> thickness [g cm-2]
1368//--------------------------------------------------
1369
1370float
1371thick( float height, float theta )
1372{
1373 float h, thk;
1374
1375 // Ibarra's modification
1376 h = sqrt( SQR(REarth) +
1377 SQR(height/cos(theta)) +
1378 (2.0*REarth*height)) - REarth;
1379
1380 if ( h < 4.e5 ) {
1381 thk = AATM[0] + BATM[0] * exp ( -h * DATM[0] );
1382 } else if ( h < 1.e6 ) {
1383 thk = AATM[1] + BATM[1] * exp ( -h * DATM[1] );
1384 } else if ( h < 4.e6 ) {
1385 thk = AATM[2] + BATM[2] * exp ( -h * DATM[2] );
1386 } else if ( h < 1.e7 ) {
1387 thk = AATM[3] + BATM[3] * exp ( -h * DATM[3] );
1388 } else {
1389 thk = AATM[4] - h * CATM[4];
1390 }
1391
1392 return ( thk );
1393
1394}
1395// @endcode
1396
1397//------------------------------------------------------------
1398// @name height
1399//
1400// @desc converts thickness [g cm-2] -> height [cm]
1401//
1402// @date Sat Jun 27 05:58:56 MET DST 1998
1403// @function @code
1404//--------------------------------------------------
1405// height
1406//
1407// converts thickness [g cm-2] -> height [cm]
1408//--------------------------------------------------
1409
1410float
1411height( float thk, float theta )
1412{
1413 float h, height;
1414
1415 if ( thk > 631.1e0 ) {
1416 h = CATM[0] * log ( BATM[0] / (thk - AATM[0]) );
1417 } else if ( thk > 271.7e0 ) {
1418 h = CATM[1] * log ( BATM[1] / (thk - AATM[1]) );
1419 } else if ( thk > 3.0395e0 ) {
1420 h = CATM[2] * log ( BATM[2] / (thk - AATM[2]) );
1421 } else if ( thk > 0.00128292e0 ) {
1422 h = CATM[3] * log ( BATM[3] / (thk - AATM[3]) );
1423 } else {
1424 h = (AATM[4] - thk) * DATM[4];
1425 }
1426
1427 // Ibarra's modification
1428 height = SQR(cos(theta)) *
1429 (-REarth + sqrt(SQR(REarth) +
1430 (SQR(h + (2.0*REarth*h))/SQR(cos(theta)))));
1431
1432 return ( height );
1433
1434}
1435// @endcode
1436
1437//------------------------------------------------------------
1438// @name monitor_event
1439//
1440// @desc compare parameter with reference and report if not equal
1441//
1442// @date Sat Jun 27 05:58:56 MET DST 1998
1443// @function @code
1444//--------------------------------------------------
1445// monitor_event
1446//
1447// compare parameter with reference
1448//--------------------------------------------------
1449void monitor_event(char diag[DIAGX][DIAGY], int *ndiag, const char* primary,
1450 const float a, const float b, const char* message)
1451{char diagt[DIAGX]; int i;
1452 if (a!=b){
1453 if (*ndiag<(DIAGY-1))
1454 {++*ndiag;
1455 sprintf(diagt,"%s%s%s%s%f%s%f\n","Primary ",primary," ",
1456 message,a," but should be ",b);}
1457 else{
1458 sprintf(diagt,"%s",
1459 "Two many errors encountered, won't print any further diagnostics\n");
1460 for (i=0;i<DIAGX;++i){
1461 diag[i][*ndiag]=diagt[i];
1462 }
1463 } // endif
1464 } // endif
1465 return;
1466}
Note: See TracBrowser for help on using the repository browser.