--- PamelaDigitizer/Digitizer.cxx 2007/10/11 11:29:21 1.3 +++ PamelaDigitizer/Digitizer.cxx 2007/10/31 18:17:58 1.4 @@ -1,2571 +1,2705 @@ -// ------ PAMELA Digitizer ------ -// -// Date, release and how-to: see file Pamelagp2Digits.cxx -// -// NB: Check length physics packet [packet type (0x10 = physics data)] -// -#include -#include -#include -#include -#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 "TMath.h" -#include "TRandom.h" -#include "TSQLServer.h" -#include "TSystem.h" -// -#include "Digitizer.h" -#include "CRC.h" -// -#include -#include -#include -#include "GLTables.h" -// -extern "C"{ - short crc(short, short); -}; -// - -Digitizer::Digitizer(TTree* tree, char* &file_raw){ - fhBookTree = tree; - fFilename = file_raw; - fCounter = 0; - fOBT = 0; - - // - // DB connections - // - TString host = "mysql://localhost/pamelaprod"; - TString user = "anonymous"; - TString psw = ""; - // - const char *pamdbhost=gSystem->Getenv("PAM_DBHOST"); - const char *pamdbuser=gSystem->Getenv("PAM_DBUSER"); - const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW"); - if ( !pamdbhost ) pamdbhost = ""; - if ( !pamdbuser ) pamdbuser = ""; - if ( !pamdbpsw ) pamdbpsw = ""; - if ( strcmp(pamdbhost,"") ) host = pamdbhost; - if ( strcmp(pamdbuser,"") ) user = pamdbuser; - if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw; - fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data()); - // - GL_TABLES *glt = new GL_TABLES(host,user,psw); - if ( glt->IsConnected(fDbc) ) printf("\n DB INFORMATION:\n SQL: %s Version: %s Host %s Port %i \n\n",fDbc->GetDBMS(),fDbc->ServerInfo(),fDbc->GetHost(),fDbc->GetPort()); - // - // Use UTC in the DB and make timeout bigger - // - stringstream myquery; - myquery.str(""); - myquery << "SET time_zone='+0:00'"; - fDbc->Query(myquery.str().c_str()); - myquery.str(""); - myquery << "SET wait_timeout=173000;"; - fDbc->Query(myquery.str().c_str()); - // - - std:: cout << "preparing tree" << endl; - - // prepare tree - fhBookTree->SetBranchAddress("Irun",&Irun); - fhBookTree->SetBranchAddress("Ievnt",&Ievnt); - fhBookTree->SetBranchAddress("Ipa",&Ipa); - fhBookTree->SetBranchAddress("X0",&X0); - fhBookTree->SetBranchAddress("Y0",&Y0); - fhBookTree->SetBranchAddress("Z0",&Z0); - fhBookTree->SetBranchAddress("Theta",&Theta); - fhBookTree->SetBranchAddress("Phi",&Phi); - fhBookTree->SetBranchAddress("P0",&P0); - fhBookTree->SetBranchAddress("Nthtof",&Nthtof); - fhBookTree->SetBranchAddress("Ipltof",Ipltof); - fhBookTree->SetBranchAddress("Ipaddle",Ipaddle); - fhBookTree->SetBranchAddress("Ipartof",Ipartof); - fhBookTree->SetBranchAddress("Xintof",Xintof); - fhBookTree->SetBranchAddress("Yintof",Yintof); - fhBookTree->SetBranchAddress("Zintof",Zintof); - fhBookTree->SetBranchAddress("Xouttof",Xouttof); - fhBookTree->SetBranchAddress("Youttof",Youttof); - fhBookTree->SetBranchAddress("Zouttof",Zouttof); - fhBookTree->SetBranchAddress("Ereltof",Ereltof); - fhBookTree->SetBranchAddress("Timetof",Timetof); - fhBookTree->SetBranchAddress("Pathtof",Pathtof); - fhBookTree->SetBranchAddress("P0tof",P0tof); - fhBookTree->SetBranchAddress("Nthcat",&Nthcat); - fhBookTree->SetBranchAddress("Iparcat",Iparcat); - fhBookTree->SetBranchAddress("Icat",Icat); - fhBookTree->SetBranchAddress("Xincat",Xincat); - fhBookTree->SetBranchAddress("Yincat",Yincat); - fhBookTree->SetBranchAddress("Zincat",Zincat); - fhBookTree->SetBranchAddress("Xoutcat",Xoutcat); - fhBookTree->SetBranchAddress("Youtcat",Youtcat); - fhBookTree->SetBranchAddress("Zoutcat",Zoutcat); - fhBookTree->SetBranchAddress("Erelcat",Erelcat); - fhBookTree->SetBranchAddress("Timecat",Timecat); - fhBookTree->SetBranchAddress("Pathcat",Pathcat); - fhBookTree->SetBranchAddress("P0cat",P0cat); - fhBookTree->SetBranchAddress("Nthcas",&Nthcas); - fhBookTree->SetBranchAddress("Iparcas",Iparcas); - fhBookTree->SetBranchAddress("Icas",Icas); - fhBookTree->SetBranchAddress("Xincas",Xincas); - fhBookTree->SetBranchAddress("Yincas",Yincas); - fhBookTree->SetBranchAddress("Zincas",Zincas); - fhBookTree->SetBranchAddress("Xoutcas",Xoutcas); - fhBookTree->SetBranchAddress("Youtcas",Youtcas); - fhBookTree->SetBranchAddress("Zoutcas",Zoutcas); - fhBookTree->SetBranchAddress("Erelcas",Erelcas); - fhBookTree->SetBranchAddress("Timecas",Timecas); - fhBookTree->SetBranchAddress("Pathcas",Pathcas); - fhBookTree->SetBranchAddress("P0cas",P0cas); - fhBookTree->SetBranchAddress("Nthspe",&Nthspe); - fhBookTree->SetBranchAddress("Iparspe",Iparspe); - fhBookTree->SetBranchAddress("Itrpb",Itrpb); - fhBookTree->SetBranchAddress("Itrsl",Itrsl); - fhBookTree->SetBranchAddress("Itspa",Itspa); - fhBookTree->SetBranchAddress("Xinspe",Xinspe); - fhBookTree->SetBranchAddress("Yinspe",Yinspe); - fhBookTree->SetBranchAddress("Zinspe",Zinspe); - fhBookTree->SetBranchAddress("Xoutspe",Xoutspe); - fhBookTree->SetBranchAddress("Youtspe",Youtspe); - fhBookTree->SetBranchAddress("Zoutspe",Zoutspe); - fhBookTree->SetBranchAddress("Xavspe",Xavspe); - fhBookTree->SetBranchAddress("Yavspe",Yavspe); - fhBookTree->SetBranchAddress("Zavspe",Zavspe); - fhBookTree->SetBranchAddress("Erelspe",Erelspe); - fhBookTree->SetBranchAddress("Pathspe",Pathspe); - fhBookTree->SetBranchAddress("P0spe",P0spe); - fhBookTree->SetBranchAddress("Nxmult",Nxmult); - fhBookTree->SetBranchAddress("Nymult",Nymult); - fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx); - fhBookTree->SetBranchAddress("Npstripx",Npstripx); - fhBookTree->SetBranchAddress("Ntstripx",Ntstripx); - fhBookTree->SetBranchAddress("Istripx",Istripx); - fhBookTree->SetBranchAddress("Qstripx",Qstripx); - fhBookTree->SetBranchAddress("Xstripx",Xstripx); - fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy); - fhBookTree->SetBranchAddress("Npstripy",Npstripy); - fhBookTree->SetBranchAddress("Ntstripy",Ntstripy); - fhBookTree->SetBranchAddress("Istripy",Istripy); - fhBookTree->SetBranchAddress("Qstripy",Qstripy); - fhBookTree->SetBranchAddress("Ystripy",Ystripy); - fhBookTree->SetBranchAddress("Nthcali",&Nthcali); - fhBookTree->SetBranchAddress("Icaplane",Icaplane); - fhBookTree->SetBranchAddress("Icastrip",Icastrip); - fhBookTree->SetBranchAddress("Icamod",Icamod); - fhBookTree->SetBranchAddress("Enestrip",Enestrip); - fhBookTree->SetBranchAddress("Nthcal",&Nthcal); - fhBookTree->SetBranchAddress("Icapl",Icapl); - fhBookTree->SetBranchAddress("Icasi",Icasi); - fhBookTree->SetBranchAddress("Icast",Icast); - fhBookTree->SetBranchAddress("Xincal",Xincal); - fhBookTree->SetBranchAddress("Yincal",Yincal); - fhBookTree->SetBranchAddress("Zincal",Zincal); - fhBookTree->SetBranchAddress("Erelcal",Erelcal); - fhBookTree->SetBranchAddress("Nthnd",&Nthnd); - fhBookTree->SetBranchAddress("Itubend",Itubend); - fhBookTree->SetBranchAddress("Iparnd",Iparnd); - fhBookTree->SetBranchAddress("Xinnd",Xinnd); - fhBookTree->SetBranchAddress("Yinnd",Yinnd); - fhBookTree->SetBranchAddress("Zinnd",Zinnd); - fhBookTree->SetBranchAddress("Xoutnd",Xoutnd); - fhBookTree->SetBranchAddress("Youtnd",Youtnd); - fhBookTree->SetBranchAddress("Zoutnd",Zoutnd); - fhBookTree->SetBranchAddress("Erelnd",Erelnd); - fhBookTree->SetBranchAddress("Timend",Timend); - fhBookTree->SetBranchAddress("Pathnd",Pathnd); - fhBookTree->SetBranchAddress("P0nd",P0nd); - fhBookTree->SetBranchAddress("Nthcard",&Nthcard); - fhBookTree->SetBranchAddress("Iparcard",Iparcard); - fhBookTree->SetBranchAddress("Icard",Icard); - fhBookTree->SetBranchAddress("Xincard",Xincard); - fhBookTree->SetBranchAddress("Yincard",Yincard); - fhBookTree->SetBranchAddress("Zincard",Zincard); - fhBookTree->SetBranchAddress("Xoutcard",Xoutcard); - fhBookTree->SetBranchAddress("Youtcard",Youtcard); - fhBookTree->SetBranchAddress("Zoutcard",Zoutcard); - fhBookTree->SetBranchAddress("Erelcard",Erelcard); - fhBookTree->SetBranchAddress("Timecard",Timecard); - fhBookTree->SetBranchAddress("Pathcard",Pathcard); - fhBookTree->SetBranchAddress("P0card",P0card); - - fhBookTree->SetBranchStatus("*",0); - -}; - - - -void Digitizer::Close(){ - - delete fhBookTree; - -}; - - - - -void Digitizer::Loop() { - // - // opens the raw output file and loops over the events - // - fOutputfile.open(fFilename, ios::out | ios::binary); - //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary); - // - // Load in memory and save at the beginning of file the calorimeter calibration - // - CaloLoadCalib(); - DigitizeCALOCALIB(); - - // load, digitize and write tracker calibration - LoadTrackCalib(); - - DigitizeTrackCalib(1); - UInt_t length=fTracklength*2; - DigitizePSCU(length,0x12); - AddPadding(); - WriteTrackCalib(); - - DigitizeTrackCalib(2); - length=fTracklength*2; - DigitizePSCU(length,0x13); - AddPadding(); - WriteTrackCalib(); - - LoadMipCor(); // some initialization of parameters -not used now- - // end loading, digitizing and writing tracker calibration - - // - // loops over the events - // - - Int_t nentries = fhBookTree->GetEntriesFast(); - Long64_t nbytes = 0; - for (Int_t i=0; iGetEntry(i); - // read detectors sequentially: - // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf - // on pamelatov: - // /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp - DigitizeTRIGGER(); - DigitizeTOF(); - DigitizeAC(); - DigitizeCALO(); - DigitizeTrack(); - DigitizeS4(); - DigitizeND(); - // - // Add padding to 64 bits - // - AddPadding(); -// - // Create CPU header, we need packet type (0x10 = physics data) and packet length. - // - UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer; - //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer; - DigitizePSCU(length,0x10); - if ( !i%100 ) std::cout << "writing event " << i << endl; - WriteData(); - }; - - fOutputfile.close(); - std::cout << "files closed" << endl << flush; - -}; - -void Digitizer::AddPadding(){ - // - Float_t pd0 = (fLen+16)/64.; - Float_t pd1 = pd0 - (Float_t)int(pd0); - Float_t padfrac = 64. - pd1 * 64.; - // - UInt_t padbytes = (UInt_t)padfrac; - if ( padbytes > 0 && padbytes < 64 ){ - // - // here the padding length - // - fPadding = padbytes+64; - // - // random padding bytes - // - for (Int_t ur=0; ur<32; ur++){ - fDataPadding[ur] = (UShort_t)rand(); - }; - }; -}; - - -void Digitizer::DigitizePSCU(UInt_t length, UChar_t type) { - // - UChar_t buff[16]; - // - // CPU signature - // - buff[0] = 0xFA; - buff[1] = 0xFE; - buff[2] = 0xDE; - // - // packet type (twice) - // - buff[3] = type; - buff[4] = type; - // - // counter - // - fCounter++; - while ( fCounter > 16777215 ){ - fCounter -= 16777215; - }; - // - buff[5] = (UChar_t)(fCounter >> 16); - buff[6] = (UChar_t)(fCounter >> 8); - buff[7] = (UChar_t)fCounter; - // - // on board time - // - ULong64_t obt = fOBT + 30LL; - // - while ( obt > 4294967295LL ){ - obt -= 4294967295LL; - }; - fOBT = (UInt_t)obt; - // - buff[8] = (UChar_t)(fOBT >> 24); - buff[9] = (UChar_t)(fOBT >> 16); - buff[10] = (UChar_t)(fOBT >> 8); - buff[11] = (UChar_t)fOBT; - // - // Packet length - // - fLen = length; - // - buff[12] = (UChar_t)(fLen >> 16); - buff[13] = (UChar_t)(fLen >> 8); - buff[14] = (UChar_t)fLen; - // - // CPU header CRC - // - buff[15] = (BYTE)CM_Compute_CRC16((UINT16)0, (BYTE*)&buff, (UINT32)15); - // - memcpy(fDataPSCU,buff,16*sizeof(UChar_t)); - // -}; - -void Digitizer::ClearCaloCalib(Int_t s){ - // - fcstwerr[s] = 0; - fcperror[s] = 0.; - for ( Int_t d=0 ; d<11 ;d++ ){ - Int_t pre = -1; - for ( Int_t j=0; j<96 ;j++){ - if ( j%16 == 0 ) pre++; - fcalped[s][d][j] = 0.; - fcstwerr[s] = 0.; - fcperror[s] = 0.; - fcalgood[s][d][j] = 0.; - fcalthr[s][d][pre] = 0.; - fcalrms[s][d][j] = 0.; - fcalbase[s][d][pre] = 0.; - fcalvar[s][d][pre] = 0.; - }; - }; - return; -} - -Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){ - // - // - UInt_t e = 0; - if ( s == 0 ) e = 0; - if ( s == 1 ) e = 2; - if ( s == 2 ) e = 3; - if ( s == 3 ) e = 1; - // - ifstream myfile; - myfile.open(fcalname.Data()); - if ( !myfile ){ - return(-107); - }; - myfile.close(); - // - TFile *File = new TFile(fcalname.Data()); - if ( !File ) return(-108); - TTree *tr = (TTree*)File->Get("CalibCalPed"); - if ( !tr ) return(-109); - // - TBranch *calo = tr->GetBranch("CalibCalPed"); - // - pamela::CalibCalPedEvent *ce = 0; - tr->SetBranchAddress("CalibCalPed", &ce); - // - Long64_t ncalibs = calo->GetEntries(); - // - if ( !ncalibs ) return(-110); - // - calo->GetEntry(calibno); - // - if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { - fcstwerr[s] = ce->cstwerr[s]; - fcperror[s] = ce->cperror[s]; - for ( Int_t d=0 ; d<11 ;d++ ){ - Int_t pre = -1; - for ( Int_t j=0; j<96 ;j++){ - if ( j%16 == 0 ) pre++; - fcalped[s][d][j] = ce->calped[e][d][j]; - fcalgood[s][d][j] = ce->calgood[e][d][j]; - fcalthr[s][d][pre] = ce->calthr[e][d][pre]; - fcalrms[s][d][j] = ce->calrms[e][d][j]; - fcalbase[s][d][pre] = ce->calbase[e][d][pre]; - fcalvar[s][d][pre] = ce->calvar[e][d][pre]; - }; - }; - } else { - printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n "); - File->Close(); - return(-111); - }; - File->Close(); - return(0); -} - - -void Digitizer::DigitizeCALOCALIB() { - // - // Header of the four sections - // - fSecCalo[0] = 0xAA00; // XE - fSecCalo[1] = 0xB100; // XO - fSecCalo[2] = 0xB600; // YE - fSecCalo[3] = 0xAD00; // YO - // - // length of the data is 0x1215 - // - fSecCALOLength[0] = 0x1215; // XE - fSecCALOLength[1] = 0x1215; // XO - fSecCALOLength[2] = 0x1215; // YE - fSecCALOLength[3] = 0x1215; // YO - // - Int_t chksum = 0; - UInt_t tstrip = 0; - UInt_t fSecPointer = 0; - // - for (Int_t sec=0; sec < 4; sec++){ - // - // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO - // - fCALOlength = 0; - memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer); - fSecPointer = fCALOlength; - // - // First of all we have section header and packet length - // - fDataCALO[fCALOlength] = fSecCalo[sec]; - fCALOlength++; - fDataCALO[fCALOlength] = fSecCALOLength[sec]; - fCALOlength++; - // - // Section XO is read in the opposite direction respect to the others - // - chksum = 0; - // - for (Int_t plane=0; plane < 11; plane++){ - // - if ( sec == 1 ) tstrip = fCALOlength + 96*2; - // - for (Int_t strip=0; strip < 96; strip++){ - // - chksum += (Int_t)fcalped[sec][plane][strip]; - // - // save value - // - if ( sec == 1 ){ - tstrip -= 2; - fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip]; - fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip]; - } else { - fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip]; - fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip]; - }; - fCALOlength +=2; - }; - // - }; - // - fDataCALO[fCALOlength] = (UShort_t)chksum; - fCALOlength++; - fDataCALO[fCALOlength] = 0; - fCALOlength++; - fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16)); - fCALOlength++; - // - // Section XO is read in the opposite direction respect to the others - // - chksum = 0; - // - for (Int_t plane=0; plane < 11; plane++){ - // - if ( sec == 1 ) tstrip = fCALOlength+6*2; - // - for (Int_t strip=0; strip < 6; strip++){ - // - chksum += (Int_t)fcalthr[sec][plane][strip]; - // - // save value - // - if ( sec == 1 ){ - tstrip -= 2; - fDataCALO[tstrip] = 0; - fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip]; - } else { - fDataCALO[fCALOlength] = 0; - fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip]; - }; - fCALOlength +=2; - }; - // - }; - // - fDataCALO[fCALOlength] = 0; - fCALOlength++; - fDataCALO[fCALOlength] = (UShort_t)chksum; - fCALOlength++; - fDataCALO[fCALOlength] = 0; - fCALOlength++; - fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16)); - fCALOlength++; - // - // Section XO is read in the opposite direction respect to the others - // - for (Int_t plane=0; plane < 11; plane++){ - // - if ( sec == 1 ) tstrip = fCALOlength+96*2; - // - for (Int_t strip=0; strip < 96; strip++){ - // - // save value - // - if ( sec == 1 ){ - tstrip -= 2; - fDataCALO[tstrip] = 0; - fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip]; - } else { - fDataCALO[fCALOlength] = 0; - fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip]; - }; - fCALOlength += 2; - }; - // - }; - // - // Section XO is read in the opposite direction respect to the others - // - for (Int_t plane=0; plane < 11; plane++){ - // - if ( sec == 1 ) tstrip = fCALOlength+6*4; - // - for (Int_t strip=0; strip < 6; strip++){ - // - // save value - // - if ( sec == 1 ){ - tstrip -= 4; - fDataCALO[tstrip] = 0; - fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip]; - fDataCALO[tstrip+2] = 0; - fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip]; - } else { - fDataCALO[fCALOlength] = 0; - fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip]; - fDataCALO[fCALOlength+2] = 0; - fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip]; - }; - fCALOlength +=4; - }; - // - }; - // - // - // here we calculate and save the CRC - // - fDataCALO[fCALOlength] = 0; - fCALOlength++; - Short_t CRC = 0; - for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){ - CRC=crc(CRC,fDataCALO[i+fSecPointer]); - }; - fDataCALO[fCALOlength] = (UShort_t)CRC; - fCALOlength++; - // - UInt_t length=fCALOlength*2; - DigitizePSCU(length,0x18); - // - // Add padding to 64 bits - // - AddPadding(); - // - fOutputfile.write(reinterpret_cast(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer); - UShort_t temp[1000000]; - memset(temp,0,sizeof(UShort_t)*1000000); - swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!! - fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fCALOlength); - // - // padding to 64 bytes - // - if ( fPadding ){ - fOutputfile.write(reinterpret_cast(fDataPadding),sizeof(UChar_t)*fPadding); - }; - // - // - }; - // -}; - -void Digitizer::CaloLoadCalib() { - // - fGivenCaloCalib = 0; // ####@@@@ should be given as input par @@@@#### - // - // first of all load the MIP to ADC conversion values - // - stringstream calfile; - Int_t error = 0; - GL_PARAM *glparam = new GL_PARAM(); - // - // determine where I can find calorimeter ADC to MIP conversion file - // - error = 0; - error = glparam->Query_GL_PARAM(3,101,fDbc); - // - calfile.str(""); - calfile << glparam->PATH.Data() << "/"; - calfile << glparam->NAME.Data(); - // - printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str()); - FILE *f; - f = fopen(calfile.str().c_str(),"rb"); - // - memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0])); - // - for (Int_t m = 0; m < 2 ; m++ ){ - for (Int_t k = 0; k < 22; k++ ){ - for (Int_t l = 0; l < 96; l++ ){ - fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f); - }; - }; - }; - fclose(f); - // - // determine which calibration has to be used and load it for each section - // - GL_CALO_CALIB *glcalo = new GL_CALO_CALIB(); - GL_ROOT *glroot = new GL_ROOT(); - TString fcalname; - UInt_t idcalib; - UInt_t calibno; - UInt_t utime = 0; - // - for (UInt_t s=0; s<4; s++){ - // - // clear calo calib variables for section s - // - ClearCaloCalib(s); - // - if ( fGivenCaloCalib ){ - // - // a time has been given as input on the command line so retrieve the calibration that preceed that time - // - glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc); - // - calibno = glcalo->EV_ROOT; - idcalib = glcalo->ID_ROOT_L0; - // - // determine path and name and entry of the calibration file - // - printf("\n"); - printf(" ** SECTION %i **\n",s); - // - glroot->Query_GL_ROOT(idcalib,fDbc); - // - stringstream name; - name.str(""); - name << glroot->PATH.Data() << "/"; - name << glroot->NAME.Data(); - // - fcalname = (TString)name.str().c_str(); - // - printf("\n Section %i : using file %s calibration at entry %i: \n",s,fcalname.Data(),calibno); - // - } else { - error = 0; - error = glparam->Query_GL_PARAM(1,104,fDbc); - // - calfile.str(""); - calfile << glparam->PATH.Data() << "/"; - calfile << glparam->NAME.Data(); - // - printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str()); - // - fcalname = (TString)calfile.str().c_str(); - calibno = s; - // - }; - // - // load calibration variables in memory - // - CaloLoadCalib(s,fcalname,calibno); - // - }; - // - // at this point we have in memory the calorimeter calibration and we can save it to disk in the correct format and use it to digitize the data - // - delete glparam; - delete glcalo; - delete glroot; -}; - -void Digitizer::DigitizeCALO() { - // - fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL ####@@@@ should be given as input par @@@@#### - // - // - // - fCALOlength = 0; // reset total dimension of calo data - // - // gpamela variables to be used - // - fhBookTree->SetBranchStatus("Nthcali",1); - fhBookTree->SetBranchStatus("Icaplane",1); - fhBookTree->SetBranchStatus("Icastrip",1); - fhBookTree->SetBranchStatus("Icamod",1); - fhBookTree->SetBranchStatus("Enestrip",1); - // - // call different routines depending on the acq mode you want to simulate - // - switch ( fModCalo ){ - case 0: - this->DigitizeCALORAW(); - break; - case 1: - this->DigitizeCALOCOMPRESS(); - break; - case 2: - this->DigitizeCALOFULL(); - break; - }; - // -}; - -Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){ - // - // determine plane and strip - // - Int_t mplane = 0; - // - // wrong! - // - // if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1; - // if ( sec == 1 ) mplane = (plane * 4) + 2 + 1; - // if ( sec == 2 ) mplane = (plane * 4) + 1 + 1; - // - if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1) from plane = 0, 11 - if ( sec == 1 ) mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1) from plane = 0, 11 - if ( sec == 2 ) mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1) from plane = 0, 11 - if ( sec == 3 ) mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1) from plane = 0, 11 - // - Int_t mstrip = strip + 1; - // - // search energy release in gpamela output - // - for (Int_t i=0; i XE 1 -> XO 2-> YE 3 -> YO - // - l = 0; // XE and XO are Y planes - if ( sec < 2 ) l = 1; // while YE and YO are X planes - // - fSecPointer = fCALOlength; - // - // First of all we have section header and packet length - // - fDataCALO[fCALOlength] = fSecCalo[sec]; - fCALOlength++; - fDataCALO[fCALOlength] = fSecCALOLength[sec]; - fCALOlength++; - // - // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers - // - for (Int_t autoplane=0; autoplane < 7; autoplane++){ - fDataCALO[fCALOlength] = 0x0000; - fCALOlength++; - }; - // - // - // here comes data - // - // - // Section XO is read in the opposite direction respect to the others - // - if ( sec == 1 ){ - tstrip = 96*11 + fCALOlength; - } else { - tstrip = 0; - }; - // - pre = -1; - // - for (Int_t strip=0; strip < 96; strip++){ - // - // which is the pre for this strip? - // - if (strip%16 == 0) { - pre++; - }; - // - if ( sec == 1 ) tstrip -= 11; - // - for (Int_t plane=0; plane < 11; plane++){ - // - // here is wrong!!!! - // - // - // if ( plane%2 == 0 && sec%2 != 0){ - // lpl = plane*2; - // } else { - // lpl = (plane*2) + 1; - // }; - // - if ( sec == 0 || sec == 3 ) lpl = plane * 2; - if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1; - // - // get the energy in GeV from the simulation for that strip - // - ens = this->GetCALOen(sec,plane,strip); - // - // convert it into ADC channels - // - adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio); - // - // sum baselines - // - adcbase = (UInt_t)fcalbase[sec][plane][pre]; - // - // add noise and pedestals - // - pedestal = fcalped[sec][plane][strip]; - rms = fcalrms[sec][plane][strip]/4.; - // - // Add random gaussian noise of RMS rms and Centered in the pedestal - // - pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms); - // - // Sum all contribution - // - adc = adcsig + adcbase + (Int_t)round(pedenoise); - // - // Signal saturation - // - if ( adc > 0x7FFF ) adc = 0x7FFF; - // - // save value - // - if ( sec == 1 ){ - fDataCALO[tstrip] = adc; - tstrip++; - } else { - fDataCALO[fCALOlength] = adc; - }; - fCALOlength++; - // - }; - // - if ( sec == 1 ) tstrip -= 11; - // - }; - // - // here we calculate and save the CRC - // - Short_t CRC = 0; - for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){ - CRC=crc(CRC,fDataCALO[i+fSecPointer]); - }; - fDataCALO[fCALOlength] = (UShort_t)CRC; - fCALOlength++; - // - }; - // - // for (Int_t i=0; iDigitizeCALORAW(); - return; - // - // - // - fSecCalo[0] = 0xEA00; - fSecCalo[1] = 0xF100; - fSecCalo[2] = 0xF600; - fSecCalo[3] = 0xED00; - // - // length of the data in DSP mode must be calculated on fly during digitization - // - memset(fSecCALOLength,0x0,4*sizeof(UShort_t)); - // - // here comes raw data - // - Int_t en = 0; - // - for (Int_t sec=0; sec < 4; sec++){ - fDataCALO[en] = fSecCalo[sec]; - en++; - fDataCALO[en] = fSecCALOLength[sec]; - en++; - for (Int_t plane=0; plane < 11; plane++){ - for (Int_t strip=0; strip < 11; strip++){ - fDataCALO[en] = 0x0; - en++; - }; - }; - }; - // -}; - -void Digitizer::DigitizeCALOFULL() { - // - printf(" FULL MODE STILL NOT IMPLEMENTED! \n"); - // - this->DigitizeCALORAW(); - return; - // - fSecCalo[0] = 0xEA00; - fSecCalo[1] = 0xF100; - fSecCalo[2] = 0xF600; - fSecCalo[3] = 0xED00; - // - // length of the data in DSP mode must be calculated on fly during digitization - // - memset(fSecCALOLength,0x0,4*sizeof(UShort_t)); - // - // here comes raw data - // - Int_t en = 0; - // - for (Int_t sec=0; sec < 4; sec++){ - fDataCALO[en] = fSecCalo[sec]; - en++; - fDataCALO[en] = fSecCALOLength[sec]; - en++; - for (Int_t plane=0; plane < 11; plane++){ - for (Int_t strip=0; strip < 11; strip++){ - fDataCALO[en] = 0x0; - en++; - }; - }; - }; - // -}; - -void Digitizer::DigitizeTRIGGER() { - //fDataTrigger: 153 bytes - for (Int_t j=0; j < 153; j++) - fDataTrigger[0]=0x00; -}; - -Int_t Digitizer::DigitizeTOF() { - //fDataTof: 12 x 23 bytes (=276 bytes) - UChar_t *pTof=fDataTof; - Bool_t DEBUG=false; - - // --- activate branches: - fhBookTree->SetBranchStatus("Nthtof",1); - fhBookTree->SetBranchStatus("Ipltof",1); - fhBookTree->SetBranchStatus("Ipaddle",1); - fhBookTree->SetBranchStatus("Xintof",1); - fhBookTree->SetBranchStatus("Yintof",1); - fhBookTree->SetBranchStatus("Xouttof",1); - fhBookTree->SetBranchStatus("Youttof",1); - fhBookTree->SetBranchStatus("Ereltof",1); - fhBookTree->SetBranchStatus("Timetof",1); - // not yet used: Zintof, Xouttof, Youttof, Zouttof - - // ------ evaluate energy in each pmt: ------ - // strip geometry (lenght/width) - Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0}; - //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0}; - - // S11 8 paddles 33.0 x 5.1 cm - // S12 6 paddles 40.8 x 5.5 cm - // S21 2 paddles 18.0 x 7.5 cm - // S22 2 paddles 15.0 x 9.0 cm - // S31 3 paddles 15.0 x 6.0 cm - // S32 3 paddles 18.0 x 5.0 cm - - // distance from the interaction point to the pmts (right,left) - Float_t xpath[2]={0., 0.}; /*path(cm) in X per S12,S21,S32 verso il pmt DX o SX*/ - Float_t ypath[2]={0., 0.}; /*path(cm) in Y per S11,S22,S31 verso il pmt DX o SX*/ - Float_t FGeo[2]={0., 0.}; /* fattore geometrico */ - - const Float_t Pho_keV = 10.; // photons per keV in scintillator - const Float_t echarge = 1.6e-19; // carica dell'elettrone - Float_t Npho=0.; - Float_t QevePmt_pC[48]; - Float_t QhitPad_pC[2]={0., 0.}; - Float_t QhitPmt_pC[2]={0., 0.}; - Float_t pmGain = 3.5e6; /* Gain: per il momento uguale per tutti */ - Float_t effi=0.21; /* Efficienza di fotocatodo */ - Float_t ADC_pC=1.666667; // ADC_ch/pC conversion = 0.6 pC/channel (+30 di offset) - Float_t ADCoffset=30.; - Int_t ADClast=4095; // no signal --> ADC ch=4095 - Int_t ADCtof[48]; - //Float_t ADCsat=3100; ci pensiamo in futuro ! - //Float_t pCsat=2500; - for(Int_t i=0; i<48; i++){ - QevePmt_pC[i] = 0; - ADCtof[i]=0; - } - - // ------ read calibration file (get A1, A2, lambda1, lambda2) - ifstream fileTriggerCalib; - TString ftrigname="TrigCalibParam.txt"; - fileTriggerCalib.open(ftrigname.Data()); - if ( !fileTriggerCalib ) { - printf("debug: no trigger calib file!\n"); - return(-117); //check output! - }; - Float_t atte1[48],atte2[48],lambda1[48],lambda2[48]; - Int_t temp=0; - // correct readout WM Oct '07 - for(Int_t i=0; i<48; i++){ - fileTriggerCalib >> temp; - fileTriggerCalib >> atte1[i]; - fileTriggerCalib >> lambda1[i]; - fileTriggerCalib >> atte2[i]; - fileTriggerCalib >> lambda2[i]; - fileTriggerCalib >> temp; - } - fileTriggerCalib.close(); - - // Read from file the 48*4 values of the attenuation fit function - // NB: lambda<0; x,y defined in gpamela (=0 in the centre of the cavity) - // Qhitpmt_pC = atte1 * exp(x/lambda1) + atte2 * exp(x/lambda2) - - // fine lettura dal file */ - - //const Int_t nmax=??; = Nthtof - Int_t ip, ipad; - //Int_t ipmt; - Int_t pmtleft=0, pmtright=0; - Int_t *pl, *pr; - pl = &pmtleft; - pr = &pmtright; - - // TDC variables: - Int_t TDClast=4095; // no signal --> ADC ch=4095 - Int_t TDCint[48]; - Float_t tdc[48],tdc1[48],tdcpmt[48]; - for(Int_t i=0; i<48; i++) - tdcpmt[i] = 1000.; - Float_t thresh=10.; // to be defined better... (Wolfgang) - - // === TDC: simulate timing for each paddle - Float_t dt1 = 285.e-12 ; // single PMT resolution - Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50]; - for(Int_t j=0;j<48;j++) tdcres[j] = 50.E-12; // TDC resolution 50 picosec - for(Int_t j=0;j<48;j++) c1_S[j] = 500.; // cable length in channels - for(Int_t j=0;j<48;j++) c2_S[j] = 0.; - for(Int_t j=0;j<48;j++) c3_S[j] = 1000.; - for(Int_t j=0;j<48;j++) c1_S[j] = c1_S[j]*tdcres[j]; // cable length in sec - for(Int_t j=0;j<48;j++) c2_S[j] = c2_S[j]*tdcres[j]; - // ih = 0 + i1; // not used?? (Silvio) - - /* ********************************** start loop over hits */ - - for(Int_t nh=0; nh half board e canale - // metodo GetPMTIndex(Int_t ipmt, Int_t &hb, Int_t &ch) // non lo usiamo x ora - - // evaluates mean position and path inside the paddle - - Float_t tpos=0.; - Float_t path[2] = {0., 0.}; - //--- Strip in Y = S11,S22,S31 ------ - if(ip==0 || ip==3 || ip==4) - tpos = (Yintof[nh]+Youttof[nh])/2.; - else - if(ip==1 || ip==2 || ip==5) //--- Strip in X per S12,S21,S32 - tpos = (Xintof[nh]+Xouttof[nh])/2.; - else //if (ip!=6) - printf("*** WARNING TOF: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh); - path[0]= tpos + dimel[ip]/2.; - path[1]= dimel[ip]/2.- tpos; - - // cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n"; - - if (DEBUG) { - cout <<" plane "< thresh) { - if (tdcpmt[pmtleft] == 1000.) { // fill for the first time - tdcpmt[pmtleft] = t1; - tdc[pmtleft] = t1 + c2_S[pmtleft] ; // Signal reaches Coincidence - } - if (tdcpmt[pmtleft] < 1000.) // is already filled! - if (t1 < tdcpmt[pmtleft]) { - tdcpmt[pmtleft] = t1; - t1 = t1 + c2_S[pmtleft] ; // Signal reaches Coincidence - tdc[pmtleft] = t1; - } - } - if (QhitPmt_pC[1] > thresh) { - if (tdcpmt[pmtright] == 1000.) { // fill for the first time - tdcpmt[pmtright] = t2; - tdc[pmtright] = t2 + c2_S[pmtright] ; // Signal reaches Coincidence - } - if (tdcpmt[pmtright] < 1000.) // is already filled! - if (t2 < tdcpmt[pmtright]) { - tdcpmt[pmtright] = t2; - t2 = t2 + c2_S[pmtright] ; - tdc[pmtright] = t2; - } - } - - if (DEBUG) - cout< ADClast) ADCtof[i]=ADClast; - } else - ADCtof[i]= ADClast; - } - - - // ====== build TDC coincidence ====== - - Float_t t_coinc = 0; - Int_t ilast = 100; - for (Int_t ii=0; ii<48;ii++) - if (tdc[ii] > t_coinc) { - t_coinc = tdc[ii]; - ilast = ii; - } - - // cout< ADClast) ADC[i]=ADClast; - } else - TDCint[i]= TDClast; - } - - if (DEBUG) - cout<<"-----------"<SetBranchAddress("Ievnt",&Ievnt); - fhBookTree->SetBranchStatus("Nthcat",1); - fhBookTree->SetBranchStatus("Iparcat",1); - fhBookTree->SetBranchStatus("Icat",1); - fhBookTree->SetBranchStatus("Xincat",1); - fhBookTree->SetBranchStatus("Yincat",1); - fhBookTree->SetBranchStatus("Zincat",1); - fhBookTree->SetBranchStatus("Xoutcat",1); - fhBookTree->SetBranchStatus("Youtcat",1); - fhBookTree->SetBranchStatus("Zoutcat",1); - fhBookTree->SetBranchStatus("Erelcat",1); - fhBookTree->SetBranchStatus("Timecat",1); - fhBookTree->SetBranchStatus("Pathcat",1); - fhBookTree->SetBranchStatus("P0cat",1); - fhBookTree->SetBranchStatus("Nthcas",1); - fhBookTree->SetBranchStatus("Iparcas",1); - fhBookTree->SetBranchStatus("Icas",1); - fhBookTree->SetBranchStatus("Xincas",1); - fhBookTree->SetBranchStatus("Yincas",1); - fhBookTree->SetBranchStatus("Zincas",1); - fhBookTree->SetBranchStatus("Xoutcas",1); - fhBookTree->SetBranchStatus("Youtcas",1); - fhBookTree->SetBranchStatus("Zoutcas",1); - fhBookTree->SetBranchStatus("Erelcas",1); - fhBookTree->SetBranchStatus("Timecas",1); - fhBookTree->SetBranchStatus("Pathcas",1); - fhBookTree->SetBranchStatus("P0cas",1); - fhBookTree->SetBranchStatus("Nthcard",1); - fhBookTree->SetBranchStatus("Iparcard",1); - fhBookTree->SetBranchStatus("Icard",1); - fhBookTree->SetBranchStatus("Xincard",1); - fhBookTree->SetBranchStatus("Yincard",1); - fhBookTree->SetBranchStatus("Zincard",1); - fhBookTree->SetBranchStatus("Xoutcard",1); - fhBookTree->SetBranchStatus("Youtcard",1); - fhBookTree->SetBranchStatus("Zoutcard",1); - fhBookTree->SetBranchStatus("Erelcard",1); - fhBookTree->SetBranchStatus("Timecard",1); - fhBookTree->SetBranchStatus("Pathcard",1); - fhBookTree->SetBranchStatus("P0card",1); - - // In this simpliefied approach we will assume that once - // a particle releases > 0.5 mip in one of the 12 AC detectors it - // will fire. We will furthermore assume that both cards read out - // identical data. - - // If you develop your digitization algorithm, you should start by - // identifying the information present in level2 (post-darth-vader) - // data. - - Float_t SumEcat[5]; - Float_t SumEcas[5]; - Float_t SumEcard[5]; - for (Int_t k= 0;k<5;k++){ - SumEcat[k]=0.; - SumEcas[k]=0.; - SumEcard[k]=0.; - }; - - if (Nthcat>50 || Nthcas>50 || Nthcard>50) - printf("*** ERROR AC! NthAC out of range!\n\n"); - - // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi) - // based on J.Lundquist's calculations (PhD thesis, page 94) - // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are: - // 8.25470e-01 +- 1.79489e-02 - // 6.41609e-01 +- 2.65846e-02 - // 9.81177e+00 +- 1.21284e+00 - // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip - // - // NB: the PMT positions are needed! - - // look in CAT - // for (UInt_t k= 0;k<50;k++){ - for (Int_t k= 0;k 0) - SumEcat[Icat[k]] += Erelcat[k]; - }; - - // look in CAS - for (Int_t k= 0;k0) - SumEcas[Icas[k]] += Erelcas[k]; - }; - - // look in CARD - for (Int_t k= 0;k0) - SumEcard[Icard[k]] += Erelcard[k]; - }; - - // channel mapping Hit Map - // 1 CARD4 0 LSB - // 2 CAT2 0 - // 3 CAS1 0 - // 4 NC 0 - // 5 CARD2 0 - // 6 CAT4 1 - // 7 CAS4 0 - // 8 NC 0 - // 9 CARD3 0 - // 10 CAT3 0 - // 11 CAS3 0 - // 12 NC 0 - // 13 CARD1 0 - // 14 CAT1 0 - // 15 CAS2 0 - // 16 NC 0 MSB - - // In the first version only the hit-map is filled, not the SR. - - // Threshold: 0.8 MeV. - - Float_t thr = 8e-4; - - fDataAC[3] = 0x0000; - - if (SumEcas[0] > thr) fDataAC[3] = 0x0004; - if (SumEcas[1] > thr) fDataAC[3] += 0x4000; - if (SumEcas[2] > thr) fDataAC[3] += 0x0400; - if (SumEcas[3] > thr) fDataAC[3] += 0x0040; - - if (SumEcat[0] > thr) fDataAC[3] += 0x2000; - if (SumEcat[1] > thr) fDataAC[3] += 0x0002; - if (SumEcat[2] > thr) fDataAC[3] += 0x0200; - if (SumEcat[3] > thr) fDataAC[3] += 0x0020; - - if (SumEcard[0] > thr) fDataAC[3] += 0x1000; - if (SumEcard[1] > thr) fDataAC[3] += 0x0010; - if (SumEcard[2] > thr) fDataAC[3] += 0x0100; - if (SumEcard[3] > thr) fDataAC[3] += 0x0001; - - fDataAC[67] = fDataAC[3]; - - // shift registers - // the central bin is equal to the hitmap, all other bins in the shift register are 0 - for (UInt_t i=0; i<=15; i++){ - fDataAC[i+11] = 0x0000; - fDataAC[i+75] = 0x0000; - } - fDataAC[18] = fDataAC[3]; - fDataAC[82] = fDataAC[3]; - - // for (Int_t i=0; iSetBranchStatus("Nthtof",1); - fhBookTree->SetBranchStatus("Ipltof",1); - fhBookTree->SetBranchStatus("Ipaddle",1); - - fhBookTree->SetBranchStatus("Xintof",1); - fhBookTree->SetBranchStatus("Yintof",1); - fhBookTree->SetBranchStatus("Xouttof",1); - fhBookTree->SetBranchStatus("Youttof",1); - - fhBookTree->SetBranchStatus("Ereltof",1); - fhBookTree->SetBranchStatus("Timetof",1); - NdFT=0; - Ert=0; - for(i=0;iXs[0] && xYs[0] && yZs[0] && zXs[0] && xYs[0] && yZs[0] && z>ciao; - } - if((xXs[1]-p)&&(y>Ys[0]+p && yZs[0]+p && z=Yp[2*t] && yXs[0]+p && xYs[1]-p)&&(z>Zs[0]+p && zXs[0]+p && xYs[0]+p && yZs[1]-p))V[2]=-V[2]; - x+=V[0]; - y+=V[1]; - z+=V[2]; - l=0; - //cout<>ciao; - } - } - } - Ert=Ert/0.002; - q=(Float_t)(random())/(Float_t)0x7fffffff; - w=0.7; - //E0=(Float_t)(4064./7.); - E0=4064./7.; - if(Ert<1) S4=0; - else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0))); - i=S4/4; - if(S4%4==0) - S4v=S4+S4p; - else if(S4%4==1){ - if(q>ciao; -} - - - -void Digitizer::DigitizeND(){ - // creato: S. Borisov, INFN Roma2 e MEPHI, Sett 2007 - Int_t i=0; - UShort_t NdN=0; - fhBookTree->SetBranchStatus("Nthnd",1); - fhBookTree->SetBranchStatus("Itubend",1); - fhBookTree->SetBranchStatus("Iparnd",1); - fhBookTree->SetBranchStatus("Xinnd",1); - fhBookTree->SetBranchStatus("Yinnd",1); - fhBookTree->SetBranchStatus("Zinnd",1); - fhBookTree->SetBranchStatus("Xoutnd",1); - fhBookTree->SetBranchStatus("Youtnd",1); - fhBookTree->SetBranchStatus("Zoutnd",1); - fhBookTree->SetBranchStatus("Erelnd",1); - fhBookTree->SetBranchStatus("Timend",1); - fhBookTree->SetBranchStatus("Pathnd",1); - fhBookTree->SetBranchStatus("P0nd",1); - //cout<<"n="<SetBranchStatus("Enestrip",1); - - // dumy header - fDataDummy[0] = 0xCAAA; - - for (Int_t i=1; i(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer); - // TRG - fOutputfile.write(reinterpret_cast(fDataTrigger),sizeof(UChar_t)*153); - // TOF - fOutputfile.write(reinterpret_cast(fDataTof),sizeof(UChar_t)*276); - // AC - UShort_t temp[1000000]; - memset(temp,0,sizeof(UShort_t)*1000000); - swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer); // WE MUST SWAP THE BYTES!!! - fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fACbuffer); - // CALO - memset(temp,0,sizeof(UShort_t)*1000000); - swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!! - fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fCALOlength); - // TRK - memset(temp,0,sizeof(UShort_t)*1000000); - swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength); // WE MUST SWAP THE BYTES!!! - fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fTracklength); - fTracklength=0; - // padding to 64 bytes - // - if ( fPadding ){ - fOutputfile.write(reinterpret_cast(fDataPadding),sizeof(UChar_t)*fPadding); - }; - // S4 - memset(temp,0,sizeof(UShort_t)*1000000); - swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer); // WE MUST SWAP THE BYTES!!! - fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fS4buffer); - // ND - memset(temp,0,sizeof(UShort_t)*1000000); - swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer); // WE MUST SWAP THE BYTES!!! - fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fNDbuffer); -}; - - -void Digitizer::ReadData(){ - - UShort_t InData[64]; - - // for debuggigng purposes only, write your own routine if you like (many - // hardwired things. - - ifstream InputFile; - - // if (!InputFile) { - - // std::cout << "ERROR" << endl; - // // An error occurred! - // // myFile.gcount() returns the number of bytes read. - // // calling myFile.clear() will reset the stream state - // // so it is usable again. - // }; - - - - //InputFile.seekg(0); - - InputFile.open(fFilename, ios::in | ios::binary); - // fOutputfile.seekg(0); - if (!InputFile.is_open()) std::cout << "ERROR" << endl; - - InputFile.seekg(0); - - for (Int_t k=0; k<=1000; k++){ - InputFile.read(reinterpret_cast(InData),384*sizeof(UShort_t)); - - std::cout << "Read back: " << endl << endl; - - for (Int_t i=0; i<=384; i++){ - printf("%4x ", InData[i]); - if ((i+1)%8 ==0) cout << endl; - } - - } - cout << endl; - InputFile.close(); - -}; - - - -void Digitizer::DigitizeTrack() { -//std:: cout << "Entering DigitizeTrack " << endl; -Float_t AdcTrack[fNviews][fNstrips_view]; // Vector of strips to be compressed - -Int_t Iview; -Int_t Nstrip; - - for (Int_t j=0; jGaus(0.,fSigmaCommon); - Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon); - for (Int_t k=0; kGaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]); - if(k<4*128) {AdcTrack[j][Nstrip] += commonN1;} // full correlation of 4 VA1 Com. Noise - else {AdcTrack[j][Nstrip] += commonN2;} // full correlation of 4 VA1 Com. Noise - - }; - - - }; - - - }; - - - fhBookTree->SetBranchStatus("Nstrpx",1); - fhBookTree->SetBranchStatus("Npstripx",1); - fhBookTree->SetBranchStatus("Ntstripx",1); - fhBookTree->SetBranchStatus("Istripx",1); - fhBookTree->SetBranchStatus("Qstripx",1); - fhBookTree->SetBranchStatus("Xstripx",1); - fhBookTree->SetBranchStatus("Nstrpy",1); - fhBookTree->SetBranchStatus("Npstripy",1); - fhBookTree->SetBranchStatus("Ntstripy",1); - fhBookTree->SetBranchStatus("Istripy",1); - fhBookTree->SetBranchStatus("Qstripy",1); - fhBookTree->SetBranchStatus("Ystripy",1); - - - - Float_t ADCfull; - Int_t iladd=0; - for (Int_t ix=0; ix=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1; - if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2; - ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview]; - AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull); - - }; - - - for (Int_t iy=0; iy=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1; - if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2; - ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview]; - AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull); - - }; - -CompressTrackData(AdcTrack); // Compress and Digitize data of one Ladder in turn for all ladders - -}; - - - -void Digitizer::DigitizeTrackCalib(Int_t ii) { - -std:: cout << "Entering DigitizeTrackCalib " << ii << endl; -if( (ii!=1)&&(ii!=2) ) { - std:: cout << "error wrong DigitizeTrackCalib argument" << endl; - return; -}; - -memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); -fTracklength=0; - -UShort_t Dato; - -Float_t dato1; -Float_t dato2; -Float_t dato3; -Float_t dato4; - -UShort_t DatoDec; -UShort_t DatoDec1; -UShort_t DatoDec2; -UShort_t DatoDec3; -UShort_t DatoDec4; - -UShort_t EVENT_CAL; -UShort_t PED_L1; -UShort_t ReLength; -UShort_t OveCheckCode; -//UShort_t PED_L2; -//UShort_t PED_L3HI; -//UShort_t PED_L3LO; -//UShort_t SIG_L1HI; -//UShort_t SIG_L1LO; -//UShort_t SIG_L2HI; -//UShort_t SIG_L2LO; -//UShort_t SIG_L3; -//UShort_t BAD_L1; -//UShort_t BAD_L2LO; -//UShort_t BAD_L3HI; -//UShort_t BAD_L3LO; -//UShort_t FLAG; - - - Int_t DSPpos; - for (Int_t j=ii-1; j> 8); - fTracklength++; - - fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) ); - fTracklength++; - - // end TRAILER - }; - - // write in buffer the DSP header - - fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) ); - - fDataTrack[DSPpos+1]=0x01A9; - - fDataTrack[DSPpos+2]=0x8740; - - EVENT_CAL=0; - fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) ); - - PED_L1=0; - fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) ); - - // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0 - - fDataTrack[DSPpos+5]=0x8014; - - fDataTrack[DSPpos+6]=0x00A0; - - fDataTrack[DSPpos+7]=0x0500; - - fDataTrack[DSPpos+8]=0x2801; - - fDataTrack[DSPpos+9]=0x400A; - - fDataTrack[DSPpos+10]=0x0050; - - CkSum=(CkSum >> 8)^(CkSum&0x00FF); - fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3)); - - fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) ); - - // end dsp header - - // write in buffer the TRAILER - - ReLength=(UShort_t)((13*2)+3); - OveCheckCode=0x0000; - fDataTrack[DSPpos+13]=0x0000; - - fDataTrack[DSPpos+14]=(ReLength >> 8); - - fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) ); - - // end TRAILER - - - - - // end DSP - }; - - - -}; - -void Digitizer::WriteTrackCalib() { - - -std:: cout << " Entering WriteTrackCalib " << endl; - -fOutputfile.write(reinterpret_cast(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer); - -UShort_t temp[1000000]; -memset(temp,0,sizeof(UShort_t)*1000000); -swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength); // WE MUST SWAP THE BYTES!!! -fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fTracklength); -fTracklength=0; -if ( fPadding ){ - fOutputfile.write(reinterpret_cast(fDataPadding),sizeof(UChar_t)*fPadding); -}; - -}; - - -void Digitizer::ClearTrackCalib() { - -std:: cout << "Entering ClearTrackCalib " << endl; - - -}; - - -void Digitizer::LoadTrackCalib() { -std:: cout << "Entering LoadTrackCalib " << endl; - -// Generate the pedestals and sigmas according to parametrization - for (Int_t j=0; jGaus(fAvePedex,fSigmaPedex); - fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax); - }; - if((j+1)%2==1) { - fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey); - fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay); - }; - - }; - }; - - - -}; - -void Digitizer::LoadMipCor() { -std:: cout << "Entering LoadMipCor" << endl; - Float_t xfactor=1./151.6*1.04; - Float_t yfactor=1./152.1; - - fMipCor[0][0]=140.02*yfactor; - fMipCor[0][1]=140.99*xfactor; - fMipCor[0][2]=134.48*yfactor; - fMipCor[0][3]=144.41*xfactor; - fMipCor[0][4]=140.74*yfactor; - fMipCor[0][5]=142.28*xfactor; - fMipCor[0][6]=134.53*yfactor; - fMipCor[0][7]=140.63*xfactor; - fMipCor[0][8]=135.55*yfactor; - fMipCor[0][9]=138.00*xfactor; - fMipCor[0][10]=154.95*yfactor; - fMipCor[0][11]=158.44*xfactor; - - - fMipCor[1][0]=136.07*yfactor; - fMipCor[1][1]=135.59*xfactor; - fMipCor[1][2]=142.69*yfactor; - fMipCor[1][3]=138.19*xfactor; - fMipCor[1][4]=137.35*yfactor; - fMipCor[1][5]=140.23*xfactor; - fMipCor[1][6]=153.15*yfactor; - fMipCor[1][7]=151.42*xfactor; - fMipCor[1][8]=129.76*yfactor; - fMipCor[1][9]=140.63*xfactor; - fMipCor[1][10]=157.87*yfactor; - fMipCor[1][11]=153.64*xfactor; - - fMipCor[2][0]=134.98*yfactor; - fMipCor[2][1]=143.95*xfactor; - fMipCor[2][2]=140.23*yfactor; - fMipCor[2][3]=138.88*xfactor; - fMipCor[2][4]=137.95*yfactor; - fMipCor[2][5]=134.87*xfactor; - fMipCor[2][6]=157.56*yfactor; - fMipCor[2][7]=157.31*xfactor; - fMipCor[2][8]=141.37*yfactor; - fMipCor[2][9]=143.39*xfactor; - fMipCor[2][10]=156.15*yfactor; - fMipCor[2][11]=158.79*xfactor; - -/* - for (Int_t j=0; j 0.5) inte=inte+1; - newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k]; - // first strip of ladder is transmitted - // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl; - DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); - DSPlength++; - ntrastot++; - trasmesso=1; - oldval=newval; - kt=k; - k++; - continue; - }; - real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte); - if (real > 0.5) inte=inte+1; - newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]); - cercacluster=1; // ????????? - if (cercacluster==1) { - - // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!! - Int_t diff=0; - - - switch ((iv+1)%2) { - case 0: diff=newval-oldval; - break; - case 1: diff=oldval-newval; - break; - }; - - if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) { - Int_t clval=newval; - Int_t klp=k; // go on to search for maximum - klp++; - - while(klp 0.5) inte=inte+1; - Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp]; - if((iv+1)%2==0) { - - if(clvalp>clval) { - clval=clvalp; - k=klp;} - else break; // max of cluster found - - } else { - - if(clvalp=kl1) kl1=kt+1; - if( (kt+1)==kl1 ) trasmesso=1; - - - - Int_t kl2=k+fNclst; - if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1; - - for(Int_t klt=kl1 ; klt<=kl2 ; klt++) { - if(trasmesso==0) { - // std:: cout << "STRIP " << klt << endl; - // std:: cout << "ADC_TOT " < 0.5) inte=inte+1; - DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); - DSPlength++; - ntrastot++; - - } - else { - // std:: cout << "ADC_TOT " < 0.5) inte=inte+1; - DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); - DSPlength++; - ntrastot++; - }; - trasmesso=1; - }; // end trasmission - kt=kl2; - k=kl2; - real=modff(AdcTrack[iv][ladder*fNstrips_ladder+kt],&inte); - if (real > 0.5) inte=inte+1; - oldval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+kt]; - k++; - continue; - - - }; // end cercacluster - }; // end cercacluster - -// start ZOP check for strips no - - if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) { - - if(trasmesso==0) { - // std:: cout << "STRIP " << k << endl; - // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl; - - DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000); - DSPlength++; - ntrastot++; - - - real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte); - if (real > 0.5) inte=inte+1; - DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); - DSPlength++; - ntrastot++; - - } - else { - // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl; - real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte); - if (real > 0.5) inte=inte+1; - DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); - DSPlength++; - ntrastot++; - }; - trasmesso=1; - oldval=newval; - kt=k; - - } - else trasmesso=0; - // end zop - - k++; - }; // end cycle inside ladder -// write here the end ladder bytes -// std:: cout << "FINE LADDER " << ladder+1 << endl; - - DataDSP[DSPlength]=( ((UShort_t)(ladder+1)) | 0x1800); - DSPlength++; - ntrastot++; - trasmesso=0; - - }; //end cycle inside dsp -// std:: cout << "FINE DSP " << iv+1 << endl; -// here put DSP header - DataDSP[0]=(0x1CA0 | ((UShort_t)(iv+1)) ); - UShort_t Nword=(DSPlength*13)/16; - if( ((DSPlength*13)%16)!=0) Nword++; - DataDSP[1]=(0x1400 | ( Nword >> 10)); - DataDSP[2]=(0x1400 | ( Nword & 0x03FF) ); - DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) ); - DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) ) & 0x03FF) ); - DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) ) - | ( (UShort_t)fCutzop ) ); - DataDSP[6]=0x1400; - DataDSP[7]=0x1400; - DataDSP[8]=0x1400; - DataDSP[9]=0x1400; - DataDSP[10]=0x1400; - DataDSP[11]=0x1400; - DataDSP[12]=0x1400; - DataDSP[13]=0x1400; - DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) ); - DataDSP[15]=0x1C00; -// end DSP header - - -// write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer - Int_t Bit16free=16; - UShort_t Dato; - for (Int_t NDSP=0; NDSP0) { - if(Bit13ToWrite<=Bit16free) { - Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite)); - fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ; - Bit16free=Bit16free-Bit13ToWrite; - Bit13ToWrite=0; - if(Bit16free==0) { - if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength]; - fTracklength++; - Bit16free=16; - }; - } - else if(Bit13ToWrite>Bit16free) { - Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) ); - fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ; - if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength]; - fTracklength++; - Bit13ToWrite=Bit13ToWrite-Bit16free; - Bit16free=16; - }; - - }; // end cycle while(Bit13ToWrite>0) - - }; // end cycle DataDSP - if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; }; - CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF); - fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3)); - fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) ); - -// end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer - -//write trailer on buffer - UShort_t ReLength=(UShort_t)((Nword+13)*2+3); - UShort_t OveCheckCode=0x0000; - - fDataTrack[fTracklength]=0x0000; - fTracklength++; - - fDataTrack[fTracklength]=(ReLength >> 8); - fTracklength++; - - fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) ); - fTracklength++; -// end trailer -// std:: cout << "DSPlength " <3000.) { SatFact=3000./ADC; }; - return SatFact; -}; - - - - - - +// ------ PAMELA Digitizer ------ +// +// Date, release and how-to: see file Pamelagp2Digits.cxx +// +// NB: Check length physics packet [packet type (0x10 = physics data)] +// +#include +#include +#include +#include +#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 "TMath.h" +#include "TRandom.h" +#include "TSQLServer.h" +#include "TSystem.h" +// +#include "Digitizer.h" +#include "CRC.h" +// +#include +#include +#include +#include "GLTables.h" +// +extern "C"{ + short crc(short, short); +}; +// + +Digitizer::Digitizer(TTree* tree, char* &file_raw){ + fhBookTree = tree; + fFilename = file_raw; + fCounter = 0; + fOBT = 0; + + // + // DB connections + // + TString host = "mysql://localhost/pamelaprod"; + TString user = "anonymous"; + TString psw = ""; + // + const char *pamdbhost=gSystem->Getenv("PAM_DBHOST"); + const char *pamdbuser=gSystem->Getenv("PAM_DBUSER"); + const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW"); + if ( !pamdbhost ) pamdbhost = ""; + if ( !pamdbuser ) pamdbuser = ""; + if ( !pamdbpsw ) pamdbpsw = ""; + if ( strcmp(pamdbhost,"") ) host = pamdbhost; + if ( strcmp(pamdbuser,"") ) user = pamdbuser; + if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw; + fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data()); + // + GL_TABLES *glt = new GL_TABLES(host,user,psw); + if ( glt->IsConnected(fDbc) ) printf("\n DB INFORMATION:\n SQL: %s Version: %s Host %s Port %i \n\n",fDbc->GetDBMS(),fDbc->ServerInfo(),fDbc->GetHost(),fDbc->GetPort()); + // + // Use UTC in the DB and make timeout bigger + // + stringstream myquery; + myquery.str(""); + myquery << "SET time_zone='+0:00'"; + fDbc->Query(myquery.str().c_str()); + myquery.str(""); + myquery << "SET wait_timeout=173000;"; + fDbc->Query(myquery.str().c_str()); + // + + std:: cout << "preparing tree" << endl; + + // prepare tree + fhBookTree->SetBranchAddress("Irun",&Irun); + fhBookTree->SetBranchAddress("Ievnt",&Ievnt); + fhBookTree->SetBranchAddress("Ipa",&Ipa); + fhBookTree->SetBranchAddress("X0",&X0); + fhBookTree->SetBranchAddress("Y0",&Y0); + fhBookTree->SetBranchAddress("Z0",&Z0); + fhBookTree->SetBranchAddress("Theta",&Theta); + fhBookTree->SetBranchAddress("Phi",&Phi); + fhBookTree->SetBranchAddress("P0",&P0); + fhBookTree->SetBranchAddress("Nthtof",&Nthtof); + fhBookTree->SetBranchAddress("Ipltof",Ipltof); + fhBookTree->SetBranchAddress("Ipaddle",Ipaddle); + fhBookTree->SetBranchAddress("Ipartof",Ipartof); + fhBookTree->SetBranchAddress("Xintof",Xintof); + fhBookTree->SetBranchAddress("Yintof",Yintof); + fhBookTree->SetBranchAddress("Zintof",Zintof); + fhBookTree->SetBranchAddress("Xouttof",Xouttof); + fhBookTree->SetBranchAddress("Youttof",Youttof); + fhBookTree->SetBranchAddress("Zouttof",Zouttof); + fhBookTree->SetBranchAddress("Ereltof",Ereltof); + fhBookTree->SetBranchAddress("Timetof",Timetof); + fhBookTree->SetBranchAddress("Pathtof",Pathtof); + fhBookTree->SetBranchAddress("P0tof",P0tof); + fhBookTree->SetBranchAddress("Nthcat",&Nthcat); + fhBookTree->SetBranchAddress("Iparcat",Iparcat); + fhBookTree->SetBranchAddress("Icat",Icat); + fhBookTree->SetBranchAddress("Xincat",Xincat); + fhBookTree->SetBranchAddress("Yincat",Yincat); + fhBookTree->SetBranchAddress("Zincat",Zincat); + fhBookTree->SetBranchAddress("Xoutcat",Xoutcat); + fhBookTree->SetBranchAddress("Youtcat",Youtcat); + fhBookTree->SetBranchAddress("Zoutcat",Zoutcat); + fhBookTree->SetBranchAddress("Erelcat",Erelcat); + fhBookTree->SetBranchAddress("Timecat",Timecat); + fhBookTree->SetBranchAddress("Pathcat",Pathcat); + fhBookTree->SetBranchAddress("P0cat",P0cat); + fhBookTree->SetBranchAddress("Nthcas",&Nthcas); + fhBookTree->SetBranchAddress("Iparcas",Iparcas); + fhBookTree->SetBranchAddress("Icas",Icas); + fhBookTree->SetBranchAddress("Xincas",Xincas); + fhBookTree->SetBranchAddress("Yincas",Yincas); + fhBookTree->SetBranchAddress("Zincas",Zincas); + fhBookTree->SetBranchAddress("Xoutcas",Xoutcas); + fhBookTree->SetBranchAddress("Youtcas",Youtcas); + fhBookTree->SetBranchAddress("Zoutcas",Zoutcas); + fhBookTree->SetBranchAddress("Erelcas",Erelcas); + fhBookTree->SetBranchAddress("Timecas",Timecas); + fhBookTree->SetBranchAddress("Pathcas",Pathcas); + fhBookTree->SetBranchAddress("P0cas",P0cas); + fhBookTree->SetBranchAddress("Nthspe",&Nthspe); + fhBookTree->SetBranchAddress("Iparspe",Iparspe); + fhBookTree->SetBranchAddress("Itrpb",Itrpb); + fhBookTree->SetBranchAddress("Itrsl",Itrsl); + fhBookTree->SetBranchAddress("Itspa",Itspa); + fhBookTree->SetBranchAddress("Xinspe",Xinspe); + fhBookTree->SetBranchAddress("Yinspe",Yinspe); + fhBookTree->SetBranchAddress("Zinspe",Zinspe); + fhBookTree->SetBranchAddress("Xoutspe",Xoutspe); + fhBookTree->SetBranchAddress("Youtspe",Youtspe); + fhBookTree->SetBranchAddress("Zoutspe",Zoutspe); + fhBookTree->SetBranchAddress("Xavspe",Xavspe); + fhBookTree->SetBranchAddress("Yavspe",Yavspe); + fhBookTree->SetBranchAddress("Zavspe",Zavspe); + fhBookTree->SetBranchAddress("Erelspe",Erelspe); + fhBookTree->SetBranchAddress("Pathspe",Pathspe); + fhBookTree->SetBranchAddress("P0spe",P0spe); + fhBookTree->SetBranchAddress("Nxmult",Nxmult); + fhBookTree->SetBranchAddress("Nymult",Nymult); + fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx); + fhBookTree->SetBranchAddress("Npstripx",Npstripx); + fhBookTree->SetBranchAddress("Ntstripx",Ntstripx); + fhBookTree->SetBranchAddress("Istripx",Istripx); + fhBookTree->SetBranchAddress("Qstripx",Qstripx); + fhBookTree->SetBranchAddress("Xstripx",Xstripx); + fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy); + fhBookTree->SetBranchAddress("Npstripy",Npstripy); + fhBookTree->SetBranchAddress("Ntstripy",Ntstripy); + fhBookTree->SetBranchAddress("Istripy",Istripy); + fhBookTree->SetBranchAddress("Qstripy",Qstripy); + fhBookTree->SetBranchAddress("Ystripy",Ystripy); + fhBookTree->SetBranchAddress("Nthcali",&Nthcali); + fhBookTree->SetBranchAddress("Icaplane",Icaplane); + fhBookTree->SetBranchAddress("Icastrip",Icastrip); + fhBookTree->SetBranchAddress("Icamod",Icamod); + fhBookTree->SetBranchAddress("Enestrip",Enestrip); + fhBookTree->SetBranchAddress("Nthcal",&Nthcal); + fhBookTree->SetBranchAddress("Icapl",Icapl); + fhBookTree->SetBranchAddress("Icasi",Icasi); + fhBookTree->SetBranchAddress("Icast",Icast); + fhBookTree->SetBranchAddress("Xincal",Xincal); + fhBookTree->SetBranchAddress("Yincal",Yincal); + fhBookTree->SetBranchAddress("Zincal",Zincal); + fhBookTree->SetBranchAddress("Erelcal",Erelcal); + fhBookTree->SetBranchAddress("Nthnd",&Nthnd); + fhBookTree->SetBranchAddress("Itubend",Itubend); + fhBookTree->SetBranchAddress("Iparnd",Iparnd); + fhBookTree->SetBranchAddress("Xinnd",Xinnd); + fhBookTree->SetBranchAddress("Yinnd",Yinnd); + fhBookTree->SetBranchAddress("Zinnd",Zinnd); + fhBookTree->SetBranchAddress("Xoutnd",Xoutnd); + fhBookTree->SetBranchAddress("Youtnd",Youtnd); + fhBookTree->SetBranchAddress("Zoutnd",Zoutnd); + fhBookTree->SetBranchAddress("Erelnd",Erelnd); + fhBookTree->SetBranchAddress("Timend",Timend); + fhBookTree->SetBranchAddress("Pathnd",Pathnd); + fhBookTree->SetBranchAddress("P0nd",P0nd); + fhBookTree->SetBranchAddress("Nthcard",&Nthcard); + fhBookTree->SetBranchAddress("Iparcard",Iparcard); + fhBookTree->SetBranchAddress("Icard",Icard); + fhBookTree->SetBranchAddress("Xincard",Xincard); + fhBookTree->SetBranchAddress("Yincard",Yincard); + fhBookTree->SetBranchAddress("Zincard",Zincard); + fhBookTree->SetBranchAddress("Xoutcard",Xoutcard); + fhBookTree->SetBranchAddress("Youtcard",Youtcard); + fhBookTree->SetBranchAddress("Zoutcard",Zoutcard); + fhBookTree->SetBranchAddress("Erelcard",Erelcard); + fhBookTree->SetBranchAddress("Timecard",Timecard); + fhBookTree->SetBranchAddress("Pathcard",Pathcard); + fhBookTree->SetBranchAddress("P0card",P0card); + + fhBookTree->SetBranchStatus("*",0); + +}; + + + +void Digitizer::Close(){ + + delete fhBookTree; + +}; + + + + +void Digitizer::Loop() { + // + // opens the raw output file and loops over the events + // + fOutputfile.open(fFilename, ios::out | ios::binary); + //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary); + // + // Load in memory and save at the beginning of file the calorimeter calibration + // + CaloLoadCalib(); + DigitizeCALOCALIB(); + + // load, digitize and write tracker calibration + LoadTrackCalib(); + + DigitizeTrackCalib(1); + UInt_t length=fTracklength*2; + DigitizePSCU(length,0x12); + AddPadding(); + WriteTrackCalib(); + + DigitizeTrackCalib(2); + length=fTracklength*2; + DigitizePSCU(length,0x13); + AddPadding(); + WriteTrackCalib(); + + LoadMipCor(); // some initialization of parameters -not used now- + // end loading, digitizing and writing tracker calibration + + // + // loops over the events + // + + Int_t nentries = fhBookTree->GetEntriesFast(); + Long64_t nbytes = 0; + for (Int_t i=0; iGetEntry(i); + // read detectors sequentially: + // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf + // on pamelatov: + // /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp + DigitizeTRIGGER(); + DigitizeTOF(); + DigitizeAC(); + DigitizeCALO(); + DigitizeTrack(); + DigitizeS4(); + DigitizeND(); + // + // Add padding to 64 bits + // + AddPadding(); +// + // Create CPU header, we need packet type (0x10 = physics data) and packet length. + // + UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer; + //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer; + DigitizePSCU(length,0x10); + if ( !i%100 ) std::cout << "writing event " << i << endl; + WriteData(); + }; + + fOutputfile.close(); + std::cout << "files closed" << endl << flush; + +}; + +void Digitizer::AddPadding(){ + // + Float_t pd0 = (fLen+16)/64.; + Float_t pd1 = pd0 - (Float_t)int(pd0); + Float_t padfrac = 64. - pd1 * 64.; + // + UInt_t padbytes = (UInt_t)padfrac; + if ( padbytes > 0 && padbytes < 64 ){ + // + // here the padding length + // + fPadding = padbytes+64; + // + // random padding bytes + // + for (Int_t ur=0; ur<32; ur++){ + fDataPadding[ur] = (UShort_t)rand(); + }; + }; +}; + + +void Digitizer::DigitizePSCU(UInt_t length, UChar_t type) { + // + UChar_t buff[16]; + // + // CPU signature + // + buff[0] = 0xFA; + buff[1] = 0xFE; + buff[2] = 0xDE; + // + // packet type (twice) + // + buff[3] = type; + buff[4] = type; + // + // counter + // + fCounter++; + while ( fCounter > 16777215 ){ + fCounter -= 16777215; + }; + // + buff[5] = (UChar_t)(fCounter >> 16); + buff[6] = (UChar_t)(fCounter >> 8); + buff[7] = (UChar_t)fCounter; + // + // on board time + // + ULong64_t obt = fOBT + 30LL; + // + while ( obt > 4294967295LL ){ + obt -= 4294967295LL; + }; + fOBT = (UInt_t)obt; + // + buff[8] = (UChar_t)(fOBT >> 24); + buff[9] = (UChar_t)(fOBT >> 16); + buff[10] = (UChar_t)(fOBT >> 8); + buff[11] = (UChar_t)fOBT; + // + // Packet length + // + fLen = length; + // + buff[12] = (UChar_t)(fLen >> 16); + buff[13] = (UChar_t)(fLen >> 8); + buff[14] = (UChar_t)fLen; + // + // CPU header CRC + // + buff[15] = (BYTE)CM_Compute_CRC16((UINT16)0, (BYTE*)&buff, (UINT32)15); + // + memcpy(fDataPSCU,buff,16*sizeof(UChar_t)); + // +}; + +void Digitizer::ClearCaloCalib(Int_t s){ + // + fcstwerr[s] = 0; + fcperror[s] = 0.; + for ( Int_t d=0 ; d<11 ;d++ ){ + Int_t pre = -1; + for ( Int_t j=0; j<96 ;j++){ + if ( j%16 == 0 ) pre++; + fcalped[s][d][j] = 0.; + fcstwerr[s] = 0.; + fcperror[s] = 0.; + fcalgood[s][d][j] = 0.; + fcalthr[s][d][pre] = 0.; + fcalrms[s][d][j] = 0.; + fcalbase[s][d][pre] = 0.; + fcalvar[s][d][pre] = 0.; + }; + }; + return; +} + +Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){ + // + // + UInt_t e = 0; + if ( s == 0 ) e = 0; + if ( s == 1 ) e = 2; + if ( s == 2 ) e = 3; + if ( s == 3 ) e = 1; + // + ifstream myfile; + myfile.open(fcalname.Data()); + if ( !myfile ){ + return(-107); + }; + myfile.close(); + // + TFile *File = new TFile(fcalname.Data()); + if ( !File ) return(-108); + TTree *tr = (TTree*)File->Get("CalibCalPed"); + if ( !tr ) return(-109); + // + TBranch *calo = tr->GetBranch("CalibCalPed"); + // + pamela::CalibCalPedEvent *ce = 0; + tr->SetBranchAddress("CalibCalPed", &ce); + // + Long64_t ncalibs = calo->GetEntries(); + // + if ( !ncalibs ) return(-110); + // + calo->GetEntry(calibno); + // + if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) { + fcstwerr[s] = ce->cstwerr[s]; + fcperror[s] = ce->cperror[s]; + for ( Int_t d=0 ; d<11 ;d++ ){ + Int_t pre = -1; + for ( Int_t j=0; j<96 ;j++){ + if ( j%16 == 0 ) pre++; + fcalped[s][d][j] = ce->calped[e][d][j]; + fcalgood[s][d][j] = ce->calgood[e][d][j]; + fcalthr[s][d][pre] = ce->calthr[e][d][pre]; + fcalrms[s][d][j] = ce->calrms[e][d][j]; + fcalbase[s][d][pre] = ce->calbase[e][d][pre]; + fcalvar[s][d][pre] = ce->calvar[e][d][pre]; + }; + }; + } else { + printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n "); + File->Close(); + return(-111); + }; + File->Close(); + return(0); +} + + +void Digitizer::DigitizeCALOCALIB() { + // + // Header of the four sections + // + fSecCalo[0] = 0xAA00; // XE + fSecCalo[1] = 0xB100; // XO + fSecCalo[2] = 0xB600; // YE + fSecCalo[3] = 0xAD00; // YO + // + // length of the data is 0x1215 + // + fSecCALOLength[0] = 0x1215; // XE + fSecCALOLength[1] = 0x1215; // XO + fSecCALOLength[2] = 0x1215; // YE + fSecCALOLength[3] = 0x1215; // YO + // + Int_t chksum = 0; + UInt_t tstrip = 0; + UInt_t fSecPointer = 0; + // + for (Int_t sec=0; sec < 4; sec++){ + // + // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO + // + fCALOlength = 0; + memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer); + fSecPointer = fCALOlength; + // + // First of all we have section header and packet length + // + fDataCALO[fCALOlength] = fSecCalo[sec]; + fCALOlength++; + fDataCALO[fCALOlength] = fSecCALOLength[sec]; + fCALOlength++; + // + // Section XO is read in the opposite direction respect to the others + // + chksum = 0; + // + for (Int_t plane=0; plane < 11; plane++){ + // + if ( sec == 1 ) tstrip = fCALOlength + 96*2; + // + for (Int_t strip=0; strip < 96; strip++){ + // + chksum += (Int_t)fcalped[sec][plane][strip]; + // + // save value + // + if ( sec == 1 ){ + tstrip -= 2; + fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip]; + fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip]; + } else { + fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip]; + fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip]; + }; + fCALOlength +=2; + }; + // + }; + // + fDataCALO[fCALOlength] = (UShort_t)chksum; + fCALOlength++; + fDataCALO[fCALOlength] = 0; + fCALOlength++; + fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16)); + fCALOlength++; + // + // Section XO is read in the opposite direction respect to the others + // + chksum = 0; + // + for (Int_t plane=0; plane < 11; plane++){ + // + if ( sec == 1 ) tstrip = fCALOlength+6*2; + // + for (Int_t strip=0; strip < 6; strip++){ + // + chksum += (Int_t)fcalthr[sec][plane][strip]; + // + // save value + // + if ( sec == 1 ){ + tstrip -= 2; + fDataCALO[tstrip] = 0; + fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip]; + } else { + fDataCALO[fCALOlength] = 0; + fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip]; + }; + fCALOlength +=2; + }; + // + }; + // + fDataCALO[fCALOlength] = 0; + fCALOlength++; + fDataCALO[fCALOlength] = (UShort_t)chksum; + fCALOlength++; + fDataCALO[fCALOlength] = 0; + fCALOlength++; + fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16)); + fCALOlength++; + // + // Section XO is read in the opposite direction respect to the others + // + for (Int_t plane=0; plane < 11; plane++){ + // + if ( sec == 1 ) tstrip = fCALOlength+96*2; + // + for (Int_t strip=0; strip < 96; strip++){ + // + // save value + // + if ( sec == 1 ){ + tstrip -= 2; + fDataCALO[tstrip] = 0; + fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip]; + } else { + fDataCALO[fCALOlength] = 0; + fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip]; + }; + fCALOlength += 2; + }; + // + }; + // + // Section XO is read in the opposite direction respect to the others + // + for (Int_t plane=0; plane < 11; plane++){ + // + if ( sec == 1 ) tstrip = fCALOlength+6*4; + // + for (Int_t strip=0; strip < 6; strip++){ + // + // save value + // + if ( sec == 1 ){ + tstrip -= 4; + fDataCALO[tstrip] = 0; + fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip]; + fDataCALO[tstrip+2] = 0; + fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip]; + } else { + fDataCALO[fCALOlength] = 0; + fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip]; + fDataCALO[fCALOlength+2] = 0; + fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip]; + }; + fCALOlength +=4; + }; + // + }; + // + // + // here we calculate and save the CRC + // + fDataCALO[fCALOlength] = 0; + fCALOlength++; + Short_t CRC = 0; + for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){ + CRC=crc(CRC,fDataCALO[i+fSecPointer]); + }; + fDataCALO[fCALOlength] = (UShort_t)CRC; + fCALOlength++; + // + UInt_t length=fCALOlength*2; + DigitizePSCU(length,0x18); + // + // Add padding to 64 bits + // + AddPadding(); + // + fOutputfile.write(reinterpret_cast(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer); + UShort_t temp[1000000]; + memset(temp,0,sizeof(UShort_t)*1000000); + swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!! + fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fCALOlength); + // + // padding to 64 bytes + // + if ( fPadding ){ + fOutputfile.write(reinterpret_cast(fDataPadding),sizeof(UChar_t)*fPadding); + }; + // + // + }; + // +}; + +void Digitizer::CaloLoadCalib() { + // + fGivenCaloCalib = 0; // ####@@@@ should be given as input par @@@@#### + // + // first of all load the MIP to ADC conversion values + // + stringstream calfile; + Int_t error = 0; + GL_PARAM *glparam = new GL_PARAM(); + // + // determine where I can find calorimeter ADC to MIP conversion file + // + error = 0; + error = glparam->Query_GL_PARAM(3,101,fDbc); + // + calfile.str(""); + calfile << glparam->PATH.Data() << "/"; + calfile << glparam->NAME.Data(); + // + printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str()); + FILE *f; + f = fopen(calfile.str().c_str(),"rb"); + // + memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0])); + // + for (Int_t m = 0; m < 2 ; m++ ){ + for (Int_t k = 0; k < 22; k++ ){ + for (Int_t l = 0; l < 96; l++ ){ + fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f); + }; + }; + }; + fclose(f); + // + // determine which calibration has to be used and load it for each section + // + GL_CALO_CALIB *glcalo = new GL_CALO_CALIB(); + GL_ROOT *glroot = new GL_ROOT(); + TString fcalname; + UInt_t idcalib; + UInt_t calibno; + UInt_t utime = 0; + // + for (UInt_t s=0; s<4; s++){ + // + // clear calo calib variables for section s + // + ClearCaloCalib(s); + // + if ( fGivenCaloCalib ){ + // + // a time has been given as input on the command line so retrieve the calibration that preceed that time + // + glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc); + // + calibno = glcalo->EV_ROOT; + idcalib = glcalo->ID_ROOT_L0; + // + // determine path and name and entry of the calibration file + // + printf("\n"); + printf(" ** SECTION %i **\n",s); + // + glroot->Query_GL_ROOT(idcalib,fDbc); + // + stringstream name; + name.str(""); + name << glroot->PATH.Data() << "/"; + name << glroot->NAME.Data(); + // + fcalname = (TString)name.str().c_str(); + // + printf("\n Section %i : using file %s calibration at entry %i: \n",s,fcalname.Data(),calibno); + // + } else { + error = 0; + error = glparam->Query_GL_PARAM(1,104,fDbc); + // + calfile.str(""); + calfile << glparam->PATH.Data() << "/"; + calfile << glparam->NAME.Data(); + // + printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str()); + // + fcalname = (TString)calfile.str().c_str(); + calibno = s; + // + }; + // + // load calibration variables in memory + // + CaloLoadCalib(s,fcalname,calibno); + // + }; + // + // at this point we have in memory the calorimeter calibration and we can save it to disk in the correct format and use it to digitize the data + // + delete glparam; + delete glcalo; + delete glroot; +}; + +void Digitizer::DigitizeCALO() { + // + fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL ####@@@@ should be given as input par @@@@#### + // + // + // + fCALOlength = 0; // reset total dimension of calo data + // + // gpamela variables to be used + // + fhBookTree->SetBranchStatus("Nthcali",1); + fhBookTree->SetBranchStatus("Icaplane",1); + fhBookTree->SetBranchStatus("Icastrip",1); + fhBookTree->SetBranchStatus("Icamod",1); + fhBookTree->SetBranchStatus("Enestrip",1); + // + // call different routines depending on the acq mode you want to simulate + // + switch ( fModCalo ){ + case 0: + this->DigitizeCALORAW(); + break; + case 1: + this->DigitizeCALOCOMPRESS(); + break; + case 2: + this->DigitizeCALOFULL(); + break; + }; + // +}; + +Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){ + // + // determine plane and strip + // + Int_t mplane = 0; + // + // wrong! + // + // if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1; + // if ( sec == 1 ) mplane = (plane * 4) + 2 + 1; + // if ( sec == 2 ) mplane = (plane * 4) + 1 + 1; + // + if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1) from plane = 0, 11 + if ( sec == 1 ) mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1) from plane = 0, 11 + if ( sec == 2 ) mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1) from plane = 0, 11 + if ( sec == 3 ) mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1) from plane = 0, 11 + // + Int_t mstrip = strip + 1; + // + // search energy release in gpamela output + // + for (Int_t i=0; i XE 1 -> XO 2-> YE 3 -> YO + // + l = 0; // XE and XO are Y planes + if ( sec < 2 ) l = 1; // while YE and YO are X planes + // + fSecPointer = fCALOlength; + // + // First of all we have section header and packet length + // + fDataCALO[fCALOlength] = fSecCalo[sec]; + fCALOlength++; + fDataCALO[fCALOlength] = fSecCALOLength[sec]; + fCALOlength++; + // + // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers + // + for (Int_t autoplane=0; autoplane < 7; autoplane++){ + fDataCALO[fCALOlength] = 0x0000; + fCALOlength++; + }; + // + // + // here comes data + // + // + // Section XO is read in the opposite direction respect to the others + // + if ( sec == 1 ){ + tstrip = 96*11 + fCALOlength; + } else { + tstrip = 0; + }; + // + pre = -1; + // + for (Int_t strip=0; strip < 96; strip++){ + // + // which is the pre for this strip? + // + if (strip%16 == 0) { + pre++; + }; + // + if ( sec == 1 ) tstrip -= 11; + // + for (Int_t plane=0; plane < 11; plane++){ + // + // here is wrong!!!! + // + // + // if ( plane%2 == 0 && sec%2 != 0){ + // lpl = plane*2; + // } else { + // lpl = (plane*2) + 1; + // }; + // + if ( sec == 0 || sec == 3 ) lpl = plane * 2; + if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1; + // + // get the energy in GeV from the simulation for that strip + // + ens = this->GetCALOen(sec,plane,strip); + // + // convert it into ADC channels + // + adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio); + // + // sum baselines + // + adcbase = (UInt_t)fcalbase[sec][plane][pre]; + // + // add noise and pedestals + // + pedestal = fcalped[sec][plane][strip]; + rms = fcalrms[sec][plane][strip]/4.; + // + // Add random gaussian noise of RMS rms and Centered in the pedestal + // + pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms); + // + // Sum all contribution + // + adc = adcsig + adcbase + (Int_t)round(pedenoise); + // + // Signal saturation + // + if ( adc > 0x7FFF ) adc = 0x7FFF; + // + // save value + // + if ( sec == 1 ){ + fDataCALO[tstrip] = adc; + tstrip++; + } else { + fDataCALO[fCALOlength] = adc; + }; + fCALOlength++; + // + }; + // + if ( sec == 1 ) tstrip -= 11; + // + }; + // + // here we calculate and save the CRC + // + Short_t CRC = 0; + for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){ + CRC=crc(CRC,fDataCALO[i+fSecPointer]); + }; + fDataCALO[fCALOlength] = (UShort_t)CRC; + fCALOlength++; + // + }; + // + // for (Int_t i=0; iDigitizeCALORAW(); + return; + // + // + // + fSecCalo[0] = 0xEA00; + fSecCalo[1] = 0xF100; + fSecCalo[2] = 0xF600; + fSecCalo[3] = 0xED00; + // + // length of the data in DSP mode must be calculated on fly during digitization + // + memset(fSecCALOLength,0x0,4*sizeof(UShort_t)); + // + // here comes raw data + // + Int_t en = 0; + // + for (Int_t sec=0; sec < 4; sec++){ + fDataCALO[en] = fSecCalo[sec]; + en++; + fDataCALO[en] = fSecCALOLength[sec]; + en++; + for (Int_t plane=0; plane < 11; plane++){ + for (Int_t strip=0; strip < 11; strip++){ + fDataCALO[en] = 0x0; + en++; + }; + }; + }; + // +}; + +void Digitizer::DigitizeCALOFULL() { + // + printf(" FULL MODE STILL NOT IMPLEMENTED! \n"); + // + this->DigitizeCALORAW(); + return; + // + fSecCalo[0] = 0xEA00; + fSecCalo[1] = 0xF100; + fSecCalo[2] = 0xF600; + fSecCalo[3] = 0xED00; + // + // length of the data in DSP mode must be calculated on fly during digitization + // + memset(fSecCALOLength,0x0,4*sizeof(UShort_t)); + // + // here comes raw data + // + Int_t en = 0; + // + for (Int_t sec=0; sec < 4; sec++){ + fDataCALO[en] = fSecCalo[sec]; + en++; + fDataCALO[en] = fSecCALOLength[sec]; + en++; + for (Int_t plane=0; plane < 11; plane++){ + for (Int_t strip=0; strip < 11; strip++){ + fDataCALO[en] = 0x0; + en++; + }; + }; + }; + // +}; + +void Digitizer::DigitizeTRIGGER() { + //fDataTrigger: 153 bytes + for (Int_t j=0; j < 153; j++) + fDataTrigger[j]=0x00; +}; + +Int_t Digitizer::DigitizeTOF() { + //fDataTof: 12 x 23 bytes (=276 bytes) + UChar_t *pTof=fDataTof; + Bool_t DEBUG=false; + + // --- activate branches: + fhBookTree->SetBranchStatus("Nthtof",1); + fhBookTree->SetBranchStatus("Ipltof",1); + fhBookTree->SetBranchStatus("Ipaddle",1); + fhBookTree->SetBranchStatus("Xintof",1); + fhBookTree->SetBranchStatus("Yintof",1); + fhBookTree->SetBranchStatus("Xouttof",1); + fhBookTree->SetBranchStatus("Youttof",1); + fhBookTree->SetBranchStatus("Ereltof",1); + fhBookTree->SetBranchStatus("Timetof",1); + // not yet used: Zintof, Xouttof, Youttof, Zouttof + + // ------ evaluate energy in each pmt: ------ + // strip geometry (lenght/width) + Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0}; + //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0}; + + // S11 8 paddles 33.0 x 5.1 cm + // S12 6 paddles 40.8 x 5.5 cm + // S21 2 paddles 18.0 x 7.5 cm + // S22 2 paddles 15.0 x 9.0 cm + // S31 3 paddles 15.0 x 6.0 cm + // S32 3 paddles 18.0 x 5.0 cm + + Float_t FGeo[2]={0., 0.}; /* geometrical factor */ + + const Float_t Pho_keV = 10.; // photons per keV in scintillator + const Float_t echarge = 1.6e-19; // electron charge + Float_t Npho=0.; + Float_t QevePmt_pC[48]; + Float_t QhitPad_pC[2]={0., 0.}; + Float_t QhitPmt_pC[2]={0., 0.}; + Float_t pmGain = 3.5e6; /* PMT Gain: the same for all PMTs */ + Float_t effi=0.21; /* Efficienza di fotocatodo */ + + Float_t ADC_pC0=-58.1; // ADC/pC conversion coefficient 0 + Float_t ADC_pC1=1.728; // ADC/pC conversion coefficient 1 + Float_t ADC_pC2=-4.063e-05; // ADC/pC conversion coefficient 2 + Float_t ADC_pC3=-5.763e-08; // ADC/pC conversion coefficient 3 + + Float_t pCthres=40.; // threshold in charge + Int_t ADClast=4095; // no signal --> ADC ch=4095 + Int_t ADCsat=3100; // saturation value for the ADCs + Int_t ADCtof[48]; + + +// ---- introduce scale factors to tune simul ADC to real data 24-oct DC + Float_t ScaleFact[48]={0.18,0.22,0.35,0.26,0.47,0.35,0.31,0.37, + 0.44,0.23,0.38,0.60,0.39,0.29,0.40,0.23, + 0.30,0.66,0.22,1.53,0.17,0.55, + 0.84,0.19,0.21,1.64,0.62,0.13, + 0.18,0.15,0.10,0.14,0.14,0.14,0.14,0.12, + 0.26,0.18,0.25,0.23,0.20,0.40, + 0.19,0.23,0.25,0.23,0.25,0.20}; + + for(Int_t i=0; i<48; i++){ + QevePmt_pC[i] = 0; + ADCtof[i]=0; + } + + // ------ read calibration file (get A1, A2, lambda1, lambda2) + ifstream fileTriggerCalib; + TString ftrigname="TrigCalibParam.txt"; + fileTriggerCalib.open(ftrigname.Data()); + if ( !fileTriggerCalib ) { + printf("debug: no trigger calib file!\n"); + return(-117); //check output! + }; + Float_t atte1[48],atte2[48],lambda1[48],lambda2[48]; + Int_t temp=0; + // correct readout WM Oct '07 + for(Int_t i=0; i<48; i++){ + fileTriggerCalib >> temp; + fileTriggerCalib >> atte1[i]; + fileTriggerCalib >> lambda1[i]; + fileTriggerCalib >> atte2[i]; + fileTriggerCalib >> lambda2[i]; + fileTriggerCalib >> temp; + } + fileTriggerCalib.close(); + + Int_t ip, ipad; + //Int_t ipmt; + Int_t pmtleft=0, pmtright=0; + Int_t *pl, *pr; + pl = &pmtleft; + pr = &pmtright; + + // TDC variables: + Int_t TDClast=4095; // no signal --> TDC ch=4095 + Int_t TDCint[48]; + Float_t tdc[48],tdc1[48],tdcpmt[48]; + for(Int_t i=0; i<48; i++) { + tdcpmt[i] = 1000.; + tdc[i] = 0.; // 18-oct WM + tdc1[i] = 0.; // 18-oct WM + } + + Float_t thresh=10.; // to be defined better... (Wolfgang) + + // === TDC: simulate timing for each paddle + Float_t dt1 = 285.e-12 ; // single PMT resolution + // Float_t dt1 = 10.e-12 ; // TEST + Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50]; + for(Int_t j=0;j<48;j++) tdcres[j] = 50.E-12; // TDC resolution 50 picosec + for(Int_t j=0;j<48;j++) c1_S[j] = 500.; // cable length in channels + for(Int_t j=0;j<48;j++) c2_S[j] = 0.; + for(Int_t j=0;j<48;j++) c3_S[j] = 1000.; + for(Int_t j=0;j<48;j++) c1_S[j] = c1_S[j]*tdcres[j]; // cable length in sec + for(Int_t j=0;j<48;j++) c2_S[j] = c2_S[j]*tdcres[j]; + // ih = 0 + i1; // not used?? (Silvio) + + /* ********************************** start loop over hits */ + + for(Int_t nh=0; nh +// t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ; // Signal reaches PMT +// t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1; + + t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1; + t2 = t2 + fabs(path[1]/veff) + s_l_g[ip]/veff1 ; // Signal reaches PMT + + Float_t t1save = t1; + Float_t t2save = t2; + +/* + TRandom r; +// This does not work... WM - but works in my simulation code ?? +// t1 = r.Gaus(t1,dt1); //apply gaussian error dt +// t2 = r.Gaus(t2,dt1); //apply gaussian error dt +*/ + t1 = gRandom->Gaus(t1,dt1); //apply gaussian error dt + t2 = gRandom->Gaus(t2,dt1); //apply gaussian error dt + +// cout<<1E12*(t1save-t1)<<" "<<1E12*(t2save-t2)< thresh) { + if (tdcpmt[pmtleft] == 1000.) { // fill for the first time + tdcpmt[pmtleft] = t1; + tdc[pmtleft] = t1 + c2_S[pmtleft] ; // Signal reaches Coincidence + } + if (tdcpmt[pmtleft] < 1000.) // is already filled! + if (t1 < tdcpmt[pmtleft]) { + tdcpmt[pmtleft] = t1; + t1 = t1 + c2_S[pmtleft] ; // Signal reaches Coincidence + tdc[pmtleft] = t1; + } + } + if (QhitPmt_pC[1] > thresh) { + if (tdcpmt[pmtright] == 1000.) { // fill for the first time + tdcpmt[pmtright] = t2; + tdc[pmtright] = t2 + c2_S[pmtright] ; // Signal reaches Coincidence + } + if (tdcpmt[pmtright] < 1000.) // is already filled! + if (t2 < tdcpmt[pmtright]) { + tdcpmt[pmtright] = t2; + t2 = t2 + c2_S[pmtright] ; + tdc[pmtright] = t2; + } + } + + if (DEBUG) + cout<= pCthres){ + ADCtof[i]= (Int_t)(ADC_pC0 + ADC_pC1*QevePmt_pC[i] + ADC_pC2*pow(QevePmt_pC[i],2) + ADC_pC3*pow(QevePmt_pC[i],3)); + } else + ADCtof[i]= ADClast; +} + +// ---- introduce scale factors to tune simul ADC to real data 24-oct DC + + for(Int_t i=0; i<48; i++){ + if(ADCtof[i] != ADClast){ + // printf("%3d, %4d, %4.2f\n",i, ADCtof[i],ScaleFact[i]); + ADCtof[i]= Int_t (ADCtof[i]*ScaleFact[i]); + // printf("%3d, %4d,\n",i, ADCtof[i]); + } + } + + for(Int_t i=0; i<48; i++){ + if(ADCtof[i] != ADClast){ + if(ADCtof[i]> ADCsat) ADCtof[i]=ADCsat; + else if(ADCtof[i]< 0) ADCtof[i]=ADClast; + } + } + // ====== build TDC coincidence ====== + + Float_t t_coinc = 0; + Int_t ilast = 100; + for (Int_t ii=0; ii<48;ii++) + if (tdc[ii] > t_coinc) { + t_coinc = tdc[ii]; + ilast = ii; + } + + // cout<4093) TDCint[i]=TDClast; // 18-oct WM + if (DEBUG) + cout< ADClast) ADC[i]=ADClast; + } else + TDCint[i]= TDClast; + } + + if (DEBUG) + cout<<"-----------"<0)&&(ADCtof[i]<4095)) || ((TDCint[i]>0)&&(TDCint[i]<4095))) cout<SetBranchAddress("Ievnt",&Ievnt); + fhBookTree->SetBranchStatus("Nthcat",1); + fhBookTree->SetBranchStatus("Iparcat",1); + fhBookTree->SetBranchStatus("Icat",1); + fhBookTree->SetBranchStatus("Xincat",1); + fhBookTree->SetBranchStatus("Yincat",1); + fhBookTree->SetBranchStatus("Zincat",1); + fhBookTree->SetBranchStatus("Xoutcat",1); + fhBookTree->SetBranchStatus("Youtcat",1); + fhBookTree->SetBranchStatus("Zoutcat",1); + fhBookTree->SetBranchStatus("Erelcat",1); + fhBookTree->SetBranchStatus("Timecat",1); + fhBookTree->SetBranchStatus("Pathcat",1); + fhBookTree->SetBranchStatus("P0cat",1); + fhBookTree->SetBranchStatus("Nthcas",1); + fhBookTree->SetBranchStatus("Iparcas",1); + fhBookTree->SetBranchStatus("Icas",1); + fhBookTree->SetBranchStatus("Xincas",1); + fhBookTree->SetBranchStatus("Yincas",1); + fhBookTree->SetBranchStatus("Zincas",1); + fhBookTree->SetBranchStatus("Xoutcas",1); + fhBookTree->SetBranchStatus("Youtcas",1); + fhBookTree->SetBranchStatus("Zoutcas",1); + fhBookTree->SetBranchStatus("Erelcas",1); + fhBookTree->SetBranchStatus("Timecas",1); + fhBookTree->SetBranchStatus("Pathcas",1); + fhBookTree->SetBranchStatus("P0cas",1); + fhBookTree->SetBranchStatus("Nthcard",1); + fhBookTree->SetBranchStatus("Iparcard",1); + fhBookTree->SetBranchStatus("Icard",1); + fhBookTree->SetBranchStatus("Xincard",1); + fhBookTree->SetBranchStatus("Yincard",1); + fhBookTree->SetBranchStatus("Zincard",1); + fhBookTree->SetBranchStatus("Xoutcard",1); + fhBookTree->SetBranchStatus("Youtcard",1); + fhBookTree->SetBranchStatus("Zoutcard",1); + fhBookTree->SetBranchStatus("Erelcard",1); + fhBookTree->SetBranchStatus("Timecard",1); + fhBookTree->SetBranchStatus("Pathcard",1); + fhBookTree->SetBranchStatus("P0card",1); + + // In this simpliefied approach we will assume that once + // a particle releases > 0.5 mip in one of the 12 AC detectors it + // will fire. We will furthermore assume that both cards read out + // identical data. + + // If you develop your digitization algorithm, you should start by + // identifying the information present in level2 (post-darth-vader) + // data. + + Float_t SumEcat[5]; + Float_t SumEcas[5]; + Float_t SumEcard[5]; + for (Int_t k= 0;k<5;k++){ + SumEcat[k]=0.; + SumEcas[k]=0.; + SumEcard[k]=0.; + }; + + if (Nthcat>50 || Nthcas>50 || Nthcard>50) + printf("*** ERROR AC! NthAC out of range!\n\n"); + + // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi) + // based on J.Lundquist's calculations (PhD thesis, page 94) + // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are: + // 8.25470e-01 +- 1.79489e-02 + // 6.41609e-01 +- 2.65846e-02 + // 9.81177e+00 +- 1.21284e+00 + // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip + // + // NB: the PMT positions are needed! + + // look in CAT + // for (UInt_t k= 0;k<50;k++){ + for (Int_t k= 0;k 0) + SumEcat[Icat[k]] += Erelcat[k]; + }; + + // look in CAS + for (Int_t k= 0;k0) + SumEcas[Icas[k]] += Erelcas[k]; + }; + + // look in CARD + for (Int_t k= 0;k0) + SumEcard[Icard[k]] += Erelcard[k]; + }; + + // channel mapping Hit Map + // 1 CARD4 0 LSB + // 2 CAT2 0 + // 3 CAS1 0 + // 4 NC 0 + // 5 CARD2 0 + // 6 CAT4 1 + // 7 CAS4 0 + // 8 NC 0 + // 9 CARD3 0 + // 10 CAT3 0 + // 11 CAS3 0 + // 12 NC 0 + // 13 CARD1 0 + // 14 CAT1 0 + // 15 CAS2 0 + // 16 NC 0 MSB + + // In the first version only the hit-map is filled, not the SR. + + // Threshold: 0.8 MeV. + + Float_t thr = 8e-4; + + fDataAC[3] = 0x0000; + + if (SumEcas[0] > thr) fDataAC[3] = 0x0004; + if (SumEcas[1] > thr) fDataAC[3] += 0x4000; + if (SumEcas[2] > thr) fDataAC[3] += 0x0400; + if (SumEcas[3] > thr) fDataAC[3] += 0x0040; + + if (SumEcat[0] > thr) fDataAC[3] += 0x2000; + if (SumEcat[1] > thr) fDataAC[3] += 0x0002; + if (SumEcat[2] > thr) fDataAC[3] += 0x0200; + if (SumEcat[3] > thr) fDataAC[3] += 0x0020; + + if (SumEcard[0] > thr) fDataAC[3] += 0x1000; + if (SumEcard[1] > thr) fDataAC[3] += 0x0010; + if (SumEcard[2] > thr) fDataAC[3] += 0x0100; + if (SumEcard[3] > thr) fDataAC[3] += 0x0001; + + fDataAC[67] = fDataAC[3]; + + // shift registers + // the central bin is equal to the hitmap, all other bins in the shift register are 0 + for (UInt_t i=0; i<=15; i++){ + fDataAC[i+11] = 0x0000; + fDataAC[i+75] = 0x0000; + } + fDataAC[18] = fDataAC[3]; + fDataAC[82] = fDataAC[3]; + + // for (Int_t i=0; iSetBranchStatus("Nthtof",1); + fhBookTree->SetBranchStatus("Ipltof",1); + fhBookTree->SetBranchStatus("Ipaddle",1); + + fhBookTree->SetBranchStatus("Xintof",1); + fhBookTree->SetBranchStatus("Yintof",1); + fhBookTree->SetBranchStatus("Xouttof",1); + fhBookTree->SetBranchStatus("Youttof",1); + + fhBookTree->SetBranchStatus("Ereltof",1); + fhBookTree->SetBranchStatus("Timetof",1); + NdFT=0; + Ert=0; + for(i=0;iXs[0] && xYs[0] && yZs[0] && zXs[0] && xYs[0] && yZs[0] && z>ciao; + } + if((xXs[1]-p)&&(y>Ys[0]+p && yZs[0]+p && z=Yp[2*t] && yXs[0]+p && xYs[1]-p)&&(z>Zs[0]+p && zXs[0]+p && xYs[0]+p && yZs[1]-p))V[2]=-V[2]; + x+=V[0]; + y+=V[1]; + z+=V[2]; + l=0; + //cout<>ciao; + } + } + } + Ert=Ert/0.002; + q=(Float_t)(random())/(Float_t)0x7fffffff; + w=0.7; + //E0=(Float_t)(4064./7.); + E0=4064./7.; + if(Ert<1) S4=0; + else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0))); + i=S4/4; + if(S4%4==0) + S4v=S4+S4p; + else if(S4%4==1){ + if(q>ciao; +} + + + +void Digitizer::DigitizeND(){ + // creato: S. Borisov, INFN Roma2 e MEPHI, Sett 2007 + Int_t i=0; + UShort_t NdN=0; + fhBookTree->SetBranchStatus("Nthnd",1); + fhBookTree->SetBranchStatus("Itubend",1); + fhBookTree->SetBranchStatus("Iparnd",1); + fhBookTree->SetBranchStatus("Xinnd",1); + fhBookTree->SetBranchStatus("Yinnd",1); + fhBookTree->SetBranchStatus("Zinnd",1); + fhBookTree->SetBranchStatus("Xoutnd",1); + fhBookTree->SetBranchStatus("Youtnd",1); + fhBookTree->SetBranchStatus("Zoutnd",1); + fhBookTree->SetBranchStatus("Erelnd",1); + fhBookTree->SetBranchStatus("Timend",1); + fhBookTree->SetBranchStatus("Pathnd",1); + fhBookTree->SetBranchStatus("P0nd",1); + //cout<<"n="<SetBranchStatus("Enestrip",1); + + // dumy header + fDataDummy[0] = 0xCAAA; + + for (Int_t i=1; i(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer); + // TRG + fOutputfile.write(reinterpret_cast(fDataTrigger),sizeof(UChar_t)*153); + // TOF + fOutputfile.write(reinterpret_cast(fDataTof),sizeof(UChar_t)*276); + // AC + UShort_t temp[1000000]; + memset(temp,0,sizeof(UShort_t)*1000000); + swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer); // WE MUST SWAP THE BYTES!!! + fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fACbuffer); + // CALO + memset(temp,0,sizeof(UShort_t)*1000000); + swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!! + fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fCALOlength); + // TRK + memset(temp,0,sizeof(UShort_t)*1000000); + swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength); // WE MUST SWAP THE BYTES!!! + fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fTracklength); + fTracklength=0; + // padding to 64 bytes + // + if ( fPadding ){ + fOutputfile.write(reinterpret_cast(fDataPadding),sizeof(UChar_t)*fPadding); + }; + // S4 + memset(temp,0,sizeof(UShort_t)*1000000); + swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer); // WE MUST SWAP THE BYTES!!! + fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fS4buffer); + // ND + memset(temp,0,sizeof(UShort_t)*1000000); + swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer); // WE MUST SWAP THE BYTES!!! + fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fNDbuffer); +}; + + +void Digitizer::ReadData(){ + + UShort_t InData[64]; + + // for debuggigng purposes only, write your own routine if you like (many + // hardwired things. + + ifstream InputFile; + + // if (!InputFile) { + + // std::cout << "ERROR" << endl; + // // An error occurred! + // // myFile.gcount() returns the number of bytes read. + // // calling myFile.clear() will reset the stream state + // // so it is usable again. + // }; + + + + //InputFile.seekg(0); + + InputFile.open(fFilename, ios::in | ios::binary); + // fOutputfile.seekg(0); + if (!InputFile.is_open()) std::cout << "ERROR" << endl; + + InputFile.seekg(0); + + for (Int_t k=0; k<=1000; k++){ + InputFile.read(reinterpret_cast(InData),384*sizeof(UShort_t)); + + std::cout << "Read back: " << endl << endl; + + for (Int_t i=0; i<=384; i++){ + printf("%4x ", InData[i]); + if ((i+1)%8 ==0) cout << endl; + } + + } + cout << endl; + InputFile.close(); + +}; + + + +void Digitizer::DigitizeTrack() { +//std:: cout << "Entering DigitizeTrack " << endl; +Float_t AdcTrack[fNviews][fNstrips_view]; // Vector of strips to be compressed + +Int_t Iview; +Int_t Nstrip; + + for (Int_t j=0; jGaus(0.,fSigmaCommon); + Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon); + for (Int_t k=0; kGaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]); + if(k<4*128) {AdcTrack[j][Nstrip] += commonN1;} // full correlation of 4 VA1 Com. Noise + else {AdcTrack[j][Nstrip] += commonN2;} // full correlation of 4 VA1 Com. Noise + + }; + + + }; + + + }; + + + fhBookTree->SetBranchStatus("Nstrpx",1); + fhBookTree->SetBranchStatus("Npstripx",1); + fhBookTree->SetBranchStatus("Ntstripx",1); + fhBookTree->SetBranchStatus("Istripx",1); + fhBookTree->SetBranchStatus("Qstripx",1); + fhBookTree->SetBranchStatus("Xstripx",1); + fhBookTree->SetBranchStatus("Nstrpy",1); + fhBookTree->SetBranchStatus("Npstripy",1); + fhBookTree->SetBranchStatus("Ntstripy",1); + fhBookTree->SetBranchStatus("Istripy",1); + fhBookTree->SetBranchStatus("Qstripy",1); + fhBookTree->SetBranchStatus("Ystripy",1); + + + + Float_t ADCfull; + Int_t iladd=0; + for (Int_t ix=0; ix=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1; + if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2; + ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview]; + AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull); + + }; + + + for (Int_t iy=0; iy=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1; + if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2; + ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview]; + AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull); + + }; + +CompressTrackData(AdcTrack); // Compress and Digitize data of one Ladder in turn for all ladders + +}; + + + +void Digitizer::DigitizeTrackCalib(Int_t ii) { + +std:: cout << "Entering DigitizeTrackCalib " << ii << endl; +if( (ii!=1)&&(ii!=2) ) { + std:: cout << "error wrong DigitizeTrackCalib argument" << endl; + return; +}; + +memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); +fTracklength=0; + +UShort_t Dato; + +Float_t dato1; +Float_t dato2; +Float_t dato3; +Float_t dato4; + +UShort_t DatoDec; +UShort_t DatoDec1; +UShort_t DatoDec2; +UShort_t DatoDec3; +UShort_t DatoDec4; + +UShort_t EVENT_CAL; +UShort_t PED_L1; +UShort_t ReLength; +UShort_t OveCheckCode; +//UShort_t PED_L2; +//UShort_t PED_L3HI; +//UShort_t PED_L3LO; +//UShort_t SIG_L1HI; +//UShort_t SIG_L1LO; +//UShort_t SIG_L2HI; +//UShort_t SIG_L2LO; +//UShort_t SIG_L3; +//UShort_t BAD_L1; +//UShort_t BAD_L2LO; +//UShort_t BAD_L3HI; +//UShort_t BAD_L3LO; +//UShort_t FLAG; + + + Int_t DSPpos; + for (Int_t j=ii-1; j> 8); + fTracklength++; + + fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) ); + fTracklength++; + + // end TRAILER + }; + + // write in buffer the DSP header + + fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) ); + + fDataTrack[DSPpos+1]=0x01A9; + + fDataTrack[DSPpos+2]=0x8740; + + EVENT_CAL=0; + fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) ); + + PED_L1=0; + fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) ); + + // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0 + + fDataTrack[DSPpos+5]=0x8014; + + fDataTrack[DSPpos+6]=0x00A0; + + fDataTrack[DSPpos+7]=0x0500; + + fDataTrack[DSPpos+8]=0x2801; + + fDataTrack[DSPpos+9]=0x400A; + + fDataTrack[DSPpos+10]=0x0050; + + CkSum=(CkSum >> 8)^(CkSum&0x00FF); + fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3)); + + fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) ); + + // end dsp header + + // write in buffer the TRAILER + + ReLength=(UShort_t)((13*2)+3); + OveCheckCode=0x0000; + fDataTrack[DSPpos+13]=0x0000; + + fDataTrack[DSPpos+14]=(ReLength >> 8); + + fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) ); + + // end TRAILER + + + + + // end DSP + }; + + + +}; + +void Digitizer::WriteTrackCalib() { + + +std:: cout << " Entering WriteTrackCalib " << endl; + +fOutputfile.write(reinterpret_cast(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer); + +UShort_t temp[1000000]; +memset(temp,0,sizeof(UShort_t)*1000000); +swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength); // WE MUST SWAP THE BYTES!!! +fOutputfile.write(reinterpret_cast(temp),sizeof(UShort_t)*fTracklength); +fTracklength=0; +if ( fPadding ){ + fOutputfile.write(reinterpret_cast(fDataPadding),sizeof(UChar_t)*fPadding); +}; + +}; + + +void Digitizer::ClearTrackCalib() { + +std:: cout << "Entering ClearTrackCalib " << endl; + + +}; + + +void Digitizer::LoadTrackCalib() { +std:: cout << "Entering LoadTrackCalib " << endl; + +// Generate the pedestals and sigmas according to parametrization + for (Int_t j=0; jGaus(fAvePedex,fSigmaPedex); + fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax); + }; + if((j+1)%2==1) { + fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey); + fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay); + }; + + }; + }; + + + +}; + +void Digitizer::LoadMipCor() { +std:: cout << "Entering LoadMipCor" << endl; + Float_t xfactor=1./151.6*1.04; + Float_t yfactor=1./152.1; + + fMipCor[0][0]=140.02*yfactor; + fMipCor[0][1]=140.99*xfactor; + fMipCor[0][2]=134.48*yfactor; + fMipCor[0][3]=144.41*xfactor; + fMipCor[0][4]=140.74*yfactor; + fMipCor[0][5]=142.28*xfactor; + fMipCor[0][6]=134.53*yfactor; + fMipCor[0][7]=140.63*xfactor; + fMipCor[0][8]=135.55*yfactor; + fMipCor[0][9]=138.00*xfactor; + fMipCor[0][10]=154.95*yfactor; + fMipCor[0][11]=158.44*xfactor; + + + fMipCor[1][0]=136.07*yfactor; + fMipCor[1][1]=135.59*xfactor; + fMipCor[1][2]=142.69*yfactor; + fMipCor[1][3]=138.19*xfactor; + fMipCor[1][4]=137.35*yfactor; + fMipCor[1][5]=140.23*xfactor; + fMipCor[1][6]=153.15*yfactor; + fMipCor[1][7]=151.42*xfactor; + fMipCor[1][8]=129.76*yfactor; + fMipCor[1][9]=140.63*xfactor; + fMipCor[1][10]=157.87*yfactor; + fMipCor[1][11]=153.64*xfactor; + + fMipCor[2][0]=134.98*yfactor; + fMipCor[2][1]=143.95*xfactor; + fMipCor[2][2]=140.23*yfactor; + fMipCor[2][3]=138.88*xfactor; + fMipCor[2][4]=137.95*yfactor; + fMipCor[2][5]=134.87*xfactor; + fMipCor[2][6]=157.56*yfactor; + fMipCor[2][7]=157.31*xfactor; + fMipCor[2][8]=141.37*yfactor; + fMipCor[2][9]=143.39*xfactor; + fMipCor[2][10]=156.15*yfactor; + fMipCor[2][11]=158.79*xfactor; + +/* + for (Int_t j=0; j 0.5) inte=inte+1; + newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k]; + // first strip of ladder is transmitted + // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl; + DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); + DSPlength++; + ntrastot++; + trasmesso=1; + oldval=newval; + kt=k; + k++; + continue; + }; + real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte); + if (real > 0.5) inte=inte+1; + newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]); + cercacluster=1; // ????????? + if (cercacluster==1) { + + // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!! + Int_t diff=0; + + + switch ((iv+1)%2) { + case 0: diff=newval-oldval; + break; + case 1: diff=oldval-newval; + break; + }; + + if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) { + Int_t clval=newval; + Int_t klp=k; // go on to search for maximum + klp++; + + while(klp 0.5) inte=inte+1; + Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp]; + if((iv+1)%2==0) { + + if(clvalp>clval) { + clval=clvalp; + k=klp;} + else break; // max of cluster found + + } else { + + if(clvalp=kl1) kl1=kt+1; + if( (kt+1)==kl1 ) trasmesso=1; + + + + Int_t kl2=k+fNclst; + if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1; + + for(Int_t klt=kl1 ; klt<=kl2 ; klt++) { + if(trasmesso==0) { + // std:: cout << "STRIP " << klt << endl; + // std:: cout << "ADC_TOT " < 0.5) inte=inte+1; + DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); + DSPlength++; + ntrastot++; + + } + else { + // std:: cout << "ADC_TOT " < 0.5) inte=inte+1; + DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); + DSPlength++; + ntrastot++; + }; + trasmesso=1; + }; // end trasmission + kt=kl2; + k=kl2; + real=modff(AdcTrack[iv][ladder*fNstrips_ladder+kt],&inte); + if (real > 0.5) inte=inte+1; + oldval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+kt]; + k++; + continue; + + + }; // end cercacluster + }; // end cercacluster + +// start ZOP check for strips no + + if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) { + + if(trasmesso==0) { + // std:: cout << "STRIP " << k << endl; + // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl; + + DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000); + DSPlength++; + ntrastot++; + + + real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte); + if (real > 0.5) inte=inte+1; + DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); + DSPlength++; + ntrastot++; + + } + else { + // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl; + real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte); + if (real > 0.5) inte=inte+1; + DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF); + DSPlength++; + ntrastot++; + }; + trasmesso=1; + oldval=newval; + kt=k; + + } + else trasmesso=0; + // end zop + + k++; + }; // end cycle inside ladder +// write here the end ladder bytes +// std:: cout << "FINE LADDER " << ladder+1 << endl; + + DataDSP[DSPlength]=( ((UShort_t)(ladder+1)) | 0x1800); + DSPlength++; + ntrastot++; + trasmesso=0; + + }; //end cycle inside dsp +// std:: cout << "FINE DSP " << iv+1 << endl; +// here put DSP header + DataDSP[0]=(0x1CA0 | ((UShort_t)(iv+1)) ); + UShort_t Nword=(DSPlength*13)/16; + if( ((DSPlength*13)%16)!=0) Nword++; + DataDSP[1]=(0x1400 | ( Nword >> 10)); + DataDSP[2]=(0x1400 | ( Nword & 0x03FF) ); + DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) ); + DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) ) & 0x03FF) ); + DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) ) + | ( (UShort_t)fCutzop ) ); + DataDSP[6]=0x1400; + DataDSP[7]=0x1400; + DataDSP[8]=0x1400; + DataDSP[9]=0x1400; + DataDSP[10]=0x1400; + DataDSP[11]=0x1400; + DataDSP[12]=0x1400; + DataDSP[13]=0x1400; + DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) ); + DataDSP[15]=0x1C00; +// end DSP header + + +// write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer + Int_t Bit16free=16; + UShort_t Dato; + for (Int_t NDSP=0; NDSP0) { + if(Bit13ToWrite<=Bit16free) { + Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite)); + fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ; + Bit16free=Bit16free-Bit13ToWrite; + Bit13ToWrite=0; + if(Bit16free==0) { + if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength]; + fTracklength++; + Bit16free=16; + }; + } + else if(Bit13ToWrite>Bit16free) { + Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) ); + fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ; + if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength]; + fTracklength++; + Bit13ToWrite=Bit13ToWrite-Bit16free; + Bit16free=16; + }; + + }; // end cycle while(Bit13ToWrite>0) + + }; // end cycle DataDSP + if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; }; + CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF); + fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3)); + fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) ); + +// end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer + +//write trailer on buffer + UShort_t ReLength=(UShort_t)((Nword+13)*2+3); + UShort_t OveCheckCode=0x0000; + + fDataTrack[fTracklength]=0x0000; + fTracklength++; + + fDataTrack[fTracklength]=(ReLength >> 8); + fTracklength++; + + fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) ); + fTracklength++; +// end trailer +// std:: cout << "DSPlength " <3000.) { SatFact=3000./ADC; }; + return SatFact; +}; + + + + + +