/* ============================================================== */ /* INCLUDES AND DECLARATIONS */ /* ============================================================== */ #include #include #include "partap.h" #include "cfortran.h" #include "mcEvKey.h" #include "mcEvent.h" #include "mcBeam.h" #include "mcTrack.h" #include "mcVert.h" #include "mcHit.h" #include "dataPuls.h" #include "dataHodo.h" #include "dataCalo.h" int i, i1, i2, j1, j2, igaf, idfl, isel, iok, ievent; int ind, indname, indname1, indname2; int inddigi, inddigi2; int decay_vertex, idet; int compare; float values[4]; int ok; /* will use as logical; T=1, F=0 */ char colnam[2][17]; char corder[2][5]; char cline[81]; char cname[33]; /* dummy string variable... */ char digitable[17]; /* another dummy string variable */ void printwin_mctrack(void); void printwin_mchit(void); void printwin_datahodo(void); /* =============================================================== */ /* START PROGRAM */ /* =============================================================== */ void main() { /* =============================================================== */ /* INITIALIZE ADAMO */ /* =============================================================== */ INITAP(); /* =============================================================== */ /* OPEN AN EXISTING GAF FILE */ /* =============================================================== */ /* The routine OPEGAF opens an existing GAF. As you can see, the syntax of this routine is very similar to that of CREGAF. */ sprintf(cline,"NAME=toyMC.events,DRIVER=FZ,FILFOR=EXCH"); OPEGAF(igaf,cline,iok); if (iok != 0) { printf("\n OPEGAF failed \n"); exit(0); } /* Let's demonstrate the command LISTAB here, which simply lists all tables defined in a given dataflow (here, our uber-dataflow mcEvents). */ printf("\n LISTAB of mcEvents \n"); idfl = GETTDF("mcEvents"); LISTAB(idfl); /* Recall that using LISTAB( GETTDF("mcEvents") ) does NOT work for some odd reason... */ /* ================================================================ */ /* READ THE FIRST RECORD FROM THE GAF */ /* ================================================================ */ /* Using NEXGAF and FETGAF ----------------------- NEXGAF and FETGAF are the two routines you need to read in each record of the GAF. */ printf("\n\n ========== NEXGAF AND FETGAF ========== \n\n"); /* As you can see in the IE version of your GAF, the key table information appears at the top of each record, in a simple row. NEXGAF retrieves this information, and loads the key table information to the key table structure */ NEXGAF(igaf, mcEvKey, iok); /* We'll check the error code variable iok now. This is more than just an error check ... it will also tell us if we have reached the end of the GAF. The end-of-file code is iok = -1. */ if (iok != 0) { CLOGAF(igaf, iok); exit(0); } /* Now we use some of the key table information from this record. We store the event number in a convenient variable ievent, and inform the user. Also, remember that the last entry in any key table MUST be a CH32 variable, and is used to hold the name of the DATAFLOW which is stored in the record. We already know that the dataflow name is mcEvents, since we WROTE this GAF. When you are reading in someone else's GAF, however, it is important to know that the dataflow name is available here. */ /* Note the use of another bizarre cfortran command, FCB2CSTR (similar to C2FCBSTR, which was used in write-gaf.c). FCB2CSTR retrieves a string from a Fortran common block. The string here is mcEvKey.cName which like all other WCB variables is internally stored by ADAMO in a Fortran common block. We now use FCB2CSTR to dump it into a C string (cname) while getting rid of the trailing blanks. */ ievent = mcEvKey.iEvent; printf("\n <<<<< EVENT NUMBER %d >>>>>", ievent); FCB2CSTR(mcEvKey.cName, cname, 0); printf("\n mcEvKey.cName = %s ", cname); printf("(should be the dataflow name, mcEvents)"); idfl = GETTDF(cname); /* Now that we know the dataflow identifier, we can call FETGAF to actually read the tables from this record into memory. Note that FETGAF can read in any record you choose. It determines which record you want from the information in the key table WIN. Since we just loaded this WIN with NEXGAF, we will be getting the NEXT record in the GAF. */ FETGAF(igaf, mcEvKey, idfl, iok); /* ================================================================= */ /* INTRODUCING INDEXES */ /* ================================================================= */ printf("\n\n ========== USING INDEXES ========== \n"); /* An index is a logical ordering of the rows in a table. You ask ADAMO to create these indexes for you, given the sorting criteria you want. Just like tables, GAFs, and dataflows, indexes are stored with a name and an identifier number. */ /* The simplest indexes are based on the values in a single table column. The simplest index of all is based on the ID column ... this is the natural sorting order of the table, by row number. This index is used so often in ADAMO that partap.inc provides you with a 'secret code' for it, the parameter ID. */ /* We saw the use of PRITAB in the write-gaf example code. However, we didn't go into much detail. Here once more is the standard way to print the entire mcTrack table: */ printf("\n\n mcTrack sorted by ID \n"); PRITAB(mcTrack, ID, MINC, MAXC, ALLCOL); /* Now look more closely at the second argument. There is the secret code ID that I just mentioned! The second argument is actually looking for an index identifier. We have used the ID index, and so the table is printed out sorted by row number. */ /* Now let's create a simple index, sorted on the column 'cName'. GETIND can create these simple indexes on any (single) column name. */ indname = GETIND(mcTrack, "cName"); printf("\n mcTrack sorted by index on cName \n"); PRITAB(mcTrack, indname, MINC, MAXC, ALLCOL); /* This time, the rows of the table are printed out in a different order: they are arranged so that the particle names are in alphabetical order. */ /* Finally, we introduce CREIND, which can create more complex indexes sorted on more than one column. You supply the column names in a vector (here colnam), by order of precedence for your sort. You also supply a vector of the same length containing either 'ASC' (ascending) or 'DESC' (descending), to specify the order in which each column should be sorted. The 3rd last argument to CREIND specifies the number of entries in your vectors. */ strcpy(colnam[0], "cName"); strcpy(colnam[1], "iCharge"); strcpy(corder[0], "ASC"); strcpy(corder[1], "DESC"); CREIND(mcTrack, indname2, "NameOrder2", 2, colnam, corder); printf("\n mcTrack sorted by ASC order on cName, "); printf("then DESC on iCharge \n"); PRITAB(mcTrack, indname2, MINC, MAXC, ALLCOL); /* ================================================================= */ /* RETRIEVING INFORMATION FROM TABLES */ /* ================================================================= */ /* We have read in the first record of the GAF, so the information is all stored in ADAMO memory. But to work with it, we need to retrieve this information row by row into our WCB variables. There are three basic routines that do this for us: FETTAB, GETTAB, and FETCOL. */ /* Introducing FETTAB ------------------ Very simple. FETTAB takes three arguments: the table identifier, an index identifier, and a 'cursor' number. This cursor specifies: which row do you want from the index you suppplied? In this first example, we use the row-sorted index ID, and request that the first row be pulled into the structure. (The subroutine printwin_mctrack is just a little utility routine whose code is at the end of this file. It prints out some of the variables from the structure to the screen, in a more compressed format than PRITAB.) */ printf("\n\n ========== USING FETTAB ========== \n"); printf("\n Retrieving mcTrack #1, from index ID"); FETTAB(mcTrack, ID, 1); printwin_mctrack(); /* Let's try that again using a different index: the index on the cName column which we defined earlier. Cursor number 1 means the first entry in this index, which is NOT row #1. */ printf("\n Retrieving mcTrack #1, from index on cName"); FETTAB(mcTrack, indname, 1); printwin_mctrack(); /* Introducing GETTAB ------------------ GETTAB retrieves one row from a table in a simpler way. It does NOT allow you to use indices. Instead, you set the row number that you want in the ID variable of the structure, and simply call GETTAB with the table name. GETTAB retrieves the row you requested. */ printf("\n\n ========== USING GETTAB ========== \n"); printf("\n Retrieving mcTrack row #2"); mcTrack.ID = 2; GETTAB(mcTrack); printwin_mctrack(); /* The principal virtue of the simpler GETTAB is that it is FASTER than FETTAB. Keep this in mind if you are writing code that does lots of table accesses and is running slowly. Of course you can only use it if you have no need of the indexing capabilities supplied by FETTAB. */ /* Introducing FETCOL ------------------ FETCOL retrieves table information by COLUMN rather than by row. Since the structure is row-oriented, the information is NOT placed there, but rather is stored in a vector that you supply. Here, we retrieve the momentum ('P') column from all rows of table mcTrack and store the values in the real vector 'values'. This vector was declared at the top of the program with 4 entries, since we already know there can be at most 4 tracks in the mcTrack table. In general: be sure your vector is LONG ENOUGH ... */ printf("\n\n ========== USING FETCOL ========== \n"); printf("\n Retrieving momenta from mcTrack, sorted by ID"); FETCOL(mcTrack, "P", ID, MINC, MAXC, values); for (i = 0; i <=3; i++) { printf("\n Momentum #%d", i+1); printf(" = %f", values[i]); } /* Arguments 3, 4, and 5 of FETCOL are the same as those from PRITAB: you can supply an index, and a first and last cursor on this index. In this example, we retrieve the momenta using an index on the momentum column itself, and retrieve only the last two rows. This should give us the highest two momenta in the table. Note that although we are retrieving the 3rd and 4th entries of the index, the values go to the first two elements of the array 'values'. To illustrate this, we will first clear the values array. */ for (i = 0; i <= 3; i++) { values[i] = 0.; } printf("\n Retrieving highest two momenta from mcTrack"); ind = GETIND(mcTrack, "P"); FETCOL(mcTrack, "P", ind, 3, 4, values); for (i = 0; i<=1; i++) { printf("\n Momentum #%d", i+1); printf(" = %f", values[i]); } /* Note that FETCOL is also faster than FETTAB. */ /* ================================================================== */ /* NAVIGATING RELATIONSHIPS */ /* ================================================================== */ /* Example 1: navigating mcHit -> mcTrack -------------------------------------- */ printf("\n\n ========== NAVIGATING RELATIONSHIPS ========== \n"); /* Remember all those nice relationships that we set up in the DDL, relating rows from one table to rows from another table? We're now going to learn how to NAVIGATE (i.e. follow) these relationships, using NATREL and NAFREL. This concept is sufficiently important (and sufficiently tricky to wrap your brain around ...) that I will supply two examples. In the first example, we will navigate the relationship mcHit -> mcTrack. Each hit (i.e. each row in the mcHit table) is associated with the particle track that caused it. If you recall, we set up our GAF so that only track #1 (the scattered positron) caused any hits, so we already know what our answers will be: all hits are linked to track #1. */ /* NATREL navigates in the TO direction, which means FORWARDS along the link. The link points from mcHit to mcTrack. So we specify a row in mcHit, and then ask NATREL to tell us which row in mcTrack it is linked TO. */ printf("\n Navigate mcHit #1 FORWARDS to mcTrack using NATREL"); mcHit.ID = 1; NATREL(mcHit, mcHit.mcTrack, mcTrack, ok); printwin_mctrack(); /* As you will see in the output, NATREL retrieved row #1 of mcTrack for us and stored it to the structure. A couple of notes: - We did not have to load up the mcHit structure with the full information from mcHit row #1. NATREL only looks at the ID column (it doesn't need to know anything else, after all, just a row number). - Please note that unlike our previous error codes, the error code here is 'ok', which is not an integer variable but a LOGICAL variable. It is perfectly possible that mcHit row #1 is not linked to ANYTHING in mcTrack (the link could be NULL). If ADAMO is unable to navigate the link, the 'ok' variable will simply return FALSE. This is not an error, really, just information you might need. */ /* OK, that was straightforward. Now we will learn about the second routine, NAFREL. This one you DO have to wrap your brain around, a bit. This routine navigates BACKWARDS along a link. Our link is mcHit -> mcTrack. NAFREL will start with a row from mcTrack, and search BACKWARDS along the link to find all rows from mcHit that are linked to it. Unlike forward navigation, this routine may well come up with more than one linked row. */ /* In the following example, we start with mcTrack row 1 (the scattered positron) and find out which hits it caused. We know already that this should give us all 5 rows of the mcHit table. */ printf("\n Navigate mcTrack #1 BACKWARDS to mcHit"); printf(" using NAFREL"); NULWIN(mcTrack); mcTrack.ID = 1; ind = GETIND(mcHit, "mcTrack"); NAFREL(mcTrack, mcHit.mcTrack, mcHit, ind, i1, i2); for (i = i1; i <= i2; i++) { FETTAB(mcHit, ind, i); printwin_mchit(); } /* Yup, you will see from the output that all 5 rows of mcHit were found. Some more notes: - The starting point this time was row 1 from mcTrack. Again, ONLY the ID column mcTrack_ID has to be set, with this information. I proved it to you this time by first clearing the mcTrack structure using a call to NULWIN. - The output of NAFREL is contained in the cursor variables i1 and i2. These specify the range of rows on mcHit that were found by the navigation. But: range of rows on what index? You have to give NAFREL exactly the right index as its fourth argument: an index sorted on the RELATIONSHIP COLUMN itself. If you think about this for a minute, it makes sense. There is no reason that the rows NAFREL finds will be contiguous within any other index ... - NAFREL found 5 rows for us this time. It did not load up any of these rows to the mcHit structure, because it does not know which one you want. We had to load them up one at a time ourselves using FETTAB. If NAFREL had found only ONE row, however, then the cursors i1 and i2 would have come back with the same values, AND the mcHit structure would have been automatically loaded up with this one row. How convenient! You will see this behaviour in operation later on, with some other (similar) commands that return cursors on an index. - Finally, what if NAFREL had found NO matching rows in mcHit? In this case, it would have returned cursor values i1 and i2 such that i1 > i2. This is how NAFREL (and all other commands that return cursors on an index) indicate that no match was found. It is a nice way of returning this information, since a loop over the range i1 --> i2 will simply not execute. But sometimes you may have reason to check explicitly for this failure condition. */ /* Example 2: navigating mcTrack -> mcTrack BY Parent -------------------------------------------------- This second example is really not any more complicated that the first one. This time we will navigate forwards and backwards along the Parent link, from one row of mcTrack to another. We will explore this relationship using track #2, the unstable K0_s particle, and tracks #3-4, the charged pions to which the K0_s decays. The K0_s is the Parent of both pions. */ /* In the forward direction, we start with one of the child particles, the pi+ (row#4) and navigate FORWARDS to its parent. */ printf("\n Navigate mcTrack #4 FORWARDS to mcTrack"); printf(" by Parent, using NATREL"); mcTrack.ID = 4; NATREL(mcTrack, mcTrack.Parent, mcTrack, ok); printwin_mctrack(); /* Let me just show you something before we continue, maybe this will help to clarify what is going on. We did not actually need to call NATREL to do this for us. We could have done it 'by hand' as follows: */ printf("\n Navigate mcTrack #4 FORWARDS to mcTrack"); printf(" by Parent, BY HAND"); mcTrack.ID = 4; GETTAB(mcTrack); mcTrack.ID = mcTrack.Parent; GETTAB(mcTrack); printwin_mctrack(); /* Makes sense? The point is this: whenever you have an EXACTLY one-to-one relationship from table1 --> table2 (they could be the same table), and you only want to navigate forwards along it (i.e. from table1 to table2), you actually don't need NATREL or NAFREL or anything else. The link variable in table1 simply gives you the corresponding row number in table2, so you can just load it up yourself. You'll often see such examples in existing codes. But note: this simple technique fails of course if the relationship is NOT one-to-one. */ /* One last forwards example illustrates how the status variable 'ok' can sometimes come back false. Let's start with the K0_s this time, and try to navigate to ITS parent particle. The K0_s does not HAVE any parent ... the link variable mcTrack_Parent is equal to INULL for the K0_s track. */ printf("\n NATREL failure example: finding Parent of K0_s"); mcTrack.ID = 2; NATREL(mcTrack, mcTrack.Parent, mcTrack, ok); printf("\n NATREL came back with ok = %d", ok); /* Finally, we go BACKWARDS. We start with the K0_s track, and navigate the Parent link backwards to its CHILD tracks. Hopefully, we will recover the two pions. */ printf("\n Navigate mcTrack #2 BACKWARDS to mcTrack"); printf(" by Parent, using NAFREL"); mcTrack.ID = 2; ind = GETIND(mcTrack, "Parent"); NAFREL(mcTrack, mcTrack.Parent, mcTrack, ind, i1, i2); for (i = i1; i <= i2; i++) { FETTAB(mcTrack, ind, i); printwin_mctrack(); } /* As before, NAFREL needs an index on the relationship column in question, and returns two cursors on that index. */ /* ================================================================== */ /* USING SELTAB */ /* ================================================================== */ printf("\n\n ========== USING SELTAB ========== \n"); /* SELTAB is a very useful function that adds a selection capability to the index mechanism. The concept is this: an index on a particular column X sorts a table for you based on the values of X. Obviously, any rows that have the SAME value of X will end up in a contiguous block in the index. SELTAB now allows you to SELECT the row or rows from the table that have some PARTICULAR value of X. The answer comes back as a cursor (range of rows) on the index. */ /* Remember our generalized relationship? mcHit -> dataPuls | dataHodo | dataCalo BY digiTable This links hits to digitizationss, where the digitizations can appear in one of three different tables. Two structure variables characterize the generalized relationship: mcHit.digiTable contains the NAME of the link table (dataPuls, ...) mcHit.digiTable_ contains the ROW NUMBER within the link table In this example, we are going to find all the rows of mcHit which are linked to digitizations in the dataHodo table. */ /* What we want to do, then, is SELECT all mcHit rows which have mcHit.digiTable = "dataHodo". SELTAB can do this easily. */ printf("\n Find all hodoscope hits using SELTAB"); inddigi = GETIND(mcHit, "digiTable"); NULWIN(mcHit); C2FCBSTR("dataHodo", mcHit.digiTable, 0); SELTAB(mcHit, inddigi, i1, i2); for (i = i1; i <= i2; i++) { FETTAB(mcHit, inddigi, i); printwin_mchit(); } /* That was straightforward. As I mentioned, you have to use an index on the selection column ... otherwise the selected rows won't necessarily appear in a contiguous block which can be identified by the two cursors. Also note: the only thing SELTAB needs to know is the value you want to select, in this case mcHit.digiTable = "dataHodo". I used NULWIN again simply to demonstrate that no other variables in the structure need to be filled. */ /* Just to show you what we did, let's do it again, but BY HAND this time: we can just cycle through the table rows, picking out the ones we want. How inelegant ... but it works too. */ printf("\n Find all hodoscope hits by BRUTE FORCE SEARCH"); for (i = 1; i <= COUTAB(mcHit); i++) { mcHit.ID = i; GETTAB(mcHit); FCB2CSTR(mcHit.digiTable, digitable, 0); compare = strcmp(digitable, "dataHodo"); if (compare == 0) printwin_mchit(); } /* Finally, a failure example. What happens if SELTAB cannot find ANY rows matching your selection criteria? Let's try it: we'll ask SELTAB to find all rows with mcHit_digiTable = 'toad' :) */ printf("\n SELTAB failure example: finding hits in toad"); C2FCBSTR("toad", mcHit.digiTable, 0); SELTAB(mcHit, inddigi, j1, j2); printf("\n SELTAB came back with"); printf(" cursor 1 = %d, cursor 2 = %d", j1, j2); /* As you will see from the output, SELTAB indicates its failure by coming back with j1 (low end cursor) > j2 (high end cursor). */ /* ================================================================== */ /* NAVIGATING GENERALIZED RELATIONSHIPS */ /* ================================================================== */ printf("\n\n ========== NAVIGATING GENERALIZED"); printf(" RELATIONSHIPS ========== \n"); /* As our last lesson in navigation, we need to learn about NATGEN and NAFGEN. These do the same as NATREL and NAFREL, but for generalized relationships which can point to rows in different tables. */ /* What we will do is navigate from mcHit to dataHodo. The syntax of NATGEN is identical to that of NATREL: it needs the identifiers of the two tables, and the relationship variable along which you want to navigate. This variable is mcHit.digiTable. (Note: it is NOT mcHit.digiTable_). */ /* At this point in the program, we already know which rows in mcHit DO in fact link to dataHodo. The cursor variables i1 and i2 currently hold the row range, on the index inddigi. So let's loop over those rows, and retrieve the corresponding rows from dataHodo. */ printf("\n Navigate mcHit #1 FORWARDS to dataHodo"); printf(" BY digiTable, using NATGEN"); for (i = i1; i <= i2; i++) { FETTAB(mcHit, inddigi, i); NATGEN(mcHit, mcHit.digiTable, dataHodo, ok); printwin_datahodo(); } /* NOTE: As usual, all NATGEN needs to know is the ID number on mcHit from which you want to start your navigation. However, in my do loop, I did NOT just set mcHit_ID = i before calling NATGEN. i1 and i2 are NOT cursors on the ID index but on a different index --> they are NOT row numbers. I used FETTAB and the proper index inddigi to pull each mcHit row (including its ID value) into the WIN before calling NATGEN. */ /* Now a failure example. What happens if I navigate from the CALORIMETER hit to dataHodo? The calorimeter hit (row #5) is linked to dataCalo. This will surely fail ... First, let's find the CALO hit. We already know it is row #5, but's let's use SELTAB, for practice. */ printf("\n NATREL failure example: navigating a CALO hit"); printf(" to dataHodo, using NATGEN"); C2FCBSTR("dataCalo", mcHit.digiTable, 0); SELTAB(mcHit, inddigi, i1, i2); printwin_mchit(); /* I didn't check, but you can be sure SELTAB found exactly one row, number 5. Cursors i1 and i2 came back equal to each other. As I mentioned before, under this circumstance, the ADAMO routine will conviently load this unique row into the structure for you. I proved that by printing out the mcHit structure. */ /* Finally, we navigate from this calorimeter hit to dataHodo. You will see in the output that ok comes back FALSE. */ NATGEN(mcHit, mcHit.digiTable, dataHodo, ok); printf("\n NATGEN came back with ok = %d", ok); /* Now, navigating generalized relationships BACKWARDS. NAFGEN is again very similar to NAFREL. Recall that when using NAFREL, you had to first set up and index on the relationship column. Same here, EXCEPT, you need a more complex index, on BOTH columns of the generalized relationship. We already learned how to do this, using CREIND. */ strcpy(colnam[0], "digiTable"); strcpy(colnam[1], "digiTable_"); strcpy(corder[0], "ASC"); strcpy(corder[1], "ASC"); CREIND(mcHit, inddigi2, "DigiOrder2", 2, colnam, corder); /* The rest is just like NAFREL. In this example, we will navigate backwards from row #1 of dataHodo to the corresponding row of mcHit. */ printf("\n Navigate dataHodo #1 BACKWARDS to mcHit"); printf(" BY digiTable, using NAFGEN"); dataHodo.ID = 1; NAFGEN(dataHodo, mcHit.digiTable, mcHit, inddigi2, i1, i2); printwin_mchit(); /* Once more, I know that only one row has been found, so the mcHit structure has been loaded up automatically. */ /* ================================================================== */ /* LOOP OVER REMAINING EVENTS */ /* ================================================================== */ /* We have already demonstrated everything we need to using the first record in the GAF, so there's really no need to read in the rest of them. But let's do it anyway. Remember, we do not usually know in advance how many records are in the GAF, so you have to watch the error code out of NEXGAF for the end-of-file signal. */ iok = 0; NEXGAF(igaf, mcEvKey, iok); while (iok == 0) { printf("\n <<<<< EVENT NUMBER"); printf(" %d >>>>>", mcEvKey.iEvent); FETGAF(igaf, mcEvKey, idfl, iok); NEXGAF(igaf, mcEvKey, iok); } /* ================================================================== */ /* CLOSE THE GAF */ /* ================================================================== */ CLOGAF(igaf, iok); printf("\n"); /* ================================================================== */ /* END PROGRAM */ /* ================================================================== */ } /* ================================================================== */ /* CODING EXAMPLES: PRINTING THE WINs of mcTrack, mcHit, and dataHodo */ /* ================================================================== */ void printwin_mctrack(void) { FCB2CSTR(mcTrack.cName, cname, 0); printf("\n mcTrack: ID = %d", mcTrack.ID); printf(" cName = %s", cname); printf(" iCharge = %d", mcTrack.iCharge); printf(" P = %f", mcTrack.P); } /* ------------------------------------------------------------------ */ void printwin_mchit(void) { FCB2CSTR(mcHit.digiTable, digitable, 0); printf("\n mcHit: ID = %d", mcHit.ID); printf(" digiTable = %s, ", digitable); printf("row %d", mcHit.digiTable_); printf(" z = %f", mcHit.Z); } /* ------------------------------------------------------------------ */ void printwin_datahodo(void) { FCB2CSTR(dataHodo.cName, cname, 0); printf("\n dataHodo: ID = %d", dataHodo.ID); printf(" cName = %s", cname); printf(" iWire = %d", dataHodo.iWire); printf(" iTDC = %d", dataHodo.iTDC); }