#include "TROOT.h"
#include "TAxis.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TObjArray.h"
#include "THnSparse.h"
#include "AliLog.h"
#include "AliTRDcluster.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackletOflHelper.h"
#include "AliTRDtrackV1.h"
#include "AliTRDtrackerV1.h"
#include "AliTRDtransform.h"
#include "AliTRDcheckTRK.h"
ClassImp(AliTRDcheckTRK)
UChar_t AliTRDcheckTRK::fgSteer= 0;
Float_t AliTRDcheckTRK::fgKalmanStep = 2.;
AliTRDcheckTRK::AliTRDcheckTRK()
: AliTRDresolution()
{
SetNameTitle("TRDtrackerSys", "TRD Tracker Systematic");
}
AliTRDcheckTRK::AliTRDcheckTRK(char* name)
: AliTRDresolution(name, kFALSE)
{
SetTitle("TRD Tracker Systematic");
MakePtCalib();
InitFunctorList();
}
AliTRDcheckTRK::~AliTRDcheckTRK()
{
}
TObjArray* AliTRDcheckTRK::Histos()
{
if(!(fContainer = AliTRDresolution::Histos())) return NULL;
THnSparse *H(NULL);
if(!(H = (THnSparseI*)gROOT->FindObject("hRoads"))){
const Char_t *title[kNdim] = {"layer", "charge", fgkTitle[kPt], fgkTitle[kYrez], fgkTitle[kPrez], "#sigma^{*}/<#sigma_{y}> [a.u.]", "n_{cl}"};
const Int_t nbins[kNdim] = {AliTRDgeometry::kNlayer, 2, kNptBins, fgkNbins[kYrez], fgkNbins[kPrez], kNSigmaBins, kNclusters};
const Double_t min[kNdim] = {-0.5, -0.5, -0.5, -1., -5., 0., 8.5},
max[kNdim] = {AliTRDgeometry::kNlayer-0.5, 1.5, kNptBins-0.5, 1., 5., 5., min[6]+kNclusters};
TString st("Tracking Roads Calib;");
for(Int_t idim(0); idim<kNdim; idim++){ st += title[idim]; st+=";";}
H = new THnSparseI("hRoads", st.Data(), kNdim, nbins, min, max);
} else H->Reset();
fContainer->AddAt(H, fContainer->GetEntries());
return fContainer;
}
TH1* AliTRDcheckTRK::PlotTrack(const AliTRDtrackV1 *track)
{
if(track) fkTrack = track;
if(!fkTrack){
AliDebug(4, "No Track defined.");
return NULL;
}
AliTRDtrackV1 lt(*fkTrack);
if(!PropagateKalman(lt, UseITS()?fkESD->GetITSoutParam():fkESD->GetTPCoutParam())) return NULL;
PlotCluster(<);
PlotTracklet(<);
PlotTrackIn(<);
DoRoads(<);
return NULL;
}
void AliTRDcheckTRK::MakePtCalib(Float_t pt0, Float_t dpt)
{
for(Int_t j(0); j<=kNptBins; j++){
pt0+=(TMath::Exp(j*j*dpt)-1.);
fPtBinCalib[j]=pt0;
}
}
Int_t AliTRDcheckTRK::GetPtBinCalib(Float_t pt)
{
Int_t ipt(-1);
while(ipt<kNptBins){
if(pt<fPtBinCalib[ipt+1]) break;
ipt++;
}
return ipt;
}
TH1* AliTRDcheckTRK::DoRoads(const AliTRDtrackV1 *track)
{
if(track) fkTrack = track;
if(!fkTrack){
AliDebug(4, "No Track defined.");
return NULL;
}
if(TMath::Abs(fkESD->GetTOFbc())>1){
AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
return NULL;
}
THnSparse *H(NULL);
if(!fContainer || !(H = (THnSparse*)fContainer->At(3))){
AliWarning("No output container defined.");
return NULL;
}
Double_t val[kNdim];
AliTRDseedV1 *fTracklet(NULL);
AliTRDtrackletOflHelper helper; TObjArray cl(AliTRDseedV1::kNclusters); cl.SetOwner(kFALSE);
for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()) continue;
Int_t ipt(GetPtBinCalib(fTracklet->GetPt()));
if(ipt<0) continue;
val[0] = il;
val[1] = fkTrack->Charge()<0?0:1;
val[2] = ipt;
Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)),
dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)),
dydx(fTracklet->GetYfit(1)),
tilt(fTracklet->GetTilt());
val[3] = dyt - dzt*tilt;
dydx+= tilt*fTracklet->GetZref(1);
val[4] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
fTracklet->ResetClusterIter(kTRUE); AliTRDcluster *c(NULL);
while((c = fTracklet->NextCluster())) cl.AddLast(c);
helper.Init(AliTRDtransform::Geometry().GetPadPlane(fTracklet->GetDetector()), &cl);
Double_t r, y, s; helper.GetRMS(r, y, s, fTracklet->GetX0());
val[5] = s/helper.GetSyMean();
val[6] = fTracklet->GetN();
H->Fill(val);
}
return NULL;
}
Bool_t AliTRDcheckTRK::PropagateKalman(AliTRDtrackV1 &t, AliExternalTrackParam *ref)
{
Int_t ntracklets(t.GetNumberOfTracklets());
if(!ntracklets){
printf("E - AliTRDcheckTRK::PropagateKalman :: No tracklets attached to track.\n");
return kFALSE;
}
AliTRDseedV1 *tr(NULL);
if(!t.GetTrackIn()){
printf("E - AliTRDcheckTRK::PropagateKalman :: Track did not entered TRD fiducial volume.\n");
return kFALSE;
}
if(!ref){
return kFALSE;
}
if(ref->Pt()<1.e-3) return kFALSE;
Float_t prod(t.GetBz()*t.Charge());
AliTRDtrackV1 tt;
tt.Set(ref->GetX(), ref->GetAlpha(), ref->GetParameter(), ref->GetCovariance());
tt.SetMass(t.GetMass());
tt.SetTrackOut(t.GetTrackOut());
if(UseITS()){
if(!tt.Rotate((Int_t(ref->GetAlpha()/AliTRDgeometry::GetAlpha()) +(ref->GetAlpha()>0.?1:-1)* 0.5)*AliTRDgeometry::GetAlpha()-ref->GetAlpha())) return kFALSE;
}
for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
if(!(tr = t.GetTracklet(ily))) continue;
Int_t det(tr->GetDetector());
if(HasClRecalibrate()){
AliTRDtransform trans(det);
AliTRDcluster *c(NULL);
Float_t exb, vd, t0, s2, dl, dt; tr->GetCalibParam(exb, vd, t0, s2, dl, dt);
tr->ResetClusterIter(kFALSE);
while((c = tr->PrevCluster())){
if(!trans.Transform(c)){
printf("W - AliTRDcheckTRK::PropagateKalman :: Transform() failed for Det[%03d %02d_%d_%d]\n", det,
AliTRDgeometry::GetSector(det), AliTRDgeometry::GetStack(det), AliTRDgeometry::GetLayer(det));
break;
}
}
}
if(HasTrkltRefit()){
TGeoHMatrix *matrix = AliTRDgeometry::GetClusterMatrix(det);
if(!tr->FitRobust(AliTRDgeometry::GetPadPlane(det), matrix, tt.GetBz(), tt.Charge())) printf("W - AliTRDcheckTRK::PropagateKalman :: FitRobust() failed for Det[%03d]\n", det);
else tr->SetXYZ(matrix);
}
if(!AliTRDtrackerV1::PropagateToX(tt, tr->GetX0(), fgKalmanStep)) continue;
if(!tt.GetTrackIn()) tt.SetTrackIn();
tr->Update(&tt);
if(HasKalmanUpdate()){
Double_t x(tr->GetX0()),
p[2] = { tr->GetYfit(0), tr->GetZfit(0)},
covTrklt[3];
tr->GetCovAt(x, covTrklt);
if(!((AliExternalTrackParam&)tt).Update(p, covTrklt)) continue;
tt.SetTracklet(tr, 0);
tt.SetNumberOfClusters();
tt.UpdateChi2(((AliExternalTrackParam)tt).GetPredictedChi2(p, covTrklt));
}
}
t.~AliTRDtrackV1();
new(&t) AliTRDtrackV1(tt);
return kTRUE;
}