///////////////////////////////////////////////////////////////////////// // Program to convert an gpamela ntuple to pamela raw data (.pam). // A root file containing the ntuple is also produced. // The conversion is done using a slightly patched version of // Rene Brun's h2root. // // Created : 2006 Jan Conrad (conrad@particle.kth.se) // Modified: 2007 Silvio Orsi (silvio.orsi@roma2.infn.it) // // // HOWTO (RECOMMENDED): // ./ Pamelagp2Digits filename.his // It generates filename.pam (PAMELA raw data) and filename.gp.root (sim info) // // NB: more details are present in the file HOW-TO-DIGIT.TXT // // HOWTO (only for testing purposes) // Pamelagp2Digits hbook_filename root_filename raw_filename compress tolower // - If the 2nd and 3rd parameters are missing the names are generated from the hbook file name. // - The following flags are from h2root and default values are used: // if compress is missing (or = 1) the ROOT file will be compressed // if compress = 0 the ROOT file will not be compressed. // if tolower is missing (or = 1) ntuple column names are converted to lower case // but the first character is converted to upper case. // if tolower = 2 same as tolower=1 except that the first character is also // convertex to lower case // // MODIFICATIONS: // // v. 1.6 (S.Orsi, Roma2, 31 October 2007) // - TOF: (Campana, Menn, Orsi) // + bugs fixed // + crc included // + introduction of 3 control bits in the ADC and TDC words before writing fDataTof // + channel map corrected (correspondence between PMT# and TDC channel) // + new conversion pC to ADC // + scaling factors to tune MC data to real ones // - HOW-TO-DIGIT.TXT file updated (E.Vannuccini) // // v. 1.5 (S.Orsi, Roma2, 11 October 2007) // - file HOW-TO-DIGIT.TXT written (mainly by S.Bottai): contains information on how to run the chain from the DB creation to DarthVader through Digitizer, yoda and YodaProfiler // // v. 1.4 (S.Orsi, Roma2, 11 October 2007) // - similar to v. 1.3 // // v. 1.3 (S.Orsi, Roma2, 11 October 2007) // - TOF: + changes and major corrections in ADC and TDC (W.Menn, D.Campana, S.Orsi) // + flag DEBUG added (default: DEBUG=false; if true, print on screem some messages to debug) // - AC: + shift register filled; // + minor bugs fixed // + counters increment by 1 at each event (they were constant before) // - S4: + minor bugs fixed (S.Borisov) // + the routine DigitizeS4() is now inserted in the code. For DV versions earlier than October 2007, the option "-S4" is still needed because S4 calibration is not present // - corrected severe bug that caused yoda to process only part of large data files (file.pam). S4 data were written to file even if the DigitizeS4() routine was not called. Previous versions of the digitizer should have crashed always, but unexpectedly worked on small files. // - file HOW-TO-DIGIT.TXT written (mainly by S.Bottai): contains information on how to run the chain from the DB creation to DarthVader through Digitizer, yoda and YodaProfiler. The file is included only from release 1.5 // - digitizer version numbers modified according to CVS repository (v1.00->v1.1, v1.01->v1.2) // - some flags (LDFLAGS) in the Makefile have been changed; they are needed on some machines (e.g. pamelatov.roma2.infn.it), but give errors on others (e.g. pamelaws02.na.infn.it). See Makefile for details if the compilation fails due to library errors) // // v. 1.2 (S.Orsi, Roma2, 28 September 2007) // - TOF TDC inserted (to be tuned) // - improved tracker dEdx, tracker saturation inserted (S.Bottai) // - S4 routine inserted in code (not in read-out) // - bugs with packet length corrected // - small bugs fixed // - NB: Run DarthVader with the options: "DarthVader -zerofill -idRun 1 -ALL +RUN +CAL [ --no-crosstalk ] +TRK +TOF +AC +TRG +ND" or "DarthVader -zerofill -idRun 1 +CAL [ --no-crosstalk ] -S4", i.e. without S4 // // v. 1.1 (S.Orsi, Roma2, 13 September 2007) // - trigger (dummy), tof (preliminary, adc only), ND (preliminary) inserted // - various changes and bug fixes // - padding increased (fPADbuffer=64, fPadding = padbytes+64) // // v. beta (J.Conrad, KTH, 2006) // - compiles; includes pscu, calo, trk, ac; not present on cvs // ///////////////////////////////////////////////////////////////////////// #include #include #include #include "Riostream.h" #include "TFile.h" #include "TDirectory.h" #include "TTree.h" #include "TLeafI.h" #include "TH1.h" #include "TH2.h" #include "TProfile.h" #include "TGraph.h" #include "TMath.h" #include "Digitizer.h" int Error; //to be removed soon // Define the names of the Fortran common blocks for the different OSs #ifndef WIN32 #define PAWC_SIZE 5000000 # define pawc pawc_ # define quest quest_ # define hcbits hcbits_ # define hcbook hcbook_ # define rzcl rzcl_ int pawc[PAWC_SIZE]; int quest[100]; int hcbits[37]; int hcbook[51]; int rzcl[11]; #else // on windows /pawc/ must have the same length as in libPacklib.a !! #define PAWC_SIZE 2000000 # define pawc PAWC # define quest QUEST # define hcbits HCBITS # define hcbook HCBOOK # define rzcl RZCL extern "C" int pawc[PAWC_SIZE]; extern "C" int quest[100]; extern "C" int hcbits[37]; extern "C" int hcbook[51]; extern "C" int rzcl[11]; #endif int *iq, *lq; float *q; char idname[128]; int nentries; char chtitl[128]; int ncx,ncy,nwt,idb; int lcont, lcid, lcdir; float xmin,xmax,ymin,ymax; const Int_t kMIN1 = 7; const Int_t kMAX1 = 8; #if defined __linux //On linux Fortran wants this, so we give to it! int xargv=0; int xargc=0; void MAIN__() {} #endif // Define the names of the Fortran subroutine and functions for the different OSs #ifndef WIN32 # define hlimit hlimit_ # define hropen hropen_ # define hrin hrin_ # define hnoent hnoent_ # define hgive hgive_ # define hgiven hgiven_ # define hprntu hprntu_ # define hgnpar hgnpar_ # define hgnf hgnf_ # define hgnt hgnt_ # define rzink rzink_ # define hdcofl hdcofl_ # define hmaxim hmaxim_ # define hminim hminim_ # define hdelet hdelet_ # define hntvar2 hntvar2_ # define hbname hbname_ # define hbnamc hbnamc_ # define hbnam hbnam_ # define hi hi_ # define hie hie_ # define hif hif_ # define hij hij_ # define hix hix_ # define hijxy hijxy_ # define hije hije_ # define hcdir hcdir_ # define zitoh zitoh_ # define uhtoc uhtoc_ # define type_of_call # define DEFCHAR const char* # define PASSCHAR(string) string #else # define hlimit HLIMIT # define hropen HROPEN # define hrin HRIN # define hnoent HNOENT # define hgive HGIVE # define hgiven HGIVEN # define hprntu HPRNTU # define hgnpar HGNPAR # define hgnf HGNF # define hgnt HGNT # define rzink RZINK # define hdcofl HDCOFL # define hmaxim HMAXIM # define hminim HMINIM # define hdelet HDELET # define hntvar2 HNTVAR2 # define hbname HBNAME # define hbnamc HBNAMC # define hbnam HBNAM # define hi HI # define hie HIE # define hif HIF # define hij HIJ # define hix HIX # define hijxy HIJXY # define hije HIJE # define hcdir HCDIR # define zitoh ZITOH # define uhtoc UHTOC # define type_of_call _stdcall # define DEFCHAR const char*, const int # define PASSCHAR(string) string, strlen(string) #endif extern "C" void type_of_call hlimit(const int&); #ifndef WIN32 extern "C" void type_of_call hropen(const int&,DEFCHAR,DEFCHAR,DEFCHAR, const int&,const int&,const int,const int,const int); #else extern "C" void type_of_call hropen(const int&,DEFCHAR,DEFCHAR,DEFCHAR, const int&,const int&); #endif extern "C" void type_of_call hrin(const int&,const int&,const int&); extern "C" void type_of_call hnoent(const int&,const int&); #ifndef WIN32 extern "C" void type_of_call hgive(const int&,DEFCHAR,const int&,const float&,const float&, const int&,const float&,const float&,const int&,const int&,const int); #else extern "C" void type_of_call hgive(const int&,DEFCHAR,const int&,const float&,const float&, const int&,const float&,const float&,const int&,const int&); #endif #ifndef WIN32 extern "C" void type_of_call hgiven(const int&,DEFCHAR,const int&,DEFCHAR, const float&,const float&,const int,const int); #else extern "C" void type_of_call hgiven(const int&,DEFCHAR,const int&,DEFCHAR, const float&,const float&); #endif #ifndef WIN32 extern "C" void type_of_call hntvar2(const int&,const int&,DEFCHAR,DEFCHAR,DEFCHAR,int&,int&,int&,int&,int&,const int,const int, const int); #else extern "C" void type_of_call hntvar2(const int&,const int&,DEFCHAR,DEFCHAR,DEFCHAR,int&,int&,int&,int&,int&); #endif #ifndef WIN32 extern "C" void type_of_call hbnam(const int&,DEFCHAR,const int&,DEFCHAR,const int&,const int, const int); #else extern "C" void type_of_call hbnam(const int&,DEFCHAR,const int&,DEFCHAR,const int&); #endif extern "C" void type_of_call hprntu(const int&); extern "C" void type_of_call hgnpar(const int&,const char *,const int); extern "C" void type_of_call hgnf(const int&,const int&,const float&,const int&); extern "C" void type_of_call hgnt(const int&,const int&,const int&); extern "C" void type_of_call rzink(const int&,const int&,const char *,const int); extern "C" void type_of_call hdcofl(); extern "C" void type_of_call hmaxim(const int&,const float&); extern "C" void type_of_call hminim(const int&,const float&); extern "C" void type_of_call hdelet(const int&); extern "C" float type_of_call hi(const int&,const int&); extern "C" float type_of_call hie(const int&,const int&); extern "C" float type_of_call hif(const int&,const int&); extern "C" float type_of_call hij(const int&,const int&,const int&); extern "C" void type_of_call hix(const int&,const int&,const float&); extern "C" void type_of_call hijxy(const int&,const int&,const int&,const float&,const float&); extern "C" float type_of_call hije(const int&,const int&,const int&); #ifndef WIN32 extern "C" void type_of_call hcdir(DEFCHAR,DEFCHAR ,const int,const int); #else extern "C" void type_of_call hcdir(DEFCHAR,DEFCHAR); #endif extern "C" void type_of_call zitoh(const int&,const int&,const int&); #ifndef WIN32 extern "C" void type_of_call uhtoc(const int&,const int&,DEFCHAR,int&,const int); #else extern "C" void type_of_call uhtoc(const int&,const int&,DEFCHAR,int&); #endif extern void convert_directory(const char*, char* file_raw); extern void convert_1d(Int_t id); extern void convert_cwn(Int_t id,char* file_raw); Int_t golower = 1; Int_t bufsize = 64000; Int_t optcwn = 1; int main(int argc, char **argv) { if (argc < 2) { printf("Error: Pamelagp2Digits \n"); printf("Pamelagp2Digits file.hbook file.root file.pam [compress] [tolower] [lrecl] [bufsize] [optcwn] \n"); printf(" if file.root is not given it will be = file.root\n"); printf(" compress = 1 by default (use 0 for no compression)\n"); printf(" tolower = 1 by default (use 0 to keep case of column names)\n"); printf(" lrecl =0 by default (must be specified if >8092)\n"); printf(" bufsize = 8000 by default (branch buffer size)\n"); printf(" for cwn ntuple only: optcwn = 1 (default) 1-byte int -> char, 2-byte int -> short, (use 0 to keep 4-byte int) \n"); return 1; } lq = &pawc[9]; iq = &pawc[17]; void *qq = iq; q = (float*)qq; char *file_in=argv[1]; char *file_out = " "; char *file_raw; Int_t compress = 1; int ier=0, record_size=0; if (argc > 8) { optcwn = atoi(argv[8]); } if (argc > 7) { bufsize = atoi(argv[7]); } if (argc > 6) { record_size = atoi(argv[6]); } if (argc > 5) { golower = atoi(argv[5]); } if (argc > 4) { compress = atoi(argv[4]); } if (argc > 3) { file_raw=argv[3]; } else { file_raw= new char[2000]; strcpy(file_raw,file_in); char *dot = strrchr(file_raw,'.'); if (dot) strcpy(dot+1,"pam"); else strcat(file_out,".pam"); } if (argc > 2) { file_out=argv[2]; } else { file_out= new char[2000]; strcpy(file_out,file_in); char *dot = strrchr(file_out,'.'); if (dot) strcpy(dot+1,"gp.root"); //modified S.Orsi 09/'07 else strcat(file_out,".gp.root"); } #if defined(_HIUX_SOURCE) && !defined(__GNUC__) hf_fint((char *)NULL); #endif int pawc_size = PAWC_SIZE; hlimit(pawc_size); int lun = 10; #ifndef WIN32 hropen(lun,PASSCHAR("example"),PASSCHAR(file_in),PASSCHAR("p"),record_size,ier,7,strlen(file_in),1); #else hropen(lun,PASSCHAR("example"),PASSCHAR(file_in),PASSCHAR("p"),record_size,ier); #endif if (ier) printf (" Error on hropen was %d \n", ier); if (quest[0]) { printf("Error cannot open input file: %s\n",file_in); return 1; } char root_file_title[2000]; TFile* hfile= TFile::Open(file_out,"RECREATE",root_file_title,compress); if (!hfile) { printf("Error: can't open output file: %s \n",file_out); return 1; } convert_directory("//example",file_raw); //hfile->Write(); //hfile->ls(); hfile->Close(); delete hfile; return(0); } //____________________________________________________________________________ void convert_directory(const char *dir, char* file_raw) { printf(" Converting directory %s\n",dir); Int_t id; // Int_t nastycase=0; // Int_t nastyprint=0; // Int_t idold = 0; for (Int_t key=1;key<1000000;key++) { int z0 = 0; rzink(key,z0,"S",1); if (quest[0]) break; if (quest[13] & 8) { continue; // if (!nastyprint) { // printf("Found nasty Hbook case!! You had an Hbook error message\n"); // printf(" when creating the file (too many records)\n"); // printf(" Hbook file should have been created with a bigger LRECL\n"); // printf(" ROOT will try to recover\n"); // nastyprint = 1; // } // nastycase = 1; } id = quest[20]; // if (id == idold && nastycase) continue; // nastycase = 0; // idold = id; int i999 = 999; hrin(id,i999,0); if (quest[0]) { printf("Error cannot read ID = %d\n",id); break; } hdcofl(); lcid = hcbook[10]; lcont = lq[lcid-1]; if (hcbits[3]) { convert_cwn(id,file_raw); hdelet(id); continue; } if (hcbits[0]) { convert_1d(id); hdelet(id); continue; } } // converting subdirectories of this directory const Int_t kKLS = 26; const Int_t kKNSD = 23; lcdir = rzcl[2]; Int_t ls = iq[lcdir+kKLS]; Int_t ndir = iq[lcdir+kKNSD]; Int_t nch=16; Int_t ihdir[4]; Int_t ncw = 4; TDirectory *cursav = gDirectory; Int_t i; char chdir[17]; char hbookdir[17]; for (Int_t k=0;k0;i--) { if (chdir[i] == 0) continue; if (chdir[i] != ' ') break; chdir[i] = 0; } #ifndef WIN32 hcdir(PASSCHAR(hbookdir),PASSCHAR(" "),16,1); #else hcdir(PASSCHAR(hbookdir),PASSCHAR(" ")); #endif TDirectory *newdir = new TDirectory(chdir,chdir); newdir->cd(); convert_directory(chdir, file_raw); #ifndef WIN32 hcdir(PASSCHAR("\\"),PASSCHAR(" "),1,1); #else hcdir(PASSCHAR("\\"),PASSCHAR(" ")); #endif newdir->Write(); cursav->cd(); } } //____________________________________________________________________________ void convert_1d(Int_t id) { if (id > 0) sprintf(idname,"h%d",id); else sprintf(idname,"h_%d",-id); hnoent(id,nentries); #ifndef WIN32 hgive(id,chtitl,ncx,xmin,xmax,ncy,ymin,ymax,nwt,idb,80); #else hgive(id,chtitl,80,ncx,xmin,xmax,ncy,ymin,ymax,nwt,idb); #endif chtitl[4*nwt] = 0; TH1F *h1; Int_t i; if (hcbits[5]) { Int_t lbins = lq[lcid-2]; Double_t *xbins = new Double_t[ncx+1]; for (i=0;i<=ncx;i++) xbins[i] = q[lbins+i+1]; h1 = new TH1F(idname,chtitl,ncx,xbins); delete [] xbins; } else { h1 = new TH1F(idname,chtitl,ncx,xmin,xmax); } if (hcbits[8]) h1->Sumw2(); TGraph *gr = 0; if (hcbits[11]) { gr = new TGraph(ncx); h1->GetListOfFunctions()->Add(gr); } Float_t x; for (i=0;i<=ncx+1;i++) { x = h1->GetBinCenter(i); h1->Fill(x,hi(id,i)); if (hcbits[8]) h1->SetBinError(i,hie(id,i)); if (gr && i>0 && i<=ncx) gr->SetPoint(i,x,hif(id,i)); } Float_t ymin, ymax; if (hcbits[19]) { ymax = q[lcid+kMAX1]; h1->SetMaximum(ymax); } if (hcbits[20]) { ymin = q[lcid+kMIN1]; h1->SetMinimum(ymin); } h1->SetEntries(nentries); h1->Write(); delete h1; } //____________________________________________________________________________ void convert_cwn(Int_t id,char* file_raw) { const int kNchar=9; int nvar; int ier=0; int i,j; int nsub,itype,isize,ielem; char *chtag_out; float *x; float rmin[1000], rmax[1000]; if (id > 0) sprintf(idname,"h%d",id); else sprintf(idname,"h_%d",-id); hnoent(id,nentries); printf(" Converting CWN with ID= %d, nentries = %d\n",id,nentries); nvar=0; #ifndef WIN32 hgiven(id,chtitl,nvar,PASSCHAR(""),rmin[0],rmax[0],80,0); #else hgiven(id,chtitl,80,nvar,PASSCHAR(""),rmin[0],rmax[0]); #endif chtag_out = new char[nvar*kNchar+1]; Int_t *charflag = new Int_t[nvar]; Int_t *lenchar = new Int_t[nvar]; Int_t *boolflag = new Int_t[nvar]; Int_t *lenbool = new Int_t[nvar]; UChar_t *boolarr = new UChar_t[10000]; x = new float[nvar]; char *bigbuf = new char[2500000]; chtag_out[nvar*kNchar]=0; for (i=0;i<80;i++)chtitl[i]=0; #ifndef WIN32 hgiven(id,chtitl,nvar,chtag_out,rmin[0],rmax[0],80,kNchar); #else hgiven(id,chtitl,80,nvar,chtag_out,kNchar,rmin[0],rmax[0]); #endif #ifndef WIN32 hbnam(id,PASSCHAR(" "),bigbuf[0],PASSCHAR("$CLEAR"),0,1,6); #else hbnam(id,PASSCHAR(" "),bigbuf[0],PASSCHAR("$CLEAR"),0); #endif Int_t bufpos = 0; Int_t isachar = 0; Int_t isabool = 0; char fullname[1024]; char name[512]; char block[512]; char oldblock[512]; Int_t nbits = 0; strcpy(oldblock,"OLDBLOCK"); Int_t oldischar = -1; for (i=80;i>0;i--) {if (chtitl[i] == ' ') chtitl[i] = 0; } TTree *tree = new TTree(idname,chtitl); for(i=0; i0;j--) { if(golower) name[j] = tolower(name[j]); if (name[j] == ' ') name[j] = 0; } if (golower == 2) name[0] = tolower(name[0]); for (j=1022;j>0;j--) { if(golower && fullname[j-1] != '[') fullname[j] = tolower(fullname[j]); // convert also character after [, if golower == 2 if (golower == 2) fullname[j] = tolower(fullname[j]); if (fullname[j] == ' ') fullname[j] = 0; } // convert also first character, if golower == 2 if (golower == 2) fullname[0] = tolower(fullname[0]); for (j=510;j>0;j--) { if (block[j] == ' ') block[j] = 0; else break; } if (itype == 1) { if( isize == 4 ) strcat(fullname,"/F"); else if( isize == 8) strcat(fullname,"/D"); } // add support for 1-byte (Char_t) and 2-byte (Short_t) integers Int_t nBytesUsed = 4; // default for integers if( itype == 2 ) { if( optcwn == 1 ) { if( nbits > 16 ) { strcat(fullname,"/I"); } else { if( nbits > 8 ) { strcat(fullname,"/S"); nBytesUsed = 2; } else { strcat(fullname,"/B"); nBytesUsed = 1; } } } else { strcat(fullname,"/I"); } } // add support for 1-byte (UChar_t) and 2-byte (UShort_t) integers if ( itype == 3 ) { if( optcwn == 1 ) { if( nbits > 16) { strcat(fullname,"/i"); } else { if( nbits > 8 ) { strcat(fullname,"/s"); nBytesUsed = 2; } else { strcat(fullname,"/b"); nBytesUsed = 1; } } } else { strcat(fullname,"/i"); } } // if (itype == 4) strcat(fullname,"/i"); if (itype == 4) strcat(fullname,"/b"); if (itype == 5) strcat(fullname,"/C"); //printf("Creating branch:%s, block:%s, fullname:%s, nsub=%d, itype=%d, isize=%d, ielem=%d\n",name,block,fullname,nsub,itype,isize,ielem); Int_t ischar; if (itype == 5) ischar = 1; else ischar = 0; if (ischar != oldischar || strcmp(oldblock,block) != 0) { strcpy(oldblock,block); oldischar = ischar; Int_t lblock = strlen(block); Long_t add= (Long_t)&bigbuf[bufpos]; #ifndef WIN32 hbnam(id,PASSCHAR(block),add,PASSCHAR("$SET"),ischar,lblock,4); #else hbnam(id,PASSCHAR(block),add,PASSCHAR("$SET"),ischar); #endif } TBranch *branch = tree->Branch(name,(void*)&bigbuf[bufpos],fullname,bufsize); boolflag[i] = -10; charflag[i] = 0; if (itype == 4) {isabool++; boolflag[i] = bufpos; lenbool[i] = ielem;} bufpos += isize*ielem; if (ischar) {isachar++; charflag[i] = bufpos-1; lenchar[i] = isize*ielem;} TObjArray *ll= branch->GetListOfLeaves(); TLeaf *leaf = (TLeaf*)ll->UncheckedAt(0); if (!leaf) continue; TLeafI *leafcount = (TLeafI*)leaf->GetLeafCount(); if (leafcount) { if (leafcount->GetMaximum() <= 0) leafcount->SetMaximum(ielem); } } Int_t cf,l; for(i=1;i<=nentries;i++) { hgnt(id,i,ier); if (isabool) { // if column is boolean for (j=0;j-1) { for (l=0;l 16) { // do nothing for integers of 4 byte } else { if( nbits > 8 ) { nBytesUsed = 2; } else { nBytesUsed = 1; } } } if(nBytesUsed == 1) { for(Int_t index = 0; index < ielem; index++) { // shift all chars with data to be one after another bigbuf[bufpos + index*nBytesUsed ] = bigbuf[bufpos + index * isize]; } } else { if(nBytesUsed == 2) { for(Int_t index = 0; index < ielem; index++) { // shift all shorts ( 2 chars) with data to be one after another bigbuf[bufpos + index*nBytesUsed ] = bigbuf[bufpos + index * isize]; bigbuf[bufpos + index*nBytesUsed+1 ] = bigbuf[bufpos + index * isize+1]; } } } bufpos += isize*ielem; } } tree->Fill(); } tree->Write(); std::cout << "Invoking Digitizer" << endl << flush; Digitizer* dig = new Digitizer(tree,file_raw); dig->Loop(); dig->Close(); std::cout << "Finished" << endl << flush; }