restructure error calculation :
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Jun 2009 14:57:27 +0000 (14:57 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Jun 2009 14:57:27 +0000 (14:57 +0000)
  - move all cluster error parameterization to AliTRDcluster
  - re-implement leading cluster method for radial calculation of pad
    row cross tracklets
  - implement pull plots for cluster calibration
fix minor errors in performance tasks
update documentation

15 files changed:
TRD/AliTRDcluster.cxx
TRD/AliTRDcluster.h
TRD/AliTRDseedV1.cxx
TRD/AliTRDseedV1.h
TRD/AliTRDtracker.cxx
TRD/AliTRDtransform.cxx
TRD/qaRec/AliTRDcheckDET.cxx
TRD/qaRec/AliTRDcheckPID.cxx
TRD/qaRec/AliTRDclusterResolution.cxx
TRD/qaRec/AliTRDclusterResolution.h
TRD/qaRec/AliTRDrecoTask.h
TRD/qaRec/AliTRDresolution.cxx
TRD/qaRec/info/AliTRDclusterInfo.cxx
TRD/qaRec/info/AliTRDclusterInfo.h
TRD/qaRec/macros/makeResults.C

index 3a693f6c6b8c0ade9c1ba3c2f131a0243eb4672f..18ad69557405694e16e0b7086d4b8685e039d05a 100644 (file)
@@ -136,7 +136,7 @@ AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
 
   SetY(c.GetY());
   SetZ(c.GetZ());
-  SetSigmaY2(c.GetSigmaY2());
+  AliCluster::SetSigmaY2(c.GetSigmaY2());
   SetSigmaZ2(c.GetSigmaZ2());  
 
   for (Int_t i = 0; i < 7; i++) {
@@ -241,11 +241,11 @@ void AliTRDcluster::Clear(Option_t *)
   fQ = 0;
   fCenter = 0;
   for (Int_t i = 0; i < 3; i++) SetLabel(0,i);
-  SetX(0);
-  SetY(0);
-  SetZ(0);
-  SetSigmaY2(0);
-  SetSigmaZ2(0);
+  SetX(0.);
+  SetY(0.);
+  SetZ(0.);
+  AliCluster::SetSigmaY2(0.);
+  SetSigmaZ2(0.);
   SetVolumeId(0);
 }
 
@@ -268,6 +268,19 @@ Float_t AliTRDcluster::GetSumS() const
 //___________________________________________________________________________
 Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z)
 {
+// Returns the error parameterization in the radial direction for TRD clusters as function of 
+// the calibrated time bin (tb) and optionally distance to anode wire (z). By default (no z information) 
+// the largest value over all cluster to wire values is chosen.
+// 
+// The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. 
+//
+//Begin_Html
+//<img src="TRD/clusterXerrorDiff2D.gif">
+//End_Html
+//
+// Author
+// A.Bercuci <A.Bercuci@gsi.de>
+
   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
   const Double_t sx[24][10]={
     {0.000e+00, 9.352e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 2.309e+00},
@@ -302,10 +315,68 @@ Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z)
 }
 
 //___________________________________________________________________________
