2) Set in OCDB the type of the fitting algorithm and if the bad channels are fitted or not
3) 2 fitters available, kStandard (used until now) and kFastFit from Aleksei, the call to this one needs checks, not working properly yet.
if(fGeoName.Contains("WSUC")) fGeoName = "EMCAL_WSUC";
//check that we have a valid geometry name
- if(!(fGeoName.Contains("EMCAL_PDC06") || fGeoName.Contains("EMCAL_COMPLETE") || fGeoName.Contains("EMCAL_WSUC") || fGeoName.Contains("EMCAL_1stYear"))) {
+ if(!(fGeoName.Contains("EMCAL_PDC06") || fGeoName.Contains("EMCAL_COMPLETE") || fGeoName.Contains("EMCAL_WSUC") || fGeoName.Contains("EMCAL_FIRSTYEAR"))) {
Fatal("Init", "%s is an undefined geometry!", fGeoName.Data()) ;
}
CheckAdditionalOptions();
}
- if(fGeoName.Contains("1stYear")){
- fNumberOfSuperModules = 2;
-
- if(fGeoName.Contains("LowerEta")) {
- fNPhiSuperModule = 1;
- }
- else if(fGeoName.Contains("LowerPhi_SideA")){
- fNPhiSuperModule = 2;
- fArm1EtaMax=0;
- }
- else if(fGeoName.Contains("LowerPhi_SideC")){
- fNPhiSuperModule = 2;
- fArm1EtaMin=0;
- }
-
- CheckAdditionalOptions();
- }
-
+ //In 2009-2010 data taking runs only 4 SM, in the upper position.
+ if(fGeoName.Contains("FIRSTYEAR")){
+ fNumberOfSuperModules = 4;
+ fArm1PhiMax = 120.0;
+ CheckAdditionalOptions();
+ }
+
// constant for transition absid <--> indexes
- fNCellsInModule = fNPHIdiv*fNETAdiv;
+ fNCellsInModule = fNPHIdiv*fNETAdiv;
fNCellsInSupMod = fNCellsInModule*fNPhi*fNZ;
fNCells = fNCellsInSupMod*fNumberOfSuperModules;
if(GetKey110DEG()) fNCells -= fNCellsInSupMod;
if(fNumberOfSuperModules > 1)
fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance);
if(fNumberOfSuperModules > 2) {
- for(int i=1; i<=4; i++) { // from 2th ro 9th
+ Int_t maxPhiBlock =fNumberOfSuperModules/2-1;
+ if(fNumberOfSuperModules > 10) maxPhiBlock = 4;
+ for(int i=1; i<=maxPhiBlock; i++) { // from 2th ro 9th
fPhiBoundariesOfSM[2*i] = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i;
fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i;
- fPhiCentersOfSM[i] = fPhiCentersOfSM[0] + 20.*TMath::DegToRad()*i;
+ fPhiCentersOfSM[i] = fPhiCentersOfSM[0] + 20.*TMath::DegToRad()*i;
}
}
if(fNumberOfSuperModules > 10) {
fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad();
fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance);
- fPhiCentersOfSM[5] = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.;
+ fPhiCentersOfSM[5] = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.;
}
//called after setting of scintillator and lead layer parameters
// Jet trigger
// 3*6*10 + 2*6*2 = 204 -> matrix (nphi(17), neta(12))
fNEtaSubOfTRU = 6;
-
fgInit = kTRUE;
}
#include "AliEMCALDigit.h"
#include "AliEMCAL.h"
#include "AliCaloCalibPedestal.h"
-
+#include "AliCaloFastAltroFitv0.h"
+
ClassImp(AliEMCALRawUtils)
// Signal shape parameters
AliEMCALRawUtils::AliEMCALRawUtils()
: fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
- fNPedSamples(0), fGeom(0), fOption("")
+ fNPedSamples(0), fGeom(0), fOption(""),
+ fRemoveBadChannels(kTRUE),fFittingAlgorithm(0)
{
//These are default parameters.
//Can be re-set from without with setter functions
+ //Already set in the OCDB and passed via setter in the AliEMCALReconstructor
fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits)
fOrder = 2; // order of gamma fn
fTau = 2.35; // in units of timebin, from CERN 2007 testbeam
fNoiseThreshold = 3; // 3 ADC counts is approx. noise level
fNPedSamples = 4; // less than this value => likely pedestal samples
-
+ fRemoveBadChannels = kTRUE; //Remove bad channels before fitting
+ fFittingAlgorithm = kFastFit;//kStandard; // Use default minuit fitter
+
//Get Mapping RCU files from the AliEMCALRecParam
const TObjArray* maps = AliEMCALRecParam::GetMappings();
if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
//____________________________________________________________________________
AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry)
: fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0),
- fNPedSamples(0), fGeom(pGeometry), fOption("")
+ fNPedSamples(0), fGeom(pGeometry), fOption(""),
+ fRemoveBadChannels(kTRUE),fFittingAlgorithm(0)
{
//
// Initialize with the given geometry - constructor required by HLT
//These are default parameters.
//Can be re-set from without with setter functions
+ //Already set in the OCDB and passed via setter in the AliEMCALReconstructor
fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits)
fOrder = 2; // order of gamma fn
fTau = 2.35; // in units of timebin, from CERN 2007 testbeam
fNoiseThreshold = 3; // 3 ADC counts is approx. noise level
fNPedSamples = 4; // less than this value => likely pedestal samples
-
+ fRemoveBadChannels = kTRUE; //Remove bad channels before fitting
+ fFittingAlgorithm = kStandard; // Use default minuit fitter
+
//Get Mapping RCU files from the AliEMCALRecParam
const TObjArray* maps = AliEMCALRecParam::GetMappings();
if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
fNoiseThreshold(rawU.fNoiseThreshold),
fNPedSamples(rawU.fNPedSamples),
fGeom(rawU.fGeom),
- fOption(rawU.fOption)
+ fOption(rawU.fOption),
+ fRemoveBadChannels(rawU.fRemoveBadChannels),
+ fFittingAlgorithm(rawU.fFittingAlgorithm)
{
//copy ctor
fMapping[0] = rawU.fMapping[0];
fNPedSamples = rawU.fNPedSamples;
fGeom = rawU.fGeom;
fOption = rawU.fOption;
+ fRemoveBadChannels = rawU.fRemoveBadChannels;
+ fFittingAlgorithm = rawU.fFittingAlgorithm;
fMapping[0] = rawU.fMapping[0];
fMapping[1] = rawU.fMapping[1];
fMapping[2] = rawU.fMapping[2];
if (caloFlag != 0 && caloFlag != 1) continue;
//Do not fit bad channels
- if(pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) {
+ if(fRemoveBadChannels && pedbadmap->IsBadChannel(in.GetModule(),in.GetColumn(),in.GetRow())) {
//printf("Tower from SM %d, column %d, row %d is BAD!!! Skip \n", in.GetModule(),in.GetColumn(),in.GetRow());
continue;
}
ampEstimate = -1 ;
timeEstimate = -1 ;
pedEstimate = -1;
- if ( (max - min) > fNoiseThreshold) {
- FitRaw(gSig, signalF, maxTimeBin, amp, time, ped,
+
+ if ( (max - min) > fNoiseThreshold) {
+ switch(fFittingAlgorithm)
+ {
+ case kStandard:
+ {
+ //printf("Standard fitter \n");
+ FitRaw(gSig, signalF, maxTimeBin, amp, time, ped,
ampEstimate, timeEstimate, pedEstimate);
- }
+ break;
+ }
+ case kFastFit:
+ {
+ //printf("FastFitter \n");
+ Double_t eSignal = 0;
+ Double_t dAmp = amp;
+ Double_t dTimeEstimate = timeEstimate;
+ Double_t eTimeEstimate = 0;
+ Double_t eAmp = 0;
+ Double_t chi2 = 0;
+
+ AliCaloFastAltroFitv0::FastFit(gSig->GetX(), gSig->GetY(), gSig->GetN(),
+ eSignal, fTau,
+ dAmp, eAmp, dTimeEstimate, eTimeEstimate, chi2);
+ amp=dAmp;
+ timeEstimate = dTimeEstimate;
+ //printf("FastFitter: Amp %f, time %f, eAmp %f, eTimeEstimate %f, chi2 %f\n",amp, timeEstimate,eAmp,eTimeEstimate,chi2);
+
+ break;
+ }
+ }
+ }
if ( amp>0 && amp<2000 && time>0 && time<(maxTimeBin*GetRawFormatTimeBinWidth()) ) { //check both high and low end of amplitude result, and time
//2000 is somewhat arbitrary - not nice with magic numbers in the code..
AliEMCALRawUtils();
AliEMCALRawUtils(AliEMCALGeometry *pGeometry);
virtual ~AliEMCALRawUtils();
-
+
+ enum fitAlgorithm {kStandard = 0, kFastFit= 1};
+
AliEMCALRawUtils(const AliEMCALRawUtils& rawUtils); //copy ctor
AliEMCALRawUtils& operator =(const AliEMCALRawUtils& rawUtils);
Double_t GetPedestalValue() const {return fgPedestalValue;}
Double_t GetFEENoise() const {return fgFEENoise;}
+ Bool_t GetRemoveBadChannels() const {return fRemoveBadChannels;}
+ Int_t GetFittingAlgorithm() const {return fFittingAlgorithm; }
+
void SetRawFormatHighLowGainFactor(Double_t val) {fHighLowGainFactor=val;}
void SetRawFormatOrder(Int_t val) {fOrder=val; }
void SetRawFormatTau(Double_t val) {fTau=val; }
void SetNoiseThreshold(Int_t val) {fNoiseThreshold=val; }
void SetNPedSamples(Int_t val) {fNPedSamples=val; }
+ void SetRemoveBadChannels(Bool_t val) {fRemoveBadChannels=val; }
+ void SetFittingAlgorithm(Int_t val) {fFittingAlgorithm=val; }
// set methods for fast fit simulation
void SetFEENoise(Double_t val) {fgFEENoise = val;}
Double_t GetRawFormatTimeTrigger() const { return fgTimeTrigger ; }
Int_t GetRawFormatThreshold() const { return fgThreshold ; }
Int_t GetRawFormatDDLPerSuperModule() const { return fgDDLPerSuperModule ; }
-
+
virtual Option_t* GetOption() const { return fOption.Data(); }
void SetOption(Option_t* opt) { fOption = opt; }
// Signal shape functions
+
void FitRaw(TGraph * gSig, TF1* signalF, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & ped, Float_t & ampEstimate, Float_t & timeEstimate, Float_t & pedEstimate, const Float_t cut = 0) const ;
static Double_t RawResponseFunction(Double_t *x, Double_t *par);
Bool_t RawSampledResponse(Double_t dtime, Double_t damp, Int_t * adcH, Int_t * adcL) const;
static Int_t fgPedestalValue; // pedestal value for Digits2Raw
static Double_t fgFEENoise; // electronics noise in ADC units
- AliEMCALGeometry* fGeom; //geometry
- AliAltroMapping* fMapping[4]; //only two for now
+ AliEMCALGeometry* fGeom; //geometry
+ AliAltroMapping* fMapping[4]; //only two for now
TString fOption; //! option passed from Reconstructor
- ClassDef(AliEMCALRawUtils,3) // utilities for raw signal fitting
+ Bool_t fRemoveBadChannels; // select if bad channels are removed before fitting
+ Int_t fFittingAlgorithm; // select the fitting algorithm
+
+ ClassDef(AliEMCALRawUtils,4) // utilities for raw signal fitting
};
#endif
fOrderParameter(2),
fTau(2.35),
fNoiseThreshold(3),
- fNPedSamples(5) //raw signal
+ fNPedSamples(5),
+ fRemoveBadChannels(kTRUE),
+ fFittingAlgorithm(0)//raw signal
{
// default reco values
fGamma[i][j] = fPiZero[i][j] = fHadron[i][j] = 0.;
fGamma1to10[i][j] = fHadron1to10[i][j]= 0.;
}
- fGammaEnergyProb[i]=0.; // not yet implemented
+ fGammaEnergyProb[i] =0.; // not yet implemented
fHadronEnergyProb[i]=0.;
fPiZeroEnergyProb[i]=0.; // not yet implemented
fOrderParameter(rp.fOrderParameter),
fTau(rp.fTau),
fNoiseThreshold(rp.fNoiseThreshold),
- fNPedSamples(rp.fNPedSamples) //raw signal
+ fNPedSamples(rp.fNPedSamples),
+ fRemoveBadChannels(rp.fRemoveBadChannels),
+ fFittingAlgorithm(rp.fFittingAlgorithm) //raw signal
{
//copy constructor
Int_t i, j;
for (i = 0; i < 6; i++) {
for (j = 0; j < 6; j++) {
- fGamma[i][j] = rp.fGamma[i][j];
- fGamma1to10[i][j] = rp.fGamma1to10[i][j];
- fHadron[i][j] = rp.fHadron[i][j];
+ fGamma[i][j] = rp.fGamma[i][j];
+ fGamma1to10[i][j] = rp.fGamma1to10[i][j];
+ fHadron[i][j] = rp.fHadron[i][j];
fHadron1to10[i][j] = rp.fHadron1to10[i][j];
- fPiZero[i][j] = rp.fPiZero[i][j];
+ fPiZero[i][j] = rp.fPiZero[i][j];
}
- fGammaEnergyProb[i] = rp.fGammaEnergyProb[i];
+ fGammaEnergyProb[i] = rp.fGammaEnergyProb[i];
fPiZeroEnergyProb[i] = rp.fPiZeroEnergyProb[i];
fHadronEnergyProb[i] = rp.fHadronEnergyProb[i];
if(this != &rp) {
fClusteringThreshold = rp.fClusteringThreshold;
- fW0 = rp.fW0;
- fMinECut = rp.fMinECut;
- fUnfold = rp.fUnfold;
+ fW0 = rp.fW0;
+ fMinECut = rp.fMinECut;
+ fUnfold = rp.fUnfold;
fLocMaxCut = rp.fLocMaxCut;
- fTimeCut = rp.fTimeCut;//clustering
- fTrkCutX = rp.fTrkCutX;
- fTrkCutY = rp.fTrkCutY;
- fTrkCutZ = rp.fTrkCutZ;
- fTrkCutR = rp.fTrkCutR;
+ fTimeCut = rp.fTimeCut;//clustering
+ fTrkCutX = rp.fTrkCutX;
+ fTrkCutY = rp.fTrkCutY;
+ fTrkCutZ = rp.fTrkCutZ;
+ fTrkCutR = rp.fTrkCutR;
fTrkCutAlphaMin = rp.fTrkCutAlphaMin;
fTrkCutAlphaMax = rp.fTrkCutAlphaMax;
- fTrkCutAngle = rp.fTrkCutAngle;
- fTrkCutNITS = rp.fTrkCutNITS;
- fTrkCutNTPC = rp.fTrkCutNTPC; //track matching
+ fTrkCutAngle = rp.fTrkCutAngle;
+ fTrkCutNITS = rp.fTrkCutNITS;
+ fTrkCutNTPC = rp.fTrkCutNTPC; //track matching
fHighLowGainFactor = rp.fHighLowGainFactor;
- fOrderParameter = rp.fOrderParameter;
- fTau = rp.fTau;
- fNoiseThreshold = rp.fNoiseThreshold;
- fNPedSamples = rp.fNPedSamples; //raw signal
-
+ fOrderParameter = rp.fOrderParameter;
+ fTau = rp.fTau;
+ fNoiseThreshold = rp.fNoiseThreshold;
+ fNPedSamples = rp.fNPedSamples;
+ fRemoveBadChannels = rp.fRemoveBadChannels;
+ fFittingAlgorithm = rp.fFittingAlgorithm;//raw signal
+
//PID values
Int_t i, j;
for (i = 0; i < 6; i++) {
for (j = 0; j < 6; j++) {
- fGamma[i][j] = rp.fGamma[i][j];
- fGamma1to10[i][j] = rp.fGamma1to10[i][j];
- fHadron[i][j] = rp.fHadron[i][j];
+ fGamma[i][j] = rp.fGamma[i][j];
+ fGamma1to10[i][j] = rp.fGamma1to10[i][j];
+ fHadron[i][j] = rp.fHadron[i][j];
fHadron1to10[i][j] = rp.fHadron1to10[i][j];
- fPiZero[i][j] = rp.fPiZero[i][j];
+ fPiZero[i][j] = rp.fPiZero[i][j];
}
- fGammaEnergyProb[i] = rp.fGammaEnergyProb[i];
+ fGammaEnergyProb[i] = rp.fGammaEnergyProb[i];
fPiZeroEnergyProb[i] = rp.fPiZeroEnergyProb[i];
fHadronEnergyProb[i] = rp.fHadronEnergyProb[i];
}
AliInfo(Form("Raw signal parameters: \n gain factor=%f, order=%d, tau=%f, noise threshold=%d, nped samples=%d \n",
fHighLowGainFactor,fOrderParameter,fTau,fNoiseThreshold,fNPedSamples));
+ AliInfo(Form("Raw signal: with bad channels? %d, \n \t with fitting algorithm %d \n",
+ fRemoveBadChannels, fFittingAlgorithm));
}
void SetOrderParameter(Int_t value) {fOrderParameter = value;}
void SetTau(Double_t value) {fTau = value;}
void SetNoiseThreshold(Int_t value) {fNoiseThreshold = value;}
- void SetNPedSamples(Int_t value) {fNPedSamples = value;}
+ void SetNPedSamples(Int_t value) {fNPedSamples = value;}
+ void SetRemoveBadChannels(Bool_t val) {fRemoveBadChannels=val; }
+ void SetFittingAlgorithm(Int_t val) {fFittingAlgorithm=val; }
+
/* raw signal getters */
Double_t GetHighLowGainFactor() const {return fHighLowGainFactor;}
Int_t GetOrderParameter() const {return fOrderParameter;}
Double_t GetTau() const {return fTau;}
Int_t GetNoiseThreshold() const {return fNoiseThreshold;}
Int_t GetNPedSamples() const {return fNPedSamples;}
-
- virtual void Print(Option_t * option="") const ;
+ Bool_t GetRemoveBadChannels() const {return fRemoveBadChannels;}
+ Int_t GetFittingAlgorithm() const {return fFittingAlgorithm; }
+
+
+ virtual void Print(Option_t * option="") const ;
static AliEMCALRecParam* GetDefaultParameters();
static AliEMCALRecParam* GetLowFluxParam();
Double_t fTrkCutNTPC; // Number of TPC hits for track matching
//Raw signal fitting parameters (Jenn)
- Double_t fHighLowGainFactor; //gain factor to convert between high and low gain
- Int_t fOrderParameter; //order parameter for raw signal fit
- Double_t fTau; //decay constant for raw signal fit
- Int_t fNoiseThreshold; //threshold to consider signal or noise
- Int_t fNPedSamples; //number of time samples to use in pedestal calculation
-
+ Double_t fHighLowGainFactor; // gain factor to convert between high and low gain
+ Int_t fOrderParameter; // order parameter for raw signal fit
+ Double_t fTau; // decay constant for raw signal fit
+ Int_t fNoiseThreshold; // threshold to consider signal or noise
+ Int_t fNPedSamples; // number of time samples to use in pedestal calculation
+ Bool_t fRemoveBadChannels; // select if bad channels are removed before fitting
+ Int_t fFittingAlgorithm; // select the fitting algorithm
+
static TObjArray* fgkMaps; // ALTRO mappings for RCU0..RCUX
- ClassDef(AliEMCALRecParam,8) // Reconstruction parameters
+ ClassDef(AliEMCALRecParam,9) // Reconstruction parameters
} ;
fgRawUtils->SetRawFormatTau(GetRecParam()->GetTau());
fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
fgRawUtils->SetNPedSamples(GetRecParam()->GetNPedSamples());
+ fgRawUtils->SetRemoveBadChannels(GetRecParam()->GetRemoveBadChannels());
+ fgRawUtils->SetFittingAlgorithm(GetRecParam()->GetFittingAlgorithm());
fgRawUtils->Raw2Digits(rawReader,digitsArr,fPedestalData);