#include "Digitizer.h" extern "C"{ short crc(short, short); }; 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,fDataPSCU); // // 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() { // // fCALOlength = 0; // reset total dimension of calo data // // gpamela variables to be used // // fhBookTree->SetBranchStatus("Nthcali",1);//modified by E.Vannuccini 03/08 // 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; 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] = 0; // Unknown: length must be calculated on fly 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++; }; // // second level trigger // fDataCALO[fCALOlength] = 0x0000; fCALOlength++; // // Nof strips transmetted: must be calculated on fly // fNofTStripsPointer = fCALOlength; fDataCALO[fCALOlength] = 0x0000; fCALOlength++; NofTransmittedStrips=0; // // Identifier of calo data // fDataCALO[fCALOlength] = 0xCA50; fCALOlength++; fDataCALO[fCALOlength] = 0xCA50; fCALOlength++; fDataCALO[fCALOlength] = 0xFFFF; // compress mode fCALOlength++; // // Pedestal threashold table checksum // fDataCALO[fCALOlength] = 0x0000; fCALOlength++; // // Calorimeter event counter // fDataCALO[fCALOlength] = fEvent; fCALOlength++; // // Start here with data // plane=-1; npre =-1; for (Int_t ipre=0; ipre< 66; ipre++){ // (11 planes*6 preampl) // // which plane if ( (ipre % 6) == 0) { plane++; } // pre=ipre; // // Adjust counter for plane X0 if (sec==1) // conto invertito { remainder = pre % 6 ; pre = ((plane+1)*6) - remainder ; } // if ( sec == 0 || sec == 3 ) lpl = plane * 2; if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1; // // initialize min_adc min_adc = 0x7FFF; for (Int_t ch=0; ch <16; ch++){ // 16 channels each pre // // strip number // strip=((pre-(6*plane))*16)+ch; if(sec==1) strip = ((pre-(6*plane))*16)+(15-ch)-16; // // calculate npre: a number between 0-5 // if( sec==1) { if ( ((95-strip) % 16) == 0) { npre++; if(npre>5) npre=0; } } else { if ( (strip % 16) == 0) { npre++; if(npre>5) npre=0; } } // // 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][npre]; // // 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[ch] = adcsig + adcbase + (Int_t)round(pedenoise); // // Signal saturation // if ( adc[ch] > 0x7FFF ) adc[ch] = 0x7FFF; // // save infos // pedround[ch] = (Int_t)round(pedestal) ; thres[ch] = ( fcalthr[sec][plane][npre] ); goodflag[ch] = ( fcalgood[sec][plane][strip] ); // if bad should be 255 // // Find minimum adc in this preamp // if ( goodflag[ch]==0 && (adc[ch]-pedround[ch])=9 // { if(sec==1) { fDataCALO[fCALOlength] = 0x0800 + ipre ; } else { fDataCALO[fCALOlength] = 0x0800 + pre; }; fCALOlength++; // // calculate baseline and save it // basenof=0; baseline=0; basesum=0; for (Int_t ch=0; ch <16; ch++){ if( goodflag[ch]==0 && ( (adc[ch]-pedround[ch])<(min_adc+thres[ch]) ) ) { basesum = basesum + (adc[ch]-pedround[ch]) ; basenof++; }; }; baseline = (Int_t)round( basesum / basenof ); fDataCALO[fCALOlength] = baseline; fCALOlength++; // // Transmit only channels > (min_adc+thres[ch]) // for (Int_t ch=0; ch <16; ch++){ if ( (adc[ch]-pedround[ch] )>(min_adc+thres[ch]) ) { fDataCALO[fCALOlength] = ch; fCALOlength++; fDataCALO[fCALOlength] = adc[ch]; fCALOlength++; NofTransmittedStrips++; }; }; }; // close if nof_chs_below }; // close preampl loop // // Write the correct length // fDataCALO[fSecPointer+1] = fCALOlength-fSecPointer+1 ; fDataCALO[fNofTStripsPointer] = NofTransmittedStrips ; // // 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++; // }; // close sec loop // The End } void Digitizer::DigitizeCALOFULL() { // printf(" FULL MODE STILL NOT IMPLEMENTED! %d\n",fEvent); // 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++; }; }; }; // }