#include <TSystem.h>
#include <TPDGCode.h>
#include <TObjArray.h>
+#include <TH3.h>
#include <TH2.h>
#include <TH1.h>
#include <TF1.h>
#include <TCanvas.h>
+#include <TGaxis.h>
#include <TBox.h>
#include <TProfile.h>
#include <TGraphErrors.h>
+#include <TGraphAsymmErrors.h>
#include <TMath.h>
+#include <TMatrixT.h>
+#include <TVectorT.h>
#include "TTreeStream.h"
#include "TGeoManager.h"
#include "AliTrackReference.h"
#include "AliTrackPointArray.h"
#include "AliCDBManager.h"
+#include "AliPID.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDCommonParam.h"
#include "AliTRDSimParam.h"
#include "AliTRDgeometry.h"
#include "AliTRDpadPlane.h"
#include "AliTRDReconstructor.h"
#include "AliTRDrecoParam.h"
-#include "AliTRDtrackInfo/AliTRDclusterInfo.h"
-#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
+#include "info/AliTRDclusterInfo.h"
+#include "info/AliTRDtrackInfo.h"
#include "AliTRDresolution.h"
ClassImp(AliTRDresolution)
+UChar_t AliTRDresolution::fNElements[kNhistos] = {
+ 2, 2, 5, 5,
+ 2, 5, 12, 2, 11
+};
+Char_t *AliTRDresolution::fAxTitle[46][3] = {
+ // Charge
+ {"x [cm]", "I_{mpv}", "x/x_{0}"}
+ ,{"x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
+ // Clusters to Kalman
+ ,{"tg(#phi)", "#mu_{y}^{cl} [#mum]", "#sigma_{y}^{cl} [#mum]"}
+ ,{"tg(#phi)", "PULL: #mu_{y}^{cl}", "PULL: #sigma_{y}^{cl}"}
+ // TRD tracklet to Kalman fit
+ ,{"tg(#phi)", "#mu_{y}^{trklt} [#mum]", "#sigma_{y}^{trklt} [#mum]"}
+ ,{"tg(#phi)", "PULL: #mu_{y}^{trklt}", "PULL: #sigma_{y}^{trklt}"}
+ ,{"tg(#theta)", "#mu_{z}^{trklt} [#mum]", "#sigma_{z}^{trklt} [#mum]"}
+ ,{"tg(#theta)", "PULL: #mu_{z}^{trklt}", "PULL: #sigma_{z}^{trklt}"}
+ ,{"tg(#phi)", "#mu_{#phi}^{trklt} [mrad]", "#sigma_{#phi}^{trklt} [mrad]"}
+ // TPC track 2 first TRD tracklet
+ ,{"tg(#phi)", "#mu_{y}^{TPC trklt} [#mum]", "#sigma_{y}^{TPC trklt} [#mum]"}
+ ,{"tg(#phi)", "PULL: #mu_{y}^{TPC trklt}", "PULL: #sigma_{y}^{TPC trklt}"}
+ ,{"tg(#theta)", "#mu_{z}^{TPC trklt} [#mum]", "#sigma_{z}^{TPC trklt} [#mum]"}
+ ,{"tg(#theta)", "PULL: #mu_{z}^{TPC trklt}", "PULL: #sigma_{z}^{TPC trklt}"}
+ ,{"tg(#phi)", "#mu_{#phi}^{TPC trklt} [mrad]", "#sigma_{#phi}^{TPC trklt} [mrad]"}
+ // MC cluster
+ ,{"tg(#phi)", "MC: #mu_{y}^{cl} [#mum]", "MC: #sigma_{y}^{cl} [#mum]"}
+ ,{"tg(#phi)", "MC PULL: #mu_{y}^{cl}", "MC PULL: #sigma_{y}^{cl}"}
+ // MC tracklet
+ ,{"tg(#phi)", "MC: #mu_{y}^{trklt} [#mum]", "MC: #sigma_{y}^{trklt} [#mum]"}
+ ,{"tg(#phi)", "MC PULL: #mu_{y}^{trklt}", "MC PULL: #sigma_{y}^{trklt}"}
+ ,{"tg(#theta)", "MC: #mu_{z}^{trklt} [#mum]", "MC: #sigma_{z}^{trklt} [#mum]"}
+ ,{"tg(#theta)", "MC PULL: #mu_{z}^{trklt}", "MC PULL: #sigma_{z}^{trklt}"}
+ ,{"tg(#phi)", "MC: #mu_{#phi}^{trklt} [mrad]", "MC: #sigma_{#phi}^{trklt} [mrad]"}
+ // MC track TPC
+ ,{"tg(#phi)", "MC: #mu_{y}^{TPC} [#mum]", "MC: #sigma_{y}^{TPC} [#mum]"}
+ ,{"tg(#phi)", "MC PULL: #mu_{y}^{TPC}", "MC PULL: #sigma_{y}^{TPC}"}
+ ,{"tg(#theta)", "MC: #mu_{z}^{TPC} [#mum]", "MC: #sigma_{z}^{TPC} [#mum]"}
+ ,{"tg(#theta)", "MC PULL: #mu_{z}^{TPC}", "MC PULL: #sigma_{z}^{TPC}"}
+ ,{"tg(#phi)", "MC: #mu_{#phi}^{TPC} [mrad]", "MC: #sigma_{#phi}^{TPC} [mrad]"}
+ ,{"tg(#phi)", "MC PULL: #mu_{snp}^{TPC}", "MC PULL: #sigma_{snp}^{TPC}"}
+ ,{"tg(#theta)", "MC: #mu_{#theta}^{TPC} [mrad]", "MC: #sigma_{#theta}^{TPC} [mrad]"}
+ ,{"tg(#theta)", "MC PULL: #mu_{tgl}^{TPC}", "MC PULL: #sigma_{tgl}^{TPC}"}
+ ,{"p_{t}^{MC} [GeV/c]", "MC: #mu^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]", "MC: #sigma^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]"}
+ ,{"1/p_{t}^{MC} [c/GeV]", "MC PULL: #mu_{1/p_{t}}^{TPC}", "MC PULL: #sigma_{1/p_{t}}^{TPC}"}
+ ,{"p^{MC} [GeV/c]", "MC: #mu^{TPC}(#Deltap/p^{MC}) [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
+ ,{"p^{MC} [GeV/c]", "MC PULL: #mu^{TPC}(#Deltap/#sigma_{p})", "MC PULL: #sigma^{TPC}(#Deltap/#sigma_{p})"}
+ // MC track TOF
+ ,{"tg(#theta)", "MC: #mu_{z}^{TOF} [#mum]", "MC: #sigma_{z}^{TOF} [#mum]"}
+ ,{"tg(#theta)", "MC PULL: #mu_{z}^{TOF}", "MC PULL: #sigma_{z}^{TOF}"}
+ // MC track in TRD
+ ,{"tg(#phi)", "MC: #mu_{y}^{Trk} [#mum]", "MC: #sigma_{y}^{Trk} [#mum]"}
+ ,{"tg(#phi)", "MC PULL: #mu_{y}^{Trk}", "MC PULL: #sigma_{y}^{Trk}"}
+ ,{"tg(#theta)", "MC: #mu_{z}^{Trk} [#mum]", "MC: #sigma_{z}^{Trk} [#mum]"}
+ ,{"tg(#theta)", "MC PULL: #mu_{z}^{Trk}", "MC PULL: #sigma_{z}^{Trk}"}
+ ,{"tg(#phi)", "MC: #mu_{#phi}^{Trk} [mrad]", "MC: #sigma_{#phi}^{Trk} [mrad]"}
+ ,{"tg(#phi)", "MC PULL: #mu_{snp}^{Trk}", "MC PULL: #sigma_{snp}^{Trk}"}
+ ,{"tg(#theta)", "MC: #mu_{#theta}^{Trk} [mrad]", "MC: #sigma_{#theta}^{Trk} [mrad]"}
+ ,{"tg(#theta)", "MC PULL: #mu_{tgl}^{Trk}", "MC PULL: #sigma_{tgl}^{Trk}"}
+ ,{"p_{t}^{MC} [GeV/c]", "MC: #mu^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]", "MC: #sigma^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]"}
+ ,{"1/p_{t}^{MC} [c/GeV]", "MC PULL: #mu_{1/p_{t}}^{Trk}", "MC PULL: #sigma_{1/p_{t}}^{Trk}"}
+ ,{"p^{MC} [GeV/c]", "MC: #mu^{Trk}(#Deltap/p^{MC}) [%]", "MC: #sigma^{Trk}(#Deltap/p^{MC}) [%]"}
+};
//________________________________________________________
AliTRDresolution::AliTRDresolution()
- :AliTRDrecoTask("Resolution", "Spatial and Momentum Resolution")
+ :AliTRDrecoTask("resolution", "Spatial and momentum TRD resolution checker")
,fStatus(0)
+ ,fIdxPlot(0)
,fReconstructor(0x0)
,fGeo(0x0)
,fGraphS(0x0)
PostData(4, fMCtrklt);
}
+//________________________________________________________
+TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
+{
+ if(track) fTrack = track;
+ if(!fTrack){
+ AliWarning("No Track defined.");
+ return 0x0;
+ }
+ TObjArray *arr = 0x0;
+ if(!(arr = ((TObjArray*)fContainer->At(kCharge)))){
+ AliWarning("No output container defined.");
+ return 0x0;
+ }
+ TH3S* h = 0x0;
+
+ AliTRDseedV1 *fTracklet = 0x0;
+ AliTRDcluster *c = 0x0;
+ for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
+ if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
+ if(!fTracklet->IsOK()) continue;
+ Float_t x0 = fTracklet->GetX0();
+ Float_t dq, 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);
+ dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
+ (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
+ }
+
+// if(!HasMCdata()) continue;
+// UChar_t s;
+// Float_t pt0, y0, z0, dydx0, dzdx0;
+// if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
+
+ }
+ return h;
+}
+
+
//________________________________________________________
TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
{
AliWarning("No Track defined.");
return 0x0;
}
- TH1 *h = 0x0;
- if(!(h = ((TH2I*)fContainer->At(kCluster)))){
- AliWarning("No output histogram defined.");
+ TObjArray *arr = 0x0;
+ if(!(arr = ((TObjArray*)fContainer->At(kCluster)))){
+ AliWarning("No output container defined.");
return 0x0;
}
- Double_t cov[3];
+ Double_t cov[7];
Float_t x0, y0, z0, dy, dydx, dzdx;
AliTRDseedV1 *fTracklet = 0x0;
for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
yc -= tilt*(zc-zt); // tilt correction
dy = yt - yc;
- h->Fill(dydx, dy/TMath::Sqrt(cov[0] + c->GetSigmaY2()));
+ //Float_t sx2 = dydx*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2;
+ Float_t sy2 = c->GetSigmaY2();
+ if(sy2<=0.) continue;
+ ((TH2I*)arr->At(0))->Fill(dydx, dy);
+ ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2));
if(fDebugLevel>=1){
// Get z-position with respect to anode wire
}
}
}
- return h;
+ return (TH2I*)arr->At(0);
}
//________________________________________________________
TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
{
+// Plot normalized residuals for tracklets to track.
+//
+// We start from the result that if X=N(|m|, |Cov|)
+// BEGIN_LATEX
+// (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
+// END_LATEX
+// in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
+// reference position.
if(track) fTrack = track;
if(!fTrack){
AliWarning("No Track defined.");
return 0x0;
}
- TH1 *h = 0x0;
- if(!(h = ((TH2I*)fContainer->At(kTrackletY)))){
- AliWarning("No output histogram defined.");
+ TObjArray *arr = 0x0;
+ if(!(arr = (TObjArray*)fContainer->At(kTrackTRD ))){
+ AliWarning("No output container defined.");
return 0x0;
}
- Double_t cov[3], covR[3];
- Float_t x, y0, dx, dy, dydx;
+ Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
+ Float_t x, dx, dy, dz;
AliTRDseedV1 *fTracklet = 0x0;
for(Int_t il=AliTRDgeometry::kNlayer; il--;){
if(!(fTracklet = fTrack->GetTracklet(il))) continue;
if(!fTracklet->IsOK()) continue;
- y0 = fTracklet->GetYref(0);
- dydx = fTracklet->GetYref(1);
x = fTracklet->GetX();
dx = fTracklet->GetX0() - x;
- dy = y0-dx*dydx - fTracklet->GetY();
+ // compute dy^2 and dz^2
+ dy = fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
+ dz = fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
+ // compute covariance matrix
fTracklet->GetCovAt(x, cov);
fTracklet->GetCovRef(covR);
- h->Fill(dydx, dy/TMath::Sqrt(cov[0] + covR[0]));
+ cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
+/* // Correct PULL calculation by considering off
+ // diagonal elements in the covariance matrix
+ // compute square root matrix
+ if(AliTRDseedV1::GetCovInv(cov, inv)==0.) continue;
+ if(AliTRDseedV1::GetCovSqrt(inv, sqr)<0.) continue;
+ Double_t y = sqr[0]*dy+sqr[1]*dz;
+ Double_t z = sqr[1]*dy+sqr[2]*dz;
+ ((TH3*)h)->Fill(y, z, fTracklet->GetYref(1));*/
+
+ ((TH2I*)arr->At(0))->Fill(fTracklet->GetYref(1), dy);
+ ((TH2I*)arr->At(1))->Fill(fTracklet->GetYref(1), dy/TMath::Sqrt(cov[0]));
+ ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), TMath::ATan((fTracklet->GetYref(1)-fTracklet->GetYfit(1))/(1-fTracklet->GetYref(1)*fTracklet->GetYfit(1))));
+ if(!fTracklet->IsRowCross()) continue;
+ ((TH2I*)arr->At(2))->Fill(fTracklet->GetZref(1), dz);
+ ((TH2I*)arr->At(3))->Fill(fTracklet->GetZref(1), dz/TMath::Sqrt(cov[2]));
}
- return h;
+
+
+ return (TH2I*)arr->At(0);
}
+
//________________________________________________________
-TH1* AliTRDresolution::PlotTrackletPhi(const AliTRDtrackV1 *track)
+TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
{
+// Store resolution/pulls of Kalman before updating with the TRD information
+// at the radial position of the first tracklet. The following points are used
+// for comparison
+// - the (y,z,snp) of the first TRD tracklet
+// - the (y, z, snp, tgl, pt) of the MC track reference
+//
+// Additionally the momentum resolution/pulls are calculated for usage in the
+// PID calculation.
+
if(track) fTrack = track;
if(!fTrack){
AliWarning("No Track defined.");
return 0x0;
}
+ AliExternalTrackParam *tin = 0x0;
+ if(!(tin = fTrack->GetTrackLow())){
+ AliWarning("Track did not entered TRD fiducial volume.");
+ return 0x0;
+ }
TH1 *h = 0x0;
- if(!(h = ((TH2I*)fContainer->At(kTrackletPhi)))){
- AliWarning("No output histogram defined.");
+
+ Double_t x = tin->GetX();
+ AliTRDseedV1 *tracklet = 0x0;
+ for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
+ if(!(tracklet = fTrack->GetTracklet(ily))) continue;
+ break;
+ }
+ if(!tracklet || TMath::Abs(x-tracklet->GetX())>1.e-3){
+ AliWarning("Tracklet did not match TRD entrance.");
return 0x0;
}
+ const Int_t kNPAR(5);
+ Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
+ Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
+ Double_t cov[3]; tracklet->GetCovAt(x, cov);
+
+ // define sum covariances
+ TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
+ Double_t *pc = &covR[0], *pp = &parR[0];
+ for(Int_t ir=0; ir<kNPAR; ir++, pp++){
+ PAR(ir) = (*pp);
+ for(Int_t ic = 0; ic<=ir; ic++,pc++){
+ COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
+ }
+ }
+ PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
+ //COV.Print(); PAR.Print();
+
+ //TODO Double_t dydx = TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
+ Double_t dy = parR[0] - tracklet->GetY();
+ TObjArray *arr = (TObjArray*)fContainer->At(kTrackTPC);
+ ((TH2I*)arr->At(0))->Fill(tracklet->GetYref(1), dy);
+ ((TH2I*)arr->At(1))->Fill(tracklet->GetYref(1), dy/TMath::Sqrt(COV(0,0)+cov[0]));
+ if(tracklet->IsRowCross()){
+ Double_t dz = parR[1] - tracklet->GetZ();
+ ((TH2I*)arr->At(2))->Fill(tracklet->GetZref(1), dz);
+ ((TH2I*)arr->At(3))->Fill(tracklet->GetZref(1), dz/TMath::Sqrt(COV(1,1)+cov[2]));
+ }
+ Double_t dphi = TMath::ASin(PAR[2])-TMath::ATan(tracklet->GetYfit(1)); ((TH2I*)arr->At(4))->Fill(tracklet->GetYref(1), dphi);
- AliTRDseedV1 *tracklet = 0x0;
- for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
- if(!(tracklet = fTrack->GetTracklet(il))) continue;
- if(!tracklet->IsOK()) continue;
- h->Fill(tracklet->GetYref(1), TMath::ATan(tracklet->GetYref(1))-TMath::ATan(tracklet->GetYfit(1)));
+
+ // register reference histo for mini-task
+ h = (TH2I*)arr->At(0);
+
+ if(fDebugLevel>=1){
+ (*fDebugStream) << "trackIn"
+ << "x=" << x
+ << "P=" << &PAR
+ << "C=" << &COV
+ << "\n";
+
+ Double_t y = tracklet->GetY();
+ Double_t z = tracklet->GetZ();
+ (*fDebugStream) << "trackletIn"
+ << "y=" << y
+ << "z=" << z
+ << "Vy=" << cov[0]
+ << "Cyz=" << cov[1]
+ << "Vz=" << cov[2]
+ << "\n";
+ }
+
+
+ if(!HasMCdata()) return h;
+ UChar_t s;
+ Float_t dx, pt0, x0=tracklet->GetX0(), y0, z0, dydx0, dzdx0;
+ if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
+ // translate to reference radial position
+ dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
+ Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
+ //Fill MC info
+ TVectorD PARMC(kNPAR);
+ PARMC[0]=y0; PARMC[1]=z0;
+ PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
+ PARMC[4]=1./pt0;
+
+// TMatrixDSymEigen eigen(COV);
+// TVectorD evals = eigen.GetEigenValues();
+// TMatrixDSym evalsm(kNPAR);
+// for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
+// TMatrixD evecs = eigen.GetEigenVectors();
+// TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
+
+ // fill histos
+ arr = (TObjArray*)fContainer->At(kMCtrackTPC);
+ // y resolution/pulls
+ ((TH2I*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0]);
+ ((TH2I*)arr->At(1))->Fill(dydx0, (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)));
+ // z resolution/pulls
+ ((TH2I*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1]);
+ ((TH2I*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
+ // phi resolution/snp pulls
+ ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
+ ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
+ // theta resolution/tgl pulls
+ ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
+ ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
+ // pt resolution\\1/pt pulls\\p resolution/pull
+ for(Int_t is=AliPID::kSPECIES; is--;){
+ if(TMath::Abs(fMC->GetPDG())!=AliPID::ParticleCode(is)) continue;
+ ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., is);
+ ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), is);
+
+ Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
+ Float_t sp;
+ p = tracklet->GetMomentum(&sp);
+ ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., is);
+ ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/sp, is);
+ break;
+ }
+
+ // fill debug for MC
+ if(fDebugLevel>=1){
+ (*fDebugStream) << "trackInMC"
+ << "P=" << &PARMC
+ << "\n";
}
return h;
}
-
//________________________________________________________
-TH1* AliTRDresolution::PlotResolution(const AliTRDtrackV1 *track)
+TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
{
if(!HasMCdata()){
AliWarning("No MC defined. Results will not be available.");
AliWarning("No Track defined.");
return 0x0;
}
+ TObjArray *arr = 0x0;
TH1 *h = 0x0;
UChar_t s;
Int_t pdg = fMC->GetPDG(), det=-1;
Int_t label = fMC->GetLabel();
- Double_t x, y, z;
- Float_t p, pt, x0, y0, z0, dx, dy, dz, dydx, dzdx;
- Double_t covR[3];
+ Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
+ Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
+ Double_t covR[7]/*, cov[3]*/;
if(fDebugLevel>=1){
Double_t DX[12], DY[12], DZ[12], DPt[12], COV[12][15];
<< "\n";
}
+ AliTRDReconstructor rec;
AliTRDseedV1 *fTracklet = 0x0;
for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
- if(!(fTracklet = fTrack->GetTracklet(ily)) ||
- !fTracklet->IsOK()) continue;
+ if(!(fTracklet = fTrack->GetTracklet(ily)))/* ||
+ !fTracklet->IsOK())*/ continue;
det = fTracklet->GetDetector();
x0 = fTracklet->GetX0();
//radial shift with respect to the MC reference (radial position of the pad plane)
x= fTracklet->GetX();
+ if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
+ xAnode = fTracklet->GetX0();
+
+ // MC track position at reference radial position
dx = x0 - x;
- if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, pt, s)) continue;
if(fDebugLevel>=1){
(*fDebugStream) << "MC"
<< "det=" << det
<< "pdg=" << pdg
- << "pt=" << pt
- << "x0=" << x0
- << "y0=" << y0
- << "z0=" << z0
- << "dydx=" << dydx
- << "dzdx=" << dzdx
+ << "pt=" << pt0
+ << "x=" << x0
+ << "y=" << y0
+ << "z=" << z0
+ << "dydx=" << dydx0
+ << "dzdx=" << dzdx0
<< "\n";
}
- // MC track position at reference radial position
- Float_t yt = y0 - dx*dydx;
- Float_t zt = z0 - dx*dzdx;
- p = pt*(1.+dzdx*dzdx); // pt -> p
-
- // add Kalman residuals for y, z and pt
- Float_t yr = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
- dy = yt - yr;
- Float_t zr = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
- dz = zt - zr;
- Float_t tgl = fTracklet->GetTgl();
- Float_t dpt = pt - fTracklet->GetMomentum()/(1.+tgl*tgl);
+ Float_t ymc = y0 - dx*dydx0;
+ Float_t zmc = z0 - dx*dzdx0;
+ //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
+
+ // Kalman position at reference radial position
+ dx = xAnode - x;
+ dydx = fTracklet->GetYref(1);
+ dzdx = fTracklet->GetZref(1);
+ dzdl = fTracklet->GetTgl();
+ y = fTracklet->GetYref(0) - dx*dydx;
+ dy = y - ymc;
+ z = fTracklet->GetZref(0) - dx*dzdx;
+ dz = z - zmc;
+ pt = TMath::Abs(fTracklet->GetPt());
fTracklet->GetCovRef(covR);
- ((TH2I*)fContainer->At(kMCtrackY))->Fill(dydx, dy);
- ((TH2I*)fContainer->At(kMCtrackYPull))->Fill(dydx, dy/TMath::Sqrt(covR[0]));
- if(ily==0){
- ((TH2I*)fContainer->At(kMCtrackZIn))->Fill(dzdx, dz);
- ((TH2I*)fContainer->At(kMCtrackZInPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
- } else if(ily==AliTRDgeometry::kNlayer-1) {
- ((TH2I*)fContainer->At(kMCtrackZOut))->Fill(dzdx, dz);
- ((TH2I*)fContainer->At(kMCtrackZOutPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
+ arr = (TObjArray*)fContainer->At(kMCtrackTRD);
+ // y resolution/pulls
+ ((TH2I*)arr->At(0))->Fill(dydx0, dy);
+ ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(covR[0]));
+ // z resolution/pulls
+ ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
+ ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
+ // phi resolution/ snp pulls
+ Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
+ ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
+ Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
+ ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
+ // theta resolution/ tgl pulls
+ Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
+ dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
+ ((TH2I*)arr->At(6))->Fill(dzdl0,
+ TMath::ATan(dtgl));
+ ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
+ // pt resolution \\ 1/pt pulls \\ p resolution for PID
+ for(Int_t is=AliPID::kSPECIES; is--;){
+ if(TMath::Abs(pdg)!=AliPID::ParticleCode(is)) continue;
+ ((TH3S*)((TObjArray*)arr->At(8))->At(ily))->Fill(pt0, pt/pt0-1., is);
+ ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), is);
+ Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
+ p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
+ ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, p/p0-1., is);
+ break;
}
- if(pdg!=kElectron && pdg!=kPositron) ((TH2I*)fContainer->At(kMCtrackPt))->Fill(1./pt, dpt);
+
// Fill Debug stream for Kalman track
if(fDebugLevel>=1){
- Float_t dydxr = fTracklet->GetYref(1);
(*fDebugStream) << "MCtrack"
+ << "pt=" << pt
<< "x=" << x
- << "dpt=" << dpt
- << "dy=" << dy
- << "dz=" << dz
- << "dydxr=" << dydxr
- << "dzdxr=" << tgl
+ << "y=" << y
+ << "z=" << z
+ << "dydx=" << dydx
+ << "dzdx=" << dzdx
<< "s2y=" << covR[0]
<< "s2z=" << covR[2]
<< "\n";
// recalculate tracklet based on the MC info
AliTRDseedV1 tt(*fTracklet);
- tt.SetZref(0, z0);
- tt.SetZref(1, dzdx);
- tt.Fit(kTRUE);
- x= tt.GetX(); // the true one
+ tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
+ tt.SetZref(1, dzdx0);
+ tt.SetReconstructor(&rec);
+ tt.Fit(kTRUE, kTRUE);
+ x= tt.GetX();y= tt.GetY();z= tt.GetZ();
+ dydx = tt.GetYfit(1);
dx = x0 - x;
- yt = y0 - dx*dydx;
+ ymc = y0 - dx*dydx0;
+ zmc = z0 - dx*dzdx0;
Bool_t rc = tt.IsRowCross();
// add tracklet residuals for y and dydx
- Float_t yf = tt.GetY();
- dy = yt - yf; dz = 100.;
- Float_t dphi = (tt.GetYfit(1) - dydx);
- dphi /= 1.- tt.GetYfit(1)*dydx;
- Double_t s2y = tt.GetS2Y(), s2z = tt.GetS2Z();
+ arr = (TObjArray*)fContainer->At(kMCtracklet);
if(!rc){
- ((TH2I*)fContainer->At(kMCtrackletY))->Fill(dydx, dy);
- ((TH2I*)fContainer->At(kMCtrackletYPull))->Fill(dydx, dy/TMath::Sqrt(s2y));
- ((TH2I*)fContainer->At(kMCtrackletPhi))->Fill(dydx, dphi*TMath::RadToDeg());
+ dy = y-ymc;
+
+ Float_t dphi = (dydx - dydx0);
+ dphi /= (1.- dydx*dydx0);
+
+ ((TH2I*)arr->At(0))->Fill(dydx0, dy);
+ if(tt.GetS2Y()>0.) ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y()));
+ ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
} else {
// add tracklet residuals for z
- dz = tt.GetZ() - (z0 - dx*dzdx) ;
- ((TH2I*)fContainer->At(kMCtrackletZ))->Fill(dzdx, dz);
- ((TH2I*)fContainer->At(kMCtrackletZPull))->Fill(dzdx, dz/TMath::Sqrt(s2z));
+ dz = z-zmc;
+ ((TH2I*)arr->At(2))->Fill(dzdl0, dz);
+ if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()));
}
// Fill Debug stream for tracklet
if(fDebugLevel>=1){
+ Float_t s2y = tt.GetS2Y();
+ Float_t s2z = tt.GetS2Z();
(*fDebugStream) << "MCtracklet"
<< "rc=" << rc
<< "x=" << x
- << "dy=" << dy
- << "dz=" << dz
- << "dphi=" << dphi
+ << "y=" << y
+ << "z=" << z
+ << "dydx=" << dydx
<< "s2y=" << s2y
<< "s2z=" << s2z
<< "\n";
AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
Float_t tilt = fTracklet->GetTilt();
+ //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
+ arr = (TObjArray*)fContainer->At(kMCcluster);
AliTRDcluster *c = 0x0;
fTracklet->ResetClusterIter(kFALSE);
while((c = fTracklet->PrevCluster())){
Float_t q = TMath::Abs(c->GetQ());
- //AliTRDseedV1::GetClusterXY(c,x,y);
- x = c->GetX(); y = c->GetY(); z = c->GetZ();
- Float_t xc = x;
- Float_t yc = y;
- Float_t zc = z;
- dx = x0 - xc;
- yt = y0 - dx*dydx;
- zt = z0 - dx*dzdx;
- dy = yt - (yc - tilt*(zc-zt));
-
+ x = c->GetX(); y = c->GetY();z = c->GetZ();
+ dx = x0 - x;
+ ymc= y0 - dx*dydx0;
+ zmc= z0 - dx*dzdx0;
+ dy = ymc - (y - tilt*(z-zmc));
// Fill Histograms
- if(q>20. && q<250.) ((TH2I*)fContainer->At(kMCcluster))->Fill(dydx, dy);
+ if(q>20. && q<250.){
+ ((TH2I*)arr->At(0))->Fill(dydx0, dy);
+ ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(c->GetSigmaY2()));
+ }
// Fill calibration container
- Float_t d = zr0 - zt;
+ Float_t d = zr0 - zmc;
d -= ((Int_t)(2 * d)) / 2.0;
if (d > 0.25) d = 0.5 - d;
AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
fMCcl->Add(clInfo);
clInfo->SetCluster(c);
clInfo->SetMC(pdg, label);
- clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
+ clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
clInfo->SetResolution(dy);
clInfo->SetAnisochronity(d);
- clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
+ clInfo->SetDriftLength(((c->GetPadTime()+0.5)*.1)*1.5);
+ //dx-.5*AliTRDgeometry::CamHght());
clInfo->SetTilt(tilt);
// Fill Debug Tree
//________________________________________________________
Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
{
- Float_t y[2] = {0., 0.};
- TBox *b = 0x0;
- TAxis *ax = 0x0;
- TGraphErrors *g = 0x0;
+ Float_t xy[4] = {0., 0., 0., 0.};
+ if(!gPad){
+ AliWarning("Please provide a canvas to draw results.");
+ return kFALSE;
+ }
+ TList *l = 0x0; TVirtualPad *pad=0x0;
switch(ifig){
+ case kCharge:
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ ((TVirtualPad*)l->At(0))->cd();
+ ((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0))->Draw("apl");
+ ((TVirtualPad*)l->At(1))->cd();
+ ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0))->Draw("apl");
+ break;
case kCluster:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- g->Draw("apl");
- ax = g->GetHistogram()->GetYaxis();
- y[0] = -0.5; y[1] = 2.5;
- ax->SetRangeUser(y[0], y[1]);
- ax->SetTitle("Cluster-Track Pools #sigma/#mu [mm]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetTitle("tg(#phi)");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
- b = new TBox(-.15, y[0], .15, y[1]);
- b->SetFillStyle(3002);b->SetFillColor(kGreen);
- b->SetLineColor(0); b->Draw();
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0] = -.3; xy[1] = -200.; xy[2] = .3; xy[3] = 1000.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kCluster, 0)) break;
+ xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kCluster, 1)) break;
+ return kTRUE;
+ case kTrackTRD :
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kTrackTRD , 0)) break;
+ xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kTrackTRD , 1)) break;
return kTRUE;
- case kTrackletY:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- g->Draw("apl");
- ax = g->GetHistogram()->GetYaxis();
- ax->SetRangeUser(-.5, 3.);
- ax->SetTitle("Tracklet-Track Y-Pools #sigma/#mu [mm]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetTitle("tg(#phi)");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
- b = new TBox(-.15, y[0], .15, y[1]);
- b->SetFillStyle(3002);b->SetFillColor(kGreen);
- b->SetLineColor(0); b->Draw();
+ case 3: // kTrackTRD z
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kTrackTRD , 2)) break;
+ xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kTrackTRD , 3)) break;
return kTRUE;
- case kTrackletPhi:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- g->Draw("apl");
- ax = g->GetHistogram()->GetYaxis();
- ax->SetRangeUser(-.5, 2.);
- ax->SetTitle("Tracklet-Track Phi-Residuals #sigma/#mu [rad]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetTitle("tg(#phi)");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
- b = new TBox(-.15, y[0], .15, y[1]);
- b->SetFillStyle(3002);b->SetFillColor(kGreen);
- b->SetLineColor(0); b->Draw();
+ case 4: // kTrackTRD phi
+ xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
+ if(GetGraphPlot(&xy[0], kTrackTRD , 4)) return kTRUE;
+ break;
+ case 5: // kTrackTPC y
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
+ pad = ((TVirtualPad*)l->At(0)); pad->cd();
+ pad->SetMargin(0.1, 0.1, 0.1, 0.01);
+ if(!GetGraphPlot(&xy[0], kTrackTPC, 0)) break;
+ xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
+ pad=((TVirtualPad*)l->At(1)); pad->cd();
+ pad->SetMargin(0.1, 0.1, 0.1, 0.01);
+ if(!GetGraphPlot(&xy[0], kTrackTPC, 1)) break;
return kTRUE;
- case kMCcluster:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- ax = g->GetHistogram()->GetYaxis();
- y[0] = -.05; y[1] = 0.6;
- ax->SetRangeUser(y[0], y[1]);
- ax->SetTitle("Y_{cluster} #sigma/#mu [mm]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetRangeUser(-.3, .3);
- ax->SetTitle("tg(#phi)");
- g->Draw("apl");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
- b = new TBox(-.15, y[0], .15, y[1]);
- b->SetFillStyle(3002);b->SetFillColor(kBlue);
- b->SetLineColor(0); b->Draw();
+ case 6: // kTrackTPC z
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
+ pad = ((TVirtualPad*)l->At(0)); pad->cd();
+ pad->SetMargin(0.1, 0.1, 0.1, 0.01);
+ if(!GetGraphPlot(&xy[0], kTrackTPC, 2)) break;
+ xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
+ pad = ((TVirtualPad*)l->At(1)); pad->cd();
+ pad->SetMargin(0.1, 0.1, 0.1, 0.01);
+ if(!GetGraphPlot(&xy[0], kTrackTPC, 3)) break;
return kTRUE;
- case kMCtrackletY:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- ax = g->GetHistogram()->GetYaxis();
- y[0] = -.05; y[1] = 0.25;
- ax->SetRangeUser(y[0], y[1]);
- ax->SetTitle("Y_{tracklet} #sigma/#mu [mm]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetRangeUser(-.2, .2);
- ax->SetTitle("tg(#phi)");
- g->Draw("apl");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
- b = new TBox(-.15, y[0], .15, y[1]);
- b->SetFillStyle(3002);b->SetFillColor(kBlue);
- b->SetLineColor(0); b->Draw();
+ case 7: // kTrackTPC phi
+ xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
+ if(GetGraphPlot(&xy[0], kTrackTPC, 4)) return kTRUE;
+ break;
+ case 8: // kMCcluster
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCcluster, 0)) break;
+ xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCcluster, 1)) break;
return kTRUE;
- case kMCtrackletZ:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- ax = g->GetHistogram()->GetYaxis();
- ax->SetRangeUser(-.5, 1.);
- ax->SetTitle("Z_{tracklet} #sigma/#mu [mm]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetTitle("tg(#theta)");
- g->Draw("apl");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
+ case 9: //kMCtracklet [y]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtracklet, 0)) break;
+ xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtracklet, 1)) break;
return kTRUE;
- case kMCtrackletPhi:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- ax = g->GetHistogram()->GetYaxis();
- y[0] = -.05; y[1] = .2;
- ax->SetRangeUser(y[0], y[1]);
- ax->SetTitle("#Phi_{tracklet} #sigma/#mu [deg]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetTitle("tg(#phi)");
- g->Draw("apl");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
+ case 10: //kMCtracklet [z]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtracklet, 2)) break;
+ xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtracklet, 3)) break;
return kTRUE;
- case kMCtrackY:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- ax = g->GetHistogram()->GetYaxis();
- y[0] = -.05; y[1] = 0.25;
- ax->SetRangeUser(y[0], y[1]);
- ax->SetTitle("Y_{track} #sigma/#mu [mm]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetRangeUser(-.2, .2);
- ax->SetTitle("tg(#phi)");
- g->Draw("apl");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
- b = new TBox(-.15, y[0], .15, y[1]);
- b->SetFillStyle(3002);b->SetFillColor(kBlue);
- b->SetLineColor(0); b->Draw();
+ case 11: //kMCtracklet [phi]
+ xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
+ if(!GetGraphPlot(&xy[0], kMCtracklet, 4)) break;
return kTRUE;
- case kMCtrackZIn:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- ax = g->GetHistogram()->GetYaxis();
- ax->SetRangeUser(-.5, 2.);
- ax->SetTitle("Z_{track} #sigma/#mu [mm]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetTitle("tg(#theta)");
- g->Draw("apl");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
+ case 12: //kMCtrackTRD [y]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =200.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 0)) break;
+ xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 1)) break;
return kTRUE;
- case kMCtrackPt:
- if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
- ax = g->GetHistogram()->GetYaxis();
- ax->SetRangeUser(-.5, 2.);
- ax->SetTitle("#epsilon_{P_{t}}^{track} / #mu [%]");
- ax = g->GetHistogram()->GetXaxis();
- ax->SetTitle("1/p_{t}");
- g->Draw("apl");
- if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
- g->Draw("pl");
+ case 13: //kMCtrackTRD [z]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-1.; xy[1]=-700.; xy[2]=1.; xy[3] =1500.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 2)) break;
+ xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 3)) break;
+ return kTRUE;
+ case 14: //kMCtrackTRD [phi/snp]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-.2; xy[1]=-0.2; xy[2]=.2; xy[3] =2.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 4)) break;
+ xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 5)) break;
+ return kTRUE;
+ case 15: //kMCtrackTRD [theta/tgl]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 6)) break;
+ xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTRD, 7)) break;
+ return kTRUE;
+ case 16: //kMCtrackTRD [pt]
+ xy[0] = 0.; xy[1] = -5.; xy[2] = 12.; xy[3] = 7.;
+ gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
+ pad = (TVirtualPad*)l->At(il); pad->cd();
+ pad->SetMargin(0.07, 0.07, 0.1, 0.);
+ if(!GetGraphTrack(&xy[0], 8, il)) break;
+ }
+ return kTRUE;
+ case 17: //kMCtrackTRD [1/pt] pulls
+ xy[0] = 0.; xy[1] = -1.5; xy[2] = 2.; xy[3] = 2.;
+ gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
+ pad = (TVirtualPad*)l->At(il); pad->cd();
+ pad->SetMargin(0.07, 0.07, 0.1, 0.);
+ if(!GetGraphTrack(&xy[0], 9, il)) break;
+ }
+ return kTRUE;
+ case 18: //kMCtrackTRD [p]
+ xy[0] = 0.; xy[1] = -7.5; xy[2] = 12.; xy[3] = 10.5;
+ gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
+ pad = (TVirtualPad*)l->At(il); pad->cd();
+ pad->SetMargin(0.07, 0.07, 0.1, 0.);
+ if(!GetGraphTrack(&xy[0], 10, il)) break;
+ }
+ return kTRUE;
+ case 19: // kMCtrackTPC [y]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-.25; xy[1]=-50.; xy[2]=.25; xy[3] =800.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 0)) break;
+ xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break;
+ return kTRUE;
+ case 20: // kMCtrackTPC [z]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 2)) break;
+ xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break;
+ return kTRUE;
+ case 21: // kMCtrackTPC [phi|snp]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 4)) break;
+ xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 5)) break;
+ return kTRUE;
+ case 22: // kMCtrackTPC [theta|tgl]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 6)) break;
+ xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphPlot(&xy[0], kMCtrackTPC, 7)) break;
+ return kTRUE;
+ case 23: // kMCtrackTPC [pt]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3;
+ ((TVirtualPad*)l->At(0))->cd();
+ if(!GetGraphTrackTPC(xy, 8)) break;
+ xy[0]=0.; xy[1]=-0.5; xy[2]=2.; xy[3] =2.5;
+ ((TVirtualPad*)l->At(1))->cd();
+ if(!GetGraphTrackTPC(xy, 9)) break;
+ return kTRUE;
+ case 24: // kMCtrackTPC [p]
+ gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
+ xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3;
+ pad = ((TVirtualPad*)l->At(0));pad->cd();
+ pad->SetMargin(0.12, 0.12, 0.1, 0.04);
+ if(!GetGraphTrackTPC(xy, 10)) break;
+ xy[0]=0.; xy[1]=-1.5; xy[2]=12.; xy[3] =2.5;
+ pad = ((TVirtualPad*)l->At(1)); pad->cd();
+ pad->SetMargin(0.12, 0.12, 0.1, 0.04);
+ if(!GetGraphTrackTPC(xy, 11)) break;
+ return kTRUE;
+ case 25: // kMCtrackTOF [z]
return kTRUE;
}
- AliInfo(Form("Reference plot [%d] missing result", ifig));
+ AliWarning(Form("Reference plot [%d] missing result", ifig));
return kFALSE;
}
Printf("ERROR: list not available");
return kFALSE;
}
- fNRefFigures = fContainer->GetEntriesFast();
- TGraphErrors *gm = 0x0, *gs = 0x0;
- if(!fGraphS){
- fGraphS = new TObjArray(fNRefFigures);
- fGraphS->SetOwner();
- for(Int_t ig=0; ig<fNRefFigures; ig++){
- gs = new TGraphErrors();
- gs->SetLineColor(kRed);
- gs->SetMarkerStyle(23);
- gs->SetMarkerColor(kRed);
- gs->SetNameTitle(Form("s_%d", ig), "");
- fGraphS->AddAt(gs, ig);
- }
- }
- if(!fGraphM){
- fGraphM = new TObjArray(fNRefFigures);
- fGraphM->SetOwner();
- for(Int_t ig=0; ig<fNRefFigures; ig++){
- gm = new TGraphErrors();
- gm->SetLineColor(kBlue);
- gm->SetMarkerStyle(7);
- gm->SetMarkerColor(kBlue);
- gm->SetNameTitle(Form("m_%d", ig), "");
- fGraphM->AddAt(gm, ig);
+ TGraph *gm= 0x0, *gs= 0x0;
+ if(!fGraphS && !fGraphM){
+ TObjArray *aM(0x0), *aS(0x0), *a(0x0);
+ Int_t n = fContainer->GetEntriesFast();
+ fGraphS = new TObjArray(n); fGraphS->SetOwner();
+ fGraphM = new TObjArray(n); fGraphM->SetOwner();
+ for(Int_t ig=0; ig<n; ig++){
+ fGraphM->AddAt(aM = new TObjArray(fNElements[ig]), ig);
+ fGraphS->AddAt(aS = new TObjArray(fNElements[ig]), ig);
+
+ for(Int_t ic=0; ic<fNElements[ig]; ic++){
+ if(ig==kMCtrackTPC&&(ic>=8&&ic<=12)){ // TPC momentum plot
+ aS->AddAt(a = new TObjArray(AliPID::kSPECIES), ic);
+ for(Int_t is=AliPID::kSPECIES; is--;){
+ a->AddAt(gs = new TGraphErrors(), is);
+ gs->SetMarkerStyle(23);
+ gs->SetMarkerColor(is ? kRed : kMagenta);
+ gs->SetLineStyle(is);
+ gs->SetLineColor(is ? kRed : kMagenta);
+ gs->SetLineWidth(is ? 1 : 3);
+ gs->SetNameTitle(Form("s_%d%02d%d", ig, ic, is), "");
+ }
+ aM->AddAt(a = new TObjArray(AliPID::kSPECIES), ic);
+ for(Int_t is=AliPID::kSPECIES; is--;){
+ a->AddAt(gm = new TGraphErrors(), is);
+ gm->SetLineColor(is ? kBlack : kBlue);
+ gm->SetLineStyle(is);
+ gm->SetMarkerStyle(7);
+ gm->SetMarkerColor(is ? kBlack : kBlue);
+ gm->SetLineWidth(is ? 1 : 3);
+ gm->SetNameTitle(Form("m_%d%02d%d", ig, ic, is), "");
+ }
+ continue;
+ } else if(ig==kMCtrackTRD&&(ic==8||ic==9||ic==10)){ // TRD momentum plot
+ TObjArray *aaS, *aaM;
+ aS->AddAt(aaS = new TObjArray(AliTRDgeometry::kNlayer), ic);
+ aM->AddAt(aaM = new TObjArray(AliTRDgeometry::kNlayer), ic);
+ for(Int_t il=AliTRDgeometry::kNlayer; il--;){
+ aaS->AddAt(a = new TObjArray(AliPID::kSPECIES), il);
+ for(Int_t is=AliPID::kSPECIES; is--;){
+ a->AddAt(gs = new TGraphErrors(), is);
+ gs->SetMarkerStyle(23);
+ gs->SetMarkerColor(is ? kRed : kMagenta);
+ gs->SetLineStyle(is);
+ gs->SetLineColor(is ? kRed : kMagenta);
+ gs->SetLineWidth(is ? 1 : 3);
+ gs->SetNameTitle(Form("s_%d%02d%d%d", ig, ic, is, il), "");
+ }
+ aaM->AddAt(a = new TObjArray(AliPID::kSPECIES), il);
+ for(Int_t is=AliPID::kSPECIES; is--;){
+ a->AddAt(gm = new TGraphErrors(), is);
+ gm->SetMarkerStyle(7);
+ gm->SetMarkerColor(is ? kBlack : kBlue);
+ gm->SetLineStyle(is);
+ gm->SetLineColor(is ? kBlack : kBlue);
+ gm->SetLineWidth(is ? 1 : 3);
+ gm->SetNameTitle(Form("m_%d%02d%d%d", ig, ic, is, il), "");
+ }
+ }
+ continue;
+ }
+
+ aS->AddAt(gs = new TGraphErrors(), ic);
+ gs->SetMarkerStyle(23);
+ gs->SetMarkerColor(kRed);
+ gs->SetLineColor(kRed);
+ gs->SetNameTitle(Form("s_%d%02d", ig, ic), "");
+
+ aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
+ gm->SetLineColor(kBlack);
+ gm->SetMarkerStyle(7);
+ gm->SetMarkerColor(kBlack);
+ gm->SetNameTitle(Form("m_%d%02d", ig, ic), "");
+ }
}
}
- TH2I *h2 = 0x0;
- TH1D *h = 0x0;
-
- // define models
- TF1 f("f1", "gaus", -.5, .5);
-
- TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
-
- TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
-
- TCanvas *c = 0x0;
- if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
- char opt[5];
- sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
-
+/* printf("\n\n\n"); fGraphS->ls();
+ printf("\n\n\n"); fGraphM->ls();*/
+
- //PROCESS RESIDUAL DISTRIBUTIONS
+ // DEFINE MODELS
+ // simple gauss
+ TF1 fg("fGauss", "gaus", -.5, .5);
+ // Landau for charge resolution
+ TF1 fl("fLandau", "landau", 0., 1000.);
+ //PROCESS EXPERIMENTAL DISTRIBUTIONS
+ // Charge resolution
+ //Process3DL(kCharge, 0, &fl);
// Clusters residuals
- h2 = (TH2I *)(fContainer->At(kCluster));
- gm = (TGraphErrors*)fGraphM->At(kCluster);
- gs = (TGraphErrors*)fGraphS->At(kCluster);
- 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);}
-
- 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 residuals
- h2 = (TH2I *)(fContainer->At(kTrackletY));
- gm = (TGraphErrors*)fGraphM->At(kTrackletY);
- gs = (TGraphErrors*)fGraphS->At(kTrackletY);
- 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);}
-
- 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 phi residuals
- h2 = (TH2I *)(fContainer->At(kTrackletPhi));
- gm = (TGraphErrors*)fGraphM->At(kTrackletPhi);
- gs = (TGraphErrors*)fGraphS->At(kTrackletPhi);
- 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);}
-
- 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));
- }
+ Process2D(kCluster, 0, &fg, 1.e4);
+ Process2D(kCluster, 1, &fg);
+ fNRefFigures = 1;
+ // Tracklet residual/pulls
+ Process2D(kTrackTRD , 0, &fg, 1.e4);
+ Process2D(kTrackTRD , 1, &fg);
+ Process2D(kTrackTRD , 2, &fg, 1.e4);
+ Process2D(kTrackTRD , 3, &fg);
+ Process2D(kTrackTRD , 4, &fg, 1.e3);
+ fNRefFigures = 4;
+ // TPC track residual/pulls
+ Process2D(kTrackTPC, 0, &fg, 1.e4);
+ Process2D(kTrackTPC, 1, &fg);
+ Process2D(kTrackTPC, 2, &fg, 1.e4);
+ Process2D(kTrackTPC, 3, &fg);
+ Process2D(kTrackTPC, 4, &fg, 1.e3);
+ fNRefFigures = 7;
+
+ if(!HasMCdata()) return kTRUE;
- if(!HasMCdata()){
- if(c) delete c;
- return kTRUE;
- }
//PROCESS MC RESIDUAL DISTRIBUTIONS
- // cluster y resolution
- h2 = (TH2I*)fContainer->At(kMCcluster);
- gm = (TGraphErrors*)fGraphM->At(kMCcluster);
- gs = (TGraphErrors*)fGraphS->At(kMCcluster);
- for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
- h = h2->ProjectionY("py", iphi, iphi);
- if(h->GetEntries()<100) continue;
- AdjustF1(h, &f);
-
- if(IsVisual()){c->cd(); c->SetLogy();}
- h->Fit(&f, opt, "", -0.5, 0.5);
- if(IsVerbose()){
- printf("phi[%d] mean[%e] sigma[%e]\n\n", iphi, 10.*f.GetParameter(1), 10.*f.GetParameter(2));
- }
- if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
-
- Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
- 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(kMCtrackletY);
- gm = (TGraphErrors*)fGraphM->At(kMCtrackletY);
- gs = (TGraphErrors*)fGraphS->At(kMCtrackletY);
- for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
- h = h2->ProjectionY("py", iphi, iphi);
- 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);}
-
- Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
- 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 z resolution
- h2 = (TH2I*)fContainer->At(kMCtrackletZ);
- gm = (TGraphErrors*)fGraphM->At(kMCtrackletZ);
- gs = (TGraphErrors*)fGraphS->At(kMCtrackletZ);
- 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();}
- h->Fit(&fb, opt, "", -0.5, 0.5);
- if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
-
- Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
- 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(kMCtrackletPhi);
- gm = (TGraphErrors*)fGraphM->At(kMCtrackletPhi);
- gs = (TGraphErrors*)fGraphS->At(kMCtrackletPhi);
- 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 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));
- }
-
- // track y resolution
- h2 = (TH2I*)fContainer->At(kMCtrackY);
- gm = (TGraphErrors*)fGraphM->At(kMCtrackY);
- gs = (TGraphErrors*)fGraphS->At(kMCtrackY);
- for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
- h = h2->ProjectionY("py", iphi, iphi);
- 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);}
-
- Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
- 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));
- }
-
- // track z resolution
- h2 = (TH2I*)fContainer->At(kMCtrackZIn);
- gm = (TGraphErrors*)fGraphM->At(kMCtrackZIn);
- gs = (TGraphErrors*)fGraphS->At(kMCtrackZIn);
- for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
- h = h2->ProjectionY("pz", iphi, iphi);
- if(h->GetEntries()<70) 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);}
-
- Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
- 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));
- }
-
- // track Pt resolution
- h2 = (TH2I*)fContainer->At(kMCtrackPt);
- TAxis *ax = h2->GetXaxis();
- gm = (TGraphErrors*)fGraphM->At(kMCtrackPt);
- gs = (TGraphErrors*)fGraphS->At(kMCtrackPt);
- TF1 fg("fg", "gaus", -1.5, 1.5);
- TF1 fl("fl", "landau", -4., 15.);
- TF1 fgl("fgl", "gaus(0)+landau(3)", -5., 20.);
- for(Int_t ip=1; ip<=ax->GetNbins(); ip++){
- h = h2->ProjectionY("ppt", ip, ip);
- if(h->GetEntries()<70) continue;
-
- h->Fit(&fg, "QN", "", -1.5, 1.5);
- fgl.SetParameter(0, fg.GetParameter(0));
- fgl.SetParameter(1, fg.GetParameter(1));
- fgl.SetParameter(2, fg.GetParameter(2));
- h->Fit(&fl, "QN", "", -4., 15.);
- fgl.SetParameter(3, fl.GetParameter(0));
- fgl.SetParameter(4, fl.GetParameter(1));
- fgl.SetParameter(5, fl.GetParameter(2));
-
- if(IsVisual()){c->cd(); c->SetLogy();}
- h->Fit(&fgl, opt, "", -5., 20.);
- if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
-
- Float_t invpt = ax->GetBinCenter(ip);
- Int_t ip = gm->GetN();
- gm->SetPoint(ip, invpt, fgl.GetParameter(1));
- gm->SetPointError(ip, 0., fgl.GetParError(1));
- gs->SetPoint(ip, invpt, fgl.GetParameter(2)*invpt);
- gs->SetPointError(ip, 0., fgl.GetParError(2));
- // fgl.GetParameter(4) // Landau MPV
- // fgl.GetParameter(5) // Landau Sigma
- }
+ // CLUSTER Y RESOLUTION/PULLS
+ Process2D(kMCcluster, 0, &fg, 1.e4);
+ Process2D(kMCcluster, 1, &fg);
+ fNRefFigures = 8;
+
+ // TRACKLET RESOLUTION/PULLS
+ Process2D(kMCtracklet, 0, &fg, 1.e4); // y
+ Process2D(kMCtracklet, 1, &fg); // y pulls
+ Process2D(kMCtracklet, 2, &fg, 1.e4); // z
+ Process2D(kMCtracklet, 3, &fg); // z pulls
+ Process2D(kMCtracklet, 4, &fg, 1.e3); // phi
+ fNRefFigures = 11;
+
+ // TRACK RESOLUTION/PULLS
+ Process2D(kMCtrackTRD, 0, &fg, 1.e4); // y
+ Process2D(kMCtrackTRD, 1, &fg); // y PULL
+ Process2D(kMCtrackTRD, 2, &fg, 1.e4); // z
+ Process2D(kMCtrackTRD, 3, &fg); // z PULL
+ Process2D(kMCtrackTRD, 4, &fg, 1.e3); // phi
+ Process2D(kMCtrackTRD, 5, &fg); // snp PULL
+ Process2D(kMCtrackTRD, 6, &fg, 1.e3); // theta
+ Process2D(kMCtrackTRD, 7, &fg); // tgl PULL
+ Process4D(kMCtrackTRD, 8, &fg, 1.e2); // pt resolution
+ Process4D(kMCtrackTRD, 9, &fg); // 1/pt pulls
+ Process4D(kMCtrackTRD, 10, &fg, 1.e2); // p resolution
+ fNRefFigures = 18;
+
+ // TRACK TPC RESOLUTION/PULLS
+ Process2D(kMCtrackTPC, 0, &fg, 1.e4);// y resolution
+ Process2D(kMCtrackTPC, 1, &fg); // y pulls
+ Process2D(kMCtrackTPC, 2, &fg, 1.e4);// z resolution
+ Process2D(kMCtrackTPC, 3, &fg); // z pulls
+ Process2D(kMCtrackTPC, 4, &fg, 1.e3);// phi resolution
+ Process2D(kMCtrackTPC, 5, &fg); // snp pulls
+ Process2D(kMCtrackTPC, 6, &fg, 1.e3);// theta resolution
+ Process2D(kMCtrackTPC, 7, &fg); // tgl pulls
+ Process3D(kMCtrackTPC, 8, &fg, 1.e2);// pt resolution
+ Process3D(kMCtrackTPC, 9, &fg); // 1/pt pulls
+ Process3D(kMCtrackTPC, 10, &fg, 1.e2);// p resolution
+ Process3D(kMCtrackTPC, 11, &fg); // p pulls
+ fNRefFigures = 24;
+
+ // TRACK HMPID RESOLUTION/PULLS
+ Process2D(kMCtrackTOF, 0, &fg, 1.e4); // z towards TOF
+ Process2D(kMCtrackTOF, 1, &fg); // z towards TOF
+ fNRefFigures = 25;
-
- if(c) delete c;
return kTRUE;
}
for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
f->SetParLimits(0, 0., 3.*sum);
f->SetParameter(0, .9*sum);
+ Double_t rms = h->GetRMS();
+ f->SetParLimits(1, -rms, rms);
+ f->SetParameter(1, h->GetMean());
- f->SetParLimits(1, -.2, .2);
- f->SetParameter(1, -0.1);
-
- f->SetParLimits(2, 0., 4.e-1);
- f->SetParameter(2, 2.e-2);
+ f->SetParLimits(2, 0., 2.*rms);
+ f->SetParameter(2, rms);
if(f->GetNpar() <= 4) return;
f->SetParLimits(3, 0., sum);
fContainer = new TObjArray(kNhistos);
//fContainer->SetOwner(kTRUE);
-
TH1 *h = 0x0;
- // cluster to tracklet residuals [2]
+ TObjArray *arr = 0x0;
+
+ // cluster to track residuals/pulls
+ fContainer->AddAt(arr = new TObjArray(fNElements[kCharge]), kCharge);
+ arr->SetName("Charge");
+ if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
+ h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
+ h->GetXaxis()->SetTitle("dx/dx_{0}");
+ h->GetYaxis()->SetTitle("x_{d} [cm]");
+ h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
+ } else h->Reset();
+ arr->AddAt(h, 0);
+
+ // cluster to track residuals/pulls
+ fContainer->AddAt(arr = new TObjArray(fNElements[kCluster]), kCluster);
+ arr->SetName("Cl");
if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
- h = new TH2I("hCl", "Clusters-Track Residuals", 21, -.33, .33, 100, -.5, .5);
+ h = new TH2I("hCl", "Cluster Residuals", 21, -.33, .33, 100, -.5, .5);
h->GetXaxis()->SetTitle("tg(#phi)");
h->GetYaxis()->SetTitle("#Delta y [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kCluster);
+ arr->AddAt(h, 0);
+ if(!(h = (TH2I*)gROOT->FindObject("hClpull"))){
+ h = new TH2I("hClpull", "Cluster Pulls", 21, -.33, .33, 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 1);
- // tracklet to track residuals [2]
+ // tracklet to track residuals/pulls in y direction
+ fContainer->AddAt(arr = new TObjArray(fNElements[kTrackTRD ]), kTrackTRD );
+ arr->SetName("Trklt");
if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
- h = new TH2I("hTrkltY", "Tracklets-Track Residuals (Y)", 21, -.33, .33, 100, -.5, .5);
- h->GetXaxis()->SetTitle("tg(#phi)");
+ h = new TH2I("hTrkltY", "Tracklet Y Residuals", 21, -.33, .33, 100, -.5, .5);
+ h->GetXaxis()->SetTitle("#tg(#phi)");
h->GetYaxis()->SetTitle("#Delta y [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kTrackletY);
-
- // tracklet to track residuals angular [2]
+ arr->AddAt(h, 0);
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkltYpull"))){
+ h = new TH2I("hTrkltYpull", "Tracklet Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("#tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 1);
+ // tracklet to track residuals/pulls in z direction
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkltZ"))){
+ h = new TH2I("hTrkltZ", "Tracklet Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
+ h->GetXaxis()->SetTitle("#tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta z [cm]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 2);
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkltZpull"))){
+ h = new TH2I("hTrkltZpull", "Tracklet Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
+ h->GetXaxis()->SetTitle("#tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 3);
+ // tracklet to track phi residuals
if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
- h = new TH2I("hTrkltPhi", "Tracklets-Track Residuals (#Phi)", 21, -.33, .33, 100, -.5, .5);
+ h = new TH2I("hTrkltPhi", "Tracklet #phi Residuals", 21, -.33, .33, 100, -.5, .5);
h->GetXaxis()->SetTitle("tg(#phi)");
- h->GetYaxis()->SetTitle("#Delta phi [#circ]");
+ h->GetYaxis()->SetTitle("#Delta phi [rad]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kTrackletPhi);
+ arr->AddAt(h, 4);
+
+
+ // tracklet to TPC track residuals/pulls in y direction
+ fContainer->AddAt(arr = new TObjArray(fNElements[kTrackTPC]), kTrackTPC);
+ arr->SetName("TrkTPC");
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCY"))){
+ h = new TH2I("hTrkTPCY", "Track[TPC] Y Residuals", 21, -.33, .33, 100, -.5, .5);
+ h->GetXaxis()->SetTitle("#tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta y [cm]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 0);
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCYpull"))){
+ h = new TH2I("hTrkTPCYpull", "Track[TPC] Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("#tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 1);
+ // tracklet to TPC track residuals/pulls in z direction
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZ"))){
+ h = new TH2I("hTrkTPCZ", "Track[TPC] Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
+ h->GetXaxis()->SetTitle("#tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta z [cm]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 2);
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZpull"))){
+ h = new TH2I("hTrkTPCZpull", "Track[TPC] Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
+ h->GetXaxis()->SetTitle("#tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 3);
+ // tracklet to TPC track phi residuals
+ if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCPhi"))){
+ h = new TH2I("hTrkTPCPhi", "Track[TPC] #phi Residuals", 21, -.33, .33, 100, -.5, .5);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta phi [rad]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 4);
// Resolution histos
if(!HasMCdata()) return fContainer;
// cluster y resolution [0]
- if(!(h = (TH2I*)gROOT->FindObject("hMCcl"))){
- h = new TH2I("hMCcl", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5);
+ fContainer->AddAt(arr = new TObjArray(fNElements[kMCcluster]), kMCcluster);
+ arr->SetName("McCl");
+ if(!(h = (TH2I*)gROOT->FindObject("hMcCl"))){
+ h = new TH2I("hMcCl", "Cluster Resolution", 48, -.48, .48, 100, -.3, .3);
h->GetXaxis()->SetTitle("tg(#phi)");
h->GetYaxis()->SetTitle("#Delta y [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCcluster);
+ arr->AddAt(h, 0);
+ if(!(h = (TH2I*)gROOT->FindObject("hMcClPull"))){
+ h = new TH2I("hMcClPull", "Cluster Pulls", 48, -.48, .48, 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Deltay/#sigma_{y}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 1);
+
- // tracklet y resolution [0]
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltY"))){
- h = new TH2I("hMCtrkltY", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.5, .5);
+ // TRACKLET RESOLUTION
+ fContainer->AddAt(arr = new TObjArray(fNElements[kMCtracklet]), kMCtracklet);
+ arr->SetName("McTrklt");
+ // tracklet y resolution
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltY"))){
+ h = new TH2I("hMcTrkltY", "Tracklet Resolution (Y)", 48, -.48, .48, 100, -.2, .2);
h->GetXaxis()->SetTitle("tg(#phi)");
h->GetYaxis()->SetTitle("#Delta y [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackletY);
-
- // tracklet y resolution [0]
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltYPull"))){
- h = new TH2I("hMCtrkltYPull", "Tracklet Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
+ arr->AddAt(h, 0);
+ // tracklet y pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltYPull"))){
+ h = new TH2I("hMcTrkltYPull", "Tracklet Pulls (Y)", 48, -.48, .48, 100, -4.5, 4.5);
h->GetXaxis()->SetTitle("tg(#phi)");
h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackletYPull);
-
- // tracklet y resolution [0]
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZ"))){
- h = new TH2I("hMCtrkltZ", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5);
+ arr->AddAt(h, 1);
+ // tracklet z resolution
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZ"))){
+ h = new TH2I("hMcTrkltZ", "Tracklet Resolution (Z)", 100, -1., 1., 100, -1., 1.);
h->GetXaxis()->SetTitle("tg(#theta)");
h->GetYaxis()->SetTitle("#Delta z [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackletZ);
-
- // tracklet y resolution [0]
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZPull"))){
- h = new TH2I("hMCtrkltZ", "Tracklet Pulls (Z)", 31, -.48, .48, 100, -3.5, 3.5);
+ arr->AddAt(h, 2);
+ // tracklet z pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZPull"))){
+ h = new TH2I("hMcTrkltZPull", "Tracklet Pulls (Z)", 100, -1., 1., 100, -3.5, 3.5);
h->GetXaxis()->SetTitle("tg(#theta)");
h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackletZPull);
-
- // tracklet angular resolution [1]
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltPhi"))){
- h = new TH2I("hMCtrkltPhi", "Tracklet Resolution (#Phi)", 31, -.48, .48, 100, -10., 10.);
+ arr->AddAt(h, 3);
+ // tracklet phi resolution
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltPhi"))){
+ h = new TH2I("hMcTrkltPhi", "Tracklet Resolution (#Phi)", 48, -.48, .48, 100, -.15, .15);
h->GetXaxis()->SetTitle("tg(#phi)");
- h->GetYaxis()->SetTitle("#Delta #phi [deg]");
+ h->GetYaxis()->SetTitle("#Delta #phi [rad]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackletPhi);
+ arr->AddAt(h, 4);
+
+ // KALMAN TRACK RESOLUTION
+ fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackTRD]), kMCtrackTRD);
+ arr->SetName("McTrkTRD");
// Kalman track y resolution
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkY"))){
- h = new TH2I("hMCtrkY", "Kalman Track Resolution (Y)", 31, -.48, .48, 100, -.5, .5);
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkY"))){
+ h = new TH2I("hMcTrkY", "Track Y Resolution", 48, -.48, .48, 100, -.2, .2);
h->GetXaxis()->SetTitle("tg(#phi)");
h->GetYaxis()->SetTitle("#Delta y [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackY);
-
- // Kalman track y resolution
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkYPull"))){
- h = new TH2I("hMCtrkYPull", "Kalman Track Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
+ arr->AddAt(h, 0);
+ // Kalman track y pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkYPull"))){
+ h = new TH2I("hMcTrkYPull", "Track Y Pulls", 48, -.48, .48, 100, -4., 4.);
h->GetXaxis()->SetTitle("tg(#phi)");
h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackYPull);
-
- // Kalman track Z resolution
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZIn"))){
- h = new TH2I("hMCtrkZIn", "Kalman Track Resolution (Zin)", 20, -1., 1., 100, -1.5, 1.5);
+ arr->AddAt(h, 1);
+ // Kalman track Z
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZ"))){
+ h = new TH2I("hMcTrkZ", "Track Z Resolution", 100, -1., 1., 100, -1., 1.);
h->GetXaxis()->SetTitle("tg(#theta)");
h->GetYaxis()->SetTitle("#Delta z [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackZIn);
+ arr->AddAt(h, 2);
+ // Kalman track Z pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZPull"))){
+ h = new TH2I("hMcTrkZPull", "Track Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 3);
+ // Kalman track SNP
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNP"))){
+ h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta #phi [rad]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 4);
+ // Kalman track SNP pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNPPull"))){
+ h = new TH2I("hMcTrkSNPPull", "Track SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 5);
+ // Kalman track TGL
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGL"))){
+ h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
+ h->GetXaxis()->SetTitle("tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta#theta [rad]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 6);
+ // Kalman track TGL pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGLPull"))){
+ h = new TH2I("hMcTrkTGLPull", "Track TGL Pulls", 100, -1., 1., 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 7);
+ // Kalman track Pt resolution
+ const Int_t n = AliPID::kSPECIES;
+ TObjArray *arr2 = 0x0; TH3S* h3=0x0;
+ arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 8);
+ arr2->SetName("Track Pt Resolution");
+ for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
+ if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPt%d", il)))){
+ h3 = new TH3S(Form("hMcTrkPt%d", il), "Track Pt Resolution", 40, 0., 20., 150, -.1, .2, n, -.5, n-.5);
+ h3->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
+ h3->GetZaxis()->SetTitle("SPECIES");
+ } else h3->Reset();
+ arr2->AddAt(h3, il);
+ }
+ // Kalman track Pt pulls
+ arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 9);
+ arr2->SetName("Track 1/Pt Pulls");
+ for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
+ if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPtPulls%d", il)))){
+ h3 = new TH3S(Form("hMcTrkPtPulls%d", il), "Track 1/Pt Pulls", 40, 0., 2., 100, -4., 4., n, -.5, n-.5);
+ h3->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
+ h3->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
+ h3->GetZaxis()->SetTitle("SPECIES");
+ } else h3->Reset();
+ arr2->AddAt(h3, il);
+ }
+ // Kalman track P resolution
+ arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 10);
+ arr2->SetName("Track P Resolution [PID]");
+ for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
+ if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkP%d", il)))){
+ h3 = new TH3S(Form("hMcTrkP%d", il), "Track P Resolution", 40, 0., 20., 150, -.15, .35, n, -.5, n-.5);
+ h3->GetXaxis()->SetTitle("p [GeV/c]");
+ h3->GetYaxis()->SetTitle("#Delta p/p^{MC}");
+ h3->GetZaxis()->SetTitle("SPECIES");
+ } else h3->Reset();
+ arr2->AddAt(h3, il);
+ }
- // Kalman track Z resolution
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOut"))){
- h = new TH2I("hMCtrkZOut", "Kalman Track Resolution (Zout)", 20, -1., 1., 100, -1.5, 1.5);
+ // TPC TRACK RESOLUTION
+ fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackTPC]), kMCtrackTPC);
+ arr->SetName("McTrkTPC");
+ // Kalman track Y
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCY"))){
+ h = new TH2I("hMcTrkTPCY", "Track[TPC] Y Resolution", 60, -.3, .3, 100, -.5, .5);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta y [cm]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 0);
+ // Kalman track Y pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCYPull"))){
+ h = new TH2I("hMcTrkTPCYPull", "Track[TPC] Y Pulls", 60, -.3, .3, 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 1);
+ // Kalman track Z
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZ"))){
+ h = new TH2I("hMcTrkTPCZ", "Track[TPC] Z Resolution", 100, -1., 1., 100, -1., 1.);
h->GetXaxis()->SetTitle("tg(#theta)");
h->GetYaxis()->SetTitle("#Delta z [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackZOut);
-
- // Kalman track Z resolution
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZInPull"))){
- h = new TH2I("hMCtrkZInPull", "Kalman Track Pulls (Zin)", 20, -1., 1., 100, -4.5, 4.5);
+ arr->AddAt(h, 2);
+ // Kalman track Z pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZPull"))){
+ h = new TH2I("hMcTrkTPCZPull", "Track[TPC] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
h->GetXaxis()->SetTitle("tg(#theta)");
h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackZInPull);
-
- // Kalman track Z resolution
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOutPull"))){
- h = new TH2I("hMCtrkZOutPull", "Kalman Track Pulls (Zout)", 20, -1., 1., 100, -4.5, 4.5);
+ arr->AddAt(h, 3);
+ // Kalman track SNP
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNP"))){
+ h = new TH2I("hMcTrkTPCSNP", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta #phi [rad]");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 4);
+ // Kalman track SNP pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNPPull"))){
+ h = new TH2I("hMcTrkTPCSNPPull", "Track[TPC] SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#phi)");
+ h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 5);
+ // Kalman track TGL
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGL"))){
+ h = new TH2I("hMcTrkTPCTGL", "Track[TPC] Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
h->GetXaxis()->SetTitle("tg(#theta)");
- h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
+ h->GetYaxis()->SetTitle("#Delta#theta [rad]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackZOutPull);
-
+ arr->AddAt(h, 6);
+ // Kalman track TGL pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGLPull"))){
+ h = new TH2I("hMcTrkTPCTGLPull", "Track[TPC] TGL Pulls", 100, -1., 1., 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 7);
// Kalman track Pt resolution
- if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPt"))){
- h = new TH2I("hMCtrkPt", "Kalman Track Resolution (Pt)", 100, 0., 2., 150, -5., 20.);
- h->GetXaxis()->SetTitle("1/p_{t}");
- h->GetYaxis()->SetTitle("#Delta p_{t} [GeV/c]");
+ if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPt"))){
+ h3 = new TH3S("hMcTrkTPCPt", "Track[TPC] Pt Resolution", 80, 0., 20., 150, -.1, .2, n, -.5, n-.5);
+ h3->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
+ h3->GetZaxis()->SetTitle("SPECIES");
+ } else h3->Reset();
+ arr->AddAt(h3, 8);
+ // Kalman track Pt pulls
+ if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPtPulls"))){
+ h3 = new TH3S("hMcTrkTPCPtPulls", "Track[TPC] 1/Pt Pulls", 80, 0., 2., 100, -4., 4., n, -.5, n-.5);
+ h3->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
+ h3->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
+ h3->GetZaxis()->SetTitle("SPECIES");
+ } else h3->Reset();
+ arr->AddAt(h3, 9);
+ // Kalman track P resolution
+ if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCP"))){
+ h3 = new TH3S("hMcTrkTPCP", "Track[TPC] P Resolution", 80, 0., 20., 150, -.15, .35, n, -.5, n-.5);
+ h3->GetXaxis()->SetTitle("p [GeV/c]");
+ h3->GetYaxis()->SetTitle("#Delta p/p^{MC}");
+ h3->GetZaxis()->SetTitle("SPECIES");
+ } else h3->Reset();
+ arr->AddAt(h3, 10);
+ // Kalman track Pt pulls
+ if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPPulls"))){
+ h3 = new TH3S("hMcTrkTPCPPulls", "Track[TPC] P Pulls", 80, 0., 20., 100, -5., 5., n, -.5, n-.5);
+ h3->GetXaxis()->SetTitle("p^{MC} [GeV/c]");
+ h3->GetYaxis()->SetTitle("#Deltap/#sigma_{p}");
+ h3->GetZaxis()->SetTitle("SPECIES");
+ } else h3->Reset();
+ arr->AddAt(h3, 11);
+
+
+
+ // Kalman track Z resolution [TOF]
+ fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackTOF]), kMCtrackTOF);
+ arr->SetName("McTrkTOF");
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZ"))){
+ h = new TH2I("hMcTrkTOFZ", "Track[TOF] Z Resolution", 100, -1., 1., 100, -1., 1.);
+ h->GetXaxis()->SetTitle("tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta z [cm]");
h->GetZaxis()->SetTitle("entries");
} else h->Reset();
- fContainer->AddAt(h, kMCtrackPt);
+ arr->AddAt(h, 0);
+ // Kalman track Z pulls
+ if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZPull"))){
+ h = new TH2I("hMcTrkTOFZPull", "Track[TOF] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
+ h->GetXaxis()->SetTitle("tg(#theta)");
+ h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
+ h->GetZaxis()->SetTitle("entries");
+ } else h->Reset();
+ arr->AddAt(h, 1);
return fContainer;
}
+//________________________________________________________
+Bool_t AliTRDresolution::Process(TH2* h2, TF1 *f, Float_t k, TGraphErrors **g)
+{
+ Char_t pn[10]; sprintf(pn, "p%02d", fIdxPlot);
+ 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);
+
+ for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
+ Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
+ TH1D *h = h2->ProjectionY(pn, ibin, ibin);
+ if(h->GetEntries()<100) continue;
+ //AdjustF1(h, f);
+
+ h->Fit(f, "QN");
+
+ Int_t ip = g[0]->GetN();
+ g[0]->SetPoint(ip, x, k*f->GetParameter(1));
+ g[0]->SetPointError(ip, 0., k*f->GetParError(1));
+ g[1]->SetPoint(ip, x, k*f->GetParameter(2));
+ g[1]->SetPointError(ip, 0., k*f->GetParError(2));
+
+/*
+ g[0]->SetPoint(ip, x, k*h->GetMean());
+ g[0]->SetPointError(ip, 0., k*h->GetMeanError());
+ g[1]->SetPoint(ip, x, k*h->GetRMS());
+ g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
+ }
+ fIdxPlot++;
+ return kTRUE;
+}
+
+//________________________________________________________
+Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
+{
+ if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
+
+ // retrive containers
+ TH2I *h2 = idx<0 ? (TH2I*)(fContainer->At(plot)) : (TH2I*)((TObjArray*)(fContainer->At(plot)))->At(idx);
+ if(!h2) return kFALSE;
+ TGraphErrors *g[2];
+ if(!(g[0] = idx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
+
+ if(!(g[1] = idx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
+
+ return Process(h2, f, k, g);
+}
+
+//________________________________________________________
+Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
+{
+ if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
+
+ // retrive containers
+ TH3S *h3 = idx<0 ? (TH3S*)(fContainer->At(plot)) : (TH3S*)((TObjArray*)(fContainer->At(plot)))->At(idx);
+ if(!h3) return kFALSE;
+
+ TObjArray *gm, *gs;
+ if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
+ if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
+ TGraphErrors *g[2];
+
+ TAxis *az = h3->GetZaxis();
+ for(Int_t iz=1; iz<=az->GetNbins(); iz++){
+ if(!(g[0] = (TGraphErrors*)gm->At(iz-1))) return kFALSE;
+ if(!(g[1] = (TGraphErrors*)gs->At(iz-1))) return kFALSE;
+ az->SetRange(iz, iz);
+ if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________
+Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
+{
+ if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
+
+ // retrive containers
+ TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
+ if(!h3) return kFALSE;
+
+ TGraphAsymmErrors *gm;
+ TGraphErrors *gs;
+ if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
+ if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
+
+ Float_t x, r, mpv, xM, xm;
+ TAxis *ay = h3->GetYaxis();
+ for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
+ ay->SetRange(iy, iy);
+ x = ay->GetBinCenter(iy);
+ TH2F *h2=(TH2F*)h3->Project3D("zx");
+ TAxis *ax = h2->GetXaxis();
+ for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
+ r = ax->GetBinCenter(ix);
+ TH1D *h1 = h2->ProjectionY("py", ix, ix);
+ if(h1->Integral()<50) continue;
+ h1->Fit(f, "QN");
+
+ GetLandauMpvFwhm(f, mpv, xm, xM);
+ Int_t ip = gm->GetN();
+ gm->SetPoint(ip, x, k*mpv);
+ gm->SetPointError(ip, 0., 0., k*xm, k*xM);
+ gs->SetPoint(ip, r, k*(xM-xm)/mpv);
+ gs->SetPointError(ip, 0., 0.);
+ }
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________
+Bool_t AliTRDresolution::Process4D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
+{
+ if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
+
+ // retrive containers
+ TObjArray *arr = (TObjArray*)((TObjArray*)(fContainer->At(plot)))->At(idx);
+ if(!arr) return kFALSE;
+
+
+ TObjArray *gm[2], *gs[2];
+ if(!(gm[0] = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
+ if(!(gs[0] = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
+
+ TGraphErrors *g[2];
+
+ TH3S *h3 = 0x0;
+ for(Int_t ix=0; ix<arr->GetEntriesFast(); ix++){
+ if(!(h3 = (TH3S*)arr->At(ix))) return kFALSE;
+ if(!(gm[1] = (TObjArray*)gm[0]->At(ix))) return kFALSE;
+ if(!(gs[1] = (TObjArray*)gs[0]->At(ix))) return kFALSE;
+ TAxis *az = h3->GetZaxis();
+ for(Int_t iz=1; iz<=az->GetNbins(); iz++){
+ if(!(g[0] = (TGraphErrors*)gm[1]->At(iz-1))) return kFALSE;
+ if(!(g[1] = (TGraphErrors*)gs[1]->At(iz-1))) return kFALSE;
+ az->SetRange(iz, iz);
+ if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
+ }
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________
+Bool_t AliTRDresolution::GetGraphPlot(Float_t *bb, ETRDresolutionPlot ip, Int_t idx)
+{
+ if(!fGraphS || !fGraphM) return kFALSE;
+ TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
+ if(!gm) return kFALSE;
+ TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
+ if(!gs) return kFALSE;
+ gs->Draw("apl"); gm->Draw("pl");
+
+ Double_t x,y;
+ gs->GetPoint(10, x, y);
+ PutTrendValue(Form("task%02d_idx%d", ip, idx), y, 0);
+
+ //printf("bb[%f %f %f %f]\n", bb[0], bb[1], bb[2], bb[3]);
+
+ // axis titles look up
+ Int_t nref = 0;
+ for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fNElements[jp];
+ UChar_t jdx = idx<0?0:idx;
+ for(Int_t jc=0; jc<TMath::Min(jdx,fNElements[ip]-1); jc++) nref++;
+ Char_t **at = fAxTitle[nref];
+
+ // axis range
+ TAxis *ax = 0x0;
+ ax = gs->GetHistogram()->GetXaxis();
+ ax->SetRangeUser(bb[0], bb[2]);
+ ax->SetTitle(*at);ax->CenterTitle();
+
+ ax = gs->GetHistogram()->GetYaxis();
+ ax->SetRangeUser(bb[1], bb[3]);
+ ax->SetTitleOffset(1.1);
+ ax->SetTitle(*(++at));ax->CenterTitle();
+
+ TGaxis *gax = 0x0;
+ gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
+ gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
+ //gax->SetVertical();
+ gax->CenterTitle(); gax->SetTitleOffset(.7);
+ gax->SetTitle(*(++at)); gax->Draw();
+
+ // bounding box
+ TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
+ b->SetFillStyle(3002);b->SetLineColor(0);
+ b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
+ b->Draw();
+ return kTRUE;
+}
+
+
+//________________________________________________________
+Bool_t AliTRDresolution::GetGraphTrack(Float_t *bb, Int_t idx, Int_t il)
+{
+ if(!fGraphS || !fGraphM) return kFALSE;
+
+ TGraphErrors *gm = 0x0, *gs = 0x0;
+ TObjArray *a0 = fGraphS, *a1 = 0x0;
+ a1 = (TObjArray*)a0->At(kMCtrackTRD); a0 = a1;
+ a1 = (TObjArray*)a0->At(idx); a0 = a1;
+ a1 = (TObjArray*)a0->At(il); a0 = a1;
+ for(Int_t is=0; is<AliPID::kSPECIES; is++){
+ if(!(gs = (TGraphErrors*)a0->At(is))) return kFALSE;
+ if(!gs->GetN()) continue;
+ gs->Draw(is ? "pl" : "apl");
+ }
+ gs = (TGraphErrors*)a0->At(0);
+ // axis titles look up
+ Int_t nref = 0;
+ for(Int_t jp=0; jp<Int_t(kMCtrackTRD); jp++) nref+=fNElements[jp];
+ for(Int_t jc=0; jc<idx; jc++) nref++;
+ Char_t **at = fAxTitle[nref];
+ // axis range
+ TAxis *ax = gs->GetHistogram()->GetXaxis();
+ ax->SetRangeUser(bb[0], bb[2]);
+ ax->SetTitle(*at);ax->CenterTitle();
+
+ ax = gs->GetHistogram()->GetYaxis();
+ ax->SetRangeUser(bb[1], bb[3]);
+ ax->SetTitleOffset(.5);ax->SetTitleSize(.06);
+ ax->SetTitle(*(++at));ax->CenterTitle();
+
+ TGaxis *gax = 0x0;
+ gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
+ gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
+ //gax->SetVertical();
+ gax->CenterTitle(); gax->SetTitleOffset(.5);gax->SetTitleSize(.06);
+ gax->SetTitle(*(++at)); gax->Draw();
+
+
+ a0 = fGraphM;
+ a1 = (TObjArray*)a0->At(kMCtrackTRD); a0 = a1;
+ a1 = (TObjArray*)a0->At(idx); a0 = a1;
+ a1 = (TObjArray*)a0->At(il); a0 = a1;
+ for(Int_t is=0; is<AliPID::kSPECIES; is++){
+ if(!(gm = (TGraphErrors*)a0->At(is))) return kFALSE;
+ if(!gm->GetN()) continue;
+ gm->Draw("pl");
+ }
+
+ return kTRUE;
+}
+
+
+//________________________________________________________
+Bool_t AliTRDresolution::GetGraphTrackTPC(Float_t *bb, Int_t sel)
+{
+ if(!fGraphS || !fGraphM) return kFALSE;
+
+ TGraphErrors *gm = 0x0, *gs = 0x0;
+ TObjArray *a0 = fGraphS, *a1 = 0x0;
+ a1 = (TObjArray*)a0->At(kMCtrackTPC); a0 = a1;
+ a1 = (TObjArray*)a0->At(sel); a0 = a1;
+ for(Int_t is=0; is<AliPID::kSPECIES; is++){
+ if(!(gs = (TGraphErrors*)a0->At(is))) return kFALSE;
+ if(!gs->GetN()) continue;
+ gs->Draw(is ? "pl" : "apl");
+ }
+ gs = (TGraphErrors*)a0->At(0);
+ // axis titles look up
+ Int_t nref = 0;
+ for(Int_t jp=0; jp<Int_t(kMCtrackTPC); jp++) nref+=fNElements[jp];
+ for(Int_t jc=0; jc<sel; jc++) nref++;
+ Char_t **at = fAxTitle[nref];
+ // axis range
+ TAxis *ax = gs->GetHistogram()->GetXaxis();
+ ax->SetRangeUser(bb[0], bb[2]);
+ ax->SetTitle(*at);ax->CenterTitle();
+
+ ax = gs->GetHistogram()->GetYaxis();
+ ax->SetRangeUser(bb[1], bb[3]);
+ ax->SetTitleOffset(1.);ax->SetTitleSize(0.05);
+ ax->SetTitle(*(++at));ax->CenterTitle();
+
+ TGaxis *gax = 0x0;
+ gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
+ gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
+ //gax->SetVertical();
+ gax->CenterTitle(); gax->SetTitleOffset(.7);gax->SetTitleSize(0.05);
+ gax->SetTitle(*(++at)); gax->Draw();
+
+
+ a0 = fGraphM;
+ a1 = (TObjArray*)a0->At(kMCtrackTPC); a0 = a1;
+ a1 = (TObjArray*)a0->At(sel); a0 = a1;
+ for(Int_t is=0; is<AliPID::kSPECIES; is++){
+ if(!(gm = (TGraphErrors*)a0->At(is))) return kFALSE;
+ if(!gm->GetN()) continue;
+ gm->Draw("pl");
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________
+void AliTRDresolution::GetLandauMpvFwhm(TF1 *f, Float_t &mpv, Float_t &xm, Float_t &xM)
+{
+ const Float_t dx = 1.;
+ mpv = f->GetParameter(1);
+ Float_t fx, max = f->Eval(mpv);
+
+ xm = mpv - dx;
+ while((fx = f->Eval(xm))>.5*max){
+ if(fx>max){
+ max = fx;
+ mpv = xm;
+ }
+ xm -= dx;
+ }
+
+ xM += 2*(mpv - xm);
+ while((fx = f->Eval(xM))>.5*max) xM += dx;
+}
+
//________________________________________________________
void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)