#include "PamVMCCaloDig.h" #include #include "CalibCalPedEvent.h" using namespace pamela; ClassImp(PamVMCCaloDig) extern "C"{ short crc(short, short); }; void PamVMCCaloDig::LoadCalib(){ cout<<"Loading Tracker Calibrations..."<Query_GL_PARAM(time,type); if(fdberr<0){ cout<<"No such record in DB for CAL: time="<GetPAR()->PATH.Data() << "/"; fquery << fsql->GetPAR()->NAME.Data(); ThrowCalFileUsage("CAL",fquery.str().c_str()); FILE *f; f = fopen(fquery.str().c_str(),"rb"); if(f==NULL) ThrowCalFileWarning("CAL"); else{ 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 UInt_t idcalib, calibno; UInt_t utime = 0; TString calname; for (UInt_t s=0; s<4; s++){ // clear calo calib variables for section s ClearCaloCalib(s); if ( GivenCaloCalib ){ // a time has been given as input on the command line // so retrieve the calibration that preceed that time fsql->Query_GL_CALO_CALIB(GivenCaloCalib,utime,s); calibno = fsql->GetCaloCalib()->EV_ROOT; idcalib = fsql->GetCaloCalib()->ID_ROOT_L0; // determine path and name and entry of the calibration file printf("\n"); printf(" ** SECTION %i **\n",s); fsql->Query_GL_ROOT(idcalib); fquery.str(""); fquery << fsql->GetROOT()->PATH.Data() << "/"; fquery << fsql->GetROOT()->NAME.Data(); calname = (TString)fquery.str().c_str(); printf("\n Section %i : using file %s calibration at entry %i: \n",s,calname.Data(),calibno); } else { fdberr = 0; fdberr = fsql->Query_GL_PARAM(1,104); fquery.str(""); fquery << fsql->GetPAR()->PATH.Data() << "/"; fquery << fsql->GetPAR()->NAME.Data(); printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,fquery.str().c_str()); calname = (TString)fquery.str().c_str(); calibno = s; }; // load calibration variables in memory CaloLoadCalib(s,calname,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 // DigitizeCALOCALIB(); } } } void PamVMCCaloDig::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 PamVMCCaloDig::CaloLoadCalib(Int_t s,TString calname, UInt_t calibno){ UInt_t e = 0; switch(s){ case 0: e = 0; break; case 1: e = 2; break; case 2: e = 3; break; case 3: e = 1; break; default: break; } fcfile.open(calname.Data()); ThrowCalFileUsage("CAL",calname.Data()); if ( !fcfile ){ ThrowCalFileWarning("CAL"); return(-107); }; fcfile.close(); fcrfile = new TFile(calname.Data()); if ( !fcrfile ){ ThrowCalFileWarning("CAL"); return(-108); } TTree *tr = (TTree*)fcrfile->Get("CalibCalPed"); if ( !tr ){ cout<<"!!!WARNING Tree CalibCalPed in file "<GetBranch("CalibCalPed"); CalibCalPedEvent *ce = 0; tr->SetBranchAddress("CalibCalPed", &ce); Long64_t ncalibs = calo->GetEntries(); if ( !ncalibs ){ cout<<"!!!WARNING No entries in calibration file!!!"<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 "); fcrfile->Close(); return(-111); }; fcrfile->Close(); delete ce; return(0); } void PamVMCCaloDig::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 Short_t CRC; for(Int_t sec=0; sec < 4; sec++){ // // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO // fData.clear(); fData.push_back(fSecCalo[sec]); fData.push_back(fSecCaloLength[sec]); if ( sec==1 ){ FillCalPedBack(&fData,sec); FillCalThrBack(&fData,sec); FillCalRmsBack(&fData,sec); FillCalBaseVarBack(&fData,sec); } else { FillCalPedNorm(&fData,sec); FillCalThrNorm(&fData,sec); FillCalRmsNorm(&fData,sec); FillCalBaseVarNorm(&fData,sec); } //CRC fData.push_back(0x0000); CRC = 0; USBuffer::const_iterator p = fData.begin(); while( p!= fData.end() ){ CRC=crc(CRC,*p); p++; } fData.push_back(CRC); //Writing calibrations DigitizePSCU(fData.size()*2,0x18); fraw->WritePSCU(&fDataPSCU); fraw->CopyUShortToBuff(&fData); //No padding for now... } } void PamVMCCaloDig::FillCalPedNorm(USBuffer *buff, Int_t sec){ Int_t chksum = 0; for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=0; strip < 96; strip++){ chksum += (Int_t)fcalped[sec][plane][strip]; buff->push_back((Int_t)fcalped[sec][plane][strip]); buff->push_back((Int_t)fcalgood[sec][plane][strip]); } }//3-2114 buff->push_back(((UShort_t)chksum)); buff->push_back(0x0000); buff->push_back(((UShort_t)((Int_t)(chksum >> 16)))); } void PamVMCCaloDig::FillCalPedBack(USBuffer *buff, Int_t sec){ Int_t chksum = 0; for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=95; strip >= 0; strip--){ chksum += (Int_t)fcalped[sec][plane][strip]; buff->push_back((Int_t)fcalped[sec][plane][strip]); buff->push_back((Int_t)fcalgood[sec][plane][strip]); } }//3-2114 buff->push_back(((UShort_t)chksum)); buff->push_back(0x0000); buff->push_back(((UShort_t)((Int_t)(chksum >> 16)))); } void PamVMCCaloDig::FillCalThrNorm(USBuffer *buff, Int_t sec){ Int_t chksum = 0; for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=0; strip < 6; strip++){ chksum += (Int_t)fcalthr[sec][plane][strip]; buff->push_back(0x0000); buff->push_back((Int_t)fcalthr[sec][plane][strip]); } }//2118-2249 buff->push_back(0x0000); buff->push_back(((UShort_t)chksum)); buff->push_back(0x0000); buff->push_back(((UShort_t)((Int_t)(chksum >> 16)))); } void PamVMCCaloDig::FillCalThrBack(USBuffer *buff, Int_t sec){ Int_t chksum = 0; for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=5; strip >= 0; strip--){ chksum += (Int_t)fcalthr[sec][plane][strip]; buff->push_back(0x0000); buff->push_back((Int_t)fcalthr[sec][plane][strip]); } }//2118-2249 buff->push_back(0x0000); buff->push_back(((UShort_t)chksum)); buff->push_back(0x0000); buff->push_back(((UShort_t)((Int_t)(chksum >> 16)))); } void PamVMCCaloDig::FillCalRmsNorm(USBuffer *buff, Int_t sec){ for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=0; strip < 96; strip++){ buff->push_back(0x0000); buff->push_back((Int_t)fcalrms[sec][plane][strip]); } }//2254-4365 } void PamVMCCaloDig::FillCalRmsBack(USBuffer *buff, Int_t sec){ for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=95; strip >= 0; strip--){ buff->push_back(0x0000); buff->push_back((Int_t)fcalrms[sec][plane][strip]); } }//2254-4365 } void PamVMCCaloDig::FillCalBaseVarNorm(USBuffer *buff, Int_t sec){ for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=0; strip < 6; strip++){ buff->push_back(0x0000); buff->push_back(((Int_t)fcalbase[sec][plane][strip])); buff->push_back(0x0000); buff->push_back(((Int_t)fcalvar[sec][plane][strip])); } }//4366-4629 } void PamVMCCaloDig::FillCalBaseVarBack(USBuffer *buff, Int_t sec){ for (Int_t plane=0; plane < 11; plane++){ for (Int_t strip=5; strip >= 0; strip--){ buff->push_back(0x0000); buff->push_back((Int_t)fcalbase[sec][plane][strip]); buff->push_back(0x0000); buff->push_back((Int_t)fcalvar[sec][plane][strip]); } }//4366-4629 } void PamVMCCaloDig::Digitize(){ cout<<"Digitizing CALO..."<DigitizeCaloRaw(); break; case 1: this->DigitizeCaloCompress(); break; case 2: this->DigitizeCaloFull(); break; }; } void PamVMCCaloDig::DigitizeCaloRaw(){ // some variables // Float_t ens = 0.; UInt_t adcsig = 0; UInt_t adcbase = 0; UInt_t adc = 0; Int_t pre = 0; UInt_t l = 0; UInt_t lpl = 0; UInt_t tstrip = 0; UInt_t fSecPointer = 0; Int_t fCALOlength; Double_t pedenoise; Float_t rms = 0.; Float_t pedestal = 0.; // // clean the data structure // UShort_t DataCALO[4264]; fData.clear(); memset(DataCALO,0,sizeof(UShort_t)*4264); static const Float_t CALOGeV2MIPratio = 0.0001059994; // Header of the four sections fSecCalo[0] = 0xEA08; // XE fSecCalo[1] = 0xF108; // XO fSecCalo[2] = 0xF608; // YE fSecCalo[3] = 0xED08; // YO // length of the data is 0x0428 in RAW mode fSecCaloLength[0] = 0x0428; // XE fSecCaloLength[1] = 0x0428; // XO fSecCaloLength[2] = 0x0428; // YE fSecCaloLength[3] = 0x0428; // YO // let's start // fCALOlength = 0; // for (Int_t sec=0; sec < 4; sec++){ // // sec = 0 -> 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 // DataCALO[fCALOlength] = fSecCalo[sec]; //1 fCALOlength++; DataCALO[fCALOlength] = fSecCaloLength[sec]; //2 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++){ DataCALO[fCALOlength] = 0x0000; fCALOlength++; }; //3-8 // // // here comes data // // // Section XO is read in the opposite direction respect to the others // if ( sec == 1 ){ tstrip = 96*11 + fCALOlength; //tstrip = 1064 } 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++){ 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->GetCaloErel(sec,plane,strip); // // convert it into ADC channels // adcsig = int(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio); // // 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 = frandom->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 ){ DataCALO[tstrip] = adc; tstrip++; } else { DataCALO[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,DataCALO[i+fSecPointer]); }; DataCALO[fCALOlength] = (UShort_t)CRC; fCALOlength++; }; for (Int_t i = 0; iAt(mplane*96+strip); if (hit) Erel=hit->GetEREL(); // if(Erel) cout<<"sec:"<GetPOS()< 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 // DataCALO[fCALOlength] = fSecCalo[sec]; fCALOlength++; DataCALO[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++){ DataCALO[fCALOlength] = 0x0000; fCALOlength++; }; // // second level trigger // DataCALO[fCALOlength] = 0x0000; fCALOlength++; // // Nof strips transmetted: must be calculated on fly // fNofTStripsPointer = fCALOlength; DataCALO[fCALOlength] = 0x0000; fCALOlength++; NofTransmittedStrips=0; // // Identifier of calo data // DataCALO[fCALOlength] = 0xCA50; fCALOlength++; DataCALO[fCALOlength] = 0xCA50; fCALOlength++; DataCALO[fCALOlength] = 0xFFFF; // compresso fCALOlength++; // // Pedestal threashold table checksum // DataCALO[fCALOlength] = 0x0000; fCALOlength++; // // Calorimeter event counter // DataCALO[fCALOlength] = Getevtcalo() ; fCALOlength++; //cout<<" evtcalo?"<5) npre=0; } } else { if ( (strip % 16) == 0) { npre++; if(npre>5) npre=0; } } // ens = this->GetCaloErel(sec,plane,strip); // // convert it into ADC channels // adcsig = int(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio); // // 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 = frandom->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) { DataCALO[fCALOlength] = 0x0800 + ipre ; } else { DataCALO[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 ); DataCALO[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]) ) { DataCALO[fCALOlength] = ch; fCALOlength++; DataCALO[fCALOlength] = adc[ch]; fCALOlength++; NofTransmittedStrips++; } }; } // close if nof_chs_below }; // close preampl loop // // Write the length // DataCALO[fSecPointer+1] = (fCALOlength-fSecPointer+1)-2 ; // total length of the packet: -2: because the words with status and length are not included DataCALO[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,DataCALO[i+fSecPointer]); }; DataCALO[fCALOlength] = (UShort_t)CRC; fCALOlength++; // }; // close sec loop Incrementevtcalo(); for (Int_t i = 0; iDigitizeCaloRaw(); return; }