#include <TMath.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TStyle.h>
#include "AliTRDCalDet.h"
#include "AliTRDgeometry.h"
#include "AliMathBase.h"
#include "AliTRDpadPlane.h"
ClassImp(AliTRDCalDet)
AliTRDCalDet::AliTRDCalDet():TNamed()
{
for (Int_t idet = 0; idet < kNdet; idet++) {
fData[idet] = 0;
}
}
AliTRDCalDet::AliTRDCalDet(const Text_t *name, const Text_t *title)
:TNamed(name,title)
{
for (Int_t idet = 0; idet < kNdet; idet++) {
fData[idet] = 0;
}
}
AliTRDCalDet::AliTRDCalDet(const AliTRDCalDet &c):TNamed(c)
{
((AliTRDCalDet &) c).Copy(*this);
}
AliTRDCalDet::~AliTRDCalDet()
{
}
AliTRDCalDet &AliTRDCalDet::operator=(const AliTRDCalDet &c)
{
if (this != &c) ((AliTRDCalDet &) c).Copy(*this);
return *this;
}
void AliTRDCalDet::Copy(TObject &c) const
{
for (Int_t idet = 0; idet < kNdet; idet++) {
((AliTRDCalDet &) c).fData[idet] = fData[idet];
}
TObject::Copy(c);
}
Double_t AliTRDCalDet::GetMean(AliTRDCalDet * const outlierDet) const
{
if (!outlierDet) return TMath::Mean(kNdet,fData);
Double_t *ddata = new Double_t[kNdet];
Int_t nPoints = 0;
for (Int_t i=0;i<kNdet;i++) {
if (!(outlierDet->GetValue(i))) {
ddata[nPoints]= fData[nPoints];
nPoints++;
}
}
Double_t mean = TMath::Mean(nPoints,ddata);
delete [] ddata;
return mean;
}
Double_t AliTRDCalDet::GetMedian(AliTRDCalDet * const outlierDet) const
{
if (!outlierDet) return (Double_t) TMath::Median(kNdet,fData);
Double_t *ddata = new Double_t[kNdet];
Int_t nPoints = 0;
for (Int_t i=0;i<kNdet;i++) {
if (!(outlierDet->GetValue(i))) {
ddata[nPoints]= fData[nPoints];
nPoints++;
}
}
Double_t mean = TMath::Median(nPoints,ddata);
delete [] ddata;
return mean;
}
Double_t AliTRDCalDet::GetRMS(AliTRDCalDet * const outlierDet) const
{
if (!outlierDet) return TMath::RMS(kNdet,fData);
Double_t *ddata = new Double_t[kNdet];
Int_t nPoints = 0;
for (Int_t i=0;i<kNdet;i++) {
if (!(outlierDet->GetValue(i))) {
ddata[nPoints]= fData[nPoints];
nPoints++;
}
}
Double_t mean = TMath::RMS(nPoints,ddata);
delete [] ddata;
return mean;
}
Double_t AliTRDCalDet::GetRMSRobust(Double_t robust) const
{
Int_t *index = new Int_t[kNdet];
TMath::Sort((Int_t)kNdet,fData,index);
Double_t reject = (Int_t) (kNdet*(1-robust)/2.0);
if(reject <= 0.0) reject = 0.0;
if(reject >= kNdet) reject = 0.0;
Double_t *ddata = new Double_t[kNdet];
Int_t nPoints = 0;
for (Int_t i=0;i<kNdet;i++) {
Bool_t rej = kFALSE;
for(Int_t k = 0; k < reject; k++){
if(i==index[k]) rej = kTRUE;
if(i==index[kNdet-(k+1)]) rej = kTRUE;
}
if(!rej){
ddata[nPoints]= fData[i];
nPoints++;
}
}
Double_t mean = TMath::RMS(nPoints,ddata);
delete [] ddata;
delete [] index;
return mean;
}
Double_t AliTRDCalDet::GetLTM(Double_t *sigma
, Double_t fraction
, AliTRDCalDet * const outlierDet)
{
Double_t *ddata = new Double_t[kNdet];
Double_t mean=0, lsigma=0;
UInt_t nPoints = 0;
for (Int_t i=0;i<kNdet;i++) {
if (!outlierDet || !(outlierDet->GetValue(i))) {
ddata[nPoints]= fData[nPoints];
nPoints++;
}
}
Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
if (sigma) *sigma=lsigma;
delete [] ddata;
return mean;
}
TH1F * AliTRDCalDet::MakeHisto1Distribution(Float_t min, Float_t max,Int_t type)
{
if (type>=0){
if (type==0){
Float_t mean = GetMean();
Float_t sigma = GetRMS();
Float_t nsigma = TMath::Abs(min);
min = mean-nsigma*sigma;
max = mean+nsigma*sigma;
}
if (type==1){
Float_t mean = GetMedian();
Float_t delta = min;
min = mean-delta;
max = mean+delta;
}
if (type==2){
Double_t sigma;
Float_t mean = GetLTM(&sigma,max);
sigma*=min;
min = mean-sigma;
max = mean+sigma;
}
}
char name[1000];
snprintf(name,1000,"%s CalDet 1Distribution",GetTitle());
TH1F * his = new TH1F(name,name,100, min,max);
for (Int_t idet=0; idet<kNdet; idet++){
his->Fill(GetValue(idet));
}
return his;
}
TH1F * AliTRDCalDet::MakeHisto1DAsFunctionOfDet(Float_t min, Float_t max,Int_t type)
{
if (type>=0){
if (type==0){
Float_t mean = GetMean();
Float_t sigma = GetRMS();
Float_t nsigma = TMath::Abs(min);
min = mean-nsigma*sigma;
max = mean+nsigma*sigma;
}
if (type==1){
Float_t mean = GetMedian();
Float_t delta = min;
min = mean-delta;
max = mean+delta;
}
if (type==2){
Double_t sigma;
Float_t mean = GetLTM(&sigma,max);
sigma*=min;
min = mean-sigma;
max = mean+sigma;
}
}
char name[1000];
snprintf(name,1000,"%s CalDet as function of det",GetTitle());
TH1F * his = new TH1F(name,name,kNdet, 0, kNdet);
for(Int_t det = 0; det< kNdet; det++){
his->Fill(det+0.5,GetValue(det));
}
his->SetMaximum(max);
his->SetMinimum(min);
return his;
}
TH2F *AliTRDCalDet::MakeHisto2DCh(Int_t ch, Float_t min, Float_t max, Int_t type)
{
gStyle->SetPalette(1);
if (type>=0){
if (type==0){
Float_t mean = GetMean();
Float_t sigma = GetRMS();
Float_t nsigma = TMath::Abs(min);
min = mean-nsigma*sigma;
max = mean+nsigma*sigma;
}
if (type==1){
Float_t mean = GetMedian();
Float_t delta = min;
min = mean-delta;
max = mean+delta;
}
if (type==2){
Double_t sigma;
Float_t mean = GetLTM(&sigma,max);
sigma*=min;
min = mean-sigma;
max = mean+sigma;
}
}
AliTRDgeometry *trdGeo = new AliTRDgeometry();
Double_t poslocal[3] = {0.0,0.0,0.0};
Double_t posglobal[3] = {0.0,0.0,0.0};
char name[1000];
snprintf(name,1000,"%s CalDet 2D ch %d",GetTitle(),ch);
TH2F * his = new TH2F(name, name, 400,-400.0,400.0,400,-400.0,400.0);
Int_t offsetch = 6*ch;
for (Int_t isec = 0; isec < kNsect; isec++){
for(Int_t ipl = 0; ipl < kNplan; ipl++){
Int_t det = offsetch+isec*30+ipl;
AliTRDpadPlane *padPlane = trdGeo->GetPadPlane(ipl,ch);
for (Int_t icol=0; icol<padPlane->GetNcols(); icol++){
poslocal[0] = trdGeo->GetTime0(ipl);
poslocal[2] = padPlane->GetRowPos(0);
poslocal[1] = padPlane->GetColPos(icol);
trdGeo->RotateBack(det,poslocal,posglobal);
Int_t binx = 1+TMath::Nint((posglobal[0]+400.0)*0.5);
Int_t biny = 1+TMath::Nint((posglobal[1]+400.0)*0.5);
his->SetBinContent(binx,biny,fData[det]);
}
}
}
his->SetXTitle("x (cm)");
his->SetYTitle("y (cm)");
his->SetStats(0);
his->SetMaximum(max);
his->SetMinimum(min);
delete trdGeo;
return his;
}
TH2F *AliTRDCalDet::MakeHisto2DSmPl(Int_t sm, Int_t pl, Float_t min, Float_t max, Int_t type)
{
gStyle->SetPalette(1);
if (type>=0){
if (type==0){
Float_t mean = GetMean();
Float_t sigma = GetRMS();
Float_t nsigma = TMath::Abs(min);
min = mean-nsigma*sigma;
max = mean+nsigma*sigma;
}
if (type==1){
Float_t mean = GetMedian();
Float_t delta = min;
min = mean-delta;
max = mean+delta;
}
if (type==2){
Double_t sigma;
Float_t mean = GetLTM(&sigma,max);
sigma*=min;
min = mean-sigma;
max = mean+sigma;
}
}
AliTRDgeometry *trdGeo = new AliTRDgeometry();
AliTRDpadPlane *padPlane0 = trdGeo->GetPadPlane(pl,0);
Double_t row0 = padPlane0->GetRow0();
Double_t col0 = padPlane0->GetCol0();
char name[1000];
snprintf(name,1000,"%s CalDet 2D sm %d and pl %d",GetTitle(),sm,pl);
TH2F * his = new TH2F( name, name, 5, -TMath::Abs(row0), TMath::Abs(row0)
, 4,-2*TMath::Abs(col0),2*TMath::Abs(col0));
Int_t offsetsmpl = 30*sm+pl;
for (Int_t k = 0; k < kNcham; k++){
Int_t det = offsetsmpl+k*6;
Int_t kb = kNcham-1-k;
his->SetBinContent(kb+1,2,fData[det]);
his->SetBinContent(kb+1,3,fData[det]);
}
his->SetXTitle("z (cm)");
his->SetYTitle("y (cm)");
his->SetStats(0);
his->SetMaximum(max);
his->SetMinimum(min);
delete trdGeo;
return his;
}
void AliTRDCalDet::Add(Float_t c1)
{
for (Int_t idet = 0; idet < kNdet; idet++) {
fData[idet] += c1;
}
}
void AliTRDCalDet::Multiply(Float_t c1)
{
for (Int_t idet = 0; idet < kNdet; idet++) {
fData[idet] *= c1;
}
}
void AliTRDCalDet::Add(const AliTRDCalDet * calDet, Double_t c1)
{
for (Int_t idet = 0; idet < kNdet; idet++) {
fData[idet] += calDet->GetValue(idet)*c1;
}
}
void AliTRDCalDet::Multiply(const AliTRDCalDet * calDet)
{
for (Int_t idet = 0; idet < kNdet; idet++) {
fData[idet] *= calDet->GetValue(idet);
}
}
void AliTRDCalDet::Divide(const AliTRDCalDet * calDet)
{
Float_t kEpsilon=0.00000000000000001;
for (Int_t idet = 0; idet < kNdet; idet++) {
if (TMath::Abs(calDet->GetValue(idet))>kEpsilon){
fData[idet] /= calDet->GetValue(idet);
}
}
}
Double_t AliTRDCalDet::CalcMean(Bool_t wghtPads, Int_t &calib)
{
Int_t iSM;
Double_t sum = 0.0;
Int_t ndet = 0;
Double_t meanALL = 0.0;
Double_t meanWP = 0.0;
Double_t pads=0.0;
Double_t padsALL=(144*16*24+144*12*6)*18;
Double_t *meanSM = new Double_t[18];
Double_t *meanSMWP = new Double_t[18];
for (Int_t i = 0; i < 18; i++) {
meanSM[i]=0.0;
meanSMWP[i]=0.0;
}
Int_t det = 0;
while(det < 540) {
Float_t val= fData[det];
iSM = (Int_t)(det / (6*5));
pads=(((Int_t) (det % (6 * 5)) / 6) == 2) ? 144*12 : 144*16;
meanALL+=val/540.;
meanSM[iSM]+=val/30.;
meanWP+=val*(pads/padsALL);
meanSMWP[iSM]+=val*(pads/(padsALL/18.));
det++;
}
det=0;
while(det < 540) {
Float_t val= fData[det];
if (( (!wghtPads) &&
(TMath::Abs(val - meanALL) > 0.0001) &&
(TMath::Abs(val - meanSM[(Int_t)(det / (6*5))]) > 0.0001) ) ||
( (wghtPads) &&
(TMath::Abs(val - meanWP) > 0.0001) &&
(TMath::Abs(val - meanSMWP[(Int_t)(det / (6*5) )]) > 0.0001) )
) {
if(val <= 50.) {
sum+=val;
ndet++;
}
}
det++;
}
delete []meanSM;
delete []meanSMWP;
calib=ndet;
return (sum!=0.0 ? sum/ndet : -1);
}
Double_t AliTRDCalDet::CalcMean(Bool_t wghtPads)
{
Int_t calib = 0;
return CalcMean(wghtPads, calib);
}
Double_t AliTRDCalDet::CalcRMS(Bool_t wghtPads, Int_t &calib)
{
Int_t iSM;
Double_t sum = 0.0;
Int_t ndet = 0;
Double_t meanALL = 0.0;
Double_t meanWP = 0.0;
Double_t pads=0.0;
Double_t padsALL=(144*16*24+144*12*6)*18;
Double_t *meanSM = new Double_t[18];
Double_t *meanSMWP = new Double_t[18];
for (Int_t i = 0; i < 18; i++) {
meanSM[i]=0.0;
meanSMWP[i]=0.0;
}
Int_t det = 0;
while(det < 540) {
Float_t val= fData[det];
iSM = (Int_t)(det / (6*5));
pads=(((Int_t) (det % (6 * 5)) / 6) == 2) ? 144*12 : 144*16;
meanALL+=val/540.;
meanSM[iSM]+=val/30.;
meanWP+=val*(pads/padsALL);
meanSMWP[iSM]+=val*(pads/(padsALL/18.));
det++;
}
Double_t mean=0.0;
if(!wghtPads) mean= meanALL;
if(wghtPads) mean= meanWP;
det=0;
while(det < 540) {
Float_t val= fData[det];
if (( (!wghtPads) &&
(TMath::Abs(val - meanALL) > 0.0001) &&
(TMath::Abs(val - meanSM[(Int_t)(det / (6*5))]) > 0.0001) ) ||
( (wghtPads) &&
(TMath::Abs(val - meanWP) > 0.0001) &&
(TMath::Abs(val - meanSMWP[(Int_t)(det / (6*5) )]) > 0.0001) )
) {
if(val <= 50.) {
sum+=(val-mean)*(val-mean);
ndet++;
}
}
det++;
}
delete []meanSM;
delete []meanSMWP;
calib=ndet;
return (sum!=0.0 ? TMath::Sqrt(sum/ndet) : -1);
}
Double_t AliTRDCalDet::CalcRMS(Bool_t wghtPads)
{
Int_t calib = 0;
return CalcRMS(wghtPads, calib);
}
Double_t AliTRDCalDet::GetMeanSM(Bool_t wghtPads, Int_t sector) const
{
Int_t iSM;
Double_t meanALL = 0.0;
Double_t meanWP = 0.0;
Double_t pads=0.0;
Double_t padsALL=(144*16*24+144*12*6)*18;
Double_t *meanSM = new Double_t[18];
Double_t *meanSMWP = new Double_t[18];
for (Int_t i = 0; i < 18; i++) {
meanSM[i]=0.0;
meanSMWP[i]=0.0;
}
Int_t det = 0;
while(det < 540) {
Float_t val= fData[det];
iSM = (Int_t)(det / (6*5));
pads=(((Int_t) (det % (6 * 5)) / 6) == 2) ? 144*12 : 144*16;
meanALL+=val/540.;
meanSM[iSM]+=val/30.;
meanWP+=val*(pads/padsALL);
meanSMWP[iSM]+=val*(pads/(padsALL/18.));
det++;
}
Double_t mean=0.0;
if(!wghtPads) mean= meanSM[sector];
if(wghtPads) mean= meanSMWP[sector];
delete []meanSM;
delete []meanSMWP;
return mean;
}