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
|
---|
29 | In this section we show the (commented) code of the program
|
---|
30 | for 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
|
---|
36 | All 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
|
---|
46 | Here we define the global variables where the values from
|
---|
47 | the parameters file are stored. See the section
|
---|
48 | \ref{paramdesc} in page \pageref{paramdesc} for a more
|
---|
49 | complete description of the available commands in the
|
---|
50 | parameters file.
|
---|
51 | @endtext */
|
---|
52 |
|
---|
53 | // @code
|
---|
54 | static char **Paths_list; // list of paths to get data from
|
---|
55 | static int Num_of_paths = 1; // number of paths
|
---|
56 | static char Output_filename[PATH_MAX_LENGTH]; // output filename
|
---|
57 | static char CT_filename[PATH_MAX_LENGTH]; // name of the CT def. file
|
---|
58 | static VerboseLevel verbose = VERBOSE_DEFAULT; // verboseness
|
---|
59 | static float fixed_Theta; // fixed target theta and phi
|
---|
60 | static float fixed_Phi;
|
---|
61 | static float low_Ecut = 0.0, high_Ecut = 100000.0;
|
---|
62 | static int is_Fixed_Target = FALSE; // are we using fixed target?
|
---|
63 | static int max_Events = 100000; // maximum number of events to read
|
---|
64 |
|
---|
65 | static float gammas = 1.;
|
---|
66 | static float protons = 14.;
|
---|
67 | static float helium = 402.;
|
---|
68 |
|
---|
69 |
|
---|
70 | inline int max(int a, int b)
|
---|
71 | {if (a>b) return a;
|
---|
72 | else return b;}
|
---|
73 |
|
---|
74 | inline int min(int a, int b)
|
---|
75 | {if (a<b) return a;
|
---|
76 | else return b;}
|
---|
77 |
|
---|
78 | // @endcode
|
---|
79 |
|
---|
80 | /* @text
|
---|
81 | Earth's radius (in cm).
|
---|
82 | @endtext */
|
---|
83 |
|
---|
84 | // @code
|
---|
85 | static const float REarth = 6348.0e5;
|
---|
86 | // @endcode
|
---|
87 |
|
---|
88 | // @subsection Main program
|
---|
89 |
|
---|
90 | /* @text
|
---|
91 | Let's start!
|
---|
92 | @endtext */
|
---|
93 |
|
---|
94 | // @code
|
---|
95 | //++++++++++++++++++++++++++++++++++++++++
|
---|
96 | // MAIN PROGRAM
|
---|
97 | //----------------------------------------
|
---|
98 |
|
---|
99 | void
|
---|
100 | main(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 |
|
---|
1060 | void
|
---|
1061 | present(void)
|
---|
1062 | {
|
---|
1063 | cout << "##################################################\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 |
|
---|
1084 | void
|
---|
1085 | log(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 |
|
---|
1112 | void
|
---|
1113 | error(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 |
|
---|
1144 | void
|
---|
1145 | readparam(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 | //--------------------------------------------------
|
---|
1267 | int
|
---|
1268 | get_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 | //--------------------------------------------------
|
---|
1287 | char *
|
---|
1288 | get_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 | //--------------------------------------------------
|
---|
1307 | char *
|
---|
1308 | get_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 | //--------------------------------------------------
|
---|
1329 | int
|
---|
1330 | get_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 | //--------------------------------------------------
|
---|
1349 | int
|
---|
1350 | get_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 |
|
---|
1370 | float
|
---|
1371 | thick( 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 |
|
---|
1410 | float
|
---|
1411 | height( 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 | //--------------------------------------------------
|
---|
1449 | void 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 | }
|
---|