-Double_t AliTRDcluster::GetSY(Int_t tb, Double_t z)
+Double_t AliTRDcluster::GetSYdrift(Int_t tb, Int_t ly, Double_t/* z*/)
 {
+// Returns the error parameterization for TRD clusters as function of the drift length (here calibrated time bin tb)
+// and optionally distance to anode wire (z) for the LUT r-phi cluster shape. By default (no z information) the largest 
+// value over all cluster to wire values is chosen.
+// 
+// For the LUT method the dependence of s_y with x and d is obtained via a fit to the cluster to MC 
+// resolution. (see class AliTRDclusterResolution for more details). A normalization to the reference radial position
+// x0 = 0.675 (tb=5 for ideal vd) is also applied (see GetSYprf()). The function is *NOT* calibration aware ! 
+// The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. A comparison 
+// with the GAUS parameterization is also given
+//
+// For the GAUS method the dependence of s_y with x is *analytic* and it is expressed by the relation.
+// BEGIN_LATEX
+// #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
+// END_LATEX
+// The result is displayed in the figure below as function of the drift time and compared with the LUT parameterization.
+//Begin_Html
+//<img src="TRD/clusterYerrorDiff2D.gif">
+//<img src="TRD/clusterYerrorDiff1D.gif">
+//End_Html
+//
+// Author
+// A.Bercuci <A.Bercuci@gsi.de>
+
+
   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
-  const Double_t sy[24][10]={
+  const Float_t lSy[6][24] = {
+    {75.7561, 0.0325, 0.0175, 0.0174, 0.0206, 0.0232,
+     0.0253, 0.0262, 0.0265, 0.0264, 0.0266, 0.0257,
+     0.0258, 0.0261, 0.0259, 0.0253, 0.0257, 0.0261,
+     0.0255, 0.0250, 0.0259, 0.0266, 0.0278, 0.0319
+    },
+    {49.2252, 0.0371, 0.0204, 0.0189, 0.0230, 0.0261,
+     0.0281, 0.0290, 0.0292, 0.0286, 0.0277, 0.0279,
+     0.0285, 0.0281, 0.0291, 0.0281, 0.0281, 0.0282,
+     0.0272, 0.0282, 0.0282, 0.0284, 0.0310, 0.0334
+    },
+    {55.1674, 0.0388, 0.0212, 0.0200, 0.0239, 0.0271,
+     0.0288, 0.0299, 0.0306, 0.0300, 0.0296, 0.0303,
+     0.0293, 0.0290, 0.0291, 0.0294, 0.0295, 0.0290,
+     0.0293, 0.0292, 0.0292, 0.0293, 0.0316, 0.0358
+    },
+    {45.1004, 0.0411, 0.0225, 0.0215, 0.0249, 0.0281,
+     0.0301, 0.0315, 0.0320, 0.0308, 0.0318, 0.0321,
+     0.0312, 0.0311, 0.0316, 0.0315, 0.0310, 0.0308,
+     0.0313, 0.0303, 0.0314, 0.0314, 0.0324, 0.0369
+    },
+    {43.8614, 0.0420, 0.0239, 0.0224, 0.0268, 0.0296,
+     0.0322, 0.0336, 0.0333, 0.0326, 0.0321, 0.0325,
+     0.0329, 0.0326, 0.0323, 0.0322, 0.0326, 0.0320,
+     0.0329, 0.0319, 0.0314, 0.0329, 0.0341, 0.0373
+    },
+    {40.5440, 0.0434, 0.0246, 0.0236, 0.0275, 0.0311,
+     0.0332, 0.0345, 0.0347, 0.0347, 0.0340, 0.0336,
+     0.0339, 0.0344, 0.0339, 0.0341, 0.0341, 0.0342,
+     0.0345, 0.0328, 0.0341, 0.0332, 0.0356, 0.0398
+    },
+  };
+  return lSy[ly][tb];
+
+/*  const Double_t sy[24][10]={
     {0.000e+00, 2.610e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 4.680e-01},
     {3.019e-02, 3.036e-02, 3.131e-02, 3.203e-02, 3.294e-02, 3.407e-02, 3.555e-02, 3.682e-02, 3.766e-02, 3.824e-02},
     {1.773e-02, 1.778e-02, 1.772e-02, 1.790e-02, 1.807e-02, 1.833e-02, 1.873e-02, 1.905e-02, 1.958e-02, 2.029e-02},
@@ -331,19 +402,83 @@ Double_t AliTRDcluster::GetSY(Int_t tb, Double_t z)
     {2.812e-02, 2.786e-02, 2.776e-02, 2.723e-02, 2.695e-02, 2.650e-02, 2.642e-02, 2.617e-02, 2.612e-02, 2.610e-02},
     {3.251e-02, 3.267e-02, 3.223e-02, 3.183e-02, 3.125e-02, 3.106e-02, 3.067e-02, 3.010e-02, 2.936e-02, 2.927e-02}
   };
-  if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)];
+  if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)] -  sy[5][Int_t(z/.025)];
 
-  Double_t m = 1.e-8; for(Int_t id=10; id--;) if(sy[tb][id]>m) m=sy[tb][id];
+  Double_t m = -1.e8; for(Int_t id=10; id--;) if((sy[tb][id] - sy[5][id])>m) m=sy[tb][id]-sy[5][id];
 
-  return m;
+  return m;*/
+}
+
+//___________________________________________________________________________
+Double_t AliTRDcluster::GetSYcharge(Float_t q)
+{
+// Parameterization of the r-phi resolution component due to cluster charge.
+// The value is the offset from the nominal cluster resolution defined as the
+// cluster resolution at average cluster charge (q0).
+// 
+// BEGIN_LATEX
+// #Delta #sigma_{y}(q) = a*(#frac{1}{q} - #frac{1}{q_{0}})
+// q_{0} #approx 50
+// END_LATEX
+// The definition is *NOT* robust against gain fluctuations and thus two approaches are possible
+// when residual miscalibration are available:
+//   - determine parameterization with a resolution matching those of the gain
+//   - define an analytic model which scales with the gain.
+// 
+// For more details please see AliTRDclusterResolution::ProcessCharge()
+//
+//Begin_Html
+//<img src="TRD/clusterQerror.gif">
+//End_Html
+//
+// Author
+// A.Bercuci <A.Bercuci@gsi.de>
+
+  const Float_t sq0inv = 0.019962; // [1/q0]
+  const Float_t sqb    = 1.0281564;// [cm]
+
+  return sqb*(1./q - sq0inv);
+}
+
+//___________________________________________________________________________
+Double_t AliTRDcluster::GetSYprf(Int_t ly, Double_t center, Double_t s2)
+{
+// Parameterization of the cluster error in the r-phi direction due to charge sharing between 
+// adiacent pads. Should be identical to what is provided in the OCDB as PRF [TODO]
+//
+// The parameterization is obtained from fitting cluster resolution at phi=exb and |x-0.675|<0.225. 
+// For more details see AliTRDclusterResolution::ProcessCenter().
+//
+//Begin_Html
+//<img src="TRD/clusterPRFerror.gif">
+//End_Html
+//
+// Author
+// A.Bercuci <A.Bercuci@gsi.de>
+
+/*  const Float_t scy[AliTRDgeometry::kNlayer][4] = {
+    {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
+    {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
+    {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
+    {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
+    {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
+    {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
+  };*/
+  const Float_t lPRF[] = {0.438, 0.403, 0.393, 0.382, 0.376, 0.345};
+
+  return s2*TMath::Gaus(center, 0., lPRF[ly]);
 }
 
+
 //___________________________________________________________________________
 Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z)
 {
-  // drift length correction [cm]
-  // TODO to be parametrized in term of drift velocity
-  // A.Bercuci (Mar 28 2009)
+// Drift length correction [cm]. Due to variation of mean drift velocity along the drift region
+// from nominal vd at xd->infinity. For drift velocity determination based on tracking information 
+// the correction should be negligible.
+//
+// TODO to be parametrized in term of drift velocity at infinite drift length
+// A.Bercuci (Mar 28 2009)
 
   if(tb<0 || tb>=24) return 0.;
   const Int_t nd = 5;
@@ -408,8 +543,8 @@ Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z)
 //___________________________________________________________________________
 Double_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y)
 {
-// PRF correction TODO to be replaced by the gaussian 
-// approximation with full error parametrization
+// PRF correction for the LUT r-phi cluster shape.
+
   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
     { 4.014e-04, 8.605e-03, -6.880e+00},
     {-3.061e-04, 9.663e-03, -6.789e+00},
@@ -465,7 +600,7 @@ Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/,
   // correction for t0
   td -= t0;
   // time bin corrected for t0
-  // BUG in TMath::Nint().root-5.23.02
+  // Bug in TMath::Nint().root-5.23.02
   // TMath::Nint(3.5) = 4 and TMath::Nint(4.5) = 4
   Double_t tmp = td*fFreq;
   fLocalTimeBin = Char_t(TMath::Floor(tmp));
@@ -483,6 +618,8 @@ Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/,
   return x;
 
 /*
+  // calculate radial posion of clusters in the drift region
+
   // invert drift time function
   Double_t xM= AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght(),
            x = vd*td + .5*AliTRDgeometry::CamHght(), 
@@ -510,7 +647,8 @@ Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W, Double_t *c
   //printf("  s[%3d %3d %3d] w[%f %f] yr[%f %f]\n", fSignals[2], fSignals[3], fSignals[4], w1/(w1+w2), w2/(w1+w2), y1r*W, y2r*W);
   if(IsRPhiMethod(kCOG)) GetDYcog();
   else if(IsRPhiMethod(kLUT)) GetDYlut();
-  if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2);
+  else if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2);
+  else return 0.;
 
   if(y1) (*y1)*=W;
   if(y2) (*y2)*=W;
@@ -518,6 +656,70 @@ Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W, Double_t *c
   return y0+fCenter*W+(IsRPhiMethod(kLUT)?GetYcorr(AliTRDgeometry::GetLayer(fDetector), fCenter):0.);
 }
 
+//___________________________________________________________________________
+void AliTRDcluster::SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x, Float_t z, Float_t tgp)
+{
+// Set variance of TRD cluster in the r-phi direction for each method.
+// Parameters :
+//   - s2  - variance due to PRF width for the case of Gauss model. Replaced by parameterization in case of LUT.
+//   - dt  - transversal diffusion coeficient 
+//   - exb - tg of lorentz angle
+//   - x   - drift length - with respect to the anode wire
+//   - z   - offset from the anode wire
+//   - tgp - local tangent of the track momentum azimuthal angle
+//
+// The ingredients from which the error is computed are:
+//   - PRF (charge sharing on adjacent pads)  - see AliTRDcluster::GetSYprf()
+//   - diffusion (dependence with drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSYdrift()
+//   - charge of the cluster (complex dependence on gain and tail cancellation) - see AliTRDcluster::GetSYcharge()
+//   - lorentz angle (dependence on the drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSX()
+//   - track angle (superposition of charges on the anode wire) - see AliTRDseedV1::Fit()
+//   - projection of radial(x) error on r-phi due to fixed value assumed in tracking for x - see AliTRDseedV1::Fit()
+//
+// The last 2 contributions to cluster error can be estimated only during tracking when the track angle 
+// is known (tgp). For this reason the errors (and optional position) of TRD clusters are recalculated during 
+// tracking and thus clusters attached to tracks might differ from bare clusters.
+// 
+// Author:
+// A.Bercuci <A.Bercuci@gsi.de>
+
+  Float_t sigmaY2 = 0.;
+  Int_t ly = AliTRDgeometry::GetLayer(fDetector);
+  if(IsRPhiMethod(kCOG)) sigmaY2 = 4.e-4;
+  else if(IsRPhiMethod(kLUT)){ 
+    Float_t sd = GetSYdrift(fLocalTimeBin, ly, z);//printf("drift[%6.2f] ", 1.e4*sd);
+    sigmaY2 = GetSYprf(ly, fCenter, sd);//printf("PRF[%6.2f] ", 1.e4*sigmaY2);
+    // add charge contribution TODO scale with respect to s2
+    sigmaY2+= GetSYcharge(TMath::Abs(fQ));//printf("Q[%6.2f] ", 1.e4*sigmaY2);
+    sigmaY2 = TMath::Max(sigmaY2, Float_t(0.)); //!! protection 
+    sigmaY2*= sigmaY2;
+  } else if(IsRPhiMethod(kGAUS)){
+    // PRF contribution
+    sigmaY2 = s2;
+    // Diffusion contribution
+    Double_t sD2 = dt/(1.+exb); sD2 *= sD2; sD2 *= x;
+    sigmaY2+= sD2; 
+    // add charge contribution TODO scale with respect to s2
+    //sigmaY2+= GetSYcharge(TMath::Abs(fQ));
+  }
+
+  // store tg^2(phi-a_L) and tg^2(a_L)
+  Double_t tgg = (tgp-exb)/(1.+tgp*exb); tgg *= tgg;
+  Double_t exb2= exb*exb;
+
+  // Lorentz angle shift contribution 
+  Float_t sx = GetSX(fLocalTimeBin, z); sx*=sx;
+  sigmaY2+= exb2*sx;//printf("Al[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
+
+  // Radial contribution due to not measuring x in Kalman model 
+  sigmaY2+= tgg*sx;//printf("x[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
+
+  // Track angle contribution
+  sigmaY2+= tgg*x*x*exb2/12.;//printf("angle[%6.2f]\n", 1.e4*TMath::Sqrt(sigmaY2));
+
+  AliCluster::SetSigmaY2(sigmaY2);
+}
+
 //_____________________________________________________________________________
 Bool_t AliTRDcluster::IsEqual(const TObject *o) const
 {
index 39c820bc4cbcb7b4c09dc11961059271a666a1df..409cf11e0f8e4572e431c1535b7e16ed69c2aca2 100644 (file)
@@ -62,7 +62,9 @@ public:
   Short_t *GetSignals()                    { return fSignals;       }
   Float_t  GetSumS() const;
   static Double_t  GetSX(Int_t tb, Double_t z=-1);
-  static Double_t  GetSY(Int_t tb, Double_t z=-1);
+  static Double_t  GetSYdrift(Int_t tb, Int_t ly=0, Double_t z=-1);
+  static Double_t  GetSYcharge(Float_t q);
+  static Double_t  GetSYprf(Int_t ly, Double_t center, Double_t s2);
   static Double_t  GetXcorr(Int_t tb, Double_t z=-1);
   static Double_t  GetYcorr(Int_t ly, Float_t y);
   Float_t  GetXloc(Double_t t0, Double_t vd, Double_t *const q=0x0, Double_t *const xq=0x0, Double_t z = 0.2);
@@ -83,6 +85,7 @@ public:
   void     SetQ(Float_t inQ){ fQ = inQ;}
   void     SetClusterMasking(UChar_t inClusterMasking){ fClusterMasking = inClusterMasking;}
   void     SetShared(Bool_t sh  = kTRUE)   { SetBit(AliCluster::kShared,sh);    }
+  void     SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x0, Float_t z=-1., Float_t tgp=0.);
   void     Use(Int_t u = 1)                  { SetBit(AliCluster::kUsed, u ? kTRUE : kFALSE);              }
   void     SetFivePad(Bool_t b = kTRUE) { SetBit(kFivePad,b);}
 
index 8cd43f2b490013f297694754106c604b2ca799cb..be2f38cfe26e67d7f0af11a715c9bb528b45ec0b 100644 (file)
@@ -293,8 +293,9 @@ void AliTRDseedV1::Update(const AliTRDtrackV1 *trk)
   Double_t fSnp = trk->GetSnp();
   Double_t fTgl = trk->GetTgl();
   fPt = trk->Pt();
-  fYref[1] = fSnp/TMath::Sqrt(1. - fSnp*fSnp);
-  fZref[1] = fTgl;
+  Double_t norm =1./TMath::Sqrt(1. - fSnp*fSnp); 
+  fYref[1] = fSnp*norm;
+  fZref[1] = fTgl*norm;
   SetCovRef(trk->GetCovariance());
 
   Double_t dx = trk->GetX() - fX0;
@@ -474,6 +475,7 @@ Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const
   }
   if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
   if(dq<1.e-3) return 0.;
+  
 
   Double_t dx = fdX;
   if(ic-1>=0 && ic+1<kNtb){
@@ -830,21 +832,28 @@ void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *p)
 //____________________________________________________________________
 Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
 {
-  //
-  // Projective algorithm to attach clusters to seeding tracklets
-  //
-  // Parameters
-  //
-  // Output
-  //
-  // Detailed description
-  // 1. Collapse x coordinate for the full detector plane
-  // 2. truncated mean on y (r-phi) direction
-  // 3. purge clusters
-  // 4. truncated mean on z direction
-  // 5. purge clusters
-  // 6. fit tracklet
-  //   
+//
+// Projective algorithm to attach clusters to seeding tracklets. The following steps are performed :
+// 1. Collapse x coordinate for the full detector plane
+// 2. truncated mean on y (r-phi) direction
+// 3. purge clusters
+// 4. truncated mean on z direction
+// 5. purge clusters
+//
+// Parameters
+//  - chamber : pointer to tracking chamber container used to search the tracklet
+//  - tilt    : switch for tilt correction during road building [default true]
+// Output
+//  - true    : if tracklet found successfully. Failure can happend because of the following:
+//      -
+// Detailed description
+//     
+// We start up by defining the track direction in the xy plane and roads. The roads are calculated based
+// on tracking information (variance in the r-phi direction) and estimated variance of clusters. For 
+// estimating the last contribution the following expression is considered :  
+// 
+// Author Alexandru Bercuci <A.Bercuci@gsi.de>
+
   Bool_t kPRINT = kFALSE;
   if(!fReconstructor->GetRecoParam() ){
     AliError("Seed can not be used without a valid RecoParam.");
@@ -1007,12 +1016,10 @@ Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
     if(!fClusters[it]) continue;
     x[irp]  = fClusters[it]->GetX();
     tb[irp] = fClusters[it]->GetLocalTimeBin();
-    printf("  x[%d]=%f t[%d]=%d\n", irp, x[irp], irp, tb[irp]);
     irp++;
   }  
   Int_t dtb = tb[1] - tb[0];
   fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
-
   return kTRUE;
 }
 
@@ -1063,7 +1070,7 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
 //   - tilt : switch for tilt pad correction of cluster y position based on 
 //            the z, dzdx info from outside [default false].
 //   - zcorr : switch for using z information to correct for anisochronity 
-//            and a finner error parametrization estimation [default false]  
+//            and a finner error parameterization estimation [default false]  
 // Output :
 //  True if successful
 //
@@ -1071,27 +1078,87 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
 //
 //            Fit in the xy plane
 // 
+// The fit is performed to estimate the y position of the tracklet and the track 
+// angle in the bending plane. The clusters are represented in the chamber coordinate 
+// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation() 
+// on how this is set). The x and y position of the cluster and also their variances 
+// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(), 
+// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()). 
+// If gaussian approximation is used to calculate y coordinate of the cluster the position 
+// is recalculated taking into account the track angle. The general formula to calculate the 
+// error of cluster position in the gaussian approximation taking into account diffusion and track
+// inclination is given for TRD by:
+// BEGIN_LATEX
+// #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}} + #frac{x^{2}tg^{2}(#phi-#alpha_{L})tg^{2}(#alpha_{L})}{12}
+// END_LATEX
+//
+// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y
+// by projection i.e.
+// BEGIN_LATEX
+// #sigma_{x|y} = tg(#phi) #sigma_{x}
+// END_LATEX
+// and also by the lorentz angle correction
+//
+//            Fit in the xz plane
+//
+// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens. 
+// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy 
+// fit (see below).
+// 
+// There are two methods to estimate the radial position of the pad row cross:
+//   1. leading cluster radial position : Here the lower part of the tracklet is considered and the last 
+// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error 
+// of the z estimate is given by :
+// BEGIN_LATEX
+// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
+// END_LATEX
+// The systematic errors for this estimation are generated by the following sources:
+//   - no charge sharing between pad rows is considered (sharp cross)
+//   - missing cluster at row cross (noise peak-up, under-threshold signal etc.).
+// 
+//   2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered 
+// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are 
+// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources:
+//   - no general model for the qx dependence
+//   - physical fluctuations of the charge deposit 
+//   - gain calibration dependence
+//
+//            Estimation of the radial position of the tracklet
 //
+// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the 
+// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error
+// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()):
+// BEGIN_LATEX
+// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
+// END_LATEX
+// and thus the radial position is:
+// BEGIN_LATEX
+// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
+// END_LATEX
+//
+//            Estimation of tracklet position error 
+//
+// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z 
+// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by:
+// BEGIN_LATEX
+// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
+// #sigma_{z} = Pad_{length}/12
+// END_LATEX
+// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error 
+// in z by the width of the crossing region - being a matter of parameterization. 
+// BEGIN_LATEX
+// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
+// END_LATEX
+// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of
+// the covariance matrix. See AliTRDseedV1::GetCovAt() for details.
+//
+// Author 
+// A.Bercuci <A.Bercuci@gsi.de>
 
   if(!IsCalibrated()) Calibrate();
 
   const Int_t kClmin = 8;
 
-
-  // cluster error parametrization parameters 
-  // 1. sy total charge
-  const Float_t sq0inv = 0.019962; // [1/q0]
-  const Float_t sqb    = 1.0281564;    //[cm]
-  // 2. sy for the PRF
-  const Float_t scy[AliTRDgeometry::kNlayer][4] = {
-    {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
-    {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
-    {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
-    {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
-    {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
-    {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
-  };
-
   // get track direction
   Double_t y0   = fYref[0];
   Double_t dydx = fYref[1]; 
@@ -1099,10 +1166,6 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
   Double_t dzdx = fZref[1];
   Double_t yt, zt;
 
-  // calculation of tg^2(phi - a_L) and tg^2(a_L)
-  Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
-  //Double_t exb2= fExB*fExB;
-
   //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
   TLinearFitter  fitterY(1, "pol1");
   TLinearFitter  fitterZ(1, "pol1");
@@ -1110,11 +1173,9 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
   // book cluster information
   Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
 
-  Int_t ily = AliTRDgeometry::GetLayer(fDet);
   Int_t n = 0;
   AliTRDcluster *c=0x0, **jc = &fClusters[0];
   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
-    //zRow[ic] = -1;
     xc[ic]  = -1.;
     yc[ic]  = 999.;
     zc[ic]  = 999.;
@@ -1125,59 +1186,31 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
     Float_t w = 1.;
     if(c->GetNPads()>4) w = .5;
     if(c->GetNPads()>5) w = .2;
-    Int_t tb = c->GetLocalTimeBin();
 
+    // cluster charge
     qc[n]   = TMath::Abs(c->GetQ());
+    // pad row of leading 
+
     // Radial cluster position
     //Int_t jc = TMath::Max(fN-3, 0);
     //xc[fN]   = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
     xc[n]   = fX0 - c->GetX();
 
-    //Double_t s2 = fS2PRF + fDiffL*fDiffL*xc[n]/(1.+2.*exb2)+tgg*xc[n]*xc[n]*exb2/12.;
-    //yc[fN]   = c->GetYloc(s2, GetPadWidth(), xc[fN], fExB);
-    yc[n]   = c->GetY()-AliTRDcluster::GetYcorr(ily, c->GetCenter());
-    zc[n]   = c->GetZ();
-
-    // extrapolated y value for the track
+    // extrapolated track to cluster position
     yt = y0 - xc[n]*dydx; 
-    // extrapolated z value for the track
     zt = z0 - xc[n]*dzdx; 
-    // tilt correction
-    if(tilt) yc[n] -= GetTilt()*(zc[n] - zt); 
-
-    // ELABORATE CLUSTER ERROR
-    // basic y error (|| to track).
-    sy[n]  = AliTRDcluster::GetSY(tb, zcorr?zt:-1.);
-    //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN,  sy[fN]*1.e4);
-    // y error due to total charge
-    sy[n] += sqb*(1./qc[n] - sq0inv);
-    //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
-    // y error due to PRF
-    sy[n] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
-    //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
-
-    sy[n] *= sy[n];
-
-    // ADD ERROR ON x
-    // error of drift length parallel to the track
-    Double_t sx = AliTRDcluster::GetSX(tb, zcorr?zt:-1.); // [cm]
-    //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
-    sx *= sx; // square sx
-
-    // add error from ExB 
-    sy[n] += fExB*fExB*sx;
-    //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
-
-    // global radial error due to misalignment/miscalibration
-    Double_t sx0  = 0.; sx0 *= sx0;
-    // add sx contribution to sy due to track angle
-    sy[n] += tgg*(sx+sx0);
-    // TODO we should add tilt pad correction here
-    //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
-    c->SetSigmaY2(sy[n]);
-
-    sy[n]  = TMath::Sqrt(sy[n]);
-    fitterY.AddPoint(&xc[n], yc[n], sy[n]);
+
+    // Recalculate cluster error based on tracking information
+    c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx);
+    sy[n]  = TMath::Sqrt(c->GetSigmaY2());
+
+    yc[n]   = fReconstructor->UseGAUS() ? 
+      c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
+    zc[n]   = c->GetZ();
+    //optional tilt correction
+    if(tilt) yc[n] -= (GetTilt()*(zc[n] - zt)); 
+
+    fitterY.AddPoint(&xc[n], yc[n], TMath::Sqrt(sy[n]));
     fitterZ.AddPoint(&xc[n], qc[n], 1.);
     n++;
   }
@@ -1199,7 +1232,24 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
 
   // fit XZ
   if(IsRowCross()){
+    // THE LEADING CLUSTER METHOD
+    Float_t xMin = fX0;
     Int_t ic=n=kNclusters-1; jc = &fClusters[ic];
+    AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1];
+    for(; ic>kNtb; ic--, --jc, --kc){
+      if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX();
+      if(!(c = (*jc))) continue;
+      if(!c->IsInChamber()) continue;
+      zc[kNclusters-1] = c->GetZ(); 
+      fX = fX0 - c->GetX();
+    }
+    fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.;
+    // Error parameterization
+    fS2Z     = fdX*fZref[1];
+    fS2Z    *= fS2Z; fS2Z    *= 0.2887; //  1/sqrt(12)
+/*
+    // THE FIT X-Q PLANE METHOD 
+    ic=n=kNclusters-1; jc = &fClusters[ic];
     for(; ic>kNtb; ic--, --jc){
       if(!(c = (*jc))) continue;
       if(!c->IsInChamber()) continue;
@@ -1224,51 +1274,13 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr)
     fS2Z     = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
     // TODO correct formula
     //fS2Z     = sigma_x*TMath::Abs(fZref[1]);
+*/
   } else {
     fZfit[0] = zc[0]; fZfit[1] = 0.;
     fS2Z     = GetPadLength()*GetPadLength()/12.;
   }
   fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
   return kTRUE;
