#include <TLegend.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
+#include <TLinearFitter.h>
#include <TMath.h>
#include <TMatrixT.h>
#include <TVectorT.h>
#include "AliPID.h"
#include "AliLog.h"
#include "AliESDtrack.h"
-#include "AliCDBManager.h"
-#include "AliCDBPath.h"
-#include "AliCDBEntry.h"
-#include "AliGeomManager.h"
#include "AliMathBase.h"
+#include "AliTrackPointArray.h"
#include "AliTRDresolution.h"
#include "AliTRDgeometry.h"
#include "AliTRDReconstructor.h"
#include "AliTRDrecoParam.h"
#include "AliTRDpidUtil.h"
+#include "AliTRDinfoGen.h"
#include "info/AliTRDclusterInfo.h"
,"TRDout2MC"
,"TRD2MC"
};
+Char_t const * AliTRDresolution::fgParticle[11]={
+ " p bar", " K -", " #pi -", " #mu -", " e -",
+ " No PID",
+ " e +", " #mu +", " #pi +", " K +", " p",
+};
+
// Configure segmentation for y resolution/residuals
Int_t const AliTRDresolution::fgkNresYsegm[3] = {
AliTRDgeometry::kNsector
,fIdxPlot(0)
,fIdxFrame(0)
,fPtThreshold(1.)
- ,fReconstructor(NULL)
- ,fGeo(NULL)
+ ,fDyRange(1.5)
,fDBPDG(NULL)
,fGraphS(NULL)
,fGraphM(NULL)
,fIdxPlot(0)
,fIdxFrame(0)
,fPtThreshold(1.)
- ,fReconstructor(NULL)
- ,fGeo(NULL)
+ ,fDyRange(1.5)
,fDBPDG(NULL)
,fGraphS(NULL)
,fGraphM(NULL)
// Default constructor
//
- fReconstructor = new AliTRDReconstructor();
- fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
- fGeo = new AliTRDgeometry();
-
InitFunctorList();
SetSegmentationLevel();
if(fGraphS){fGraphS->Delete(); delete fGraphS;}
if(fGraphM){fGraphM->Delete(); delete fGraphM;}
- delete fGeo;
- delete fReconstructor;
- if(gGeoManager) delete gGeoManager;
if(fCl){fCl->Delete(); delete fCl;}
if(fMCcl){fMCcl->Delete(); delete fMCcl;}
/* if(fTrklt){fTrklt->Delete(); delete fTrklt;}
{
// spatial resolution
- if(!fReconstructor){
- fReconstructor = new AliTRDReconstructor();
- fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
- }
-
AliTRDrecoTask::UserCreateOutputObjects();
InitExchangeContainers();
}
//________________________________________________________
void AliTRDresolution::InitExchangeContainers()
{
- fCl = new TObjArray();
+// Init containers for subsequent tasks (AliTRDclusterResolution)
+
+ fCl = new TObjArray(200);
fCl->SetOwner(kTRUE);
fMCcl = new TObjArray();
fMCcl->SetOwner(kTRUE);
// Execution part
//
- if(!IsInitGeom()){
- // create geometry
- AliCDBManager* ocdb = AliCDBManager::Instance();
- AliCDBEntry* obj = ocdb->Get(AliCDBPath("GRP", "Geometry", "Data"));
- AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
- AliGeomManager::GetNalignable("TRD");
- AliGeomManager::ApplyAlignObjsFromCDB("TRD");
- fGeo = new AliTRDgeometry();
- SetInitGeom();
- }
-
fCl->Delete();
fMCcl->Delete();
-/* fTrklt->Delete();
- fMCtrklt->Delete();*/
AliTRDrecoTask::UserExec(opt);
}
//________________________________________________________
-Bool_t AliTRDresolution::Pulls(Double_t dyz[2], Double_t cov[3], Double_t tilt)
+Bool_t AliTRDresolution::Pulls(Double_t dyz[2], Double_t cov[3], Double_t tilt) const
{
// Helper function to calculate pulls in the yz plane
// using proper tilt rotation
if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
if(!fTracklet->IsOK()) continue;
Float_t x0 = fTracklet->GetX0();
- Float_t dq, dl;
+ Float_t dqdl, dl;
for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
if(!(c = fTracklet->GetClusters(itb))){
if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
}
- dq = fTracklet->GetdQdl(itb, &dl);
+ dqdl = fTracklet->GetdQdl(itb, &dl);
+ if(dqdl<1.e-5) continue;
dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
- (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
+ (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dqdl);
}
// if(!HasMCdata()) continue;
Int_t sgm[3];
Double_t covR[7], cov[3], dy[2], dz[2];
- Float_t pt, x0, y0, z0, dydx, dzdx;
+ Float_t pt, x0, y0, z0, dydx, dzdx, tilt(0.);
+ const AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
+ // do LINEAR track refit if asked by the user
+ // it is the user responsibility to check if B=0T
+ Float_t param[10]; memset(param, 0, 10*sizeof(Float_t));
+ Int_t np(0), nrc(0); AliTrackPoint clusters[300];
+ if(HasTrackRefit()){
+ Bool_t kPrimary(kFALSE);
+ for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
+ if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
+ if(!fTracklet->IsOK()) continue;
+ x0 = fTracklet->GetX0();
+ tilt = fTracklet->GetTilt();
+ AliTRDcluster *c = NULL;
+ fTracklet->ResetClusterIter(kFALSE);
+ while((c = fTracklet->PrevCluster())){
+ Float_t xyz[3] = {c->GetX(), c->GetY(), c->GetZ()};
+ clusters[np].SetCharge(tilt);
+ clusters[np].SetClusterType(0);
+ clusters[np].SetVolumeID(ily);
+ clusters[np].SetXYZ(xyz);
+ np++;
+ }
+ if(fTracklet->IsRowCross()){
+ Float_t xcross(0.), zcross(0.);
+ if(fTracklet->GetEstimatedCrossPoint(xcross, zcross)){
+ clusters[np].SetCharge(tilt);
+ clusters[np].SetClusterType(1);
+ clusters[np].SetVolumeID(ily);
+ clusters[np].SetXYZ(xcross, 0., zcross);
+ np++;
+ nrc++;
+ }
+ }
+ if(fTracklet->IsPrimary()) kPrimary = kTRUE;
+ }
+ if(kPrimary){
+ clusters[np].SetCharge(tilt);
+ clusters[np].SetClusterType(1);
+ clusters[np].SetVolumeID(-1);
+ clusters[np].SetXYZ(0., 0., 0.);
+ np++;
+ }
+ if(!FitTrack(np, clusters, param)) {
+ AliDebug(1, "Linear track Fit failed.");
+ return NULL;
+ }
+ if(HasTrackSelection()){
+ Bool_t kReject(kFALSE);
+ if(fkTrack->GetNumberOfTracklets() != AliTRDgeometry::kNlayer) kReject = kTRUE;
+ if(!kReject && !UseTrack(np, clusters, param)) kReject = kTRUE;
+ if(kReject){
+ AliDebug(1, "Reject track for residuals analysis.");
+ return NULL;
+ }
+ }
+ }
for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
if(!fTracklet->IsOK()) continue;
sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
// retrive the track angle with the chamber
- y0 = fTracklet->GetYref(0);
- z0 = fTracklet->GetZref(0);
- dydx = fTracklet->GetYref(1);
- dzdx = fTracklet->GetZref(1);
+ if(HasTrackRefit()){
+ Float_t par[3];
+ if(!FitTracklet(ily, np, clusters, param, par)) continue;
+ dydx = par[2];//param[3];
+ dzdx = param[4];
+ y0 = par[1] + dydx * (x0 - par[0]);//param[1] + dydx * (x0 - param[0]);
+ z0 = param[2] + dzdx * (x0 - param[0]);
+ } else {
+ y0 = fTracklet->GetYref(0);
+ z0 = fTracklet->GetZref(0);
+ dydx = fTracklet->GetYref(1);
+ dzdx = fTracklet->GetZref(1);
+ }
+ /*printf("RC[%c] Primary[%c]\n"
+ " Fit : y0[%f] z0[%f] dydx[%f] dzdx[%f]\n"
+ " Ref: y0[%f] z0[%f] dydx[%f] dzdx[%f]\n", fTracklet->IsRowCross()?'y':'n', fTracklet->IsPrimary()?'y':'n', y0, z0, dydx, dzdx
+ ,fTracklet->GetYref(0),fTracklet->GetZref(0),fTracklet->GetYref(1),fTracklet->GetZref(1));*/
+ tilt = fTracklet->GetTilt();
fTracklet->GetCovRef(covR);
- Double_t tilt(fTracklet->GetTilt())
- ,t2(tilt*tilt)
+ Double_t t2(tilt*tilt)
,corr(1./(1. + t2))
,cost(TMath::Sqrt(corr));
AliTRDcluster *c = NULL;
((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
// Get z-position with respect to anode wire
- Int_t istk = fGeo->GetStack(c->GetDetector());
- AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
+ Int_t istk = geo->GetStack(c->GetDetector());
+ AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
Float_t row0 = pp->GetRow0();
Float_t d = row0 - zt + pp->GetAnodeWireOffset();
d -= ((Int_t)(2 * d)) / 2.0;
if (d > 0.25) d = 0.5 - d;
- AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
+ AliTRDclusterInfo *clInfo(NULL);
+ clInfo = new AliTRDclusterInfo;
clInfo->SetCluster(c);
Float_t covF[] = {cov[0], cov[1], cov[2]};
clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covF);
}
}
if(clInfoArr) delete clInfoArr;
+
return (TH3S*)arr->At(0);
}
<< "\n";
}
}
-
-
return (TH2I*)arr->At(0);
}
Double_t covR[7]/*, cov[3]*/;
if(DebugLevel()>=3){
- TVectorD dX(12), dY(12), dZ(12), Pt(12), dPt(12), cCOV(12*15);
- fkMC->PropagateKalman(&dX, &dY, &dZ, &Pt, &dPt, &cCOV);
+ TVectorD dX(12), dY(12), dZ(12), vPt(12), dPt(12), cCOV(12*15);
+ fkMC->PropagateKalman(&dX, &dY, &dZ, &vPt, &dPt, &cCOV);
(*DebugStream()) << "MCkalman"
<< "pdg=" << pdg
<< "dx=" << &dX
<< "dy=" << &dY
<< "dz=" << &dZ
- << "pt=" << &Pt
+ << "pt=" << &vPt
<< "dpt=" << &dPt
<< "cov=" << &cCOV
<< "\n";
}
-
- AliTRDReconstructor rec;
+ AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
AliTRDseedV1 tt(*fTracklet);
tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
tt.SetZref(1, dzdx0);
- tt.SetReconstructor(&rec);
+ tt.SetReconstructor(AliTRDinfoGen::Reconstructor());
tt.Fit(1);
x= tt.GetX();y= tt.GetY();z= tt.GetZ();
dydx = tt.GetYfit(1);
<< "\n";
}
- AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
+ AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
//Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
cOut->SaveAs(Form("%s.gif", cOut->GetName()));
- if(!HasMCdata()){
+ if(!HasMCdata() ||
+ (!fGraphS->At(kMCcluster) || !fGraphM->At(kMCcluster) ||
+ !fGraphS->At(kMCtracklet) || !fGraphM->At(kMCtracklet))){
delete cOut;
return;
}
Int_t bmax(h1.GetMaximumBin());
Int_t jBinMin(1), jBinMax(100);
for(Int_t ibin(bmax); ibin--;){
- if(h1.GetBinContent(ibin)==0){
+ if(h1.GetBinContent(ibin)<1.){
jBinMin=ibin; break;
}
}
for(Int_t ibin(bmax); ibin++;){
- if(h1.GetBinContent(ibin)==0){
+ if(h1.GetBinContent(ibin)<1.){
jBinMax=ibin; break;
}
}
//________________________________________________________
void AliTRDresolution::MakeSummaryPlot(TObjArray *a, TH2 *h2)
{
+// Core functionality for MakeSummary function.
+
h2->Reset();
Double_t x,y;
TGraphErrors *g(NULL); TAxis *ax(h2->GetXaxis());
//________________________________________________________
-Char_t const *fgParticle[11]={
- " p bar", " K -", " #pi -", " #mu -", " e -",
- " No PID",
- " e +", " #mu +", " #pi +", " K +", " p",
-};
-const Color_t fgColorS[11]={
-kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
-kGray,
-kRed, kViolet, kMagenta+1, kOrange-3, kOrange
-};
-const Color_t fgColorM[11]={
-kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
-kBlack,
-kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
-};
-const Marker_t fgMarker[11]={
-30, 30, 26, 25, 24,
-28,
-20, 21, 22, 29, 29
-};
Bool_t AliTRDresolution::PostProcess()
{
- //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
+// Fit, Project, Combine, Extract values from the containers filled during execution
+
+ /*fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));*/
if (!fContainer) {
AliError("ERROR: list not available");
return kFALSE;
}
+
+ // define general behavior parameters
+ const Color_t fgColorS[11]={
+ kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
+ kGray,
+ kRed, kViolet, kMagenta+1, kOrange-3, kOrange
+ };
+ const Color_t fgColorM[11]={
+ kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
+ kBlack,
+ kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
+ };
+ const Marker_t fgMarker[11]={
+ 30, 30, 26, 25, 24,
+ 28,
+ 20, 21, 22, 29, 29
+ };
+
TGraph *gm= NULL, *gs= NULL;
if(!fGraphS && !fGraphM){
TObjArray *aM(NULL), *aS(NULL);
Int_t nybins=fgkNresYsegm[fSegmentLevel];
if(expand) nybins*=2;
h = new TH3S(hname, htitle,
- 48, -.48, .48, 60, -1.5, 1.5, nybins, -0.5, nybins-0.5);
+ 48, -.48, .48, // phi
+ 60, -fDyRange, fDyRange, // dy
+ nybins, -0.5, nybins-0.5);// segment
} else h->Reset();
arr->AddAt(h, 0);
sprintf(hname, "%s_%s_YZpull", GetNameId(), name);
const Int_t kNpt(14);
const Int_t kNdpt(150);
const Int_t kNspc = 2*AliPID::kSPECIES+1;
- Float_t Pt=0.1, DPt=-.1, Spc=-5.5;
+ Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
- for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
- for(Int_t i=0; i<kNspc+1; i++,Spc+=1.) binsSpc[i]=Spc;
- for(Int_t i=0; i<kNdpt+1; i++,DPt+=2.e-3) binsDPt[i]=DPt;
+ for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
+ for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
+ for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
// Pt resolution
sprintf(hname, "%s_%s_Pt", GetNameId(), name);
// binnings for plots containing momentum or pt
const Int_t kNpt(14), kNphi(48), kNdy(60);
- Float_t Phi=-.48, Dy=-.3, Pt=0.1;
+ Float_t lPhi=-.48, lDy=-.3, lPt=0.1;
Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
- for(Int_t i=0; i<kNphi+1; i++,Phi+=.02) binsPhi[i]=Phi;
- for(Int_t i=0; i<kNdy+1; i++,Dy+=.01) binsDy[i]=Dy;
- for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
+ for(Int_t i=0; i<kNphi+1; i++,lPhi+=.02) binsPhi[i]=lPhi;
+ for(Int_t i=0; i<kNdy+1; i++,lDy+=.01) binsDy[i]=lDy;
+ for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
// cluster to track residuals/pulls
fContainer->AddAt(arr = new TObjArray(2), kCharge);
return kTRUE;
}
+//________________________________________________________
+Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
+{
+// Generic function to process sigma/mean for 2D plot dy(x)
+
+ if(!h2) {
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
+ return kFALSE;
+ }
+ if(!Int_t(h2->GetEntries())){
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
+ return kFALSE;
+ }
+ if(!g || !g[0]|| !g[1]) {
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
+ return kFALSE;
+ }
+ // prepare
+ TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
+ TF1 f("f", "gaus", ay->GetXmin(), ay->GetXmax());
+ Int_t n(0);
+ if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
+ if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
+ TH1D *h(NULL);
+ if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
+
+ // do actual loop
+ for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
+ Double_t x = ax->GetBinCenter(ix);
+ Double_t ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
+ h = h2->ProjectionY("py", ix, ix);
+ if((n=(Int_t)h->GetEntries())<stat){
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
+ continue;
+ }
+ f.SetParameter(1, 0.);
+ f.SetParameter(2, 3.e-2);
+ h->Fit(&f, "QN");
+ g[0]->SetPoint(np, x, f.GetParameter(1));
+ g[0]->SetPointError(np, ex, f.GetParError(1));
+ g[1]->SetPoint(np, x, f.GetParameter(2));
+ g[1]->SetPointError(np, ex, f.GetParError(2));
+ np++;
+ }
+ return kTRUE;
+}
+
//________________________________________________________
Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
return kTRUE;
}
+//____________________________________________________________________
+Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
+{
+//
+// Fit track with a staight line using the "np" clusters stored in the array "points".
+// The following particularities are stored in the clusters from points:
+// 1. pad tilt as cluster charge
+// 2. pad row cross or vertex constrain as fake cluster with cluster type 1
+// The parameters of the straight line fit are stored in the array "param" in the following order :
+// param[0] - x0 reference radial position
+// param[1] - y0 reference r-phi position @ x0
+// param[2] - z0 reference z position @ x0
+// param[3] - slope dy/dx
+// param[4] - slope dz/dx
+//
+// Attention :
+// Function should be used to refit tracks for B=0T
+//
+
+ if(np<40){
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].", np);
+ return kFALSE;
+ }
+ TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
+
+ Double_t x0(0.);
+ for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
+ x0/=Float_t(np);
+
+ Double_t x, y, z, dx, tilt(0.);
+ for(Int_t ip(0); ip<np; ip++){
+ x = points[ip].GetX(); z = points[ip].GetZ();
+ dx = x - x0;
+ zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
+ }
+ if(zfitter.Eval() != 0) return kFALSE;
+
+ Double_t z0 = zfitter.GetParameter(0);
+ Double_t dzdx = zfitter.GetParameter(1);
+ for(Int_t ip(0); ip<np; ip++){
+ if(points[ip].GetClusterType()) continue;
+ x = points[ip].GetX();
+ dx = x - x0;
+ y = points[ip].GetY();
+ z = points[ip].GetZ();
+ tilt = points[ip].GetCharge();
+ y -= tilt*(-dzdx*dx + z - z0);
+ Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
+ yfitter.AddPoint(&dx, y, 1.);
+ }
+ if(yfitter.Eval() != 0) return kFALSE;
+ Double_t y0 = yfitter.GetParameter(0);
+ Double_t dydx = yfitter.GetParameter(1);
+
+ param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f]\n", x0, y0, z0, dydx, dzdx);
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, AliTrackPoint *points, const Float_t param[10], Float_t par[3])
+{
+//
+// Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
+// See function FitTrack for the data stored in the "clusters" array
+
+// The parameters of the straight line fit are stored in the array "param" in the following order :
+// par[0] - x0 reference radial position
+// par[1] - y0 reference r-phi position @ x0
+// par[2] - slope dy/dx
+//
+// Attention :
+// Function should be used to refit tracks for B=0T
+//
+
+ TLinearFitter yfitter(2, "pol1");
+
+ // grep data for tracklet
+ Double_t x0(0.), x[60], y[60], dy[60];
+ Int_t nly(0);
+ for(Int_t ip(0); ip<np; ip++){
+ if(points[ip].GetClusterType()) continue;
+ if(points[ip].GetVolumeID() != ly) continue;
+ Float_t xt(points[ip].GetX())
+ ,yt(param[1] + param[3] * (xt - param[0]));
+ x[nly] = xt;
+ y[nly] = points[ip].GetY();
+ dy[nly]= y[nly]-yt;
+ x0 += xt;
+ nly++;
+ }
+ if(nly<10){
+ if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].", nly);
+ return kFALSE;
+ }
+ // set radial reference for fit
+ x0 /= Float_t(nly);
+
+ // find tracklet core
+ Double_t mean(0.), sig(1.e3);
+ AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
+
+ // simple cluster error parameterization
+ Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
+
+ // fit tracklet core
+ for(Int_t jly(0); jly<nly; jly++){
+ if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
+ Double_t dx(x[jly]-x0);
+ yfitter.AddPoint(&dx, y[jly], 1.);
+ }
+ if(yfitter.Eval() != 0) return kFALSE;
+ par[0] = x0;
+ par[1] = yfitter.GetParameter(0);
+ par[2] = yfitter.GetParameter(1);
+ return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDresolution::UseTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
+{
+//
+// Global selection mechanism of tracksbased on cluster to fit residuals
+// The parameters are the same as used ni function FitTrack().
+
+ const Float_t kS(0.6), kM(0.2);
+ TH1S h("h1", "", 100, -5.*kS, 5.*kS);
+ Float_t dy, dz, s, m;
+ for(Int_t ip(0); ip<np; ip++){
+ if(points[ip].GetClusterType()) continue;
+ Float_t x0(points[ip].GetX())
+ ,y0(param[1] + param[3] * (x0 - param[0]))
+ ,z0(param[2] + param[4] * (x0 - param[0]));
+ dy=points[ip].GetY() - y0; h.Fill(dy);
+ dz=points[ip].GetZ() - z0;
+ }
+ TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
+ fg.SetParameter(1, 0.);
+ fg.SetParameter(2, 2.e-2);
+ h.Fit(&fg, "QN");
+ m=fg.GetParameter(1); s=fg.GetParameter(2);
+ if(s>kS || TMath::Abs(m)>kM) return kFALSE;
+ return kTRUE;
+}
+
//________________________________________________________
void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
{
}
-//________________________________________________________
-void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
-{
-
- fReconstructor->SetRecoParam(r);
-}
-
-
//________________________________________________________
void AliTRDresolution::SetSegmentationLevel(Int_t l)
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff