From b2dc316dab8aa5f556fb488dff607e39ed411c2b Mon Sep 17 00:00:00 2001 From: abercuci Date: Tue, 25 Nov 2008 22:57:19 +0000 Subject: [PATCH] inserted the first calibration task for reconstruction algorithms - task for cluster error parameterization (first version) - helper class for matching TRD clusters with global tracking info - updates in the tracking resolution task --- TRD/TRDqaRecLinkDef.h | 4 + TRD/libTRDqaRec.pkg | 14 +- TRD/qaRec/AliTRDcheckDetector.h | 6 +- TRD/qaRec/AliTRDclusterResolution.cxx | 113 ++++++ TRD/qaRec/AliTRDclusterResolution.h | 39 ++ .../AliTRDtrackInfo/AliTRDclusterInfo.cxx | 47 +++ TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.h | 81 +++++ TRD/qaRec/AliTRDtrackingResolution.cxx | 340 +++++------------- TRD/qaRec/AliTRDtrackingResolution.h | 18 +- TRD/qaRec/run.C | 14 + 10 files changed, 410 insertions(+), 266 deletions(-) create mode 100644 TRD/qaRec/AliTRDclusterResolution.cxx create mode 100644 TRD/qaRec/AliTRDclusterResolution.h create mode 100644 TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.cxx create mode 100644 TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.h diff --git a/TRD/TRDqaRecLinkDef.h b/TRD/TRDqaRecLinkDef.h index 039196d50b2..d1aa14512eb 100644 --- a/TRD/TRDqaRecLinkDef.h +++ b/TRD/TRDqaRecLinkDef.h @@ -6,6 +6,7 @@ #pragma link off all classes; #pragma link off all functions; +#pragma link C++ class AliTRDclusterInfo+; #pragma link C++ class AliTRDtrackInfo+; #pragma link C++ class AliTRDeventInfo+; #pragma link C++ class AliTRDtrackInfo::AliESDinfo+; @@ -20,4 +21,7 @@ #pragma link C++ class AliTRDpidChecker+; #pragma link C++ class AliTRDpidRefMaker+; +// reconstruction calibration tasks +#pragma link C++ class AliTRDclusterResolution+; + #endif diff --git a/TRD/libTRDqaRec.pkg b/TRD/libTRDqaRec.pkg index 858d4b47843..35b77b75af5 100644 --- a/TRD/libTRDqaRec.pkg +++ b/TRD/libTRDqaRec.pkg @@ -1,16 +1,8 @@ #-*- Mode: Makefile -*- -SRCS= qaRec/AliTRDtrackInfoGen.cxx \ - qaRec/AliTRDrecoTask.cxx \ - qaRec/AliTRDcheckDetector.cxx \ - qaRec/AliTRDpidChecker.cxx \ - qaRec/AliTRDpidRefMaker.cxx \ - qaRec/AliTRDcalibration.cxx \ - qaRec/AliTRDtrackingResolution.cxx \ - qaRec/AliTRDtrackingEfficiency.cxx \ - qaRec/AliTRDtrackingEfficiencyCombined.cxx \ - qaRec/AliTRDtrackInfo/AliTRDtrackInfo.cxx \ - qaRec/AliTRDtrackInfo/AliTRDeventInfo.cxx +ORGSRCS := $(wildcard TRD/qaRec/*.cxx) +ORGSRCS += $(wildcard TRD/qaRec/AliTRDtrackInfo/*.cxx) +SRCS := $(patsubst TRD/%, %, ${ORGSRCS}) HDRS= $(SRCS:.cxx=.h) diff --git a/TRD/qaRec/AliTRDcheckDetector.h b/TRD/qaRec/AliTRDcheckDetector.h index 842041b5b18..9a4223c05e1 100644 --- a/TRD/qaRec/AliTRDcheckDetector.h +++ b/TRD/qaRec/AliTRDcheckDetector.h @@ -1,5 +1,5 @@ -#ifndef __ALITRDCHECKDETECTOR_H__ -#define __ALITRDCHECKDETECTOR_H__ +#ifndef ALITRDCHECKDETECTOR_H +#define ALITRDCHECKDETECTOR_H #ifndef ALITRDRECOTASK_H #include "AliTRDrecoTask.h" @@ -12,7 +12,7 @@ class AliESDHeader; class AliTRDgeometry; class AliTRDReconstructor; class AliTRDrecoParam; - +class AliTRDeventInfo; class AliTRDcheckDetector : public AliTRDrecoTask{ // The Histogram number typedef enum{ diff --git a/TRD/qaRec/AliTRDclusterResolution.cxx b/TRD/qaRec/AliTRDclusterResolution.cxx new file mode 100644 index 00000000000..9a87b92397d --- /dev/null +++ b/TRD/qaRec/AliTRDclusterResolution.cxx @@ -0,0 +1,113 @@ +#include "AliTRDclusterResolution.h" +#include "AliTRDtrackInfo/AliTRDclusterInfo.h" + +#include "AliLog.h" + +#include "TObjArray.h" +#include "TAxis.h" +#include "TH2I.h" +#include "TMath.h" + +ClassImp(AliTRDclusterResolution) + + +//_______________________________________________________ +AliTRDclusterResolution::AliTRDclusterResolution() + : AliTRDrecoTask("CalibClRes", "Calibrate Cluster Resolution for Tracking") + ,fInfo(0x0) +{ +} + +//_______________________________________________________ +AliTRDclusterResolution::~AliTRDclusterResolution() +{ + +} + +//_______________________________________________________ +void AliTRDclusterResolution::ConnectInputData(Option_t *) +{ + fInfo = dynamic_cast(GetInputData(0)); +} + +//_______________________________________________________ +void AliTRDclusterResolution::CreateOutputObjects() +{ + OpenFile(0, "RECREATE"); + fContainer = Histos(); +} + +//_______________________________________________________ +void AliTRDclusterResolution::GetRefFigure(Int_t /*ifig*/) +{ + +} + +//_______________________________________________________ +TObjArray* AliTRDclusterResolution::Histos() +{ + if(fContainer) return fContainer; + fContainer = new TObjArray(kN+1); + //fContainer->SetOwner(kTRUE); + + TAxis at(kNTB, -0.075, (kNTB-.5)*.15); + TAxis ad(kND, 0., .25); + Int_t ih = 0; + for(Int_t id=1; id<=ad.GetNbins(); id++){ + for(Int_t it=1; it<=at.GetNbins(); it++){ + fContainer->AddAt(new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%5.2f-%5.2f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*ad.GetBinCenter(id)- 5.*ad.GetBinWidth(id), 10.*ad.GetBinCenter(id)+ 5.*ad.GetBinWidth(id), 10.*at.GetBinCenter(it)- 5.*at.GetBinWidth(it), 10.*at.GetBinCenter(it)+ 5.*at.GetBinWidth(it)), 30, -.15, .15, 100, -.5, .5), ih++); + } + } + fContainer->AddAt(new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), ih++); + return fContainer; +} + +//_______________________________________________________ +void AliTRDclusterResolution::Exec(Option_t *) +{ + Int_t det; + Float_t x, y, z, q, dy, dydx, dzdx, cov[3]; + TAxis at(kNTB, -0.075, (kNTB-.5)*.15); Int_t it = 0; + TAxis ad(kND, 0., .25); Int_t id = 0; + TH2I *h2 = 0x0; + const AliTRDclusterInfo *cli = 0x0; + TIterator *iter=fInfo->MakeIterator(); + while((cli=dynamic_cast((*iter)()))){ + dy = cli->GetResolution(); + it = at.FindBin(cli->GetDriftLength()); + if(it==0 || it == at.GetNbins()+1){ + AliWarning(Form("Drift length %f outside allowed range", cli->GetDriftLength())); + continue; + } + id = ad.FindBin(cli->GetAnisochronity()); + if(id==0 || id == ad.GetNbins()+1){ + AliWarning(Form("Distance to anode %f outside allowed range", cli->GetAnisochronity())); + continue; + } + if(!(h2 = (TH2I*)fContainer->At((id-1)*kNTB+it-1))){ + AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", (id-1)*kNTB+it-1, id, it, cli->GetDriftLength(), cli->GetAnisochronity())); + continue; + } + + cli->GetGlobalPosition(y, z, dydx, dzdx, &cov[0]); + h2->Fill(dydx, dy); + + // resolution as a function of cluster charge + // only for phi equal exB + if(TMath::Abs(dydx)<.01){ + cli->GetCluster(det, x, y, z, q); + h2 = (TH2I*)fContainer->At(kN); + h2->Fill(TMath::Log(q), dy); + } + } + PostData(0, fContainer); +} + + +//_______________________________________________________ +Bool_t AliTRDclusterResolution::PostProcess() +{ + return kFALSE; +} + + diff --git a/TRD/qaRec/AliTRDclusterResolution.h b/TRD/qaRec/AliTRDclusterResolution.h new file mode 100644 index 00000000000..20ea52bf928 --- /dev/null +++ b/TRD/qaRec/AliTRDclusterResolution.h @@ -0,0 +1,39 @@ +#ifndef ALITRDCLUSTERRESOLUTION_H +#define ALITRDCLUSTERRESOLUTION_H + + +#ifndef ALITRDRECOTASK_H +#include "AliTRDrecoTask.h" +#endif + +class TObjArray; +class AliTRDclusterResolution : public AliTRDrecoTask +{ +public: + enum { + kNTB = 24 + ,kND = 5 + ,kN = 120 + }; + AliTRDclusterResolution(); + virtual ~AliTRDclusterResolution(); + + void ConnectInputData(Option_t *); + void CreateOutputObjects(); + void GetRefFigure(Int_t ifig); + TObjArray* Histos(); + + void Exec(Option_t *); + Bool_t PostProcess(); + void Terminate(Option_t *){}; +private: + AliTRDclusterResolution(const AliTRDclusterResolution&); + AliTRDclusterResolution& operator=(const AliTRDclusterResolution&); + + TObjArray *fInfo; + + ClassDef(AliTRDclusterResolution, 0) // cluster resolution +}; + +#endif + diff --git a/TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.cxx b/TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.cxx new file mode 100644 index 00000000000..1fa98e3032e --- /dev/null +++ b/TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.cxx @@ -0,0 +1,47 @@ +#include "TMath.h" + +#include "AliLog.h" +#include "AliTRDcluster.h" + +#include "AliTRDclusterInfo.h" + +ClassImp(AliTRDclusterInfo) + +//_________________________________________________ +AliTRDclusterInfo::AliTRDclusterInfo() + : TObject() + ,fDet(0xffff) + ,fPdg(0) + ,fLbl(-1) + ,fQ(0.) + ,fX(0.) + ,fY(0.) + ,fZ(0.) + ,fdydx(0.) + ,fdzdx(0.) + ,fXd(0.) + ,fYt(0.) + ,fZt(0.) + ,fdy(0.) + ,fD(0.) +{ + memset(fCov, 0, 3*sizeof(Float_t)); +} + +//_________________________________________________ +void AliTRDclusterInfo::SetCluster(const AliTRDcluster *c) +{ + if(!c) return; + fDet = c->GetDetector(); + fX = c->GetX(); + fY = c->GetY(); + fZ = c->GetZ(); + fQ = TMath::Abs(c->GetQ()); +} + +//_________________________________________________ +void AliTRDclusterInfo::Print(Option_t */*opt*/) const +{ + printf("Det[%3d] X[%7.2f] Y[%7.2f] Z[%7.2f] Q[%7.2f]\n", fDet==0xffff ? -1 : fDet, fX, fY, fZ, fQ); + printf("\tPdg[%d] Lbl[%d] Yt[%7.2f] Zt[%7.2f]\n", fPdg, fLbl, fYt, fZt); +} diff --git a/TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.h b/TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.h new file mode 100644 index 00000000000..95ffb89e439 --- /dev/null +++ b/TRD/qaRec/AliTRDtrackInfo/AliTRDclusterInfo.h @@ -0,0 +1,81 @@ +#ifndef ALITRDCLUSTERINFO_H +#define ALITRDCLUSTERINFO_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + +#ifndef Root_TObject +#include "TObject.h" +#endif + +class AliTRDcluster; +class AliTRDclusterInfo : public TObject +{ +public: + AliTRDclusterInfo(); + + Float_t GetAnisochronity() const {return fD;} + inline void GetCluster(Int_t &det, Float_t &x, Float_t &y, Float_t &z, Float_t &q) const; + void GetMC(Int_t &pdg, Int_t &label) const{ + pdg = fPdg; + label= fLbl; } + void GetGlobalPosition(Float_t &yt, Float_t &zt, Float_t &dydx, Float_t &dzdx, Float_t *cov) const { + dydx = fdydx; + dzdx = fdzdx; + yt = fYt; + zt = fZt; + memcpy(cov, fCov, 3*sizeof(Float_t));} + Float_t GetResolution() const {return fdy;} + Float_t GetDriftLength() const {return fXd;} + + void Print(Option_t *opt="") const; + + void SetAnisochronity(Float_t d) {fD = d;} + void SetCluster(const AliTRDcluster *c=0x0); + void SetMC(Int_t pdg, Int_t label){ + fPdg = pdg; + fLbl = label;} + void SetGlobalPosition(Float_t yt, Float_t zt, Float_t dydx, Float_t dzdx, Float_t *cov=0x0) { + fdydx = dydx; + fdzdx = dzdx; + fYt = yt; + fZt = zt; + if(cov) memcpy(fCov, cov, 3*sizeof(Float_t));} + void SetResolution(Float_t dy) {fdy = dy;} + void SetDriftLength(Float_t d) {fXd = d;} + +private: + UShort_t fDet; // detector + Short_t fPdg; // particle code + Short_t fLbl; // track label (MC) + + Float_t fQ; // cluster charge (REC) + Float_t fX; // x coordinate (REC) + Float_t fY; // y coordinate (REC) + Float_t fZ; // z coordinate (REC) + Float_t fdydx; // slope in phi (MC) + Float_t fdzdx; // slope in theta (MC) + Float_t fXd; // drift length + Float_t fYt; // y coordinate (MC) + Float_t fZt; // z coordinate (MC) + Float_t fCov[3];// covariance matrix in the yz plane (global) + Float_t fdy; // difference in y after tilt correction + Float_t fD; // distance to the anode wire + + ClassDef(AliTRDclusterInfo, 1) // extracted cluster2MC information +}; + + +//_________________________________________________ +inline void AliTRDclusterInfo::GetCluster(Int_t &det, Float_t &x, Float_t &y, Float_t &z, Float_t &q) const +{ + det = fDet; + x = fX; + y = fY; + z = fZ; + q = fQ; +} + +#endif + diff --git a/TRD/qaRec/AliTRDtrackingResolution.cxx b/TRD/qaRec/AliTRDtrackingResolution.cxx index b36f829853d..fcac3329d94 100644 --- a/TRD/qaRec/AliTRDtrackingResolution.cxx +++ b/TRD/qaRec/AliTRDtrackingResolution.cxx @@ -56,6 +56,7 @@ #include #include #include +#include #include #include #include @@ -77,6 +78,7 @@ #include "AliTRDReconstructor.h" #include "AliTRDrecoParam.h" +#include "AliTRDtrackInfo/AliTRDclusterInfo.h" #include "AliTRDtrackInfo/AliTRDtrackInfo.h" #include "AliTRDtrackingResolution.h" @@ -90,11 +92,19 @@ AliTRDtrackingResolution::AliTRDtrackingResolution() ,fGeo(0x0) ,fGraphS(0x0) ,fGraphM(0x0) + ,fClResiduals(0x0) + ,fClResolution(0x0) + ,fTrkltResolution(0x0) { fReconstructor = new AliTRDReconstructor(); fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam()); fGeo = new AliTRDgeometry(); + InitFunctorList(); + + DefineOutput(1+kClusterResidual, TObjArray::Class()); + DefineOutput(1+kClusterResolution, TObjArray::Class()); + DefineOutput(1+kTrackletYResolution, TObjArray::Class()); } //________________________________________________________ @@ -105,6 +115,9 @@ AliTRDtrackingResolution::~AliTRDtrackingResolution() delete fGeo; delete fReconstructor; if(gGeoManager) delete gGeoManager; + if(fClResiduals){fClResiduals->Delete(); delete fClResiduals;} + if(fClResolution){fClResolution->Delete(); delete fClResolution;} + if(fTrkltResolution){fTrkltResolution->Delete(); delete fTrkltResolution;} } @@ -115,205 +128,28 @@ void AliTRDtrackingResolution::CreateOutputObjects() OpenFile(0, "RECREATE"); fContainer = Histos(); + + fClResiduals = new TObjArray(); + fClResiduals->SetOwner(kTRUE); + fClResolution = new TObjArray(); + fClResolution->SetOwner(kTRUE); + fTrkltResolution = new TObjArray(); + fTrkltResolution->SetOwner(kTRUE); } -// //________________________________________________________ -// void AliTRDtrackingResolution::Exec(Option_t *) -// { -// // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire}) -// // angular Resolution: res = Tracklet angle - TrackRef Angle -// -// Int_t nTrackInfos = fTracks->GetEntriesFast(); -// if(fDebugLevel>=2 && nTrackInfos){ -// printf("Event[%d] TrackInfos[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTrackInfos); -// } -// const Int_t kNLayers = AliTRDgeometry::kNlayer; -// Int_t pdg, nly, ncrs; -// Double_t p, dy, theta/*, dphi, dymc, dzmc, dphimc*/; -// Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer -// Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map -// -// AliTRDpadPlane *pp = 0x0; -// AliTrackPoint tr[kNLayers], tk[kNLayers]; -// AliExternalTrackParam *fOp = 0x0; -// AliTRDtrackV1 *fTrack = 0x0; -// AliTRDtrackInfo *fInfo = 0x0; -// for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){ -// // check if ESD and MC-Information are available -// if(!(fInfo = dynamic_cast(fTracks->UncheckedAt(iTI)))) continue; -// if(!(fTrack = fInfo->GetTrack())) continue; -// if(!(fOp = fInfo->GetOuterParam())) continue; -// pdg = fInfo->GetPDG(); -// nly = 0; ncrs = 0; theta = 0.; -// -// if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs()); -// -// p = fOp->P(); -// Int_t npts = 0; -// memset(fP, 0, kNLayers*sizeof(Float_t)); -// memset(fX, 0, kNLayers*sizeof(Float_t)); -// memset(fY, 0, kNLayers*sizeof(Float_t)); -// memset(fZ, 0, kNLayers*sizeof(Float_t)); -// memset(fPhi, 0, kNLayers*sizeof(Float_t)); -// memset(fTheta, 0, kNLayers*sizeof(Float_t)); -// memset(fLayerMap, 0, kNLayers*sizeof(Bool_t)); -// memset(fMCMap, 0, kNLayers*sizeof(Bool_t)); -// AliTRDseedV1 *fTracklet = 0x0; -// for(Int_t iplane = 0; iplane < kNLayers; iplane++){ -// if(!(fTracklet = fTrack->GetTracklet(iplane))) continue; -// if(!fTracklet->IsOK()) continue; -// -// // Book point arrays -// fLayerMap[iplane] = kTRUE; -// tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.); -// tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0)); -// npts++; -// -// if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN()); -// -// // define reference values -// fP[iplane] = p; -// fX[iplane] = fTracklet->GetX0(); -// fY[iplane] = fTracklet->GetYref(0); -// fZ[iplane] = fTracklet->GetZref(0); -// fPhi[iplane] = TMath::ATan(fTracklet->GetYref(1)); -// fTheta[iplane] = TMath::ATan(fTracklet->GetZref(1)); -// -// -// // RESOLUTION (compare to MC) -// if(HasMCdata()){ -// if(fInfo->GetNTrackRefs() >= 2){ -// Double_t pmc, ymc, zmc, phiMC, thetaMC; -// if(Resolution(fTracklet, fInfo, pmc, ymc, zmc, phiMC, thetaMC)){ -// fMCMap[iplane] = kTRUE; -// fP[iplane] = pmc; -// fY[iplane] = ymc; -// fZ[iplane] = zmc; -// fPhi[iplane] = phiMC; -// fTheta[iplane] = thetaMC; -// } -// } -// } -// Float_t phi = fPhi[iplane]*TMath::RadToDeg(); -// theta += fTheta[iplane]; nly++; -// if(fTracklet->GetNChange()) ncrs++; -// -// // Do clusters residuals -// if(!fTracklet->Fit(kFALSE)) continue; -// AliTRDcluster *c = 0x0; -// for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){ -// if(!(c = fTracklet->GetClusters(ic))) continue; -// -// dy = fTracklet->GetYat(c->GetX()) - c->GetY(); -// ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy); -// -// if(fDebugLevel>=2){ -// Float_t q = c->GetQ(); -// // Get z-position with respect to anode wire -// AliTRDSimParam *simParam = AliTRDSimParam::Instance(); -// Int_t det = c->GetDetector(); -// Float_t x = c->GetX(); -// Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]); -// Int_t stack = fGeo->GetStack(det); -// pp = fGeo->GetPadPlane(iplane, stack); -// Float_t row0 = pp->GetRow0(); -// Float_t d = row0 - z + simParam->GetAnodeWireOffset(); -// d -= ((Int_t)(2 * d)) / 2.0; -// if (d > 0.25) d = 0.5 - d; -// -// (*fDebugStream) << "ResidualClusters" -// << "ly=" << iplane -// << "stk=" << stack -// << "pdg=" << pdg -// << "phi=" << fPhi[iplane] -// << "tht=" << fTheta[iplane] -// << "q=" << q -// << "x=" << x -// << "z=" << z -// << "d=" << d -// << "dy=" << dy -// << "\n"; -// } -// } -// pp = 0x0; -// } -// if(nly) theta /= nly; -// if(fDebugLevel>=1){ -// (*fDebugStream) << "TrackStatistics" -// << "nly=" << nly -// << "ncrs=" << ncrs -// << "tht=" << theta -// << "\n"; -// } -// -// -// // // this protection we might drop TODO -// // if(fTrack->GetNumberOfTracklets() < 6) continue; -// // -// // AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr); -// // Int_t iref = 0; -// // for(Int_t ip=0; ipGetTracklet(ip); -// // // recalculate fit based on the new tilt correction -// // fTracklet->Fit(); -// // -// // dy = fTracklet->GetYfit(0) - tr[iref].GetY(); -// // ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy); -// // -// // dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1); -// // ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi); -// // -// // if(HasMCdata()){ -// // dymc = fY[ip] - tr[iref].GetY(); -// // ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc); -// // -// // dzmc = fZ[ip] - tr[iref].GetZ(); -// // ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc); -// // -// // dphimc = fPhi[ip] - fTracklet->GetYfit(1); -// // ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc); -// // } -// // -// // iref++; -// // -// // if(fDebugLevel>=1){ -// // (*fDebugStream) << "RiemannTrack" -// // << "ly=" << ip -// // << "mc=" << fMCMap[ip] -// // << "p=" << fP[ip] -// // << "phi=" << fPhi[ip] -// // << "tht=" << fTheta[ip] -// // << "dy=" << dy -// // << "dphi=" << dphi -// // << "dymc=" << dymc -// // << "dzmc=" << dzmc -// // << "dphimc="<< dphimc -// // << "\n"; -// // } -// // } -// -// // if(!gGeoManager) TGeoManager::Import("geometry.root"); -// // AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr); -// // for(Int_t ip=0; ipAt(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy); -// // dz = cl[ip].GetZ() - tr[ip].GetZ(); -// // if(fDebugLevel>=1){ -// // (*fDebugStream) << "KalmanTrack" -// // << "dy=" << dy -// // << "dz=" << dz -// // /* << "phi=" << phi -// // << "theta=" << theta -// // << "dphi=" << dphi*/ -// // << "\n"; -// // } -// // } -// -// -// } -// PostData(0, fContainer); -// } +//________________________________________________________ +void AliTRDtrackingResolution::Exec(Option_t *opt) +{ + fClResiduals->Delete(); + fClResolution->Delete(); + fTrkltResolution->Delete(); + + AliTRDrecoTask::Exec(opt); + + PostData(1+kClusterResidual, fClResiduals); + PostData(1+kClusterResolution, fClResolution); + PostData(1+kTrackletYResolution, fTrkltResolution); +} //________________________________________________________ TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track) @@ -370,7 +206,7 @@ TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track) Float_t d = row0 - z + simParam->GetAnodeWireOffset(); d -= ((Int_t)(2 * d)) / 2.0; if (d > 0.25) d = 0.5 - d; - + (*fDebugStream) << "ClusterResidual" << "ly=" << ily << "stk=" << istk @@ -383,6 +219,8 @@ TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track) << "d=" << d << "dy=" << dy << "\n"; + + } } } @@ -496,20 +334,18 @@ TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track) d -= ((Int_t)(2 * d)) / 2.0; if (d > 0.25) d = 0.5 - d; Int_t label = fMC->GetLabel(); + + AliTRDclusterInfo *clInfo = new AliTRDclusterInfo; + fClResolution->Add(clInfo); + clInfo->SetCluster(c); + clInfo->SetMC(pdg, label); + clInfo->SetGlobalPosition(yt, zt, dydx, dzdx); + clInfo->SetResolution(dy); + clInfo->SetAnisochronity(d); + clInfo->SetDriftLength(x0-xc); + //clInfo->Print(); (*fDebugStream) << "ClusterResolution" - << "det=" << det - << "pdg=" << pdg - << "dydx=" << dydx - << "dzdx=" << dzdx - << "q=" << q - << "d=" << d - << "dy=" << dy - << "xc=" << xc - << "yc=" << yc - << "zc=" << zc - << "yt=" << yt - << "zt=" << zt - << "lbl=" << label + <<"clInfo.=" << clInfo << "\n"; } } @@ -521,6 +357,7 @@ TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track) //________________________________________________________ void AliTRDtrackingResolution::GetRefFigure(Int_t ifig) { + TBox *b = 0x0; TAxis *ax = 0x0; TGraphErrors *g = 0x0; switch(ifig){ @@ -534,6 +371,9 @@ void AliTRDtrackingResolution::GetRefFigure(Int_t ifig) ax->SetTitle("tg(#phi)"); if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break; g->Draw("pl"); + b = new TBox(-.15, -.5, .15, 1.); + b->SetFillStyle(3002);b->SetFillColor(kGreen); + b->SetLineColor(0); b->Draw(); return; case kClusterResolution: if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break; @@ -545,6 +385,9 @@ void AliTRDtrackingResolution::GetRefFigure(Int_t ifig) g->Draw("apl"); if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break; g->Draw("pl"); + b = new TBox(-.15, -.5, .15, 1.); + b->SetFillStyle(3002);b->SetFillColor(kGreen); + b->SetLineColor(0); b->Draw(); return; case kTrackletYResolution: if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break; @@ -626,13 +469,13 @@ Bool_t AliTRDtrackingResolution::PostProcess() // Clusters residuals h2 = (TH2I *)(fContainer->At(kClusterResidual)); - gm = new TGraphErrors(h2->GetNbinsX()); + gm = new TGraphErrors(); gm->SetLineColor(kBlue); gm->SetMarkerStyle(7); gm->SetMarkerColor(kBlue); gm->SetNameTitle("clm", ""); fGraphM->AddAt(gm, kClusterResidual); - gs = new TGraphErrors(h2->GetNbinsX()); + gs = new TGraphErrors(); gs->SetLineColor(kRed); gs->SetMarkerStyle(23); gs->SetMarkerColor(kRed); @@ -641,16 +484,18 @@ Bool_t AliTRDtrackingResolution::PostProcess() for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){ Double_t phi = h2->GetXaxis()->GetBinCenter(ibin); h = h2->ProjectionY("py", ibin, ibin); + if(h->GetEntries()<100) continue; AdjustF1(h, &f); if(IsVisual()){c->cd(); c->SetLogy();} h->Fit(&f, opt, "", -0.5, 0.5); if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);} - gm->SetPoint(ibin - 1, phi, 10.*f.GetParameter(1)); - gm->SetPointError(ibin - 1, 0., 10.*f.GetParError(1)); - gs->SetPoint(ibin - 1, phi, 10.*f.GetParameter(2)); - gs->SetPointError(ibin - 1, 0., 10.*f.GetParError(2)); + Int_t ip = gm->GetN(); + gm->SetPoint(ip, phi, 10.*f.GetParameter(1)); + gm->SetPointError(ip, 0., 10.*f.GetParError(1)); + gs->SetPoint(ip, phi, 10.*f.GetParameter(2)); + gs->SetPointError(ip, 0., 10.*f.GetParError(2)); } @@ -659,13 +504,13 @@ Bool_t AliTRDtrackingResolution::PostProcess() if(HasMCdata()){ // cluster y resolution h2 = (TH2I*)fContainer->At(kClusterResolution); - gm = new TGraphErrors(h2->GetNbinsX()); + gm = new TGraphErrors(); gm->SetLineColor(kBlue); gm->SetMarkerStyle(7); gm->SetMarkerColor(kBlue); gm->SetNameTitle("clym", ""); fGraphM->AddAt(gm, kClusterResolution); - gs = new TGraphErrors(h2->GetNbinsX()); + gs = new TGraphErrors(); gs->SetLineColor(kRed); gs->SetMarkerStyle(23); gs->SetMarkerColor(kRed); @@ -673,7 +518,7 @@ Bool_t AliTRDtrackingResolution::PostProcess() fGraphS->AddAt(gs, kClusterResolution); for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){ h = h2->ProjectionY("py", iphi, iphi); - if(h->GetEntries()<100.) continue; + if(h->GetEntries()<100) continue; AdjustF1(h, &f); if(IsVisual()){c->cd(); c->SetLogy();} @@ -684,22 +529,22 @@ Bool_t AliTRDtrackingResolution::PostProcess() if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);} Double_t phi = h2->GetXaxis()->GetBinCenter(iphi); - Int_t jphi = iphi -1; - gm->SetPoint(jphi, phi, 10.*f.GetParameter(1)); - gm->SetPointError(jphi, 0., 10.*f.GetParError(1)); - gs->SetPoint(jphi, phi, 10.*f.GetParameter(2)); - gs->SetPointError(jphi, 0., 10.*f.GetParError(2)); + Int_t ip = gm->GetN(); + gm->SetPoint(ip, phi, 10.*f.GetParameter(1)); + gm->SetPointError(ip, 0., 10.*f.GetParError(1)); + gs->SetPoint(ip, phi, 10.*f.GetParameter(2)); + gs->SetPointError(ip, 0., 10.*f.GetParError(2)); } // tracklet y resolution h2 = (TH2I*)fContainer->At(kTrackletYResolution); - gm = new TGraphErrors(h2->GetNbinsX()); + gm = new TGraphErrors(); gm->SetLineColor(kBlue); gm->SetMarkerStyle(7); gm->SetMarkerColor(kBlue); gm->SetNameTitle("trkltym", ""); fGraphM->AddAt(gm, kTrackletYResolution); - gs = new TGraphErrors(h2->GetNbinsX()); + gs = new TGraphErrors(); gs->SetLineColor(kRed); gs->SetMarkerStyle(23); gs->SetMarkerColor(kRed); @@ -707,6 +552,7 @@ Bool_t AliTRDtrackingResolution::PostProcess() fGraphS->AddAt(gs, kTrackletYResolution); for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){ h = h2->ProjectionY("py", iphi, iphi); + if(h->GetEntries()<100) continue; AdjustF1(h, &fb); if(IsVisual()){c->cd(); c->SetLogy();} @@ -714,22 +560,22 @@ Bool_t AliTRDtrackingResolution::PostProcess() if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);} Double_t phi = h2->GetXaxis()->GetBinCenter(iphi); - Int_t jphi = iphi -1; - gm->SetPoint(jphi, phi, 10.*fb.GetParameter(1)); - gm->SetPointError(jphi, 0., 10.*fb.GetParError(1)); - gs->SetPoint(jphi, phi, 10.*fb.GetParameter(2)); - gs->SetPointError(jphi, 0., 10.*fb.GetParError(2)); + Int_t ip = gm->GetN(); + gm->SetPoint(ip, phi, 10.*fb.GetParameter(1)); + gm->SetPointError(ip, 0., 10.*fb.GetParError(1)); + gs->SetPoint(ip, phi, 10.*fb.GetParameter(2)); + gs->SetPointError(ip, 0., 10.*fb.GetParError(2)); } // tracklet z resolution h2 = (TH2I*)fContainer->At(kTrackletZResolution); - gm = new TGraphErrors(h2->GetNbinsX()); + gm = new TGraphErrors(); gm->SetLineColor(kBlue); gm->SetMarkerStyle(7); gm->SetMarkerColor(kBlue); gm->SetNameTitle("trkltzm", ""); fGraphM->AddAt(gm, kTrackletZResolution); - gs = new TGraphErrors(h2->GetNbinsX()); + gs = new TGraphErrors(); gs->SetLineColor(kRed); gs->SetMarkerStyle(23); gs->SetMarkerColor(kRed); @@ -737,6 +583,7 @@ Bool_t AliTRDtrackingResolution::PostProcess() fGraphS->AddAt(gs, kTrackletZResolution); for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){ h = h2->ProjectionY("py", iphi, iphi); + if(h->GetEntries()<100) continue; AdjustF1(h, &fb); if(IsVisual()){c->cd(); c->SetLogy();} @@ -744,22 +591,22 @@ Bool_t AliTRDtrackingResolution::PostProcess() if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);} Double_t phi = h2->GetXaxis()->GetBinCenter(iphi); - Int_t jphi = iphi -1; - gm->SetPoint(jphi, phi, 10.*fb.GetParameter(1)); - gm->SetPointError(jphi, 0., 10.*fb.GetParError(1)); - gs->SetPoint(jphi, phi, 10.*fb.GetParameter(2)); - gs->SetPointError(jphi, 0., 10.*fb.GetParError(2)); + Int_t ip = gm->GetN(); + gm->SetPoint(ip, phi, 10.*fb.GetParameter(1)); + gm->SetPointError(ip, 0., 10.*fb.GetParError(1)); + gs->SetPoint(ip, phi, 10.*fb.GetParameter(2)); + gs->SetPointError(ip, 0., 10.*fb.GetParError(2)); } // tracklet phi resolution h2 = (TH2I*)fContainer->At(kTrackletAngleResolution); - gm = new TGraphErrors(h2->GetNbinsX()); + gm = new TGraphErrors(); gm->SetLineColor(kBlue); gm->SetMarkerStyle(7); gm->SetMarkerColor(kBlue); gm->SetNameTitle("trkltym", ""); fGraphM->AddAt(gm, kTrackletAngleResolution); - gs = new TGraphErrors(h2->GetNbinsX()); + gs = new TGraphErrors(); gs->SetLineColor(kRed); gs->SetMarkerStyle(23); gs->SetMarkerColor(kRed); @@ -767,17 +614,18 @@ Bool_t AliTRDtrackingResolution::PostProcess() fGraphS->AddAt(gs, kTrackletAngleResolution); for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){ h = h2->ProjectionY("py", iphi, iphi); + if(h->GetEntries()<100) continue; if(IsVisual()){c->cd(); c->SetLogy();} h->Fit(&f, opt, "", -0.5, 0.5); if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);} Double_t phi = h2->GetXaxis()->GetBinCenter(iphi); - Int_t jphi = iphi -1; - gm->SetPoint(jphi, phi, f.GetParameter(1)); - gm->SetPointError(jphi, 0., f.GetParError(1)); - gs->SetPoint(jphi, phi, f.GetParameter(2)); - gs->SetPointError(jphi, 0., f.GetParError(2)); + Int_t ip = gm->GetN(); + gm->SetPoint(ip, phi, f.GetParameter(1)); + gm->SetPointError(ip, 0., f.GetParError(1)); + gs->SetPoint(ip, phi, f.GetParameter(2)); + gs->SetPointError(ip, 0., f.GetParError(2)); } } if(c) delete c; diff --git a/TRD/qaRec/AliTRDtrackingResolution.h b/TRD/qaRec/AliTRDtrackingResolution.h index 7c14dd1b449..46200c38d1c 100644 --- a/TRD/qaRec/AliTRDtrackingResolution.h +++ b/TRD/qaRec/AliTRDtrackingResolution.h @@ -52,9 +52,9 @@ public: virtual ~AliTRDtrackingResolution(); void CreateOutputObjects(); - //void Exec(Option_t *); void GetRefFigure(Int_t ifig); TObjArray* Histos(); + void Exec(Option_t *); Bool_t IsVerbose() const {return TESTBIT(fStatus, kVerbose);} Bool_t IsVisual() const {return TESTBIT(fStatus, kVisual);} Bool_t PostProcess(); @@ -74,11 +74,17 @@ private: void AdjustF1(TH1 *h, TF1 *f); private: - UChar_t fStatus; // steer parameter of the task - AliTRDReconstructor *fReconstructor; //! local reconstructor - AliTRDgeometry *fGeo; //! TRD geometry - TObjArray *fGraphS; //! result holder - sigma values - TObjArray *fGraphM; //! result holder - mean values + UChar_t fStatus; // steer parameter of the task + AliTRDReconstructor *fReconstructor; //! local reconstructor + AliTRDgeometry *fGeo; //! TRD geometry + TObjArray *fGraphS; //! result holder - sigma values + TObjArray *fGraphM; //! result holder - mean values + + // calibration containers + TObjArray *fClResiduals; //! + TObjArray *fClResolution; //! + TObjArray *fTrkltResolution;//! + ClassDef(AliTRDtrackingResolution, 1) // tracking resolution task }; #endif diff --git a/TRD/qaRec/run.C b/TRD/qaRec/run.C index 7327cc1a676..50ce720bca5 100644 --- a/TRD/qaRec/run.C +++ b/TRD/qaRec/run.C @@ -196,6 +196,20 @@ void run(Char_t *tasks="ALL", const Char_t *files=0x0, Int_t nmax=-1) // Create containers for input/output mgr->ConnectInput( task, 0, coutput1); mgr->ConnectOutput(task, 0, mgr->CreateContainer(task->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%s.root", task->GetName()))); + + AliAnalysisDataContainer *co1 = mgr->CreateContainer(Form("%sClRez", task->GetName()), TObjArray::Class(), AliAnalysisManager::kExchangeContainer); + mgr->ConnectOutput(task, 1, co1); + + AliAnalysisDataContainer *co2 = mgr->CreateContainer(Form("%sClRes", task->GetName()), TObjArray::Class(), AliAnalysisManager::kExchangeContainer); + mgr->ConnectOutput(task, 2, co2); + + AliAnalysisDataContainer *co3 = mgr->CreateContainer(Form("%sTrkltRes", task->GetName()), TObjArray::Class(), AliAnalysisManager::kExchangeContainer); + mgr->ConnectOutput(task, 3, co3); + + // test reconstruction calibration plugin + mgr->AddTask(task = new AliTRDclusterResolution()); + mgr->ConnectInput(task, 0, co2); + mgr->ConnectOutput(task, 0, mgr->CreateContainer(task->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%s.root", task->GetName()))); } //____________________________________________ -- 2.39.3