X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSRawDigiProducer.cxx;h=97ecbed7f410d62107d7930c7c07703c67f79fa0;hb=bb6fc63e31ec91b569c16d02befa2046c5f3e138;hp=582b539cb3fc7ba64c351f095c4d7e6a734e2b5b;hpb=77ea1c6f5e1848e31e148c393a9d38cac3957d0b;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSRawDigiProducer.cxx b/PHOS/AliPHOSRawDigiProducer.cxx index 582b539cb3f..97ecbed7f41 100644 --- a/PHOS/AliPHOSRawDigiProducer.cxx +++ b/PHOS/AliPHOSRawDigiProducer.cxx @@ -37,9 +37,70 @@ #include "AliPHOSRawDecoder.h" #include "AliPHOSGeometry.h" #include "AliPHOSDigit.h" +#include "AliPHOSRecoParam.h" +#include "AliPHOSCalibData.h" +#include "AliPHOSPulseGenerator.h" +#include "AliLog.h" ClassImp(AliPHOSRawDigiProducer) +AliPHOSCalibData * AliPHOSRawDigiProducer::fgCalibData = 0 ; + +//-------------------------------------------------------------------------------------- +AliPHOSRawDigiProducer::AliPHOSRawDigiProducer():TObject(), + fEmcMinE(0.),fCpvMinE(0.),fSampleQualityCut(1.), + fEmcCrystals(0),fGeom(0),fPulseGenerator(0){ + +} +//-------------------------------------------------------------------------------------- +AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRecoParam* parEmc, const AliPHOSRecoParam* parCpv):TObject(), + fEmcMinE(0.),fCpvMinE(0.),fSampleQualityCut(1.),fEmcCrystals(0),fGeom(0),fPulseGenerator(0){ + + if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!"); + if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!"); + + fEmcMinE = parEmc->GetMinE(); + fCpvMinE = parCpv->GetMinE(); + + fSampleQualityCut = parEmc->GetSampleQualityCut() ; + + fGeom=AliPHOSGeometry::GetInstance() ; + if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP"); + + fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ; + + fPulseGenerator = new AliPHOSPulseGenerator(); + + GetCalibrationParameters() ; +} +//-------------------------------------------------------------------------------------- +AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp):TObject(), + fEmcMinE(0.),fCpvMinE(0.),fSampleQualityCut(1.),fEmcCrystals(0),fGeom(0),fPulseGenerator(0){ + + fEmcMinE = dp.fEmcMinE ; + fCpvMinE = dp.fCpvMinE ; + fEmcCrystals = dp.fEmcCrystals ; + fPulseGenerator = new AliPHOSPulseGenerator(); + fGeom = dp.fGeom ; +} +//-------------------------------------------------------------------------------------- +AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiProducer &dp){ + if(&dp == this) return *this; + + fEmcMinE = dp.fEmcMinE ; + fCpvMinE = dp.fCpvMinE ; + fSampleQualityCut = dp.fSampleQualityCut ; + fEmcCrystals = dp.fEmcCrystals ; + fGeom = dp.fGeom ; + if(fPulseGenerator) delete fPulseGenerator ; + fPulseGenerator = new AliPHOSPulseGenerator(); + return *this; +} +//-------------------------------------------------------------------------------------- +AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer(){ + if(fPulseGenerator) delete fPulseGenerator ; + fPulseGenerator=0 ; +} //-------------------------------------------------------------------------------------- void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* decoder) { @@ -48,67 +109,209 @@ void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* digits->Clear(); - AliPHOSGeometry* geo = AliPHOSGeometry::GetInstance(); - if(!geo) geo = AliPHOSGeometry::GetInstance("IHEP"); - Int_t iDigit = 0 ; - Double_t time = 0. ; - Int_t iOldDigit; - Bool_t seen,lowGainFlag; Int_t relId[4], absId =0; + const Double_t baseLine=1. ; //Minimal energy of digit in ADC ch. + const Double_t highLowDiff=2.; //Maximal difference between High and Low channels in LG adc channels + + //Temporary array for LowGain digits + TClonesArray tmpLG("AliPHOSDigit",10000) ; + Int_t ilgDigit=0 ; + while (decoder->NextDigit()) { - if (decoder->GetEnergy() <= 0.) continue; + Double_t energy=decoder->GetEnergy() ; + if(energy<=baseLine) //in ADC channels + continue ; + + //remove digits with bas shape. Decoder should calculate quality so that + //in default case quality [0,1], while larger values of quality mean somehow + //corrupted samples, 999 means obviously corrupted sample. + //It is difficult to fit samples with overflow (even setting cut on overflow values) + //because too few points are left to fit. So we do not evaluate samples with overflow + if(decoder->GetSampleQuality() > fSampleQualityCut && !(decoder->IsOverflow())) + continue ; - lowGainFlag = decoder->IsLowGain(); - time = decoder->GetTime(); + Bool_t lowGainFlag = decoder->IsLowGain(); relId[0] = decoder->GetModule(); relId[1] = 0; relId[2] = decoder->GetRow(); relId[3] = decoder->GetColumn(); - geo->RelToAbsNumbering(relId, absId); + fGeom->RelToAbsNumbering(relId, absId); - // Add low gain digit only - //if the high gain digit does not exist in the digits array + Double_t time = decoder->GetTime() ; + time = CalibrateT(time,relId,lowGainFlag) ; - seen = kFALSE; + energy = CalibrateE(energy,relId,lowGainFlag) ; - if(lowGainFlag) { - for (iOldDigit=iDigit-1; iOldDigit>=0; iOldDigit--) { - if ((dynamic_cast(digits->At(iOldDigit)))->GetId() == absId) { - seen = kTRUE; - break; + if(energy <= 0.) + continue; + + if(lowGainFlag){ + new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time); + ilgDigit++ ; + } + else{ + if(decoder->IsOverflow()) //Keep this digit to replace it by Low Gain later. + //If there is no LogGain it wil be removed by cut on Min E + new((*digits)[iDigit]) AliPHOSDigit(-1,absId,-1.f,(Float_t)time); + else + new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time); + iDigit++; + } + } + + //Now scan created LG and HG digits and keep only those which appeared in both lists + //replace energy of HighGain digits only if there is overflow + //negative energy (overflow) + digits->Sort() ; + tmpLG.Sort() ; + Int_t iLG = 0; + Int_t nLG1 = tmpLG.GetEntriesFast()-1 ; + + for(Int_t iDig=0 ; iDig < digits->GetEntriesFast() ; iDig++) { + AliPHOSDigit * digHG = dynamic_cast(digits->At(iDig)) ; + if (!digHG) continue; + AliPHOSDigit * digLG = dynamic_cast(tmpLG.At(iLG)) ; + while(digLG && iLGGetId()> digLG->GetId()){ + iLG++ ; + digLG = dynamic_cast(tmpLG.At(iLG)) ; + } + Int_t absId=digHG->GetId() ; + Int_t relId[4] ; + fGeom->AbsToRelNumbering(absId,relId) ; + + if(digLG && digHG->GetId() == digLG->GetId()){ //we found pair + if(digHG->GetEnergy()<0.){ //This is overflow in HG + digHG->SetTime(digLG->GetTime()) ; + digHG->SetEnergy(digLG->GetEnergy()) ; + } + else{ //Make approximate comparison of HG and LG energies + Double_t de = (digHG->GetEnergy()-digLG->GetEnergy()) ; + if(TMath::Abs(de)>CalibrateE(double(highLowDiff),relId,1)){ //too strong difference, remove digit + digits->RemoveAt(iDig) ; } } - if (!seen) { - new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)decoder->GetEnergy(),time); - iDigit++; - } } + else{ //no pair - remove + // temporary fix for dead LG channels + if(relId[2]%2==1 && relId[3]%16==4) + continue ; + if(digHG->GetEnergy()>CalibrateE(double(5.),relId,1)) //One can not always find LG with Amp<5 ADC ch. + digits->RemoveAt(iDig) ; + } + } - // Add high gain digit only if it is not saturated; - // replace low gain digit by a high gain one - else { - if (decoder->GetEnergy() >= 1023) continue; - for (iOldDigit=iDigit-1; iOldDigit>=0; iOldDigit--) { - AliPHOSDigit * dig = dynamic_cast(digits->At(iOldDigit)) ; - if (dig->GetId() == absId) { - dig->SetEnergy((Float_t)decoder->GetEnergy()) ; - dig->SetTime(time) ; - seen = kTRUE; - break; - } + CleanDigits(digits) ; + +} +//____________________________________________________________________________ +Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain) +{ + // Convert EMC measured amplitude into real energy. + // Calibration parameters are taken from calibration data base for raw data, + // or from digitizer parameters for simulated data. + if(fgCalibData){ + Int_t module = relId[0]; + Int_t column = relId[3]; + Int_t row = relId[2]; + if(relId[1]==0) { // this is EMC + if(isLowGain){ + amp*= fgCalibData->GetHighLowRatioEmc(module,column,row); } - if (!seen) { - new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)decoder->GetEnergy(),time); - iDigit++; + amp *= fgCalibData->GetADCchannelEmc(module,column,row); + return amp ; + } + } + return 0; +} +//____________________________________________________________________________ +Double_t AliPHOSRawDigiProducer::CalibrateT(Double_t time, Int_t * relId, Bool_t /* isLowGain */) +{ + //Calibrate time + time*=fPulseGenerator->GetRawFormatTimeTrigger() ; + if(fgCalibData){ + Int_t module = relId[0]; + Int_t column = relId[3]; + Int_t row = relId[2]; + if(relId[1]==0) { // this is EMC + time += fgCalibData->GetTimeShiftEmc(module,column,row); + return time ; + } + } + + return -999.; +} +//____________________________________________________________________________ +void AliPHOSRawDigiProducer::CleanDigits(TClonesArray * digits) +{ + // remove digits with amplitudes below threshold. + // remove digits in bad channels + // sort digits with icreasing AbsId + + //remove digits in bad map and below threshold + Bool_t isBadMap = 0 ; + if(fgCalibData->GetNumOfEmcBadChannels()){ + isBadMap=1 ; + } + + for(Int_t i=0; iGetEntriesFast(); i++){ + AliPHOSDigit * digit = static_cast(digits->At(i)) ; + if(!digit) + continue ; + if ( (IsInEMC(digit) && digit->GetEnergy() < fEmcMinE) || + (IsInCPV(digit) && digit->GetEnergy() < fCpvMinE) ){ + digits->RemoveAt(i) ; + continue ; + } + if(isBadMap){ //check bad map now + Int_t relid[4] ; + fGeom->AbsToRelNumbering(digit->GetId(), relid) ; + if(fgCalibData->IsBadChannelEmc(relid[0],relid[3],relid[2])){ + digits->RemoveAt(i) ; } } - } - digits->Compress(); - digits->Sort(); + //Compress, sort and set indexes + digits->Compress() ; +// digits->Sort(); already sorted earlier + for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) { + AliPHOSDigit *digit = static_cast( digits->At(i) ) ; + digit->SetIndexInList(i) ; + } +} +//____________________________________________________________________________ +Bool_t AliPHOSRawDigiProducer::IsInEMC(AliPHOSDigit * digit) const +{ + // Tells if (true) or not (false) the digit is in a PHOS-EMC module + return digit->GetId() <= fEmcCrystals ; + +} + +//____________________________________________________________________________ +Bool_t AliPHOSRawDigiProducer::IsInCPV(AliPHOSDigit * digit) const +{ + // Tells if (true) or not (false) the digit is in a PHOS-CPV module + return digit->GetId() > fEmcCrystals ; +} +//____________________________________________________________________________ +void AliPHOSRawDigiProducer::GetCalibrationParameters() +{ + // Set calibration parameters: + // if calibration database exists, they are read from database, + // otherwise, reconstruction stops in the constructor of AliPHOSCalibData + // + // It is a user responsilibity to open CDB before reconstruction, for example: + // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); + + if (!fgCalibData){ + fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number + } + if (fgCalibData->GetCalibDataEmc() == 0) + AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n"); + if (fgCalibData->GetCalibDataCpv() == 0) + AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n"); }