From 379c5c09ee1298dca59230016610e34ff4fe2cb2 Mon Sep 17 00:00:00 2001 From: kharlov Date: Fri, 3 Jul 2009 17:03:50 +0000 Subject: [PATCH] Raw data decoder is migrated from AliCaloRawStream to AliCaloRawStreamV3 --- PHOS/AliPHOS.cxx | 15 +- PHOS/AliPHOSQADataMakerRec.cxx | 244 ++++++++++------- PHOS/AliPHOSRawDigiProducer.cxx | 233 ++++++++++------- PHOS/AliPHOSRawDigiProducer.h | 19 +- PHOS/AliPHOSRawFitterv0.cxx | 220 ++++++++++++++++ PHOS/AliPHOSRawFitterv0.h | 69 +++++ PHOS/AliPHOSRawFitterv1.cxx | 445 ++++++++++++++++++++++++++++++++ PHOS/AliPHOSRawFitterv1.h | 42 +++ PHOS/AliPHOSRawFitterv2.cxx | 251 ++++++++++++++++++ PHOS/AliPHOSRawFitterv2.h | 40 +++ PHOS/AliPHOSRecoParam.cxx | 10 +- PHOS/AliPHOSRecoParam.h | 6 +- PHOS/AliPHOSReconstructor.cxx | 77 +++--- PHOS/CMake_libPHOSbase.txt | 6 +- PHOS/PHOSbaseLinkDef.h | 6 +- PHOS/libPHOSbase.pkg | 6 +- 16 files changed, 1425 insertions(+), 264 deletions(-) create mode 100644 PHOS/AliPHOSRawFitterv0.cxx create mode 100644 PHOS/AliPHOSRawFitterv0.h create mode 100644 PHOS/AliPHOSRawFitterv1.cxx create mode 100644 PHOS/AliPHOSRawFitterv1.h create mode 100644 PHOS/AliPHOSRawFitterv2.cxx create mode 100644 PHOS/AliPHOSRawFitterv2.h diff --git a/PHOS/AliPHOS.cxx b/PHOS/AliPHOS.cxx index 7cbced51f14..efe95165f46 100644 --- a/PHOS/AliPHOS.cxx +++ b/PHOS/AliPHOS.cxx @@ -94,8 +94,7 @@ class TFile; #include "AliPHOSCalibData.h" #include "AliPHOSPulseGenerator.h" #include "AliDAQ.h" -#include "AliPHOSRawDecoder.h" -#include "AliPHOSRawDecoderv1.h" +#include "AliPHOSRawFitterv0.h" #include "AliPHOSCalibData.h" #include "AliPHOSRawDigiProducer.h" #include "AliPHOSQAChecker.h" @@ -658,15 +657,15 @@ Bool_t AliPHOS::Raw2SDigits(AliRawReader* rawReader) mapping[i] = (AliAltroMapping*)maps->At(i); } - AliPHOSRawDecoderv1 dc(rawReader,mapping); + AliPHOSRawFitterv0 fitter; - dc.SubtractPedestals(AliPHOSSimParam::GetInstance()->EMCSubtractPedestals()); - dc.SetAmpOffset(AliPHOSSimParam::GetInstance()->GetGlobalAltroOffset()); - dc.SetAmpThreshold(AliPHOSSimParam::GetInstance()->GetGlobalAltroThreshold()); + fitter.SubtractPedestals(AliPHOSSimParam::GetInstance()->EMCSubtractPedestals()); + fitter.SetAmpOffset(AliPHOSSimParam::GetInstance()->GetGlobalAltroOffset()); + fitter.SetAmpThreshold(AliPHOSSimParam::GetInstance()->GetGlobalAltroThreshold()); - AliPHOSRawDigiProducer pr; + AliPHOSRawDigiProducer pr(rawReader,mapping); pr.SetSampleQualityCut(AliPHOSSimParam::GetInstance()->GetEMCSampleQualityCut()); - pr.MakeDigits(sdigits,&dc); + pr.MakeDigits(sdigits,&fitter); Int_t bufferSize = 32000 ; // TBranch * sdigitsBranch = tree->Branch("PHOS",&sdigits,bufferSize); diff --git a/PHOS/AliPHOSQADataMakerRec.cxx b/PHOS/AliPHOSQADataMakerRec.cxx index 934f8af0129..54c821eff73 100644 --- a/PHOS/AliPHOSQADataMakerRec.cxx +++ b/PHOS/AliPHOSQADataMakerRec.cxx @@ -33,20 +33,22 @@ // --- Standard library --- // --- AliRoot header files --- +#include "AliCaloRawStreamV3.h" #include "AliESDCaloCluster.h" #include "AliESDEvent.h" #include "AliLog.h" #include "AliPHOSQADataMakerRec.h" -#include "AliQAChecker.h" #include "AliPHOSDigit.h" #include "AliPHOSCpvRecPoint.h" #include "AliPHOSEmcRecPoint.h" +#include "AliPHOSRawFitterv0.h" +#include "AliPHOSRawFitterv1.h" +#include "AliPHOSRawFitterv2.h" +#include "AliPHOSReconstructor.h" #include "AliPHOSRecParticle.h" #include "AliPHOSTrackSegment.h" -#include "AliPHOSRawDecoder.h" -#include "AliPHOSRawDecoderv1.h" -#include "AliPHOSRawDecoderv2.h" -#include "AliPHOSReconstructor.h" +#include "AliQAChecker.h" +#include "AliRawReader.h" ClassImp(AliPHOSQADataMakerRec) @@ -374,111 +376,153 @@ void AliPHOSQADataMakerRec::MakeRaws(AliRawReader* rawReader) //Fill prepared histograms with Raw digit properties rawReader->Reset() ; - AliPHOSRawDecoder * decoder ; - if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0) - decoder=new AliPHOSRawDecoderv1(rawReader); - else - if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v2")==0) - decoder=new AliPHOSRawDecoderv2(rawReader); + + const TObjArray* maps = AliPHOSRecoParam::GetMappings(); + if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); + + AliAltroMapping *mapping[4]; + for(Int_t i = 0; i < 4; i++) { + mapping[i] = (AliAltroMapping*)maps->At(i); + } + + AliCaloRawStreamV3 *fRawStream = new AliCaloRawStreamV3(rawReader,"PHOS",mapping); + + AliPHOSRawFitterv0 * fitter ; + if (strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0) + fitter=new AliPHOSRawFitterv1(); + else if(strcmp(GetRecoParam()->EMCFitterVersion(),"v2")==0) + fitter=new AliPHOSRawFitterv2(); else - decoder=new AliPHOSRawDecoder(rawReader); - //decoder->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals()); - decoder->SubtractPedestals(kTRUE); + fitter=new AliPHOSRawFitterv0(); + fitter->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals()); +// fitter->SubtractPedestals(kTRUE); Double_t lgEtot=0. ; Double_t hgEtot=0. ; - Int_t lgNtot=0 ; - Int_t hgNtot=0 ; - - while (decoder->NextDigit()) { - Int_t module = decoder->GetModule() ; - Int_t row = decoder->GetRow() ; - Int_t col = decoder->GetColumn() ; - Double_t time = decoder->GetTime() ; - Double_t energy = decoder->GetEnergy() ; - Bool_t lowGain = decoder->IsLowGain(); - if (lowGain) { - //if(GetRecoParam()->EMCSubtractPedestals()) - GetRawsData(kLGpedRMS)->Fill(decoder->GetPedestalRMS()) ; - if(energy<2.) - continue ; - switch(module){ - case 1: GetRawsData(kLGmod1)->Fill(row-0.5,col-0.5) ; break ; - case 2: GetRawsData(kLGmod2)->Fill(row-0.5,col-0.5) ; break ; - case 3: GetRawsData(kLGmod3)->Fill(row-0.5,col-0.5) ; break ; - case 4: GetRawsData(kLGmod4)->Fill(row-0.5,col-0.5) ; break ; - case 5: GetRawsData(kLGmod5)->Fill(row-0.5,col-0.5) ; break ; - } - switch (module){ - case 1: ((TH2F*)GetRawsData(kLGpedRMSMod1))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 2: ((TH2F*)GetRawsData(kLGpedRMSMod2))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 3: ((TH2F*)GetRawsData(kLGpedRMSMod3))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 4: ((TH2F*)GetRawsData(kLGpedRMSMod4))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 5: ((TH2F*)GetRawsData(kLGpedRMSMod5))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - } - //if quality was evaluated, fill histo - if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0){ - switch (module){ - case 1: ((TH2F*)GetRawsData(kLGqualMod1))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 2: ((TH2F*)GetRawsData(kLGqualMod2))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 3: ((TH2F*)GetRawsData(kLGqualMod3))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 4: ((TH2F*)GetRawsData(kLGqualMod4))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 5: ((TH2F*)GetRawsData(kLGqualMod5))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - } - } - GetRawsData(kNmodLG)->Fill(module) ; - GetRawsData(kLGtime)->Fill(time) ; - GetRawsData(kSpecLG)->Fill(energy) ; - lgEtot+=energy ; - lgNtot++ ; - } else { - //if this isnon-ZS run - fill pedestal RMS - //if(GetRecoParam()->EMCSubtractPedestals()) - GetRawsData(kHGpedRMS)->Fill(decoder->GetPedestalRMS()) ; - if(energy<8.) - continue ; - switch (module){ - case 1: GetRawsData(kHGmod1)->Fill(row-0.5,col-0.5) ; break ; - case 2: GetRawsData(kHGmod2)->Fill(row-0.5,col-0.5) ; break ; - case 3: GetRawsData(kHGmod3)->Fill(row-0.5,col-0.5) ; break ; - case 4: GetRawsData(kHGmod4)->Fill(row-0.5,col-0.5) ; break ; - case 5: GetRawsData(kHGmod5)->Fill(row-0.5,col-0.5) ; break ; - } - switch (module){ - case 1: ((TH2F*)GetRawsData(kHGpedRMSMod1))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 2: ((TH2F*)GetRawsData(kHGpedRMSMod2))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 3: ((TH2F*)GetRawsData(kHGpedRMSMod3))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 4: ((TH2F*)GetRawsData(kHGpedRMSMod4))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - case 5: ((TH2F*)GetRawsData(kHGpedRMSMod5))->Fill(row-0.5,col-0.5,decoder->GetPedestalRMS()) ; break ; - } - //if quality was evaluated, fill histo - if(strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0){ - switch (module){ - case 1: ((TH2F*)GetRawsData(kHGqualMod1))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 2: ((TH2F*)GetRawsData(kHGqualMod2))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 3: ((TH2F*)GetRawsData(kHGqualMod3))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 4: ((TH2F*)GetRawsData(kHGqualMod4))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - case 5: ((TH2F*)GetRawsData(kHGqualMod5))->Fill(row-0.5,col-0.5,decoder->GetSampleQuality()) ; break ; - } - } - GetRawsData(kNmodHG)->Fill(module) ; - GetRawsData(kHGtime)->Fill(time) ; - GetRawsData(kSpecHG)->Fill(energy) ; - hgEtot+=energy ; - hgNtot++ ; - } - } - delete decoder; + Int_t lgNtot=0 ; + Int_t hgNtot=0 ; + + while (fRawStream->NextDDL()) { // !!!!!!!!!!! YK + while (fRawStream->NextChannel()) { + Int_t module = fRawStream->GetModule(); + Int_t cellX = fRawStream->GetCellX(); + Int_t cellZ = fRawStream->GetCellZ(); + Int_t caloFlag = fRawStream->GetCaloFlag(); // 0=LG, 1=HG, 2=TRU + + Int_t nBunches = 0; + while (fRawStream->NextBunch()) { + nBunches++; + if (nBunches > 1) continue; + const UShort_t *sig = fRawStream->GetSignals(); + Int_t sigStart = fRawStream->GetStartTimeBin(); + Int_t sigLength = fRawStream->GetBunchLength(); + fitter->SetSamples(sig,sigStart,sigLength); + } // End of NextBunch() + + fitter->SetNBunches(nBunches); + fitter->SetChannelGeo(module,cellX,cellZ,caloFlag); + fitter->Eval(); + + Double_t energy = fitter->GetEnergy() ; + Double_t time = fitter->GetTime() ; + + if (caloFlag == 0) { // LG + //if(GetRecoParam()->EMCSubtractPedestals()) + GetRawsData(kLGpedRMS)->Fill(fitter->GetPedestalRMS()) ; + if(energy<2.) + continue ; + switch(module){ + case 1: GetRawsData(kLGmod1)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 2: GetRawsData(kLGmod2)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 3: GetRawsData(kLGmod3)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 4: GetRawsData(kLGmod4)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 5: GetRawsData(kLGmod5)->Fill(cellX-0.5,cellZ-0.5) ; break ; + } + switch (module){ + case 1: ((TH2F*)GetRawsData(kLGpedRMSMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 2: ((TH2F*)GetRawsData(kLGpedRMSMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 3: ((TH2F*)GetRawsData(kLGpedRMSMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 4: ((TH2F*)GetRawsData(kLGpedRMSMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 5: ((TH2F*)GetRawsData(kLGpedRMSMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + } + //if quality was evaluated, fill histo + if(strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0){ + switch (module){ + case 1: ((TH2F*)GetRawsData(kLGqualMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 2: ((TH2F*)GetRawsData(kLGqualMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 3: ((TH2F*)GetRawsData(kLGqualMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 4: ((TH2F*)GetRawsData(kLGqualMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 5: ((TH2F*)GetRawsData(kLGqualMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + } + } + GetRawsData(kNmodLG)->Fill(module) ; + GetRawsData(kLGtime)->Fill(time) ; + GetRawsData(kSpecLG)->Fill(energy) ; + lgEtot+=energy ; + lgNtot++ ; + } else if (caloFlag == 1) { // HG + //if this isnon-ZS run - fill pedestal RMS + //if(GetRecoParam()->EMCSubtractPedestals()) + GetRawsData(kHGpedRMS)->Fill(fitter->GetPedestalRMS()) ; + if(energy<8.) + continue ; + switch (module){ + case 1: GetRawsData(kHGmod1)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 2: GetRawsData(kHGmod2)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 3: GetRawsData(kHGmod3)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 4: GetRawsData(kHGmod4)->Fill(cellX-0.5,cellZ-0.5) ; break ; + case 5: GetRawsData(kHGmod5)->Fill(cellX-0.5,cellZ-0.5) ; break ; + } + switch (module){ + case 1: ((TH2F*)GetRawsData(kHGpedRMSMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 2: ((TH2F*)GetRawsData(kHGpedRMSMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 3: ((TH2F*)GetRawsData(kHGpedRMSMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 4: ((TH2F*)GetRawsData(kHGpedRMSMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + case 5: ((TH2F*)GetRawsData(kHGpedRMSMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetPedestalRMS()) ; break ; + } + //if quality was evaluated, fill histo + if(strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0){ + switch (module){ + case 1: ((TH2F*)GetRawsData(kHGqualMod1))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 2: ((TH2F*)GetRawsData(kHGqualMod2))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 3: ((TH2F*)GetRawsData(kHGqualMod3))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 4: ((TH2F*)GetRawsData(kHGqualMod4))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + case 5: ((TH2F*)GetRawsData(kHGqualMod5))->Fill(cellX-0.5,cellZ-0.5,fitter->GetSignalQuality()) ; break ; + } + } + GetRawsData(kNmodHG)->Fill(module) ; + GetRawsData(kHGtime)->Fill(time) ; + GetRawsData(kSpecHG)->Fill(energy) ; + hgEtot+=energy ; + hgNtot++ ; + } + } // End of NextChannel + } // End of NextDDL + delete fitter; + GetRawsData(kEtotLG)->Fill(lgEtot) ; - TParameter * p = dynamic_cast*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kEtotLG)->GetName()))) ; + TParameter * p; + p = dynamic_cast*>(GetParameterList()-> + FindObject(Form("%s_%s_%s", GetName(), + AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), + GetRawsData(kEtotLG)->GetName()))) ; if (p) p->SetVal(lgEtot) ; GetRawsData(kEtotHG)->Fill(hgEtot) ; - p = dynamic_cast*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kEtotHG)->GetName()))) ; + p = dynamic_cast*>(GetParameterList()-> + FindObject(Form("%s_%s_%s", GetName(), + AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), + GetRawsData(kEtotHG)->GetName()))) ; if (p) p->SetVal(hgEtot) ; GetRawsData(kNtotLG)->Fill(lgNtot) ; - p = dynamic_cast*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kNtotLG)->GetName()))) ; + p = dynamic_cast*>(GetParameterList()-> + FindObject(Form("%s_%s_%s", GetName(), + AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), + GetRawsData(kNtotLG)->GetName()))) ; if (p) p->SetVal(lgNtot) ; GetRawsData(kNtotHG)->Fill(hgNtot) ; - p = dynamic_cast*>(GetParameterList()->FindObject(Form("%s_%s_%s", GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), GetRawsData(kNtotHG)->GetName()))) ; + p = dynamic_cast*>(GetParameterList()-> + FindObject(Form("%s_%s_%s", + GetName(), AliQAv1::GetTaskName(AliQAv1::kRAWS).Data(), + GetRawsData(kNtotHG)->GetName()))) ; if (p) p->SetVal(hgNtot) ; } diff --git a/PHOS/AliPHOSRawDigiProducer.cxx b/PHOS/AliPHOSRawDigiProducer.cxx index ea39f00512c..5ee95aa6015 100644 --- a/PHOS/AliPHOSRawDigiProducer.cxx +++ b/PHOS/AliPHOSRawDigiProducer.cxx @@ -16,7 +16,7 @@ /* $Id$ */ //This class produces PHOS digits of one event -//using AliPHOSRawDecoder. +//using AliPHOSRawFitter. // // For example: // TClonesArray *digits = new TClonesArray("AliPHOSDigit",100); @@ -34,11 +34,12 @@ // --- AliRoot header files --- #include "AliPHOSRawDigiProducer.h" -#include "AliPHOSRawDecoder.h" +#include "AliPHOSRawFitterv0.h" #include "AliPHOSGeometry.h" #include "AliPHOSDigit.h" #include "AliPHOSCalibData.h" #include "AliPHOSPulseGenerator.h" +#include "AliCaloRawStreamV3.h" #include "AliLog.h" ClassImp(AliPHOSRawDigiProducer) @@ -53,7 +54,9 @@ AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(): fSampleQualityCut(1.), fEmcCrystals(0), fGeom(0), - fPulseGenerator(0) + fPulseGenerator(0), + fRawReader(0), + fRawStream(0) { // Default constructor @@ -65,7 +68,32 @@ AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(): GetCalibrationParameters() ; } -//-------------------------------------------------------------------------------------- +//-------------------------------------------------------------------------------------- +AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(AliRawReader *rawReader, + AliAltroMapping **mapping): + TObject(), + fEmcMinE(0.), + fCpvMinE(0.), + fSampleQualityCut(1.), + fEmcCrystals(0), + fGeom(0), + fPulseGenerator(0), + fRawReader(rawReader), + fRawStream(0) +{ + // Default constructor + + fGeom=AliPHOSGeometry::GetInstance() ; + if(!fGeom) fGeom = AliPHOSGeometry::GetInstance("IHEP"); + + fEmcCrystals=fGeom->GetNCristalsInModule()*fGeom->GetNModules() ; + fPulseGenerator = new AliPHOSPulseGenerator(); + GetCalibrationParameters() ; + + fRawStream = new AliCaloRawStreamV3(rawReader,"PHOS",mapping); + +} +//-------------------------------------------------------------------------------------- AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp): TObject(), fEmcMinE(0.), @@ -73,7 +101,9 @@ AliPHOSRawDigiProducer::AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp) fSampleQualityCut(1.), fEmcCrystals(0), fGeom(0), - fPulseGenerator(0) + fPulseGenerator(0), + fRawReader(0), + fRawStream(0) { // Copy constructor @@ -100,23 +130,24 @@ AliPHOSRawDigiProducer& AliPHOSRawDigiProducer::operator= (const AliPHOSRawDigiP fPulseGenerator = new AliPHOSPulseGenerator(); return *this; } -//-------------------------------------------------------------------------------------- +//-------------------------------------------------------------------------------------- AliPHOSRawDigiProducer::~AliPHOSRawDigiProducer() { - // Desctructor + // Destructor if(fPulseGenerator) delete fPulseGenerator ; fPulseGenerator=0 ; + delete fRawStream; } //-------------------------------------------------------------------------------------- -void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* decoder) +void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawFitterv0* fitter) { //Makes the job. - //TClonesArray *digits and raw data decoder should be provided by calling function. + //TClonesArray *digits and raw data fitter should be provided by calling function. digits->Clear(); - Int_t iDigit = 0 ; - Int_t relId[4], absId =0; + Int_t iDigit=0 ; + Int_t relId[4], absId=-1, caloFlag=-1; 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 @@ -125,95 +156,109 @@ void AliPHOSRawDigiProducer::MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* TClonesArray tmpLG("AliPHOSDigit",10000) ; Int_t ilgDigit=0 ; - //Let decoder subtract pedestals in case of ZS - decoder->SetCalibData(fgCalibData) ; + //Let fitter subtract pedestals in case of ZS + fitter->SetCalibData(fgCalibData) ; - while (decoder->NextDigit()) { - - Double_t energy=decoder->GetEnergy() ; - if(energy<=baseLine) //in ADC channels - continue ; - - //remove digits with bad 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 ; - - Bool_t lowGainFlag = decoder->IsLowGain(); - - relId[0] = decoder->GetModule(); - relId[1] = 0; - relId[2] = decoder->GetRow(); - relId[3] = decoder->GetColumn(); - fGeom->RelToAbsNumbering(relId, absId); - - Double_t time = decoder->GetTime() ; - time = CalibrateT(time,relId,lowGainFlag) ; - - energy = CalibrateE(energy,relId,lowGainFlag) ; - - 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)) ; - } - absId=digHG->GetId() ; - fGeom->AbsToRelNumbering(absId,relId) ; + while (fRawStream->NextDDL()) { + while (fRawStream->NextChannel()) { + relId[0] = fRawStream->GetModule(); + relId[1] = 0; + relId[2] = fRawStream->GetCellX(); + relId[3] = fRawStream->GetCellZ(); + caloFlag = fRawStream->GetCaloFlag(); // 0=LG, 1=HG, 2=TRU + fGeom->RelToAbsNumbering(relId, absId); + + Int_t nBunches = 0; + while (fRawStream->NextBunch()) { + nBunches++; + if (nBunches > 1) continue; + const UShort_t *sig = fRawStream->GetSignals(); + Int_t sigStart = fRawStream->GetStartTimeBin(); + Int_t sigLength = fRawStream->GetBunchLength(); + fitter->SetSamples(sig,sigStart,sigLength); + } // End of NextBunch() + + fitter->SetNBunches(nBunches); + fitter->SetChannelGeo(relId[0],relId[2],relId[3],caloFlag); + fitter->Eval(); + + Double_t energy = fitter->GetEnergy() ; + Double_t time = fitter->GetTime() ; + if(energy<=baseLine) //in ADC channels + continue ; + + //remove digits with bad shape. Fitter 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(fitter->GetSignalQuality() > fSampleQualityCut && !(fitter->IsOverflow())) + continue ; + +// energy = CalibrateE(energy,relId,lowGainFlag) ; +// time = CalibrateT(time,relId,lowGainFlag) ; + + if(energy <= 0.) + continue; + + if (caloFlag == AliCaloRawStreamV3::kLowGain) { + new(tmpLG[ilgDigit]) AliPHOSDigit(-1,absId,(Float_t)energy,(Float_t)time); + ilgDigit++ ; + } + else if (caloFlag == AliCaloRawStreamV3::kHighGain) { + if(fitter->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++; + } + } // End of NextChannel() + + //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)) ; + } + absId=digHG->GetId() ; + 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(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) ; + } + } + } + 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) ; } } - 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) ; - } - } + } // End of NextDDL() CleanDigits(digits) ; - + } //____________________________________________________________________________ Double_t AliPHOSRawDigiProducer::CalibrateE(Double_t amp, Int_t* relId, Bool_t isLowGain) diff --git a/PHOS/AliPHOSRawDigiProducer.h b/PHOS/AliPHOSRawDigiProducer.h index 183c4db87c2..352641ce10d 100644 --- a/PHOS/AliPHOSRawDigiProducer.h +++ b/PHOS/AliPHOSRawDigiProducer.h @@ -6,13 +6,17 @@ /* $Id$ */ //This class produces PHOS digits of one event -//using AliPHOSRawDecoder. See cxx source for use case. +//using AliPHOSRawFitter. See cxx source for use case. -class AliPHOSRawDecoder; class AliPHOSCalibData ; class AliPHOSDigit ; class AliPHOSGeometry ; class AliPHOSPulseGenerator; +class AliRawReader; +class AliCaloRawStreamV3; +class AliPHOSRawFitterv0; + +#include "AliAltroMapping.h" #include "TObject.h" class AliPHOSRawDigiProducer: public TObject { @@ -20,12 +24,13 @@ class AliPHOSRawDigiProducer: public TObject { public: AliPHOSRawDigiProducer() ; + AliPHOSRawDigiProducer(AliRawReader *rawReader, AliAltroMapping **mapping = NULL); AliPHOSRawDigiProducer(const AliPHOSRawDigiProducer &dp); AliPHOSRawDigiProducer& operator= (const AliPHOSRawDigiProducer &dp); virtual ~AliPHOSRawDigiProducer(); - void MakeDigits(TClonesArray *digits, AliPHOSRawDecoder* decoder); + void MakeDigits(TClonesArray *digits, AliPHOSRawFitterv0* fitter); void SetEmcMinAmp(Float_t emcMin) { fEmcMinE=emcMin; } void SetCpvMinAmp(Float_t cpvMin) { fCpvMinE=cpvMin; } @@ -48,10 +53,12 @@ private: Float_t fSampleQualityCut; // Cut on sample shapes: 0: no samples; 1: default parameterization; 999: accept even obviously bad Int_t fEmcCrystals ; // number of EMC crystals AliPHOSGeometry * fGeom ; //! PHOS geometry - static AliPHOSCalibData * fgCalibData ; //! Calibration database if avalable - AliPHOSPulseGenerator * fPulseGenerator ; //! Class with pulse shape parameters + static AliPHOSCalibData * fgCalibData ; //! Calibration database if avalable + AliPHOSPulseGenerator * fPulseGenerator ; //! Class with pulse shape parameters + AliRawReader * fRawReader; //! Raw data reader + AliCaloRawStreamV3 * fRawStream; //! Calorimeter decoder of ALTRO format - ClassDef(AliPHOSRawDigiProducer,5) + ClassDef(AliPHOSRawDigiProducer,6) }; #endif diff --git a/PHOS/AliPHOSRawFitterv0.cxx b/PHOS/AliPHOSRawFitterv0.cxx new file mode 100644 index 00000000000..e1d112544ee --- /dev/null +++ b/PHOS/AliPHOSRawFitterv0.cxx @@ -0,0 +1,220 @@ +/************************************************************************** + * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id: $ */ + +// This class extracts the signal parameters (energy, time, quality) +// from ALTRO samples. Energy is in ADC counts, time is in time bin units. +// A coarse algorithm takes the energy as the maximum +// sample, time is the first sample index, and the quality is the number +// of bunches in the signal. +// +// AliPHOSRawFitterv0 *fitterv0=new AliPHOSRawFitterv0(); +// fitterv0->SetSamples(sig,sigStart,sigLength); +// fitterv0->SetNBunches(nBunches); +// fitterv0->SetChannelGeo(module,cellX,cellZ,caloFlag); +// fitterv0->SetCalibData(fgCalibData) ; +// fitterv0->Eval(); +// Double_t amplitude = fitterv0.GetEnergy(); +// Double_t time = fitterv0.GetTime(); +// Bool_t isLowGain = fitterv0.GetCaloFlag()==0; + +// Author: Yuri Kharlov + +// --- ROOT system --- +#include "TArrayI.h" +#include "TMath.h" +#include "TObject.h" + +// --- AliRoot header files --- +#include "AliPHOSRawFitterv0.h" +#include "AliPHOSCalibData.h" +#include "AliLog.h" + +ClassImp(AliPHOSRawFitterv0) + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv0::AliPHOSRawFitterv0(): + TObject(), + fSignal(0), + fModule(0), + fCellX(0), + fCellZ(0), + fCaloFlag(0), + fStart(0), + fLength(0), + fNBunches(0), + fPedSubtract(kFALSE), + fEnergy(-111), + fTime(-111), + fQuality(0.), + fPedestalRMS(0.), + fAmpOffset(0), + fAmpThreshold(0), + fOverflow(kFALSE), + fCalibData(0) +{ + //Default constructor +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv0::~AliPHOSRawFitterv0() +{ + //Destructor + delete [] fSignal; +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv0::AliPHOSRawFitterv0(const AliPHOSRawFitterv0 &phosFitter ): + TObject(), + fSignal (phosFitter.fSignal), + fModule (phosFitter.fModule), + fCellX (phosFitter.fCellX), + fCellZ (phosFitter.fCellZ), + fCaloFlag (phosFitter.fCaloFlag), + fStart (phosFitter.fStart), + fLength (phosFitter.fLength), + fNBunches (phosFitter.fNBunches), + fPedSubtract (phosFitter.fPedSubtract), + fEnergy (phosFitter.fEnergy), + fTime (phosFitter.fTime), + fQuality (phosFitter.fQuality), + fPedestalRMS (phosFitter.fPedestalRMS), + fAmpOffset (phosFitter.fAmpOffset), + fAmpThreshold(phosFitter.fAmpThreshold), + fOverflow (phosFitter.fOverflow), + fCalibData (phosFitter.fCalibData) +{ + //Copy constructor +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv0& AliPHOSRawFitterv0::operator = (const AliPHOSRawFitterv0 &phosFitter) +{ + //Assignment operator. + + if(this != &phosFitter) { + fSignal = phosFitter.fSignal; + fModule = phosFitter.fModule; + fCellX = phosFitter.fCellX; + fCellZ = phosFitter.fCellZ; + fCaloFlag = phosFitter.fCaloFlag; + fStart = phosFitter.fStart; + fLength = phosFitter.fLength; + fNBunches = phosFitter.fNBunches; + fPedSubtract = phosFitter.fPedSubtract; + fEnergy = phosFitter.fEnergy; + fTime = phosFitter.fTime; + fQuality = phosFitter.fQuality; + fPedestalRMS = phosFitter.fPedestalRMS; + fAmpOffset = phosFitter.fAmpOffset; + fAmpThreshold = phosFitter.fAmpThreshold; + fOverflow = phosFitter.fOverflow; + fCalibData = phosFitter.fCalibData; + } + + return *this; +} + +//----------------------------------------------------------------------------- + +void AliPHOSRawFitterv0::SetSamples(const UShort_t *sig, Int_t sigStart, Int_t sigLength) +{ + // Set the sample array, its start and length in time bin units + + fStart = sigStart; + fLength = sigLength; + fSignal = new UShort_t[fLength]; + for (Int_t i=0; i 1) { + fQuality = 1000; + return kTRUE; + } + + const Float_t kBaseLine = 1.0; + const Int_t kPreSamples = 10; + + Float_t pedMean = 0; + Float_t pedRMS = 0; + Int_t nPed = 0; + UShort_t maxSample = 0; + Int_t nMax = 0; + + for (Int_t i=0; i maxSample) maxSample = fSignal[i]; + if(fSignal[i] == maxSample) nMax++; + } + fEnergy = (Double_t)maxSample; + if (maxSample > 900 && nMax > 2) fOverflow = kTRUE; + + if (fPedSubtract) { + if (nPed > 0) { + fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ; + if(fPedestalRMS > 0.) + fPedestalRMS = TMath::Sqrt(fPedestalRMS) ; + fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction + } + else + return kFALSE; + } + else { + //take pedestals from DB + Double_t pedestal = (Double_t) fAmpOffset ; + if (fCalibData) { + Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ; + Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ; + pedestal += truePed - altroSettings ; + } + else{ + AliWarning(Form("Can not read data from OCDB")) ; + } + fEnergy-=pedestal ; + } + if (fEnergy < kBaseLine) fEnergy = 0; + + return kTRUE; + +} diff --git a/PHOS/AliPHOSRawFitterv0.h b/PHOS/AliPHOSRawFitterv0.h new file mode 100644 index 00000000000..631d1e73378 --- /dev/null +++ b/PHOS/AliPHOSRawFitterv0.h @@ -0,0 +1,69 @@ +#ifndef ALIPHOSRAWFITTERV0_H +#define ALIPHOSRAWFITTERV0_H +/* Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id: $ */ + +// This class extracts the signal parameters (energy, time, quality) +// from ALTRO samples. Energy is in ADC counts, time is in time bin units. +// A coarse algorithm is applied. + +class TArrayI; +class AliPHOSCalibData; + +class AliPHOSRawFitterv0 : public TObject +{ + +public: + + AliPHOSRawFitterv0(); + AliPHOSRawFitterv0(const AliPHOSRawFitterv0& rawFitterv0); + AliPHOSRawFitterv0& operator = (const AliPHOSRawFitterv0& rawFitterv0); + virtual ~AliPHOSRawFitterv0(); + + void SubtractPedestals(Bool_t subtract) {fPedSubtract = subtract;} + void SetAmpOffset (Int_t extPed=5) {fAmpOffset = extPed ;} + void SetAmpThreshold (Int_t thr=5) {fAmpThreshold = thr ;} + void SetSamples(const UShort_t *sig, Int_t sigStart, Int_t sigLength); + void SetNBunches(const Int_t nBunches) { fNBunches = nBunches; } + void SetChannelGeo(const Int_t module, const Int_t cellX, + const Int_t cellZ, const Int_t caloFlag); + + virtual Bool_t Eval(); + Double_t GetEnergy() const { return fEnergy; } + Double_t GetTime() const { return fTime; } + Double_t GetSignalQuality() const { return fQuality; } + Double_t GetPedestalRMS() const { return fPedestalRMS; } + Int_t GetModule() const { return fModule; } + Int_t GetCellX() const { return fCellX; } + Int_t GetCellZ() const { return fCellZ; } + Int_t GetCaloFlag() const { return fCaloFlag; } + Bool_t IsOverflow() const { return fOverflow; } + + void SetCalibData(AliPHOSCalibData * cdata){ fCalibData=cdata ;} + +protected: + + UShort_t *fSignal; // array of samples + Int_t fModule; // PHOS module number + Int_t fCellX; // cell number along X-axis + Int_t fCellZ; // cell number along Z-axis + Int_t fCaloFlag; // 0=LG, 1=HG, 2=TRU + Int_t fStart; // time bin of start signal + Int_t fLength; // signal length in time bins + Int_t fNBunches; // number of bunches in a signal + Bool_t fPedSubtract; // pedestals subtraction (kTRUE="yes") + Double_t fEnergy; // "digit" energy + Double_t fTime; // "digit" time + Double_t fQuality ; // sample quality + Double_t fPedestalRMS; // calciulated RMS of pedestal (non-ZS runs) + Int_t fAmpOffset ; // pedestal offset from ALTRO chips + Int_t fAmpThreshold ; // zero suppression threshold from ALTRO chips + Bool_t fOverflow ; // kTRUE is the signal overflows + AliPHOSCalibData * fCalibData ; //! Calibration database if avalable + + ClassDef(AliPHOSRawFitterv0,1) +}; + +#endif diff --git a/PHOS/AliPHOSRawFitterv1.cxx b/PHOS/AliPHOSRawFitterv1.cxx new file mode 100644 index 00000000000..7359f799284 --- /dev/null +++ b/PHOS/AliPHOSRawFitterv1.cxx @@ -0,0 +1,445 @@ +/************************************************************************** + * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id: $ */ + +// This class extracts the signal parameters (energy, time, quality) +// from ALTRO samples. Energy is in ADC counts, time is in time bin units. +// A fitting algorithm evaluates the energy and the time from Minuit minimization +// +// Typical use case: +// AliPHOSRawFitter *fitter=new AliPHOSRawFitter(); +// fitter->SetSamples(sig,sigStart,sigLength); +// fitter->SetNBunches(nBunches); +// fitter->SetChannelGeo(module,cellX,cellZ,caloFlag); +// fitter->SetCalibData(fgCalibData) ; +// fitter->Eval(); +// Double_t amplitude = fitter.GetEnergy(); +// Double_t time = fitter.GetTime(); +// Bool_t isLowGain = fitter.GetCaloFlag()==0; + +// Author: Dmitri Peressounko (Oct.2008) +// Modified: Yuri Kharlov (Jul.2009) + +// --- ROOT system --- +#include "TArrayD.h" +#include "TList.h" +#include "TMath.h" +#include "TMinuit.h" +#include "TCanvas.h" +#include "TH1.h" +#include "TH2.h" +#include "TF1.h" +#include "TROOT.h" + +// --- AliRoot header files --- +#include "AliLog.h" +#include "AliPHOSCalibData.h" +#include "AliPHOSRawFitterv1.h" +#include "AliPHOSPulseGenerator.h" + +ClassImp(AliPHOSRawFitterv1) + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv1::AliPHOSRawFitterv1(): + AliPHOSRawFitterv0(), + fSampleParamsLow(0x0), + fSampleParamsHigh(0x0), + fToFit(0x0) +{ + //Default constructor. + if(!gMinuit) + gMinuit = new TMinuit(100); + fSampleParamsHigh =new TArrayD(7) ; + fSampleParamsHigh->AddAt(2.174,0) ; + fSampleParamsHigh->AddAt(0.106,1) ; + fSampleParamsHigh->AddAt(0.173,2) ; + fSampleParamsHigh->AddAt(0.06106,3) ; + //last two parameters are pedestal and overflow + fSampleParamsLow=new TArrayD(7) ; + fSampleParamsLow->AddAt(2.456,0) ; + fSampleParamsLow->AddAt(0.137,1) ; + fSampleParamsLow->AddAt(2.276,2) ; + fSampleParamsLow->AddAt(0.08246,3) ; + fToFit = new TList() ; +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv1::~AliPHOSRawFitterv1() +{ + //Destructor. + if(fSampleParamsLow){ + delete fSampleParamsLow ; + fSampleParamsLow=0 ; + } + if(fSampleParamsHigh){ + delete fSampleParamsHigh ; + fSampleParamsHigh=0; + } + if(fToFit){ + delete fToFit ; + fToFit=0 ; + } +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv1::AliPHOSRawFitterv1(const AliPHOSRawFitterv1 &phosFitter ): + AliPHOSRawFitterv0(phosFitter), + fSampleParamsLow(0x0), + fSampleParamsHigh(0x0), + fToFit(0x0) +{ + //Copy constructor. + fToFit = new TList() ; + fSampleParamsLow =new TArrayD(*(phosFitter.fSampleParamsLow)) ; + fSampleParamsHigh=new TArrayD(*(phosFitter.fSampleParamsHigh)) ; +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv1& AliPHOSRawFitterv1::operator = (const AliPHOSRawFitterv1 &phosFitter) +{ + //Assignment operator. + + fToFit = new TList() ; + if(fSampleParamsLow){ + fSampleParamsLow = phosFitter.fSampleParamsLow ; + fSampleParamsHigh= phosFitter.fSampleParamsHigh ; + } + else{ + fSampleParamsLow =new TArrayD(*(phosFitter.fSampleParamsLow)) ; + fSampleParamsHigh=new TArrayD(*(phosFitter.fSampleParamsHigh)) ; + } + return *this; +} + +//----------------------------------------------------------------------------- +Bool_t AliPHOSRawFitterv1::Eval() +{ + //Extract an energy deposited in the crystal, + //crystal' position (module,column,row), + //time and gain (high or low). + //First collects sample, then evaluates it and if it has + //reasonable shape, fits it with Gamma2 function and extracts + //energy and time. + + if (fNBunches > 1) { + fQuality = 1000; + return kTRUE; + } + + Float_t pedMean = 0; + Float_t pedRMS = 0; + Int_t nPed = 0; + const Float_t kBaseLine = 1.0; + const Int_t kPreSamples = 10; + + TArrayI *fSamples = new TArrayI(fLength); // array of samples + TArrayI *fTimes = new TArrayI(fLength); // array of times corresponding to samples + for (Int_t i=0; iAddAt(fSignal[i],i); + fTimes ->AddAt( i ,i); + } + + Int_t iBin = fSamples->GetSize() ; + fEnergy = -111; + fQuality= 999. ; + const Float_t sampleMaxHG=102.332 ; //maximal height of HG sample with given parameterization + const Float_t sampleMaxLG=277.196 ; //maximal height of LG sample with given parameterization + const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it + Double_t pedestal = 0; + + if (fPedSubtract) { + if (nPed > 0) { + fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ; + if(fPedestalRMS > 0.) + fPedestalRMS = TMath::Sqrt(fPedestalRMS) ; + fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction + } + else + return kFALSE; + } + else { + //take pedestals from DB + pedestal = (Double_t) fAmpOffset ; + if (fCalibData) { + Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ; + Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ; + pedestal += truePed - altroSettings ; + } + else{ + AliWarning(Form("Can not read data from OCDB")) ; + } + fEnergy-=pedestal ; + } + + if (fEnergy < kBaseLine) fEnergy = 0; + + //calculate time and energy + Int_t maxBin=0 ; + Int_t maxAmp=0 ; + Double_t aMean=0. ; + Double_t aRMS=0. ; + Double_t wts=0 ; + Int_t tStart = 0 ; + for(Int_t i=0; iGetSize(); i++){ + if(fSamples->At(i) > pedestal){ + Double_t de=fSamples->At(i)-pedestal ; + if(de>1.){ + aMean+=de*i ; + aRMS+=de*i*i ; + wts+=de; + } + if(de>2 && tStart==0) + tStart=i ; + if(maxAmpAt(i)){ + maxBin=i ; + maxAmp=fSamples->At(i) ; + } + } + } + if(maxBin==fSamples->GetSize()-1){//bad "rising" sample + fEnergy=0. ; + fTime=-999.; + fQuality= 999. ; + return kTRUE ; + } + fEnergy=Double_t(maxAmp)-pedestal ; + fOverflow =0 ; //look for plato on the top of sample + if(fEnergy>500 && //this is not fluctuation of soft sample + maxBinGetSize()-1 && fSamples->At(maxBin+1)==maxAmp){ //and there is a plato + fOverflow = kTRUE ; + } + + if(wts>0){ + aMean/=wts; + aRMS=aRMS/wts-aMean*aMean; + } + + //do not take too small energies + if(fEnergy < kBaseLine) + fEnergy = 0; + + //do not test quality of too soft samples + if(fEnergyAt(tStart); + if(aRMS<2.) //sigle peak + fQuality=999. ; + else + fQuality= 0. ; + return kTRUE ; + } + + //IF sample has reasonable mean and RMS, try to fit it with gamma2 + + gMinuit->mncler(); // Reset Minuit's list of paramters + gMinuit->SetPrintLevel(-1) ; // No Printout + gMinuit->SetFCN(AliPHOSRawFitterv1::UnfoldingChiSquare) ; + // To set the address of the minimization function + + fToFit->Clear("nodelete") ; + Double_t b=0,bmin=0,bmax=0 ; + if(fCaloFlag == 0){ // Low gain + fSampleParamsLow->AddAt(pedestal,4) ; + if(fOverflow) + fSampleParamsLow->AddAt(double(maxAmp),5) ; + else + fSampleParamsLow->AddAt(double(1023),5) ; + fSampleParamsLow->AddAt(double(iBin),6) ; + fToFit->AddFirst((TObject*)fSampleParamsLow) ; + b=fSampleParamsLow->At(2) ; + bmin=0.5 ; + bmax=10. ; + } + else if(fCaloFlag == 1){ // High gain + fSampleParamsHigh->AddAt(pedestal,4) ; + if(fOverflow) + fSampleParamsHigh->AddAt(double(maxAmp),5) ; + else + fSampleParamsHigh->AddAt(double(1023),5); + fSampleParamsHigh->AddAt(double(iBin),6); + fToFit->AddFirst((TObject*)fSampleParamsHigh) ; + b=fSampleParamsHigh->At(2) ; + bmin=0.05 ; + bmax=0.4 ; + } + fToFit->AddLast((TObject*)fSamples) ; + fToFit->AddLast((TObject*)fTimes) ; + + gMinuit->SetObjectFit((TObject*)fToFit) ; // To tranfer pointer to UnfoldingChiSquare + Int_t ierflg ; + gMinuit->mnparm(0, "t0", 1.*tStart, 0.01, -500., 500., ierflg) ; + if(ierflg != 0){ + // AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ; + fEnergy=0. ; + fTime=-999. ; + fQuality=999 ; + return kTRUE ; //will scan further + } + Double_t amp0=0; + if(fCaloFlag == 0) // Low gain + amp0=fEnergy/sampleMaxLG; + else if(fCaloFlag == 1) // High gain + amp0=fEnergy/sampleMaxHG; + + gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ; + if(ierflg != 0){ + // AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ; + fEnergy=0. ; + fTime=-999. ; + fQuality=999 ; + return kTRUE ; //will scan further + } + + gMinuit->mnparm(2, "p2", b, 0.01*b, bmin, bmax, ierflg) ; + if(ierflg != 0){ + // AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ; + fEnergy=0. ; + fTime=-999. ; + fQuality=999 ; + return kTRUE ; //will scan further + } + + + Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly + // depends on it. + Double_t p1 = 1.0 ; + Double_t p2 = 0.0 ; + gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls + gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient + // gMinuit->SetMaxIterations(100); + gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings + + gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize + + Double_t err,t0err ; + Double_t t0,efit ; + gMinuit->GetParameter(0,t0, t0err) ; + gMinuit->GetParameter(1,efit, err) ; + + Double_t bfit, berr ; + gMinuit->GetParameter(2,bfit,berr) ; + + //Calculate total energy + //this isparameterization of depetendence of pulse height on parameter b + if(fCaloFlag == 0) // Low gain + efit*=99.54910 + 78.65038*bfit ; + else if(fCaloFlag == 1) // High gain + efit*=80.33109+128.6433*bfit ; + + if(efit<0. || efit > 10000.){ + //set energy to previously found max + fTime=-999.; + fQuality=999 ; + return kTRUE; + } + + //evaluate fit quality + Double_t fmin,fedm,errdef ; + Int_t npari,nparx,istat; + gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ; + fQuality=fmin/(fSamples->GetSize()-iBin) ; + //compare quality with some parameterization + if(fCaloFlag == 0) // Low gain + fQuality/=2.+0.002*fEnergy ; + else if(fCaloFlag == 1) // High gain + fQuality/=0.75+0.0025*fEnergy ; + + fEnergy=efit ; + fTime=t0-4.024*bfit ; //-10.402*bfit+4.669*bfit*bfit ; //Correction for 70 samples + + delete fSamples ; + delete fTimes ; + return kTRUE; +} +//_____________________________________________________________________________ +void AliPHOSRawFitterv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) +{ + // Number of parameters, Gradient, Chi squared, parameters, what to do + + TList * toFit= (TList*)gMinuit->GetObjectFit() ; + TArrayD * params=(TArrayD*)toFit->At(0) ; + TArrayI * samples = (TArrayI*)toFit->At(1) ; + TArrayI * times = (TArrayI*)toFit->At(2) ; + + fret = 0. ; + if(iflag == 2) + for(Int_t iparam = 0 ; iparam < 3 ; iparam++) + Grad[iparam] = 0 ; // Will evaluate gradient + + Double_t t0=x[0] ; + Double_t en=x[1] ; + Double_t b=x[2] ; + Double_t n=params->At(0) ; + Double_t alpha=params->At(1) ; + Double_t beta=params->At(3) ; + Double_t ped=params->At(4) ; + + Double_t overflow=params->At(5) ; + Int_t iBin = (Int_t) params->At(6) ; + Int_t nSamples=TMath::Min(iBin+70,samples->GetSize()) ; //Here we set number of points to fit (70) + // iBin - first non-zero sample + Int_t tStep=times->At(iBin+1)-times->At(iBin) ; + Double_t ddt=times->At(iBin)-t0-tStep ; + Double_t exp1=TMath::Exp(-alpha*ddt) ; + Double_t exp2=TMath::Exp(-beta*ddt) ; + Double_t dexp1=TMath::Exp(-alpha*tStep) ; + Double_t dexp2=TMath::Exp(-beta*tStep) ; + for(Int_t i = iBin; iAt(i))-t0 ; + Double_t fsample = double(samples->At(i)) ; + if(fsample>=overflow) + continue ; + Double_t diff ; + exp1*=dexp1 ; + exp2*=dexp2 ; + if(dt<=0.){ + diff=fsample - ped ; + fret += diff*diff ; + continue ; + } + Double_t dtn=TMath::Power(dt,n) ; + Double_t dtnE=dtn*exp1 ; + Double_t dt2E=dt*dt*exp2 ; + Double_t fit=ped+en*(dtnE + b*dt2E) ; + diff = fsample - fit ; + fret += diff*diff ; + if(iflag == 2){ // calculate gradient + Grad[0] += en*diff*(dtnE*(n/dt-alpha)+b*dt2E*(2./dt-beta)) ; //derivative over t0 + Grad[1] -= diff*(dtnE+b*dt2E) ; + Grad[2] -= en*diff*dt2E ; + } + } + if(iflag == 2) + for(Int_t iparam = 0 ; iparam < 3 ; iparam++) + Grad[iparam] *= 2. ; +} +//----------------------------------------------------------------------------- +Double_t AliPHOSRawFitterv1::Gamma2(Double_t dt,Double_t en,Double_t b,TArrayD * params){ //Function for fitting samples + //parameters: + //dt-time after start + //en-amplutude + //function parameters + + Double_t ped=params->At(4) ; + if(dt<0.) + return ped ; //pedestal + else + return ped+en*(TMath::Power(dt,params->At(0))*TMath::Exp(-dt*params->At(1))+b*dt*dt*TMath::Exp(-dt*params->At(3))) ; +} diff --git a/PHOS/AliPHOSRawFitterv1.h b/PHOS/AliPHOSRawFitterv1.h new file mode 100644 index 00000000000..f02e41d42f8 --- /dev/null +++ b/PHOS/AliPHOSRawFitterv1.h @@ -0,0 +1,42 @@ +#ifndef ALIPHOSRAWFITTERV1_H +#define ALIPHOSRAWFITTERV1_H +/* Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id: $ */ + +// This class extracts the PHOS "digits" of current event +// (amplitude,time, position,gain) from the raw stream +// provided by AliRawReader. See cxx source for use case. + +#include "AliPHOSRawFitterv0.h" +class TArrayD; +class TList; + +class AliPHOSRawFitterv1 : public AliPHOSRawFitterv0 { + +public: + + AliPHOSRawFitterv1(); + AliPHOSRawFitterv1(const AliPHOSRawFitterv1& rawFitter); + AliPHOSRawFitterv1& operator = (const AliPHOSRawFitterv1& rawFitter); + virtual ~AliPHOSRawFitterv1(); + + virtual Bool_t Eval(); + + static Double_t Gamma2(Double_t dt,Double_t en,Double_t b,TArrayD * fitparams) ; // Shape of correct sample + //class member function (not object member function) + static void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) ; + // Chi^2 of the fit. Should be static to be passed to MINUIT + void SetLowGainParams(Int_t n, Double_t * params){fSampleParamsLow->Set(n,params) ;} //fixed parameters of fit function + void SetHighGainParams(Int_t n,Double_t * params){fSampleParamsHigh->Set(n,params) ;} //fixed parameters of fit function + +private: + TArrayD *fSampleParamsLow ; //Fixed params of sample parameterization for Low gain + TArrayD *fSampleParamsHigh; //Fixed params of sample parameterization for High gain + TList *fToFit ; //! container to transfer parameters and data to fit + + ClassDef(AliPHOSRawFitterv1,1) +}; + +#endif diff --git a/PHOS/AliPHOSRawFitterv2.cxx b/PHOS/AliPHOSRawFitterv2.cxx new file mode 100644 index 00000000000..f38114c0783 --- /dev/null +++ b/PHOS/AliPHOSRawFitterv2.cxx @@ -0,0 +1,251 @@ +/************************************************************************** + * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id: $ */ + +// This class plots samples and qualify their quality according to their shape +// +// Typical use case: +// AliPHOSRawFitter *fitter=new AliPHOSRawFitter(); +// fitter->SetSamples(sig,sigStart,sigLength); +// fitter->SetNBunches(nBunches); +// fitter->SetChannelGeo(module,cellX,cellZ,caloFlag); +// fitter->SetCalibData(fgCalibData) ; +// fitter->Eval(); +// Double_t amplitude = fitter.GetEnergy(); +// Double_t time = fitter.GetTime(); +// Bool_t isLowGain = fitter.GetCaloFlag()==0; + +// Author: Dmitri Peressounko (Oct.2008) +// Modified: Yuri Kharlov (Jul.2009) + +// --- ROOT system --- +#include "TList.h" +#include "TMath.h" +#include "TMinuit.h" +#include "TCanvas.h" +#include "TH1.h" +#include "TH2.h" +#include "TF1.h" +#include "TROOT.h" + +// --- AliRoot header files --- +#include "AliLog.h" +#include "AliPHOSCalibData.h" +#include "AliPHOSRawFitterv2.h" +#include "AliPHOSPulseGenerator.h" + +ClassImp(AliPHOSRawFitterv2) + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv2::AliPHOSRawFitterv2(): + AliPHOSRawFitterv0(), + fNtimeSamples(25), + fRMScut(11.) +{ + //Default constructor. + fLGpar[0]=0.971 ; + fLGpar[1]=0.0465; + fLGpar[2]=1.56 ; + fHGpar[0]=0.941 ; + fHGpar[1]=0.0436; + fHGpar[2]=1.65 ; +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv2::~AliPHOSRawFitterv2() +{ + //Destructor. +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ): + AliPHOSRawFitterv0(phosFitter), + fNtimeSamples(25), + fRMScut(11.) +{ + //Copy constructor. + fNtimeSamples=phosFitter.fNtimeSamples ; + for(Int_t i=0; i<3;i++){ + fLGpar[i]=phosFitter.fLGpar[i] ; + fHGpar[i]=phosFitter.fHGpar[i] ; + } + fRMScut=phosFitter.fRMScut ; +} + +//----------------------------------------------------------------------------- +AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 &phosFitter) +{ + //Assignment operator. + + fNtimeSamples=phosFitter.fNtimeSamples ; + for(Int_t i=0; i<3;i++){ + fLGpar[i]=phosFitter.fLGpar[i] ; + fHGpar[i]=phosFitter.fHGpar[i] ; + } + fRMScut=phosFitter.fRMScut ; + return *this; +} + +//----------------------------------------------------------------------------- +Bool_t AliPHOSRawFitterv2::Eval() +{ + //Extract an energy deposited in the crystal, + //crystal' position (module,column,row), + //time and gain (high or low). + //First collects sample, then evaluates it and if it has + //reasonable shape, fits it with Gamma2 function and extracts + //energy and time. + + if (fNBunches > 1) { + fQuality = 1000; + return kTRUE; + } + + const Float_t kBaseLine = 1.0; + const Int_t kPreSamples = 10; + + Float_t pedMean = 0; + Float_t pedRMS = 0; + Int_t nPed = 0; + UShort_t maxSample = 0; + Int_t nMax = 0; + + TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ; + if(!cs) + cs = new TCanvas("CSample","CSample") ; + + TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ; + if(!h) h = new TH1D("hSample","",200,0.,200.) ; + + Double_t pedestal = 0; + for (Int_t i=0; i maxSample) maxSample = fSignal[i]; + if(fSignal[i] == maxSample) nMax++; + h->SetBinContent(i+1,fSignal[i]) ; + } + fEnergy = (Double_t)maxSample; + if (maxSample > 900 && nMax > 2) fOverflow = kTRUE; + + if (fPedSubtract) { + if (nPed > 0) { + fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ; + if(fPedestalRMS > 0.) + fPedestalRMS = TMath::Sqrt(fPedestalRMS) ; + fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction + } + else + return kFALSE; + } + else { + //take pedestals from DB + pedestal = (Double_t) fAmpOffset ; + if (fCalibData) { + Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ; + Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ; + pedestal += truePed - altroSettings ; + } + else{ + AliWarning(Form("Can not read data from OCDB")) ; + } + fEnergy-=pedestal ; + } + + if (fEnergy < kBaseLine) { + fEnergy = 0; + return kTRUE; + } + + // calculate time + fTime=0. ; + Double_t tRMS = 0. ; + Double_t tW = 0. ; + Int_t cnts=0 ; + Double_t a=0,b=0,c=0 ; + if(fCaloFlag == 0){ // Low gain + a=fLGpar[0] ; + b=fLGpar[1] ; + c=fLGpar[2] ; + } + else if(fCaloFlag == 1){ // High gain + a=fHGpar[0] ; + b=fHGpar[1] ; + c=fHGpar[2] ; + } + + + fQuality = 0. ; + + for(Int_t i=1; i0.){ + fTime/=tW ; + fQuality = tRMS/tW-fTime*fTime ; + } + else{ + fTime=-999. ; + fQuality=999. ; + } + + Bool_t isBad = 0 ; + for(Int_t i=1; i fSignal[i-1]+5 && fSignal[i] > fSignal[i+1]+5) { //single jump + isBad=1 ; + } + } + if(pedestal < 10.) + isBad=1 ; + + if(fPedestalRMS > 0.1) + isBad=1 ; + + for(Int_t i=1; i10. && !isBad ){ + printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ; + if(fOverflow)printf(" Overflow \n") ; + if(isBad)printf("bad") ; + cs->cd() ; + h->Draw() ; + cs->Update() ; + getchar() ; + } + + return kTRUE; +} diff --git a/PHOS/AliPHOSRawFitterv2.h b/PHOS/AliPHOSRawFitterv2.h new file mode 100644 index 00000000000..9d223622b84 --- /dev/null +++ b/PHOS/AliPHOSRawFitterv2.h @@ -0,0 +1,40 @@ +#ifndef ALIPHOSRAWFITTERV2_H +#define ALIPHOSRAWFITTERV2_H +/* Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id: $ */ + +// This class extracts the PHOS "digits" of current event +// (amplitude,time, position,gain) from the raw stream +// provided by AliRawReader. See cxx source for use case. + +#include "AliPHOSRawFitterv0.h" +class TList; + +class AliPHOSRawFitterv2 : public AliPHOSRawFitterv0 { + +public: + + AliPHOSRawFitterv2(); + AliPHOSRawFitterv2(const AliPHOSRawFitterv2& rawFitter); + AliPHOSRawFitterv2& operator = (const AliPHOSRawFitterv2& rawFitter); + virtual ~AliPHOSRawFitterv2(); + + virtual Bool_t Eval(); + + void SetNTimeSamples(Short_t n=25) { fNtimeSamples=n ;} + void SetLowGainTParams (Double_t *pars){ for(Int_t i=0;i<3;i++) fLGpar[i]=pars[i] ;} + void SetHighGainTParams(Double_t *pars){ for(Int_t i=0;i<3;i++) fHGpar[i]=pars[i] ;} + void SetRMScut(Double_t cut=2.) { fRMScut = cut ;} + +private: + Short_t fNtimeSamples ; //Number of samples (after start) used to extract time + Double_t fLGpar[3] ; //parameters for shape parameterization + Double_t fHGpar[3] ; //parameters for shape parameterization + Double_t fRMScut ; //cut to estmate goodness of sample + + ClassDef(AliPHOSRawFitterv2,1) +}; + +#endif diff --git a/PHOS/AliPHOSRecoParam.cxx b/PHOS/AliPHOSRecoParam.cxx index c1f47db1577..ce326e2719f 100644 --- a/PHOS/AliPHOSRecoParam.cxx +++ b/PHOS/AliPHOSRecoParam.cxx @@ -43,7 +43,7 @@ AliPHOSRecoParam::AliPHOSRecoParam() : fEMCSubtractPedestals(kTRUE), fEMCUnfold(kTRUE), fEMCEnergyCorrectionOn(kTRUE), - fEMCDecoderVersion(""), + fEMCFitterVersion(""), fGlobalAltroOffset(0), fGlobalAltroThreshold(0), fCPVClusteringThreshold(0.0), @@ -69,7 +69,7 @@ AliPHOSRecoParam::AliPHOSRecoParam(const AliPHOSRecoParam& ): fEMCSubtractPedestals(kTRUE), fEMCUnfold(kTRUE), fEMCEnergyCorrectionOn(kTRUE), - fEMCDecoderVersion(""), + fEMCFitterVersion(""), fGlobalAltroOffset(0), fGlobalAltroThreshold(0), fCPVClusteringThreshold(0.0), @@ -98,7 +98,7 @@ AliPHOSRecoParam& AliPHOSRecoParam::operator = (const AliPHOSRecoParam& recoPara fEMCSubtractPedestals = recoParam.fEMCSubtractPedestals; fEMCUnfold = recoParam.fEMCUnfold; fEMCEnergyCorrectionOn = recoParam.fEMCEnergyCorrectionOn; - fEMCDecoderVersion = recoParam.fEMCDecoderVersion; + fEMCFitterVersion = recoParam.fEMCFitterVersion; fGlobalAltroOffset = recoParam.fGlobalAltroOffset; fGlobalAltroThreshold = recoParam.fGlobalAltroThreshold; fCPVClusteringThreshold = recoParam.fCPVClusteringThreshold; @@ -126,7 +126,7 @@ void AliPHOSRecoParam::Print(Option_t * /*option*/) const "\tEMCSubtractPedestals = %d\n" "\tEMCUnfold = %d\n" "\tEMCEnergyCorrectionOn = %d\n" - "\tEMCDecoderVersion = \"%s\"\n" + "\tEMCFitterVersion = \"%s\"\n" "\tGlobalAltroOffset = %d", "\tGlobalAltroThreshold = %d", fEMCClusteringThreshold, @@ -140,7 +140,7 @@ void AliPHOSRecoParam::Print(Option_t * /*option*/) const fEMCSubtractPedestals, fEMCUnfold, fEMCEnergyCorrectionOn, - fEMCDecoderVersion.Data(), + fEMCFitterVersion.Data(), fGlobalAltroOffset, fGlobalAltroThreshold)); diff --git a/PHOS/AliPHOSRecoParam.h b/PHOS/AliPHOSRecoParam.h index fb939582c0b..bd2037073c4 100644 --- a/PHOS/AliPHOSRecoParam.h +++ b/PHOS/AliPHOSRecoParam.h @@ -29,7 +29,7 @@ public: Bool_t EMCEcore2ESD() const { return fEMCEcore2ESD; } Bool_t EMCSubtractPedestals() const { return fEMCSubtractPedestals; } Bool_t EMCToUnfold() const { return fEMCUnfold; } - const char* EMCDecoderVersion() const { return fEMCDecoderVersion.Data();} + const char* EMCFitterVersion() const { return fEMCFitterVersion.Data(); } Bool_t GetEMCEnergyCorrectionOn() const { return fEMCEnergyCorrectionOn; } Int_t GetGlobalAltroOffset() const { return fGlobalAltroOffset ; } Int_t GetGlobalAltroThreshold() const { return fGlobalAltroThreshold ; } @@ -49,7 +49,7 @@ public: void SetEMCEcoreRadius(Float_t rCore) { fEMCEcoreRadius =rCore; } void SetEMCEcore2ESD(Bool_t ecore) { fEMCEcore2ESD =ecore; } void SetEMCSubtractPedestals(Bool_t subtract) { fEMCSubtractPedestals =subtract;} - void SetEMCDecoderVersion(const char* version="v1"){ fEMCDecoderVersion =version ;} + void SetEMCFitterVersion(const char* version="v1") { fEMCFitterVersion =version ; } void SetEMCUnfolding(Bool_t toUnfold=kFALSE) { fEMCUnfold =toUnfold;} void SetEMCEnergyCorrectionOn(Bool_t on=kTRUE) { fEMCEnergyCorrectionOn =on; } void SetGlobalAltroOffset(Int_t offset=5) { fGlobalAltroOffset =offset ; } @@ -79,7 +79,7 @@ protected: Bool_t fEMCSubtractPedestals; // EMC: true if pedestal should be subtracted (in non-ZS) Bool_t fEMCUnfold; // EMC: true if overlapped clusters should be unfolded Bool_t fEMCEnergyCorrectionOn; // EMC: if true do non-linear correction of cluster energy - TString fEMCDecoderVersion ; // EMC: AliPHOSRawDecoder version + TString fEMCFitterVersion ; // EMC: AliPHOSRawFitter version Int_t fGlobalAltroOffset ; // Offset used in ALTRO chips in SZ runs Int_t fGlobalAltroThreshold ; // Threshold used in ALTRO chips in SZ runs diff --git a/PHOS/AliPHOSReconstructor.cxx b/PHOS/AliPHOSReconstructor.cxx index c3a1f6de4af..fbaab3edae7 100644 --- a/PHOS/AliPHOSReconstructor.cxx +++ b/PHOS/AliPHOSReconstructor.cxx @@ -46,9 +46,9 @@ #include "AliPHOSTrackSegment.h" #include "AliPHOSEmcRecPoint.h" #include "AliPHOSRecParticle.h" -#include "AliPHOSRawDecoder.h" -#include "AliPHOSRawDecoderv1.h" -#include "AliPHOSRawDecoderv2.h" +#include "AliPHOSRawFitterv0.h" +#include "AliPHOSRawFitterv1.h" +#include "AliPHOSRawFitterv2.h" #include "AliPHOSRawDigiProducer.h" #include "AliPHOSPulseGenerator.h" @@ -343,7 +343,7 @@ void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits // Works on a single-event basis rawReader->Reset() ; - AliPHOSRawDecoder * dc ; + AliPHOSRawFitterv0 * fitter ; const TObjArray* maps = AliPHOSRecoParam::GetMappings(); if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); @@ -353,55 +353,54 @@ void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits mapping[i] = (AliAltroMapping*)maps->At(i); } - if (strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0) - dc=new AliPHOSRawDecoderv1(rawReader,mapping); - else if (strcmp(GetRecoParam()->EMCDecoderVersion(),"v2")==0) - dc=new AliPHOSRawDecoderv2(rawReader,mapping); + if (strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0) + fitter=new AliPHOSRawFitterv1(); + else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v2")==0) + fitter=new AliPHOSRawFitterv2(); else - dc=new AliPHOSRawDecoder(rawReader,mapping); + fitter=new AliPHOSRawFitterv0(); - dc->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals()); - dc->SetAmpOffset(GetRecoParam()->GetGlobalAltroOffset()); - dc->SetAmpThreshold(GetRecoParam()->GetGlobalAltroThreshold()); + fitter->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals()); + fitter->SetAmpOffset (GetRecoParam()->GetGlobalAltroOffset()); + fitter->SetAmpThreshold (GetRecoParam()->GetGlobalAltroThreshold()); TClonesArray *digits = new TClonesArray("AliPHOSDigit",1); digits->SetName("DIGITS"); Int_t bufsize = 32000; digitsTree->Branch("PHOS", &digits, bufsize); - //AliPHOSRawDigiProducer pr(GetRecoParam()); - AliPHOSRawDigiProducer pr; + AliPHOSRawDigiProducer rdp(rawReader,mapping); - pr.SetEmcMinAmp(GetRecoParam()->GetEMCRawDigitThreshold()); // in ADC - pr.SetCpvMinAmp(GetRecoParam()->GetCPVMinE()); - pr.SetSampleQualityCut(GetRecoParam()->GetEMCSampleQualityCut()); - pr.MakeDigits(digits,dc); + rdp.SetEmcMinAmp(GetRecoParam()->GetEMCRawDigitThreshold()); // in ADC + rdp.SetCpvMinAmp(GetRecoParam()->GetCPVMinE()); + rdp.SetSampleQualityCut(GetRecoParam()->GetEMCSampleQualityCut()); + rdp.MakeDigits(digits,fitter); - delete dc ; + delete fitter ; - //!!!!for debug!!! -/* - Int_t modMax=-111; - Int_t colMax=-111; - Int_t rowMax=-111; - Float_t eMax=-333; - //!!!for debug!!! - - Int_t relId[4]; - for(Int_t iDigit=0; iDigitGetEntries(); iDigit++) { - AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit); - if(digit->GetEnergy()>eMax) { - fGeom->AbsToRelNumbering(digit->GetId(),relId); - eMax=digit->GetEnergy(); - modMax=relId[0]; - rowMax=relId[2]; - colMax=relId[3]; + if (AliLog::GetGlobalDebugLevel() == 1) { + Int_t modMax=-111; + Int_t colMax=-111; + Int_t rowMax=-111; + Float_t eMax=-333; + //!!!for debug!!! + + Int_t relId[4]; + for(Int_t iDigit=0; iDigitGetEntries(); iDigit++) { + AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit); + if(digit->GetEnergy()>eMax) { + fGeom->AbsToRelNumbering(digit->GetId(),relId); + eMax=digit->GetEnergy(); + modMax=relId[0]; + rowMax=relId[2]; + colMax=relId[3]; + } } + + AliDebug(1,Form("Digit with max. energy: modMax %d colMax %d rowMax %d eMax %f\n\n", + modMax,colMax,rowMax,eMax)); } - AliDebug(1,Form("Digit with max. energy: modMax %d colMax %d rowMax %d eMax %f\n\n", - modMax,colMax,rowMax,eMax)); -*/ digitsTree->Fill(); digits->Delete(); delete digits; diff --git a/PHOS/CMake_libPHOSbase.txt b/PHOS/CMake_libPHOSbase.txt index 6aae51f6869..8b7f60b535c 100644 --- a/PHOS/CMake_libPHOSbase.txt +++ b/PHOS/CMake_libPHOSbase.txt @@ -12,9 +12,9 @@ set(SRCS AliPHOSCpvCalibData.cxx AliPHOSEmcCalibData.cxx AliPHOSPulseGenerator.cxx - AliPHOSRawDecoder.cxx - AliPHOSRawDecoderv1.cxx - AliPHOSRawDecoderv2.cxx + AliPHOSRawFitter.cxx + AliPHOSRawFitterv1.cxx + AliPHOSRawFitterv2.cxx AliPHOSRawDigiProducer.cxx AliPHOSEmcBadChannelsMap.cxx AliPHOSSurvey.cxx diff --git a/PHOS/PHOSbaseLinkDef.h b/PHOS/PHOSbaseLinkDef.h index 28f68889e6e..b930df031c9 100644 --- a/PHOS/PHOSbaseLinkDef.h +++ b/PHOS/PHOSbaseLinkDef.h @@ -18,9 +18,9 @@ #pragma link C++ class AliPHOSCpvCalibData+; #pragma link C++ class AliPHOSEmcCalibData+; #pragma link C++ class AliPHOSPulseGenerator+; -#pragma link C++ class AliPHOSRawDecoder+; -#pragma link C++ class AliPHOSRawDecoderv1+; -#pragma link C++ class AliPHOSRawDecoderv2+; +#pragma link C++ class AliPHOSRawFitterv0+; +#pragma link C++ class AliPHOSRawFitterv1+; +#pragma link C++ class AliPHOSRawFitterv2+; #pragma link C++ class AliPHOSRawDigiProducer+; #pragma link C++ class AliPHOSEmcBadChannelsMap+; #pragma link C++ class AliPHOSSurvey+; diff --git a/PHOS/libPHOSbase.pkg b/PHOS/libPHOSbase.pkg index af5a4365591..0c463c41035 100644 --- a/PHOS/libPHOSbase.pkg +++ b/PHOS/libPHOSbase.pkg @@ -12,9 +12,9 @@ SRCS = \ AliPHOSCpvCalibData.cxx \ AliPHOSEmcCalibData.cxx \ AliPHOSPulseGenerator.cxx \ - AliPHOSRawDecoder.cxx \ - AliPHOSRawDecoderv1.cxx \ - AliPHOSRawDecoderv2.cxx \ + AliPHOSRawFitterv0.cxx \ + AliPHOSRawFitterv1.cxx \ + AliPHOSRawFitterv2.cxx \ AliPHOSRawDigiProducer.cxx \ AliPHOSEmcBadChannelsMap.cxx \ AliPHOSSurvey.cxx \ -- 2.43.0