-//   // determine z offset of the fit
-//   Float_t zslope = 0.;
-//   Int_t nchanges = 0, nCross = 0;
-//   if(nz==2){ // tracklet is crossing pad row
-//     // Find the break time allowing one chage on pad-rows
-//     // with maximal number of accepted clusters
-//     Int_t padRef = zRow[0];
-//     for (Int_t ic=1; ic<fN; ic++) {
-//       if(zRow[ic] == padRef) continue;
-//       
-//       // debug
-//       if(zRow[ic-1] == zRow[ic]){
-//         printf("ERROR in pad row change!!!\n");
-//       }
-//     
-//       // evaluate parameters of the crossing point
-//       Float_t sx = (xc[ic-1] - xc[ic])*convert;
-//       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
-//       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
-//       fCross[3] = TMath::Max(dzdx * sx, .01);
-//       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
-//       padRef    = zRow[ic];
-//       nCross    = ic;
-//       nchanges++;
-//     }
-//   }
-// 
-//   // condition on nCross and reset nchanges TODO
-// 
-//   if(nchanges==1){
-//     if(dzdx * zslope < 0.){
-//       AliInfo("Tracklet-Track mismatch in dzdx. TODO.");
-//     }
-// 
-// 
-//     //zc[nc] = fitterZ.GetFunctionParameter(0); 
-//     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
-//     fCross[0] = fX0 - fCross[0];
-//   }
 }
 
 
index b7d3d127ef5c27ab610bfc37583500d3c146a055..6429aa1819762224761ebde32c5f3bafb7741c5f 100644 (file)
@@ -59,7 +59,7 @@ public:
    ,kPID        = BIT(16) // PID contributor
    ,kCalib      = BIT(17) // calibrated tracklet
    ,kKink       = BIT(18) // kink prolongation tracklet
-   ,kStandAlone = BIT(19) // stand alone build tracklet
+   ,kStandAlone = BIT(19) // tracklet build during stand alone track finding
   };
 
   AliTRDseedV1(Int_t det = -1);
@@ -67,11 +67,7 @@ public:
   AliTRDseedV1(const AliTRDseedV1 &ref);
   AliTRDseedV1& operator=(const AliTRDseedV1 &ref);
 
-/*  Bool_t       AttachClustersIter(
-              AliTRDtrackingChamber *chamber, Float_t quality, 
-              Bool_t kZcorr = kFALSE, AliTRDcluster *c=0x0);*/
-  Bool_t         AttachClusters(
-              AliTRDtrackingChamber *chamber, Bool_t tilt = kFALSE);
+  Bool_t    AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt = kFALSE);
   void      Bootstrap(const AliTRDReconstructor *rec);
   void      Calibrate();
   void      CookdEdx(Int_t nslices);
@@ -129,7 +125,7 @@ public:
   Float_t   GetS2Z() const           { return fS2Z;}
   Float_t   GetSigmaY() const        { return fS2Y > 0. ? TMath::Sqrt(fS2Y) : 0.2;}
   Float_t   GetSnp() const           { return fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);}
-  Float_t   GetTgl() const           { return fZref[1];}
+  Float_t   GetTgl() const           { return fZref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);}
   Float_t   GetTilt() const          { return fPad[2];}
   UInt_t    GetTrackletWord() const  { return 0;}
   UShort_t  GetVolumeId() const;
@@ -137,14 +133,14 @@ public:
   Float_t   GetX() const             { return fX0 - fX;}
   Float_t   GetY() const             { return fYfit[0] - fYfit[1] * fX;}
   Double_t  GetYat(Double_t x) const { return fYfit[0] - fYfit[1] * (fX0-x);}
-  Float_t   GetYfit(Int_t id) const { return fYfit[id];}
-  Float_t   GetYref(Int_t id) const { return fYref[id];}
-  Float_t   GetZ() const            { return fZfit[0] - fZfit[1] * fX;}
+  Float_t   GetYfit(Int_t id) const  { return fYfit[id];}
+  Float_t   GetYref(Int_t id) const  { return fYref[id];}
+  Float_t   GetZ() const             { return fZfit[0] - fZfit[1] * fX;}
   Double_t  GetZat(Double_t x) const { return fZfit[0] - fZfit[1] * (fX0-x);}
-  Float_t   GetZfit(Int_t id) const { return fZfit[id];}
-  Float_t   GetZref(Int_t id) const { return fZref[id];}
-  Int_t     GetYbin() const         { return Int_t(GetY()/0.016);}
-  Int_t     GetZbin() const         { return Int_t(GetZ()/fPad[0]);}
+  Float_t   GetZfit(Int_t id) const  { return fZfit[id];}
+  Float_t   GetZref(Int_t id) const  { return fZref[id];}
+  Int_t     GetYbin() const          { return Int_t(GetY()/0.016);}
+  Int_t     GetZbin() const          { return Int_t(GetZ()/fPad[0]);}
 
   inline AliTRDcluster* NextCluster();
   inline AliTRDcluster* PrevCluster();
@@ -152,8 +148,8 @@ public:
   inline void ResetClusterIter(Bool_t forward = kTRUE);
   void      Reset();
 
-  void      SetC(Float_t c)         { fC = c;}
-  void      SetChi2(Float_t chi2)   { fChi2 = chi2;}
+  void      SetC(Float_t c)          { fC = c;}
+  void      SetChi2(Float_t chi2)    { fChi2 = chi2;}
   inline void SetCovRef(const Double_t *cov);
   void      SetIndexes(Int_t i, Int_t idx) { fIndexes[i]  = idx; }
   void      SetLabels(Int_t *lbls)   { memcpy(fLabels, lbls, 3*sizeof(Int_t)); }
@@ -199,10 +195,10 @@ private:
   Short_t          fDet;                    // TRD detector
   AliTRDcluster   *fClusters[kNclusters];   // Clusters
   Float_t          fPad[3];                 // local pad definition : length/width/tilt 
-  Float_t          fYref[2];                //  Reference y
-  Float_t          fZref[2];                //  Reference z
-  Float_t          fYfit[2];                //  Y fit position +derivation
-  Float_t          fZfit[2];                //  Z fit position
+  Float_t          fYref[2];                //  Reference y, dydx
+  Float_t          fZref[2];                //  Reference z, dz/dx
+  Float_t          fYfit[2];                //  Fit y, dy/dx
+  Float_t          fZfit[2];                //  Fit z
   Float_t          fPt;                     //  Pt estimate @ tracklet [GeV/c]
   Float_t          fdX;                     // length of time bin
   Float_t          fX0;                     // anode wire position
@@ -214,7 +210,7 @@ private:
   Float_t          fC;                      // Curvature
   Float_t          fChi2;                   // Global chi2  
   Float_t          fdEdx[kNslices];         // dE/dx measurements for tracklet
-  Float_t          fProb[AliPID::kSPECIES]; //  PID probabilities
+  Float_t          fProb[AliPID::kSPECIES]; // PID probabilities
   Int_t            fLabels[3];              // most frequent MC labels and total number of different labels
   Double_t         fRefCov[7];              // covariance matrix of the track in the yz plane + the rest of the diagonal elements
   Double_t         fCov[3];                 // covariance matrix of the tracklet in the xy plane
index 2ff6df2e36570cee14019158f0b828d85b1fea6c..918f316cfb09d8c026f8518b595fcfbfd1578b5b 100644 (file)
@@ -3555,7 +3555,7 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
     }
 
     // Set cluster error
-    cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); 
+    ((AliCluster*)cl[best[bestiter][it]][it])->SetSigmaY2(expectederr); 
     if (!cl[best[bestiter][it]][it]->IsUsed()) {
       cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY()); 
     }
index ca00a70919dd2c5d586c3c22c764b1e89b1b9184..f5d12ba2573aca710b5970b5628d6edf33c78a46 100644 (file)
@@ -234,16 +234,16 @@ Bool_t AliTRDtransform::Transform(AliTRDcluster *c)
 
   Float_t x = c->GetXloc(t0, vd);
 
-  // pad response width with diffusion corrections
-  Double_t s2  = fCalPRFROC->GetValue(col, row); s2 *= s2; 
-  Float_t dl, dt;
-  AliTRDCommonParam::Instance()->GetDiffCoeff(dl, dt, vd);
-  s2 += dl*dl*x/(1.+2.*exb*exb);
-  s2 -= - 1.5e-1;
-
   // Pad dimensions
   Double_t rs = fPadPlane->GetRowSize(row);
   Double_t cs = fPadPlane->GetColSize(col);
+
+  // cluster error with diffusion corrections
+  Double_t s2  = cs*fCalPRFROC->GetValue(col, row); 
+  s2 *= s2; 
+  Float_t dl, dt;
+  AliTRDCommonParam::Instance()->GetDiffCoeff(dl, dt, vd);
+
   Double_t y0 = fPadPlane->GetColPos(col) + .5*cs;
   Double_t loc[] = {
     kX0shift-x,                    // Invert the X-position,
@@ -257,7 +257,7 @@ Bool_t AliTRDtransform::Transform(AliTRDcluster *c)
 
   // store tracking values
   c->SetX(trk[0]);c->SetY(trk[1]);c->SetZ(trk[2]);
-  c->SetSigmaY2(s2);
+  c->SetSigmaY2(s2, dt, exb, x);
   c->SetSigmaZ2(fPadPlane->GetRowSize(row)*fPadPlane->GetRowSize(row)/12.);
   
   return kTRUE;
index 8caef1f4a631f6048313862e6c5fadeb2a864bac..3789e55b525bb2c094dab43bcdc946c62310d66c 100644 (file)
@@ -183,6 +183,8 @@ Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
   //
   // Setting Reference Figures
   //
+  gPad->SetLogy(0);
+  gPad->SetLogx(0);
   switch(ifig){
   case kNclustersTrack:
     ((TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
@@ -213,6 +215,7 @@ Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
     return kTRUE;
   case kChargeCluster:
     ((TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
+    gPad->SetLogy(1);
     return kTRUE;
   case kChargeTracklet:
     ((TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
@@ -904,17 +907,25 @@ void AliTRDcheckDET::MakePlotChi2()
 
   TH2S *h2 = (TH2S*)fContainer->At(kChi2);
   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
-
+  TLegend *leg = new TLegend(.7,.7,.95,.95);
+  leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
   Bool_t kFIRST = kTRUE;
   for(Int_t il=1; il<=h2->GetNbinsX(); il++){
-    TH1D *h1 = h2->ProjectionY(Form("py%d", il), il, il);
+    TH1D *h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
     if(h1->Integral()<50) continue;
     h1->Scale(1./h1->Integral());
-    h1->SetMarkerStyle(7);h1->SetMarkerColor(il);h1->SetLineColor(il);
+    h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
+    h1->SetLineColor(il);h1->SetLineStyle(2);
     f.SetParameter(1, .5*il);f.SetLineColor(il);
-    h1->Fit(&f, "Q+", kFIRST ? "e1": "e1 same");
+    h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
+    leg->AddEntry(h1, Form("%d", il), "l");
+    if(kFIRST){
+      h1->GetXaxis()->SetRangeUser(0., 25.);
+    }
     kFIRST = kFALSE;
   }
+  leg->Draw();
+  gPad->SetLogy();
 }
 
 
@@ -926,41 +937,37 @@ void AliTRDcheckDET::MakePlotNTracklets(){
   TH1F *hBAR = (TH1F *)fContainer->FindObject("hNtlsBAR");
   TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
   TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
+  TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89);
+
+  Float_t scale = hCON->Integral();
+  hCON->Scale(100./scale);
+  hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
+  hCON->SetBarWidth(0.2);
+  hCON->SetBarOffset(0.6);
+  hCON->SetStats(kFALSE);
+  hCON->GetYaxis()->SetRangeUser(0.,40.);
+  hCON->GetYaxis()->SetTitleOffset(1.2);
+  hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
 
-  hBAR->Scale(100./hCON->Integral());
-  hBAR->SetFillColor(kRed);
+  hBAR->Scale(100./scale);
+  hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
   hBAR->SetBarWidth(0.2);
   hBAR->SetBarOffset(0.2);
   hBAR->SetTitle("");
   hBAR->SetStats(kFALSE);
   hBAR->GetYaxis()->SetRangeUser(0.,40.);
   hBAR->GetYaxis()->SetTitleOffset(1.2);
-  hBAR->Draw("bar1");
+  hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
 
-  hSTA->Scale(100./hCON->Integral());
-  hSTA->SetFillColor(kBlue);
+  hSTA->Scale(100./scale);
+  hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
   hSTA->SetBarWidth(0.2);
   hSTA->SetBarOffset(0.4);
   hSTA->SetTitle("");
   hSTA->SetStats(kFALSE);
   hSTA->GetYaxis()->SetRangeUser(0.,40.);
   hSTA->GetYaxis()->SetTitleOffset(1.2);
-  hSTA->Draw("bar1same");
-
-  hCON->Scale(100./hCON->Integral());
-  hCON->SetFillColor(kGreen);
-  hCON->SetBarWidth(0.2);
-  hCON->SetBarOffset(0.6);
-  hCON->SetStats(kFALSE);
-  hCON->GetYaxis()->SetRangeUser(0.,40.);
-  hCON->GetYaxis()->SetTitleOffset(1.2);
-  hCON->Draw("bar1same");
-
-  TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89);
-  leg->AddEntry(hBAR, "Barrel", "f");
-  leg->AddEntry(hSTA, "Stand Alone", "f");
-  leg->AddEntry(hCON, "Convoluted", "f");
-
+  hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
   leg->Draw();
   gPad->Update();
 }
index 8a2c1f0133560db32e0513064fa9a8b39f714cb7..dca4d304d078fecb1d72f6963a6d66c96d3c0efe 100644 (file)
@@ -735,10 +735,10 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
   TList *content = 0x0;
   switch(ifig){
   case kEfficiency:{
-    gPad->Divide(2, 1, 0.,0.);
+    gPad->Divide(2, 1, 1.e-5, 1.e-5);
     TList *l=gPad->GetListOfPrimitives();
     TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd();
-    pad->SetMargin(0.07, 0.07, 0.1, 0.001);
+    pad->SetMargin(0.1, 0.01, 0.1, 0.01);
 
     TLegend *legEff = new TLegend(.7, .7, .98, .98);
     legEff->SetBorderSize(1);
@@ -768,7 +768,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
     gPad->SetGridx();
 
     pad = ((TVirtualPad*)l->At(1));pad->cd();
-    pad->SetMargin(0.07, 0.07, 0.1, 0.001);
+    pad->SetMargin(0.1, 0.01, 0.1, 0.01);
     content = (TList *)fGraph->FindObject("Thresholds");
     if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break;
     if(!g->GetN()) break;
@@ -796,6 +796,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
     FIRST = kTRUE;
     if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break;
     legdEdx->SetHeader("Particle Species");
+    gPad->SetMargin(0.1, 0.01, 0.1, 0.01);
     for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){
       Int_t bin = FindBin(is, 2.);
       h1 = h2->ProjectionY("px", bin, bin);
@@ -827,6 +828,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
     FIRST = kTRUE;
     if(!(h2 = (TH2F*)(fContainer->At(kPH)))) break;;
     legPH->SetHeader("Particle Species");
+    gPad->SetMargin(0.1, 0.01, 0.1, 0.01);
     for(Int_t is=0; is<AliPID::kSPECIES; is++){
       Int_t bin = FindBin(is, 2.);
       h1 = h2->ProjectionY("py", bin, bin);
@@ -835,8 +837,8 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
       h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
       if(FIRST){
-        h1->GetXaxis()->SetTitle("tb(1/100 ns^{-1})");
-        h1->GetYaxis()->SetTitle("<PH> (a.u.)");
+        h1->GetXaxis()->SetTitle("t_{drift} [100*ns]");
+        h1->GetYaxis()->SetTitle("<dQ/dl> [a.u./cm]");
       }
       h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
       legPH->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl");
@@ -859,21 +861,21 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
     legNClus->SetHeader("Particle Species");
     for(Int_t is=0; is<AliPID::kSPECIES; is++){
       Int_t bin = FindBin(is, 2.);
-      h1 = h2->ProjectionY("py", bin, bin);
+      h1 = h2->ProjectionY("pyNClus", bin, bin);
       if(!h1->GetEntries()) continue;
-      h1->Scale(1./h1->Integral());
+      h1->Scale(100./h1->Integral());
       //h1->SetMarkerStyle(24);
       //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
       if(FIRST) h1->GetXaxis()->SetTitle("N^{cl}/tracklet");
-      if(FIRST) h1->GetYaxis()->SetTitle("<Entries>");
+      if(FIRST) h1->GetYaxis()->SetTitle("Prob. [%]");
       h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
       legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
       FIRST = kFALSE;
     }
     if(FIRST) break;
     legNClus->Draw();
-    gPad->SetLogy();
+    gPad->SetLogy(0);
     gPad->SetLogx(0);
     gPad->SetGridy();
     gPad->SetGridx();
@@ -890,21 +892,24 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig)
     legNClus->SetHeader("Particle Species");
     for(Int_t is=0; is<AliPID::kSPECIES; is++){
       Int_t bin = FindBin(is, 2.);
-      h1 = h2->ProjectionY("py", bin, bin);
+      h1 = h2->ProjectionY("pyNTracklets", bin, bin);
       if(!h1->GetEntries()) continue;
-      h1->Scale(1./h1->Integral());
+      h1->Scale(100./h1->Integral());
       //h1->SetMarkerStyle(24);
       //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is));
       h1->SetLineColor(AliTRDCalPID::GetPartColor(is));
-      if(FIRST) h1->GetXaxis()->SetTitle("N^{tl}/track");
-      if(FIRST) h1->GetYaxis()->SetTitle("<Entries>");
+      if(FIRST){ 
+        h1->GetXaxis()->SetTitle("N^{trklt}/track");
+        h1->GetXaxis()->SetRangeUser(1.,6.);
+        h1->GetYaxis()->SetTitle("Prob. [%]");
+      }
       h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec");
       legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l");
       FIRST = kFALSE;
     }
     if(FIRST) break;
     legNClus->Draw();
-    gPad->SetLogy();
+    gPad->SetLogy(0);
     gPad->SetLogx(0);
     gPad->SetGridy();
     gPad->SetGridx();
index 0ff75f0244796b5d8b0a0df83f731ddb10354f64..ef264300f92a24fa707b59c86fad5ed61fb8a4c0 100644 (file)
 #include "TObjArray.h"
 #include "TAxis.h"
 #include "TF1.h"
+#include "TLegend.h"
 #include "TGraphErrors.h"
 #include "TLine.h"
 #include "TH2I.h"
@@ -215,6 +216,7 @@ AliTRDclusterResolution::AliTRDclusterResolution(const char *name, const char *t
   ,fZ(0.)
 {
   memset(fR, 0, 4*sizeof(Float_t));
+  memset(fP, 0, 4*sizeof(Float_t));
   // time drift axis
   fAt = new TAxis(kNTB, 0., kNTB*fgkTimeBinLength);
 
@@ -254,7 +256,7 @@ void AliTRDclusterResolution::CreateOutputObjects()
 Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig)
 {
   if(!fResults) return kFALSE;
-  
+  TLegend *leg = 0x0;
   TList *l = 0x0;
   TObjArray *arr = 0x0;
   TH2 *h2 = 0x0;TH1 *h1 = 0x0;
@@ -266,20 +268,31 @@ Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig)
     if(!(gs = (TGraphErrors*)arr->At(1))) break;
     if(!(gp = (TGraphErrors*)arr->At(2))) break;
     gs->Draw("apl");
-    gs->GetHistogram()->GetYaxis()->SetRangeUser(-.01, .6);
+    gs->GetHistogram()->GetYaxis()->SetRangeUser(-50., 700.);
     gs->GetHistogram()->SetXTitle("Q [a.u.]");
-    gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [mm] / freq");
+    gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [#mum] / freq");
     gm->Draw("pl");
     gp->Draw("pl");
     return kTRUE;
   case kCenter:
     if(!(arr = (TObjArray*)fResults->At(kCenter))) break;
-    gPad->Divide(3, 1); l = gPad->GetListOfPrimitives();
-    for(Int_t ipad = 3; ipad--;){
-      if(!(h2 = (TH2F*)arr->At(ipad))) return kFALSE;
-      ((TVirtualPad*)l->At(ipad))->cd();
-      h2->Draw("lego2fb");
+    gPad->Divide(2, 1); l = gPad->GetListOfPrimitives();
+    ((TVirtualPad*)l->At(0))->cd();
+    ((TTree*)arr->At(0))->Draw("y:x>>h(23, 0.1, 2.4, 51, -.51, .51)",
+            "m[0]*(ly==0&&abs(m[0])<1.e-1)", "colz");
+    ((TVirtualPad*)l->At(1))->cd();
+    leg= new TLegend(.7, .7, .9, .95);
+    leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetFillStyle(0);
+    leg->SetHeader("TRD Plane"); 
+    for(Int_t il = 1; il<=AliTRDgeometry::kNlayer; il++){
+      if(!(gm = (TGraphErrors*)arr->At(il))) return kFALSE;
+      gm->Draw(il>1?"pc":"apc"); leg->AddEntry(gm, Form("%d", il-1), "pl");
+      if(il>1) continue;
+      gm->GetHistogram()->SetXTitle("t_{drift} [#mus]");
+      gm->GetHistogram()->SetYTitle("#sigma_{y}(x|cen=0) [#mum]");
+      gm->GetHistogram()->GetYaxis()->SetRangeUser(150., 500.);
     }
+    leg->Draw();
     return kTRUE;
   case kSigm:
     if(!(arr = (TObjArray*)fResults->At(kSigm))) break;
@@ -326,21 +339,16 @@ TObjArray* AliTRDclusterResolution::Histos()
   fContainer = new TObjArray(kNtasks);
   //fContainer->SetOwner(kTRUE);
 
-  TH2I *h2 = 0x0;
   TH3S *h3 = 0x0;
   TObjArray *arr = 0x0;
 
-  fContainer->AddAt(h2 = new TH2I("Charge", "dy=f(q)", 50, 2.2, 7.5, 60, -.3, .3), kQRes);
-  h2->SetXTitle("log(q) [a.u.]");
-  h2->SetYTitle("#Delta y[cm]");
-  h2->SetZTitle("entries");
-
-  fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kCenter);
+  fContainer->AddAt(arr = new TObjArray(2*AliTRDgeometry::kNlayer), kCenter);
   arr->SetName("Center");
   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
-    if(!(h3=(TH3S*)gROOT->FindObject(Form("hc_l%1d", il)))){ 
+    // add resolution plot for each layer
+    if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenResLy%d", il)))){ 
       h3 = new TH3S(
-        Form("hc_l%1d", il), 
+        Form("hCenResLy%d", il), 
         Form(" ly [%d]", il), 
         kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB),   // x
         51, -.51, .51, // y 
@@ -350,8 +358,29 @@ TObjArray* AliTRDclusterResolution::Histos()
       h3->SetZTitle("#Delta y[cm]");
     } h3->Reset();
     arr->AddAt(h3, il);
+    // add Pull plot for each layer
+    if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenPullLy%d", il)))){ 
+      h3 = new TH3S(
+        Form("hCenPullLy%d", il), 
+        Form(" ly [%d]", il), 
+        kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB),   // x
+        51, -.51, .51, // y 
+        60, -4., 4.); // dy
+      h3->SetXTitle("x [#mus]");
+      h3->SetYTitle("y [pw]");
+      h3->SetZTitle("#Delta y/#sigma_{y}");
+    } h3->Reset();
+    arr->AddAt(h3, AliTRDgeometry::kNlayer+il);
   }
 
+  if(!(h3 = (TH3S*)gROOT->FindObject("Charge"))){
+    h3 = new TH3S("Charge", "dy=f(q)", 50, 2.2, 7.5, 60, -.3, .3, 60, -4., 4.);
+    h3->SetXTitle("log(q) [a.u.]");
+    h3->SetYTitle("#Delta y[cm]");
+    h3->SetZTitle("#Delta y/#sigma_{y}");
+  }
+  fContainer->AddAt(h3, kQRes);
+
   fContainer->AddAt(arr = new TObjArray(kNTB), kSigm);
   arr->SetName("Resolution");
   for(Int_t ix=0; ix<kNTB; ix++){
@@ -396,7 +425,6 @@ void AliTRDclusterResolution::Exec(Option_t *)
 
   Int_t det, t;
   Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
-  TH2I *h2 = 0x0;
   TH3S *h3 = 0x0;
 
   // define limits around ExB for which x contribution is negligible
@@ -418,8 +446,8 @@ void AliTRDclusterResolution::Exec(Option_t *)
     // resolution as a function of cluster charge
     // only for phi equal exB 
     if(TMath::Abs(dydx-fExB) < kDtgPhi){
-      h2 = (TH2I*)fContainer->At(kQRes);
-      h2->Fill(TMath::Log(q), dy);
+      h3 = (TH3S*)fContainer->At(kQRes);
+      h3->Fill(TMath::Log(q), dy, dy/TMath::Sqrt(covcl[0]));
     }
 
     // do not use problematic clusters in resolution analysis
@@ -429,12 +457,14 @@ void AliTRDclusterResolution::Exec(Option_t *)
     x = (t+.5)*fgkTimeBinLength; // conservative approach !!
 
     // resolution as a function of y displacement from pad center
-    // only for phi equal exB and clusters close to cathode wire plane
-    // for ideal simulations time bins 4,5 and 6.
-    if(TMath::Abs(dydx-fExB) < kDtgPhi &&
-       TMath::Abs(x-0.675)<0.225){
-      h3 = (TH3S*)arr0->At(AliTRDgeometry::GetLayer(det));
+    // only for phi equal exB
+    if(TMath::Abs(dydx-fExB) < kDtgPhi/* &&
+       TMath::Abs(x-0.675)<0.225*/){
+      Int_t ly(AliTRDgeometry::GetLayer(det));
+      h3 = (TH3S*)arr0->At(ly);
       h3->Fill(x, cli->GetYDisplacement(), dy);
+      h3 = (TH3S*)arr0->At(AliTRDgeometry::kNlayer+ly);
+      h3->Fill(x, cli->GetYDisplacement(), dy/TMath::Sqrt(covcl[0]));
     }
 
     Int_t ix = fAt->FindBin(x);
@@ -478,12 +508,22 @@ Bool_t AliTRDclusterResolution::PostProcess()
     g->SetMarkerStyle(7); 
 
     // pad center dependence
-    fResults->AddAt(t = new TTree("cent", "dy=f(y,x,ly)"), kCenter);
+    fResults->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer+1), kCenter);
+    arr->SetOwner();
+    arr->AddAt(
+    t = new TTree("cent", "dy=f(y,x,ly)"), 0);
     t->Branch("ly", &fLy, "ly/B");
     t->Branch("x", &fX, "x/F");
     t->Branch("y", &fY, "y/F");
     t->Branch("m", &fR[0], "m[2]/F");
     t->Branch("s", &fR[2], "s[2]/F");
+    t->Branch("pm", &fP[0], "pm[2]/F");
+    t->Branch("ps", &fP[2], "ps[2]/F");
+    for(Int_t il=1; il<=AliTRDgeometry::kNlayer; il++){
+      arr->AddAt(g = new TGraphErrors(), il);
+      g->SetLineColor(il); g->SetLineStyle(il);
+      g->SetMarkerColor(il);g->SetMarkerStyle(4); 
+    }
 
 
     fResults->AddAt(t = new TTree("sigm", "dy=f(dw,x,dydx)"), kSigm);
@@ -567,7 +607,7 @@ void AliTRDclusterResolution::SetVisual()
 void AliTRDclusterResolution::ProcessCharge()
 {
   TH2I *h2 = 0x0;
-  if((h2 = (TH2I*)fContainer->At(kQRes))) {
+  if(!(h2 = (TH2I*)fContainer->At(kQRes))) {
     AliWarning("Missing dy=f(Q) histo");
     return;
   }
@@ -593,12 +633,12 @@ void AliTRDclusterResolution::ProcessCharge()
 
     // Fill sy^2 = f(q)
     Int_t ip = gqm->GetN();
-    gqm->SetPoint(ip, q, 10.*f.GetParameter(1));
-    gqm->SetPointError(ip, 0., 10.*f.GetParError(1));
+    gqm->SetPoint(ip, q, 1.e4*f.GetParameter(1));
+    gqm->SetPointError(ip, 0., 1.e4*f.GetParError(1));
 
     // correct sigma for ExB effect
-    gqs->SetPoint(ip, q, 1.e1*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/);
-    gqs->SetPointError(ip, 0., 1.e1*f.GetParError(2)/**f.GetParameter(2)*/);
+    gqs->SetPoint(ip, q, 1.e4*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/);
+    gqs->SetPointError(ip, 0., 1.e4*f.GetParError(2)/**f.GetParameter(2)*/);
 
     // save probability
     n += entries;
@@ -611,7 +651,7 @@ void AliTRDclusterResolution::ProcessCharge()
   for(Int_t ip=gqp->GetN(); ip--;){
     gqp->GetPoint(ip, q, entries);
     entries/=n;
-    gqp->SetPoint(ip, q, entries);
+    gqp->SetPoint(ip, q, 1.e3*entries);
     gqs->GetPoint(ip, q, sy);
     sm += entries*sy;
   }
@@ -619,61 +659,108 @@ void AliTRDclusterResolution::ProcessCharge()
   // error parametrization s(q) = <sy> + b(1/q-1/q0)
   TF1 fq("fq", "[0] + [1]/x", 20., 250.);
   gqs->Fit(&fq);
-  printf("sy(Q) :: sm[%f] b[%f] 1/q0[%f]\n", sm, fq.GetParameter(1), (sm-fq.GetParameter(0))/fq.GetParameter(1));
+  //printf("sm=%f [0]=%f [1]=%f\n", 1.e-4*sm, fq.GetParameter(0), fq.GetParameter(1));
+  printf("  const Float_t sq0inv = %f; // [1/q0]\n", (sm-fq.GetParameter(0))/fq.GetParameter(1));
+  printf("  const Float_t sqb    = %f; // [cm]\n", 1.e-4*fq.GetParameter(1));
 }
 
 //_______________________________________________________
 void AliTRDclusterResolution::ProcessCenterPad()
 {
+// Resolution as a function of y displacement from pad center
+// only for phi equal exB and clusters close to cathode wire plane
+// for ideal simulations time bins 4,5 and 6.
+
   TObjArray *arr = (TObjArray*)fContainer->At(kCenter);
   if(!arr) {
     AliWarning("Missing dy=f(y | x, ly) container");
     return;
   }
+  Float_t s[AliTRDgeometry::kNlayer];
   TF1 f("f", "gaus", -.5, .5);
-  TH1D *h1 = 0x0; TH3S *h3=0x0;
-  TTree *t = (TTree*)fResults->At(kCenter);
+  TF1 fp("fp", "gaus", -3.5, 3.5);
+
+  TH1D *h1 = 0x0; TH2F *h2 = 0x0; TH3S *h3r=0x0, *h3p=0x0;
+  TObjArray *arrRes = (TObjArray*)fResults->At(kCenter);
+  TTree *t = (TTree*)arrRes->At(0);
+  TGraphErrors *gs = 0x0;
   TAxis *ax = 0x0;
-  for(Int_t il=0; il<arr->GetEntriesFast(); il++){
-    if(!(h3 = (TH3S*)arr->At(il))) continue;
 
+  printf("  const Float_t lSy[6][24] = {\n      {");
+  const Int_t nl = AliTRDgeometry::kNlayer;
+  for(Int_t il=0; il<nl; il++){
+    if(!(h3r = (TH3S*)arr->At(il))) continue;
+    if(!(h3p = (TH3S*)arr->At(nl+il))) continue;
+    gs = (TGraphErrors*)arrRes->At(il+1);
     fLy = il;
-    for(Int_t ix=1; ix<=h3->GetXaxis()->GetNbins(); ix++){
-      ax = h3->GetXaxis();
-      ax->SetRange(ix, ix);     
+    for(Int_t ix=1; ix<=h3r->GetXaxis()->GetNbins(); ix++){
+      ax = h3r->GetXaxis(); ax->SetRange(ix, ix);
+      ax = h3p->GetXaxis(); ax->SetRange(ix, ix);
       fX  = ax->GetBinCenter(ix);
       //printf("  x[%2d]=%4.2f\n", ix, fX);
-      for(Int_t iy=1; iy<=h3->GetYaxis()->GetNbins(); iy++){
-        ax = h3->GetYaxis();
-        ax->SetRange(iy, iy);
+      for(Int_t iy=1; iy<=h3r->GetYaxis()->GetNbins(); iy++){
+        ax = h3r->GetYaxis(); ax->SetRange(iy, iy);
+        ax = h3p->GetYaxis(); ax->SetRange(iy, iy);
         fY  = ax->GetBinCenter(iy);
         //printf("    y[%2d]=%5.2f\n", iy, fY);
         // finish navigation in the HnSparse
 
-        h1 = (TH1D*)h3->Project3D("z");
+        h1 = (TH1D*)h3r->Project3D("z");
         Int_t entries = (Int_t)h1->Integral();
         if(entries < 50) continue;
         //Adjust(&f, h1);
         h1->Fit(&f, "QN");
-
     
         // Fill sy,my=f(y_w,x,ly)
         fR[0] = f.GetParameter(1); fR[1] = f.GetParError(1);
         fR[2] = f.GetParameter(2); fR[3] = f.GetParError(2);
 
+        h1 = (TH1D*)h3p->Project3D("z");
+        h1->Fit(&fp, "QN");
+        fP[0] = fp.GetParameter(1); fP[1] = fp.GetParError(1);
+        fP[2] = fp.GetParameter(2); fP[3] = fp.GetParError(2);
+
         //printf("ly[%d] x[%3.1f] y[%+5.2f] m[%5.3f] s[%5.3f] \n", fLy, fX, fY, fR[0], fR[2]);
         t->Fill();
+
+
       }
     }
-    if(!fCanvas) continue;
+    t->Draw("y:x>>h(24, 0., 2.4, 51, -.51, .51)",
+            Form("s[0]*(ly==%d&&abs(m[0])<1.e-1)", fLy),
+            "goff");
+    h2=(TH2F*)gROOT->FindObject("h");
+    f.FixParameter(1, 0.);
+    Int_t n = h2->GetXaxis()->GetNbins(); s[il]=0.;
+    printf("    {");
+    for(Int_t ix=1; ix<=n; ix++){
+      ax = h2->GetXaxis();
+      fX  = ax->GetBinCenter(ix);
+      h1 = h2->ProjectionY("hCenPy", ix, ix);
+      //if((Int_t)h1->Integral() < 1.e-10) continue; 
+      h1->Fit(&f, "QN");
+      s[il]+=f.GetParameter(2);
+      printf("%6.4f,%s", f.GetParameter(0), ix%6?" ":"\n     ");
+      Int_t jx = gs->GetN();
+      gs->SetPoint(jx, fX, 1.e4*f.GetParameter(0));
+      gs->SetPointError(jx, 0., 0./*f.GetParError(0)*/);
+    }
+    printf("\b},\n");
+    s[il]/=n;
+
+    f.ReleaseParameter(2);
 
-    t->Draw("y:x>>h(23, 0.1, 2.4, 51, -.51, .51)",
-            Form("m[0]*(ly==%d&&abs(m[0])<1.e-1)", fLy),
-            "lego2fb");
+
+    if(!fCanvas) continue;
+    h2->Draw("lego2fb");
     fCanvas->Modified(); fCanvas->Update();
     if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessCenter_ly[%d].gif", fLy));
     else gSystem->Sleep(100);
   }
+  printf("  };\n");
+  printf("  const Float_t lPRF[] = {"
+    "%5.3f, %5.3f, %5.3f, %5.3f, %5.3f, %5.3f};\n",
+    s[0], s[1], s[2], s[3], s[4], s[5]);
 }
 
 //_______________________________________________________
index caaa7e7d3ed32be109ec01e4dde9bb97fddba3d3..2c923f3ade62e9c84242c7fb448ed0b269ca1854 100644 (file)
@@ -18,11 +18,11 @@ public:
    ,kN   = kND*kNTB
   };
   enum EResultContainer { // results container type
-    kQRes   = 0
-   ,kCenter = 1
-   ,kSigm   = 2
-   ,kMean   = 3
-   ,kNtasks = 4
+    kCenter = 0   // cluster2center pad calibration
+   ,kQRes   = 1   // resolution on charge dependence
+   ,kSigm   = 2   // sigma cluster as func of x and z
+   ,kMean   = 3   // shift cluster as func of x and z
+   ,kNtasks = 4   // total number os sub tasks
   };
   enum ECheckBits { // force setting the ExB
     kSaveAs    = BIT(22)
@@ -81,9 +81,10 @@ private:
   Float_t    fX;       //! local drift length 
   Float_t    fY;       //! local rphi offset 
   Float_t    fZ;       //! local anode wire offset 
-  Float_t    fR[4];    //! val/error result list
+  Float_t    fR[4];    //! mean/sgm resolution
+  Float_t    fP[4];    //! mean/sgm pulls
   
-  ClassDef(AliTRDclusterResolution, 2)  // cluster resolution
+  ClassDef(AliTRDclusterResolution, 3)  // cluster resolution
 };
 
 //___________________________________________________
index 656730b85cb1c46ee24d8f6e4b2d3d8b706a3398..fa3378a5a7f5ee305ab589e99bd9de3dd75fef2b 100644 (file)
@@ -45,7 +45,7 @@ public:
   Bool_t         HasPostProcess() const {return TestBit(kPostProcess);};
   virtual TObjArray* Histos() {return fContainer;}
 
-  virtual Bool_t Load(const Char_t *filename);
+  virtual Bool_t Load(const Char_t *filename = "TRD.Performance.root");
   virtual Bool_t Save(TObjArray *res);
   virtual Bool_t PostProcess();
   virtual Bool_t PutTrendValue(Char_t *name, Double_t val, Double_t err);
index f9420244ec6f3c9de63d75397a4117477ed31dcb..232344f9ff8bbb0703ef7c49c8849d90f1bf0a90 100644 (file)
@@ -312,10 +312,11 @@ TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
       yc -= tilt*(zc-zt); // tilt correction
       dy = yt - yc;
 
-      Float_t sx2 = dydx*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2;
-      Float_t sy2 = c->GetSY(c->GetLocalTimeBin()); sy2*=sy2;
+      //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));
+      ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2));
   
       if(fDebugLevel>=1){
         // Get z-position with respect to anode wire
@@ -495,10 +496,11 @@ TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
   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/TMath::Sqrt(1.+dydx0*dydx0); PARMC[3]=dzdx0;
+  PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
   PARMC[4]=1./pt0;
 
 //   TMatrixDSymEigen eigen(COV);
@@ -525,14 +527,13 @@ TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
   // 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, 1.-PARMC[4]/PAR[4], is);
+    ((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.+ dzdx0*dzdx0)*pt0, p;
+    Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
     Float_t sp;
     p = tracklet->GetMomentum(&sp);
-    //printf("p[%f] dp[%f] s[%d]\n", p0, 1.-p/p0, is);
-    ((TH3S*)arr->At(10))->Fill(p0, 1.-p/p0, is);
+    ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., is);
     ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/sp, is);
     break;
   }
@@ -563,7 +564,7 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
   UChar_t s;
   Int_t pdg = fMC->GetPDG(), det=-1;
   Int_t label = fMC->GetLabel();
-  Double_t xAnode, x, y, z, pt, dydx, dzdx;
+  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]*/;
 
@@ -587,6 +588,7 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
       << "\n";
   }
 
+  AliTRDReconstructor rec;
   AliTRDseedV1 *fTracklet = 0x0;  
   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
     if(!(fTracklet = fTrack->GetTracklet(ily)))/* ||
@@ -613,18 +615,19 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
         << "dzdx="    << dzdx0
         << "\n";
     }
-    Float_t yt = y0 - dx*dydx0;
-    Float_t zt = z0 - dx*dzdx0;
+    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;
-    y  = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
-    dy = yt - y;
-    z  = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
-    dz = zt - z;
     dydx = fTracklet->GetYref(1);
-    dzdx = fTracklet->GetTgl();
+    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);
 
@@ -638,20 +641,22 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
     // phi resolution/ snp pulls
     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
-    //TODO ((TH2I*)arr->At(5))->Fill(dydx0, );
+    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 dtgl = (dzdx - dzdx0)/(1.- dzdx*dzdx0);
-    ((TH2I*)arr->At(6))->Fill(dzdx0, 
+    Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
+              dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
+    ((TH2I*)arr->At(6))->Fill(dzdl0, 
     TMath::ATan(dtgl));
-    //TODO ((TH2I*)arr->At(7))->Fill(dydx0, );
+    ((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, 1.-pt/pt0, is);
-      //TODO ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (pt0/pt-1.)/TMath::Sqrt(covR[4]), is);
-      Double_t p0 = TMath::Sqrt(1.+ dzdx0*dzdx0)*pt0,
-               p  = TMath::Sqrt(1.+ dzdx*dzdx)*pt;
-     ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, 1-p/p0, is);
+      ((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;
     }
 
@@ -673,30 +678,31 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
     AliTRDseedV1 tt(*fTracklet);
     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*dydx0;
-    zt = z0 - dx*dzdx0;
+    ymc = y0 - dx*dydx0;
+    zmc = z0 - dx*dzdx0;
     Bool_t rc = tt.IsRowCross(); 
     
     // add tracklet residuals for y and dydx
     arr = (TObjArray*)fContainer->At(kMCtracklet);
     if(!rc){
-      dy    = yt-y;
+      dy    = y-ymc;
 
       Float_t dphi  = (dydx - dydx0);
-      dphi /= 1.- 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 = zt-z;
-      ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
-      if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(tt.GetS2Z()));
+      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
@@ -725,37 +731,26 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
     fTracklet->ResetClusterIter(kFALSE);
     while((c = fTracklet->PrevCluster())){
       Float_t  q = TMath::Abs(c->GetQ());
-      x = c->GetX(); y = c->GetY();
-//       Int_t col = c->GetPadCol();
-//       Int_t row = c->GetPadRow();
-//       Double_t cw = pp->GetColSize(col);
-//       Double_t y0 = pp->GetColPos(col) + 0.5 * cw;
-//       Double_t s2 = AliTRDcalibDB::Instance()->GetPRFWidth(det, col, row); s2 *= s2; s2 -= - 1.5e-1;
-//       y = c->GetYloc(y0, s2, cw); y-=(xAnode-x)*exb;
-
-      z = c->GetZ();
+      x = c->GetX(); y = c->GetY();z = c->GetZ();
       dx = x0 - x; 
-      yt = y0 - dx*dydx0;
-      zt = z0 - dx*dzdx0;
-      dy = yt - (y - tilt*(z-zt));
-      Float_t sx2 = dydx0*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2;
-      Float_t sy2 = c->GetSY(c->GetLocalTimeBin()); sy2*=sy2;
-
+      ymc= y0 - dx*dydx0;
+      zmc= z0 - dx*dzdx0;
+      dy = ymc - (y - tilt*(z-zmc));
       // Fill Histograms
       if(q>20. && q<250.){ 
         ((TH2I*)arr->At(0))->Fill(dydx0, dy);
-        ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(sx2+sy2));
+        ((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, dydx0, dzdx0);
+      clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
       clInfo->SetResolution(dy);
       clInfo->SetAnisochronity(d);
       clInfo->SetDriftLength(((c->GetPadTime()+0.5)*.1)*1.5);
@@ -891,7 +886,7 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
     return kTRUE;
   case 13: //kMCtrackTRD [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] =2000.;
+    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;
@@ -899,22 +894,22 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
     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]=-2.; xy[2]=.2; xy[3] =5.;
-    //((TVirtualPad*)l->At(0))->cd();
+    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], kMCtrack, 3)) 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]=-50.; xy[2]=1.; xy[3] =70.;
-    //((TVirtualPad*)l->At(0))->cd();
+    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], kMCtrack, 3)) 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.;
@@ -925,7 +920,16 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
       if(!GetGraphTrack(&xy[0], 8, il)) break;
     }
     return kTRUE;
-  case 17: //kMCtrackTRD [p]
+  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++){
@@ -934,7 +938,7 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
       if(!GetGraphTrack(&xy[0], 10, il)) break;
     }
     return kTRUE;
-  case 18: // kMCtrackTPC [y]
+  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();
@@ -943,7 +947,7 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
     ((TVirtualPad*)l->At(1))->cd();
     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break;
     return kTRUE;
-  case 19: // kMCtrackTPC [z]
+  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();
@@ -952,25 +956,25 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
     ((TVirtualPad*)l->At(1))->cd();
     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break;
     return kTRUE;
-  case 20: // kMCtrackTPC [phi|snp]
+  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] = 2.5;
+    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 21: // kMCtrackTPC [theta|tgl]
+  case 22: // kMCtrackTPC [theta|tgl]
     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
-    xy[0]=-1.; xy[1]=-10.5; xy[2]=1.; xy[3] =20.5;
+    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] = 2.5;
+    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 22: // kMCtrackTPC [pt]
+  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();
@@ -979,7 +983,7 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
     ((TVirtualPad*)l->At(1))->cd();
     if(!GetGraphTrackTPC(xy, 9)) break;
     return kTRUE;
-  case 23: // kMCtrackTPC [p]
+  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();
@@ -990,10 +994,10 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
     pad->SetMargin(0.12, 0.12, 0.1, 0.04);
     if(!GetGraphTrackTPC(xy, 11)) break;
     return kTRUE;
-  case 24:  // kMCtrackTOF [z]
+  case 25:  // kMCtrackTOF [z]
     return kTRUE;
   }
-  AliInfo(Form("Reference plot [%d] missing result", ifig));
+  AliWarning(Form("Reference plot [%d] missing result", ifig));
   return kFALSE;
 }
 
@@ -1039,7 +1043,7 @@ Bool_t AliTRDresolution::PostProcess()
             gm->SetNameTitle(Form("m_%d%02d%d", ig, ic, is), "");
           }
           continue;
-        } else if(ig==kMCtrackTRD&&(ic==8||ic==10)){ // TRD momentum plot
+        } 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);
@@ -1139,13 +1143,13 @@ Bool_t AliTRDresolution::PostProcess()
   Process2D(kMCtrackTRD, 2, &fg, 1.e4);   // z
   Process2D(kMCtrackTRD, 3, &fg);         // z PULL
   Process2D(kMCtrackTRD, 4, &fg, 1.e3);   // phi
-  //Process(kMCtrack, 5, &fg);         // snp PULL
+  Process2D(kMCtrackTRD, 5, &fg);         // snp PULL
   Process2D(kMCtrackTRD, 6, &fg, 1.e3);   // theta
-  //Process(kMCtrack, 7, &fg);         // tgl PULL
+  Process2D(kMCtrackTRD, 7, &fg);         // tgl PULL
   Process4D(kMCtrackTRD, 8, &fg, 1.e2);   // pt resolution
-  //Process4D(kMCtrack, 9, &fg);         // 1/pt pulls
+  Process4D(kMCtrackTRD, 9, &fg);         // 1/pt pulls
   Process4D(kMCtrackTRD, 10, &fg, 1.e2);  // p resolution
-  fNRefFigures = 17;
+  fNRefFigures = 18;
 
   // TRACK TPC RESOLUTION/PULLS
   Process2D(kMCtrackTPC, 0, &fg, 1.e4);// y resolution
@@ -1160,12 +1164,12 @@ Bool_t AliTRDresolution::PostProcess()
   Process3D(kMCtrackTPC, 9, &fg);      // 1/pt pulls
   Process3D(kMCtrackTPC, 10, &fg, 1.e2);// p resolution
   Process3D(kMCtrackTPC, 11, &fg);      // p pulls
-  fNRefFigures = 23;
+  fNRefFigures = 24;
 
   // TRACK HMPID RESOLUTION/PULLS
   Process2D(kMCtrackTOF, 0, &fg, 1.e4); // z towards TOF
   Process2D(kMCtrackTOF, 1, &fg);       // z towards TOF
-  fNRefFigures = 24;
+  fNRefFigures = 25;
 
   return kTRUE;
 }
@@ -1439,7 +1443,7 @@ TObjArray* AliTRDresolution::Histos()
   arr->AddAt(h, 3);
   // Kalman track SNP
   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNP"))){
-    h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -.02, .02);
+    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");
@@ -1455,7 +1459,7 @@ TObjArray* AliTRDresolution::Histos()
   arr->AddAt(h, 5);
   // Kalman track TGL
   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGL"))){
-    h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -.1, .1);
+    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");
@@ -1476,7 +1480,7 @@ TObjArray* AliTRDresolution::Histos()
   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, -.15, .15, n, -.5, n-.5);
+      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");
@@ -1493,13 +1497,14 @@ TObjArray* AliTRDresolution::Histos()
       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, -.25, .25, n, -.5, n-.5);
+      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");
@@ -1544,7 +1549,7 @@ TObjArray* AliTRDresolution::Histos()
   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, -.02, .02);
+    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");
@@ -1560,7 +1565,7 @@ TObjArray* AliTRDresolution::Histos()
   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, -.1, .1);
+    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#theta [rad]");
     h->GetZaxis()->SetTitle("entries");
@@ -1576,7 +1581,7 @@ TObjArray* AliTRDresolution::Histos()
   arr->AddAt(h, 7);
   // Kalman track Pt resolution
   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPt"))){
-    h3 = new TH3S("hMcTrkTPCPt", "Track[TPC] Pt Resolution", 80, 0., 20., 150, -.15, .15, n, -.5, n-.5);
+    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");
@@ -1592,7 +1597,7 @@ TObjArray* AliTRDresolution::Histos()
   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, -.25, .25, n, -.5, n-.5);
+    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");
index 93af9963ae5440c9c0b16dc8a3432493fb0ff92b..dbd8d2f2931125908deeafb8b33769ae3892f32c 100644 (file)
@@ -35,7 +35,7 @@ AliTRDclusterInfo::AliTRDclusterInfo()
 }
 
 //_________________________________________________
-void AliTRDclusterInfo::SetCluster(const AliTRDcluster *c, Float_t *cov)
+void AliTRDclusterInfo::SetCluster(const AliTRDcluster *c)
 {
   if(!c) return;
   fDet = c->GetDetector();
@@ -45,8 +45,9 @@ void AliTRDclusterInfo::SetCluster(const AliTRDcluster *c, Float_t *cov)
   fQ   = TMath::Abs(c->GetQ());
   fLocalTime = c->GetLocalTimeBin();
   fYd  = c->GetCenter();
-
-  if(cov) memcpy(cov, fCovCl, 3*sizeof(Float_t));
+  fCovCl[0] = c->GetSigmaY2();
+  fCovCl[1] = 0.;
+  fCovCl[2] = c->GetSigmaZ2();
 }
 
 //_________________________________________________
index 8b1db24ed1f07b5a0ae63b6033c3361e4181fbc9..b946fa38e88ba264a59136da29863f8dbe527020 100644 (file)
@@ -34,7 +34,7 @@ public:
   void      Print(Option_t *opt="") const;
 
   void      SetAnisochronity(Float_t d) {fD = d;}
-  void      SetCluster(const AliTRDcluster *c, Float_t *cov=0x0);
+  void      SetCluster(const AliTRDcluster *c);
   void      SetMC(Int_t pdg, Int_t label){
       fPdg  = pdg;
       fLbl  = label;}
index 59e357036ea9589c2ef2c1532e607cef087cb6b9..a75a336d677fc1f664433e556e80f2ffbd70eefd 100644 (file)
@@ -68,7 +68,6 @@
 
 Char_t *libs[] = {"libProofPlayer.so", "libANALYSIS.so", "libTRDqaRec.so"};
 // define setup
-gStyle->SetOptStat(0);
 Bool_t mc      = kTRUE;
 Bool_t friends = kTRUE;
 TCanvas *c = 0x0;
@@ -93,6 +92,8 @@ void makeResults(Char_t *opt = "ALL", const Char_t *files=0x0, Bool_t kGRID=kFAL
     return;
   }
 
+  gStyle->SetOptStat(0);
+  gStyle->SetOptFit(0);
   if(files) mergeProd("TRD.Performance.root", files);
   Int_t fSteerTask = ParseOptions(opt);