d72a0e1b27ce47593e6a85f000d860263591aa01
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDresolution.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercialf purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 /* $Id: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  TRD tracking resolution                                               //
21 //
22 // The class performs resolution and residual studies 
23 // of the TRD tracks for the following quantities :
24 //   - spatial position (y, [z])
25 //   - angular (phi) tracklet
26 //   - momentum at the track level
27 // 
28 // The class has to be used for regular detector performance checks using the official macros:
29 //   - $ALICE_ROOT/TRD/qaRec/run.C
30 //   - $ALICE_ROOT/TRD/qaRec/makeResults.C
31 // 
32 // For stand alone usage please refer to the following example: 
33 // {  
34 //   gSystem->Load("libANALYSIS.so");
35 //   gSystem->Load("libTRDqaRec.so");
36 //   AliTRDresolution *res = new AliTRDresolution();
37 //   //res->SetMCdata();
38 //   //res->SetVerbose();
39 //   //res->SetVisual();
40 //   res->Load();
41 //   if(!res->PostProcess()) return;
42 //   res->GetRefFigure(0);
43 // }  
44 //
45 //  Authors:                                                              //
46 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
47 //    Markus Fasel <M.Fasel@gsi.de>                                       //
48 //                                                                        //
49 ////////////////////////////////////////////////////////////////////////////
50
51 #include <TSystem.h>
52 #include <TStyle.h>
53 #include <TROOT.h>
54 #include <TObjArray.h>
55 #include <TH3.h>
56 #include <TH2.h>
57 #include <TH1.h>
58 #include <THnSparse.h>
59 #include <TF1.h>
60 #include <TCanvas.h>
61 #include <TGaxis.h>
62 #include <TBox.h>
63 #include <TLegend.h>
64 #include <TGraphErrors.h>
65 #include <TGraphAsymmErrors.h>
66 #include <TLinearFitter.h>
67 #include <TMath.h>
68 #include <TMatrixT.h>
69 #include <TVectorT.h>
70 #include <TTreeStream.h>
71 #include <TGeoManager.h>
72 #include <TDatabasePDG.h>
73
74 #include "AliPID.h"
75 #include "AliLog.h"
76 #include "AliESDtrack.h"
77 #include "AliMathBase.h"
78 #include "AliTrackPointArray.h"
79
80 #include "AliTRDresolution.h"
81 #include "AliTRDgeometry.h"
82 #include "AliTRDtransform.h"
83 #include "AliTRDpadPlane.h"
84 #include "AliTRDcluster.h"
85 #include "AliTRDseedV1.h"
86 #include "AliTRDtrackV1.h"
87 #include "AliTRDReconstructor.h"
88 #include "AliTRDrecoParam.h"
89 #include "AliTRDpidUtil.h"
90 #include "AliTRDinfoGen.h"
91
92 #include "info/AliTRDclusterInfo.h"
93
94 ClassImp(AliTRDresolution)
95
96 Int_t const   AliTRDresolution::fgkNbins[kNdim] = {
97   Int_t(kNbunchCross)/*bc*/,
98   180/*phi*/,
99   50/*eta*/,
100   Int_t(kNcharge)*AliPID::kSPECIES+1/*chg*species*/,
101   50/*dy*/,
102   50/*dz*/,
103   40/*dphi*/
104 //  kNpt/*pt*/,
105 };  //! no of bins/projection
106 Double_t const AliTRDresolution::fgkMin[kNdim] = {
107   -0.5,
108   -TMath::Pi(),
109   -1.,
110   -AliPID::kSPECIES-0.5,
111   -1.5,
112   -2.5,
113   -10.
114 //  -0.5,
115 };    //! low limits for projections
116 Double_t const AliTRDresolution::fgkMax[kNdim] = {
117   Int_t(kNbunchCross)-0.5,
118   TMath::Pi(),
119   1.,
120   AliPID::kSPECIES+0.5,
121   1.5,
122   2.5,
123   10.
124 //  kNpt-0.5,
125 };    //! high limits for projections
126 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
127   "bunch cross",
128   "#phi [rad]",
129   "#eta",
130   "chg*spec*rc",
131   "#Deltay [cm]",
132   "#Deltaz [cm]",
133   "#Delta#phi [deg]"
134 //  "bin_p_{t}",
135 };  //! title of projection
136
137 UChar_t const AliTRDresolution::fgNproj[kNclasses] = {
138   6, 5, 5, 5,
139   2, 5, 11, 11, 11
140 };
141 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
142     "Cluster2Track"
143     ,"Tracklet2Track"
144     ,"Tracklet2TRDin" 
145     ,"Tracklet2TRDout" 
146     ,"Cluster2MC"
147     ,"Tracklet2MC"
148     ,"TRDin2MC"
149     ,"TRDout2MC"
150     ,"TRD2MC"
151 };
152 Float_t AliTRDresolution::fgPtBin[kNpt+1] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
153
154 //________________________________________________________
155 AliTRDresolution::AliTRDresolution()
156   :AliTRDrecoTask()
157   ,fIdxPlot(0)
158   ,fIdxFrame(0)
159   ,fPtThreshold(1.)
160   ,fDyRange(0.75)
161   ,fProj(NULL)
162   ,fDBPDG(NULL)
163   ,fCl(NULL)
164   ,fMCcl(NULL)
165 {
166   //
167   // Default constructor
168   //
169   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
170   MakePtSegmentation();
171 }
172
173 //________________________________________________________
174 AliTRDresolution::AliTRDresolution(char* name)
175   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
176   ,fIdxPlot(0)
177   ,fIdxFrame(0)
178   ,fPtThreshold(1.)
179   ,fDyRange(0.75)
180   ,fProj(NULL)
181   ,fDBPDG(NULL)
182   ,fCl(NULL)
183   ,fMCcl(NULL)
184 {
185   //
186   // Default constructor
187   //
188
189   InitFunctorList();
190   MakePtSegmentation();
191
192   DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
193   DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
194 }
195
196 //________________________________________________________
197 AliTRDresolution::~AliTRDresolution()
198 {
199   //
200   // Destructor
201   //
202
203   if(fProj){fProj->Delete(); delete fProj;}
204   if(fCl){fCl->Delete(); delete fCl;}
205   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
206 }
207
208
209 //________________________________________________________
210 void AliTRDresolution::UserCreateOutputObjects()
211 {
212   // spatial resolution
213
214   AliTRDrecoTask::UserCreateOutputObjects();
215   InitExchangeContainers();
216 }
217
218 //________________________________________________________
219 void AliTRDresolution::InitExchangeContainers()
220 {
221 // Init containers for subsequent tasks (AliTRDclusterResolution)
222
223   fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
224   fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
225   PostData(kClToTrk, fCl);
226   PostData(kClToMC, fMCcl);
227 }
228
229 //________________________________________________________
230 void AliTRDresolution::UserExec(Option_t *opt)
231 {
232   //
233   // Execution part
234   //
235
236   fCl->Delete();
237   fMCcl->Delete();
238   AliTRDrecoTask::UserExec(opt);
239 }
240
241 //________________________________________________________
242 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
243 {
244 // Helper function to calculate pulls in the yz plane 
245 // using proper tilt rotation
246 // Uses functionality defined by AliTRDseedV1.
247
248   return kTRUE;
249 /*
250   Double_t t2(tilt*tilt);
251   // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
252
253   // rotate along pad
254   Double_t cc[3];
255   cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2]; 
256   cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
257   cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
258   // do sqrt
259   Double_t sqr[3]={0., 0., 0.}; 
260   if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
261   Double_t invsqr[3]={0., 0., 0.}; 
262   if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
263   Double_t tmp(dyz[0]);
264   dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
265   dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
266   return kTRUE;
267 */
268 }
269
270 //________________________________________________________
271 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
272 {
273   //
274   // Plot the cluster distributions
275   //
276
277   if(track) fkTrack = track;
278   if(!fkTrack){
279     AliDebug(4, "No Track defined.");
280     return NULL;
281   }
282   if(fkESD->GetTOFbc() > 1){
283     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
284     return NULL;
285   }
286   if(fPt<fPtThreshold){
287     AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
288     return NULL;
289   }
290   THnSparse *H(NULL);
291   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
292     AliWarning("No output container defined.");
293     return NULL;
294   }
295
296   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
297   Double_t val[kNdim]; Float_t exb, vd, t0, s2, dl, dt, corr;
298   TObjArray     *clInfoArr(NULL);
299   AliTRDseedV1  *fTracklet(NULL);
300   AliTRDcluster *c(NULL)/*, *cc(NULL)*/;
301   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
302     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
303     if(!fTracklet->IsOK()) continue;
304     fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
305     val[kBC]  = ily;
306     val[kPhi] = fPhi;
307     val[kEta] = fEta;
308     //val[kPt]  = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb))*TMath::RadToDeg();
309     corr = fkTrack->Charge()/TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1))/vd;
310     fTracklet->ResetClusterIter(kTRUE);
311     while((c = fTracklet->NextCluster())){
312       Float_t xc(c->GetX());  //Int_t tb(c->GetLocalTimeBin());
313
314       val[kYrez] = c->GetY()-fTracklet->GetYat(xc);
315       //val[kZrez] = fTracklet->GetX0()-xc;
316       //val[kPrez] = 0.; Int_t ic(0);
317       //if((cc = fTracklet->GetClusters(tb-1))) {val[kPrez] += cc->GetQ(); ic++;}
318       //if((cc = fTracklet->GetClusters(tb-2))) {val[kPrez] += cc->GetQ(); ic++;}
319       //if(ic) val[kPrez] /= (ic*c->GetQ());
320       val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:(c->GetQ()*corr);
321       H->Fill(val);
322 /*      // tilt rotation of covariance for clusters
323       Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
324       cov[0] = (sy2+t2*sz2)*corr;
325       cov[1] = tilt*(sz2 - sy2)*corr;
326       cov[2] = (t2*sy2+sz2)*corr;
327       // sum with track covariance
328       cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
329       Double_t dyz[2]= {dy[1], dz[1]};
330       Pulls(dyz, cov, tilt);*/
331   
332       // Get z-position with respect to anode wire
333       Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
334               zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
335       Int_t istk = geo->GetStack(c->GetDetector());
336       AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
337       Float_t row0 = pp->GetRow0();
338       Float_t d  = row0 - zt + pp->GetAnodeWireOffset();
339       d -= ((Int_t)(2 * d)) / 2.0;
340       if (d > 0.25) d  = 0.5 - d;
341
342       AliTRDclusterInfo *clInfo(NULL);
343       clInfo = new AliTRDclusterInfo;
344       clInfo->SetCluster(c);
345       //Float_t covF[] = {cov[0], cov[1], cov[2]};
346       clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
347       clInfo->SetResolution(val[kYrez]);
348       clInfo->SetAnisochronity(d);
349       clInfo->SetDriftLength(val[kZrez]);
350       clInfo->SetTilt(fTracklet->GetTilt());
351       if(fCl) fCl->Add(clInfo);
352       else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
353
354       if(DebugLevel()>=1){
355         if(!clInfoArr){ 
356           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
357           clInfoArr->SetOwner(kFALSE);
358         }
359         clInfoArr->Add(clInfo);
360       }
361     }
362     if(DebugLevel()>=1 && clInfoArr){
363       ULong_t status = fkESD->GetStatus();
364       (*DebugStream()) << "cluster"
365         <<"status="  << status
366         <<"clInfo.=" << clInfoArr
367         << "\n";
368       clInfoArr->Clear();
369     }
370   }
371   if(clInfoArr) delete clInfoArr;
372
373   return NULL;//H->Projection(kEta, kPhi);
374 }
375
376
377 //________________________________________________________
378 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
379 {
380 // Plot normalized residuals for tracklets to track. 
381 // 
382 // We start from the result that if X=N(|m|, |Cov|)
383 // BEGIN_LATEX
384 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
385 // END_LATEX
386 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial 
387 // reference position. 
388   if(track) fkTrack = track;
389   if(!fkTrack){
390     AliDebug(4, "No Track defined.");
391     return NULL;
392   }
393   if(fkESD->GetTOFbc()>1){
394     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
395     return NULL;
396   }
397   THnSparse *H(NULL);
398   if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
399     AliWarning("No output container defined.");
400     return NULL;
401   }
402 //  return NULL;
403   Double_t val[kNdim+1];
404   AliTRDseedV1 *fTracklet(NULL);
405   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
406     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
407     if(!fTracklet->IsOK()) continue;
408     val [kBC] = il; 
409     val[kPhi] = fPhi;
410     val[kEta] = fEta;
411     val[kSpeciesChgRC]= fSpecies;
412     //val[kPt]  = GetPtBin(fTracklet->GetPt());
413     Double_t dyt(fTracklet->GetYref(0) - fTracklet->GetYfit(0)),
414              dzt(fTracklet->GetZref(0) - fTracklet->GetZfit(0)),
415              dydx(fTracklet->GetYfit(1)),
416              tilt(fTracklet->GetTilt());
417     // correct for tilt rotation
418     val[kYrez] = dyt - dzt*tilt;
419     val[kZrez] = dzt + dyt*tilt;
420     dydx+= tilt*fTracklet->GetZref(1);
421     val[kPrez] = TMath::ATan((fTracklet->GetYref(1) - dydx)/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
422     if(fTracklet->IsRowCross()){
423       val[kSpeciesChgRC]= 0.;
424       val[kPrez] = fkTrack->Charge(); // may be better defined
425     } else {
426       Float_t exb, vd, t0, s2, dl, dt;
427       fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
428       val[kZrez] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb));
429     }
430     val[kNdim] = fTracklet->GetdQdl();
431     if(DebugLevel()>=1) H->Fill(val);
432
433 //     // compute covariance matrix
434 //     fTracklet->GetCovAt(x, cov);
435 //     fTracklet->GetCovRef(covR);
436 //     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2]; 
437 //     Double_t dyz[2]= {dy[1], dz[1]};
438 //     Pulls(dyz, cov, tilt);
439 //     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
440 //     ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
441
442     if(DebugLevel()>=2){
443       Bool_t rc(fTracklet->IsRowCross());
444       UChar_t err(fTracklet->GetErrorMsg());
445       Double_t x(fTracklet->GetX()),
446                pt(fTracklet->GetPt()),
447                yt(fTracklet->GetYref(0)),
448                zt(fTracklet->GetZref(0)),
449                phi(fTracklet->GetYref(1)),
450                tht(fTracklet->GetZref(1));
451       Int_t ncl(fTracklet->GetN()),
452             det(fTracklet->GetDetector());
453       (*DebugStream()) << "tracklet"
454         <<"pt="  << pt
455         <<"x="   << x
456         <<"yt="  << yt
457         <<"zt="  << zt
458         <<"phi=" << phi
459         <<"tht=" << tht
460         <<"det=" << det
461         <<"n="   << ncl
462         <<"dy0=" << dyt
463         <<"dz0=" << dzt
464         <<"dy="  << val[kYrez]
465         <<"dz="  << val[kZrez]
466         <<"dphi="<< val[kPrez]
467         <<"dQ  ="<< val[kNdim]
468         <<"rc="  << rc
469         <<"err=" << err
470         << "\n";
471     }
472   }
473   return NULL;//H->Projection(kEta, kPhi);
474 }
475
476
477 //________________________________________________________
478 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
479 {
480 // Store resolution/pulls of Kalman before updating with the TRD information 
481 // at the radial position of the first tracklet. The following points are used 
482 // for comparison  
483 //  - the (y,z,snp) of the first TRD tracklet
484 //  - the (y, z, snp, tgl, pt) of the MC track reference
485 // 
486 // Additionally the momentum resolution/pulls are calculated for usage in the 
487 // PID calculation. 
488
489   if(track) fkTrack = track;
490   if(!fkTrack){
491     AliDebug(4, "No Track defined.");
492     return NULL;
493   }
494   // check container
495   THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
496   if(!H){
497     AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
498     return NULL;
499   }
500   // check input track status
501   AliExternalTrackParam *tin(NULL);
502   if(!(tin = fkTrack->GetTrackIn())){
503     AliError("Track did not entered TRD fiducial volume.");
504     return NULL;
505   }
506   // check first tracklet
507   AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
508   if(!fTracklet){
509     AliDebug(3, "No Tracklet in ly[0]. Skip track.");
510     return NULL;
511   }
512   // check radial position
513   Double_t x = tin->GetX();
514   if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
515     AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
516     return NULL;
517   }
518
519   Int_t bc(fkESD->GetTOFbc()%2);
520   const Double_t *parR(tin->GetParameter());
521   Double_t dyt(parR[0] - fTracklet->GetYfit(0)), dzt(parR[1] - fTracklet->GetZfit(0)),
522             phit(fTracklet->GetYfit(1)),
523             tilt(fTracklet->GetTilt());
524
525   // correct for tilt rotation
526   Double_t dy  = dyt - dzt*tilt,
527            dz  = dzt + dyt*tilt;
528   phit       += tilt*parR[3];
529   Double_t dphi = TMath::ASin(parR[2])-TMath::ATan(phit);
530
531   Double_t val[kNdim];
532   val[kBC]          = bc;
533   val[kPhi]         = fPhi;
534   val[kEta]         = fEta;
535   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fSpecies;
536 //  val[kPt]          = GetPtBin(fPt);
537   val[kYrez]        = dy;
538   val[kZrez]        = dz;
539   val[kPrez]        = dphi*TMath::RadToDeg();
540   H->Fill(val);
541
542   if(!HasMCdata()) return NULL; // H->Projection(kEta, kPhi);
543   if(!(H = (THnSparseI*)fContainer->At(kMCtrackIn))) {
544     AliError(Form("Missing container @ %d", Int_t(kMCtrackIn)));
545     return NULL;
546   }
547
548   // get MC info
549   UChar_t s;
550   Float_t pt0, eta, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
551   if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, s)) return NULL;
552   dyt = y0 - fTracklet->GetYfit(0);
553   dzt = z0 - fTracklet->GetZfit(0);
554   phit= fTracklet->GetYfit(1) + tilt*dzdx0;
555   Float_t phi = TMath::ATan2(y0, x0);
556   dy  = dyt - dzt*tilt;
557   dz  = dzt + dyt*tilt;
558   dphi= TMath::ASin(dydx0)-TMath::ATan(phit);
559
560   Int_t pdg = fkMC->GetPDG(),
561         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
562         sign(0);
563   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
564   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
565   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
566
567
568   val[kBC]          = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
569   val[kPhi]         = phi;
570   val[kEta]         = eta;
571   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:sign*(sIdx+1);
572 //  val[kPt]          = GetPtBin(pt0);
573   val[kYrez]        = dy;
574   val[kZrez]        = dz;
575   val[kPrez]        = dphi*TMath::RadToDeg();
576   H->Fill(val);
577
578   return NULL; //H->Projection(kEta, kPhi);
579 }
580
581 //________________________________________________________
582 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
583 {
584 // Store resolution/pulls of Kalman after last update with the TRD information 
585 // at the radial position of the first tracklet. The following points are used 
586 // for comparison  
587 //  - the (y,z,snp) of the first TRD tracklet
588 //  - the (y, z, snp, tgl, pt) of the MC track reference
589 // 
590 // Additionally the momentum resolution/pulls are calculated for usage in the 
591 // PID calculation. 
592
593   if(track) fkTrack = track;
594   return NULL;
595 }
596
597 //________________________________________________________
598 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
599 {
600   //
601   // Plot MC distributions
602   //
603
604   if(!HasMCdata()){ 
605     AliDebug(2, "No MC defined. Results will not be available.");
606     return NULL;
607   }
608   if(track) fkTrack = track;
609   if(!fkTrack){
610     AliDebug(4, "No Track defined.");
611     return NULL;
612   }
613   if(!fContainer){
614     AliWarning("No output container defined.");
615     return NULL;
616   }
617   // retriev track characteristics
618   Int_t pdg = fkMC->GetPDG(),
619         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
620         sign(0),
621         sgm[3],
622         label(fkMC->GetLabel()),
623         fSegmentLevel(0);
624   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
625   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
626   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
627
628   TObjArray *arr(NULL);TH1 *h(NULL);
629   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
630   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
631   UChar_t s;
632   Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
633   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
634   Double_t covR[7]/*, cov[3]*/;
635   
636   if(DebugLevel()>=3){
637     // get first detector
638     Int_t det = -1;
639     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
640       if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
641       det = fTracklet->GetDetector();
642       break;
643     }
644     if(det>=0){
645       TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
646       Double_t m(-1.);
647       m = fkTrack->GetMass();
648       if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
649         (*DebugStream()) << "MCkalman"
650           << "pdg=" << pdg
651           << "det=" << det
652           << "x="   << &X
653           << "y="   << &Y
654           << "z="   << &Z
655           << "dx="  << &dX
656           << "dy="  << &dY
657           << "dz="  << &dZ
658           << "pt="  << &vPt
659           << "dpt=" << &dPt
660           << "bgt=" << &budget
661           << "cov=" << &cCOV
662           << "\n";
663       }
664     }
665   }
666   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
667     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
668        !fTracklet->IsOK())*/ continue;
669
670     sgm[2] = fTracklet->GetDetector();
671     sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
672     sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
673     Double_t tilt(fTracklet->GetTilt())
674             ,t2(tilt*tilt)
675             ,corr(1./(1. + t2))
676             ,cost(TMath::Sqrt(corr));
677     x0  = fTracklet->GetX0();
678     //radial shift with respect to the MC reference (radial position of the pad plane)
679     x= fTracklet->GetX();
680     Bool_t rc(fTracklet->IsRowCross()); Float_t eta;
681     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, s)) continue;
682     xAnode  = fTracklet->GetX0();
683
684     // MC track position at reference radial position
685     dx  = x0 - x;
686     if(DebugLevel()>=4){
687       (*DebugStream()) << "MC"
688         << "det="     << sgm[2]
689         << "pdg="     << pdg
690         << "sgn="     << sign
691         << "pt="      << pt0
692         << "x="       << x0
693         << "y="       << y0
694         << "z="       << z0
695         << "dydx="    << dydx0
696         << "dzdx="    << dzdx0
697         << "\n";
698     }
699     Float_t ymc = y0 - dx*dydx0;
700     Float_t zmc = z0 - dx*dzdx0;
701     //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
702
703     // Kalman position at reference radial position
704     dx = xAnode - x;
705     dydx = fTracklet->GetYref(1);
706     dzdx = fTracklet->GetZref(1);
707     dzdl = fTracklet->GetTgl();
708     y  = fTracklet->GetYref(0) - dx*dydx;
709     dy = y - ymc;
710     z  = fTracklet->GetZref(0) - dx*dzdx;
711     dz = z - zmc;
712     pt = TMath::Abs(fTracklet->GetPt());
713     fTracklet->GetCovRef(covR);
714
715     arr = (TObjArray*)((TObjArray*)fContainer->At(kMCtrack))->At(ily);
716     // y resolution/pulls
717     if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
718     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(covR[0]), dz/TMath::Sqrt(covR[2]));
719     // z resolution/pulls
720     ((TH2S*)arr->At(2))->Fill(dzdx0, dz);
721     ((TH3S*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]), 0);
722     // phi resolution/ snp pulls
723     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
724     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
725     Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
726     ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
727     // theta resolution/ tgl pulls
728     Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
729               dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
730     ((TH2I*)arr->At(6))->Fill(dzdl0, 
731     TMath::ATan(dtgl));
732     ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
733     // pt resolution  \\ 1/pt pulls \\ p resolution for PID
734     Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
735              p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
736     ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
737     ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
738     ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);
739
740     // Fill Debug stream for Kalman track
741     if(DebugLevel()>=4){
742       (*DebugStream()) << "MCtrack"
743         << "pt="      << pt
744         << "x="       << x
745         << "y="       << y
746         << "z="       << z
747         << "dydx="    << dydx
748         << "dzdx="    << dzdx
749         << "s2y="     << covR[0]
750         << "s2z="     << covR[2]
751         << "\n";
752     }
753
754     // recalculate tracklet based on the MC info
755     AliTRDseedV1 tt(*fTracklet);
756     tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
757     tt.SetZref(1, dzdx0); 
758     tt.SetReconstructor(AliTRDinfoGen::Reconstructor());
759     tt.Fit(1);
760     x= tt.GetX();y= tt.GetY();z= tt.GetZ();
761     dydx = tt.GetYfit(1);
762     dx = x0 - x;
763     ymc = y0 - dx*dydx0;
764     zmc = z0 - dx*dzdx0;
765     dy = y-ymc;
766     dz = z-zmc;
767     Float_t dphi  = (dydx - dydx0);
768     dphi /= (1.- dydx*dydx0);
769
770     // add tracklet residuals for y and dydx
771     arr = (TObjArray*)fContainer->At(kMCtracklet);
772
773     if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
774     if(tt.GetS2Y()>0. && tt.GetS2Z()>0.) ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(tt.GetS2Y()), dz/TMath::Sqrt(tt.GetS2Z()));
775     ((TH3S*)arr->At(2))->Fill(dzdl0, dz, rc);
776     if(tt.GetS2Z()>0.) ((TH3S*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()), rc);
777     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
778   
779     // Fill Debug stream for tracklet
780     if(DebugLevel()>=4){
781       Float_t s2y = tt.GetS2Y();
782       Float_t s2z = tt.GetS2Z();
783       (*DebugStream()) << "MCtracklet"
784         << "rc="    << rc
785         << "x="     << x
786         << "y="     << y
787         << "z="     << z
788         << "dydx="  << dydx
789         << "s2y="   << s2y
790         << "s2z="   << s2z
791         << "\n";
792     }
793
794     AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
795     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
796     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
797
798     arr = (TObjArray*)fContainer->At(kMCcluster);
799     AliTRDcluster *c = NULL;
800     tt.ResetClusterIter(kFALSE);
801     while((c = tt.PrevCluster())){
802       Float_t  q = TMath::Abs(c->GetQ());
803       x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()]; y = c->GetY();z = c->GetZ();
804       dx = x0 - x; 
805       ymc= y0 - dx*dydx0;
806       zmc= z0 - dx*dzdx0;
807       dy = cost*(y - ymc - tilt*(z-zmc));
808       dz = cost*(z - zmc + tilt*(y-ymc));
809       
810       // Fill Histograms
811       if(q>20. && q<250. && pt0>fPtThreshold && c->IsInChamber()){ 
812         ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
813         ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(c->GetSigmaY2()), dz/TMath::Sqrt(c->GetSigmaZ2()));
814       }
815
816       // Fill calibration container
817       Float_t d = zr0 - zmc;
818       d -= ((Int_t)(2 * d)) / 2.0;
819       if (d > 0.25) d  = 0.5 - d;
820       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
821       clInfo->SetCluster(c);
822       clInfo->SetMC(pdg, label);
823       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
824       clInfo->SetResolution(dy);
825       clInfo->SetAnisochronity(d);
826       clInfo->SetDriftLength(dx);
827       clInfo->SetTilt(tilt);
828       if(fMCcl) fMCcl->Add(clInfo);
829       else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
830       if(DebugLevel()>=5){ 
831         if(!clInfoArr){ 
832           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
833           clInfoArr->SetOwner(kFALSE);
834         }
835         clInfoArr->Add(clInfo);
836       }
837     }
838     // Fill Debug Tree
839     if(DebugLevel()>=5 && clInfoArr){
840       (*DebugStream()) << "MCcluster"
841         <<"clInfo.=" << clInfoArr
842         << "\n";
843       clInfoArr->Clear();
844     }
845   }
846   if(clInfoArr) delete clInfoArr;
847   return h;
848 }
849
850
851 //__________________________________________________________________________
852 Int_t AliTRDresolution::GetPtBin(Float_t pt)
853 {
854 // Find pt bin according to local pt segmentation
855   Int_t ipt(0);
856   while(ipt<AliTRDresolution::kNpt){
857     if(pt<fgPtBin[ipt]) break;
858     ipt++;
859   }
860   return TMath::Max(0,ipt);
861 }
862
863 //________________________________________________________
864 Float_t AliTRDresolution::GetMeanWithBoundary(TH1 *h, Float_t zm, Float_t zM, Float_t dz) const
865 {
866 // return mean of histogram "h"
867 // if histo is empty returns -infinity
868 // if few entries returns zM+epsilon
869 // if mean less than zm returns zm-epsilon
870
871   Int_t ne(Int_t(h->GetEntries()));
872   if(ne==0) return -1.e+5;
873   else if(ne<20) return zM+0.5*dz;
874   else{
875     Float_t val(h->GetMean());
876     if(val<zm) return zm-0.5*dz;
877     else if(val>zM) return zM-0.5*dz;
878     else return val;
879   }
880 }
881
882 //________________________________________________________
883 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
884 {
885   //
886   // Get the reference figures
887   //
888
889   if(!gPad){
890     AliWarning("Please provide a canvas to draw results.");
891     return kFALSE;
892   }
893 /*  Int_t selection[100], n(0), selStart(0); // 
894   Int_t ly0(0), dly(5);
895   TList *l(NULL); TVirtualPad *pad(NULL); */
896   switch(ifig){
897   case 0:
898     break;
899   }
900   AliWarning(Form("Reference plot [%d] missing result", ifig));
901   return kFALSE;
902 }
903
904
905 //________________________________________________________
906 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
907 {
908 // Build pt segments
909   for(Int_t j(0); j<=kNpt; j++){
910     pt0+=(TMath::Exp(j*j*dpt)-1.);
911     fgPtBin[j]=pt0;
912   }
913 }
914
915 //________________________________________________________
916 void AliTRDresolution::MakeSummary()
917 {
918 // Build summary plots
919
920   if(!fProj){
921     AliError("Missing results");
922     return;
923   }  
924   Int_t iSumPlot(0);
925   TVirtualPad *p(NULL); TCanvas *cOut(NULL);
926   TObjArray *arr(NULL);
927
928   // cluster resolution
929   // define palette
930   gStyle->SetPalette(1);
931   cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot++), "Cluster Resolution", 1024, 768);
932   cOut->Divide(3,2, 2.e-3, 2.e-3);
933   arr = (TObjArray*)fProj->At(kCluster);
934   for(Int_t iplot(0); iplot<fgNproj[kCluster]; iplot++){
935     p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
936     ((TH2*)arr->At(iplot))->Draw("colz");
937   }
938   cOut->SaveAs(Form("%s.gif", cOut->GetName()));
939
940
941   // trackIn systematic
942   // define palette
943   Int_t palette[50];
944   for (int i=1;i<49;i++) palette[i] = 50+i; palette[0]=kMagenta; palette[49]=kBlack;
945   gStyle->SetPalette(50, palette);
946   cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot++), "Track IN Resolution", 1024, 768);
947   cOut->Divide(3,2, 2.e-3, 2.e-3);
948   arr = (TObjArray*)fProj->At(kTrackIn); 
949   for(Int_t iplot(0); iplot<fgNproj[kTrackIn]; iplot++){
950     p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
951     ((TH2*)arr->At(iplot))->Draw("colz");
952   }
953   cOut->SaveAs(Form("%s.gif", cOut->GetName()));
954
955   delete cOut;
956   gStyle->SetPalette(1);
957 }
958
959 //________________________________________________________
960 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
961 {
962 // Returns the range of the bulk of data in histogram h2. Removes outliers. 
963 // The "range" vector should be initialized with 2 elements
964 // Option "mod" can be any of
965 //   - 0 : gaussian like distribution 
966 //   - 1 : tailed distribution 
967
968   Int_t nx(h2->GetNbinsX())
969       , ny(h2->GetNbinsY())
970       , n(nx*ny);
971   Double_t *data=new Double_t[n];
972   for(Int_t ix(1), in(0); ix<=nx; ix++){
973     for(Int_t iy(1); iy<=ny; iy++)
974       data[in++] = h2->GetBinContent(ix, iy);
975   }
976   Double_t mean, sigm;
977   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
978
979   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
980   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]); 
981   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
982   TH1S h1("h1SF0", "", 100, range[0], range[1]);
983   h1.FillN(n,data,0);
984   delete [] data;
985  
986   switch(mod){
987   case 0:// gaussian distribution  
988   {
989     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
990     h1.Fit(&fg, "QN");
991     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
992     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
993     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
994     break;
995   }
996   case 1:// tailed distribution  
997   {  
998     Int_t bmax(h1.GetMaximumBin());
999     Int_t jBinMin(1), jBinMax(100);
1000     for(Int_t ibin(bmax); ibin--;){
1001       if(h1.GetBinContent(ibin)<1.){
1002         jBinMin=ibin; break;
1003       }
1004     }
1005     for(Int_t ibin(bmax); ibin++;){
1006       if(h1.GetBinContent(ibin)<1.){
1007         jBinMax=ibin; break;
1008       }
1009     }
1010     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1011     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
1012     break;
1013   }
1014   }
1015
1016   return;
1017 }
1018
1019
1020 //________________________________________________________
1021 Bool_t AliTRDresolution::MakeProjectionCluster()
1022 {
1023 // Analyse cluster
1024   Int_t cidx = kCluster;
1025   if(fProj && fProj->At(cidx)) return kTRUE;
1026   if(!fContainer){
1027     AliError("Missing data container.");
1028     return kFALSE;
1029   }
1030   THnSparse *H(NULL);
1031   if(!(H = (THnSparse*)fContainer->At(cidx))){
1032     AliError(Form("Missing/Wrong data @ %d.", cidx));
1033     return kFALSE;
1034   }
1035
1036   TAxis *aphi(H->GetAxis(kPhi)),
1037         *aeta(H->GetAxis(kEta)),
1038         *as(H->GetAxis(kSpeciesChgRC)),
1039         //*apt(H->GetAxis(kPt)),
1040         *ay(H->GetAxis(kYrez));
1041         //*az(H->GetAxis(kZrez)),
1042         //*ap(H->GetAxis(kPrez));
1043   Int_t neta(aeta->GetNbins()), nphi(aphi->GetNbins()), rcBin(as->GetNbins()/2 + 1);
1044   TH3I *h3[fgNproj[cidx]];
1045   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1046     h3[ily] = new TH3I(Form("h3CL%d", ily), Form("r-#phi residuals @ ly[%d];%s;%s;%s", ily, aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
1047            neta, aeta->GetXmin(), aeta->GetXmax(),
1048            nphi, aphi->GetXmin(), aphi->GetXmax(),
1049            ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1050   }
1051   Int_t  coord[AliTRDresolution::kNdim]; memset(coord, 0, sizeof(Int_t) * AliTRDresolution::kNdim); Double_t v = 0.;
1052   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1053     v = H->GetBinContent(ib, coord);
1054     if(v<1.) continue;
1055     if(coord[kSpeciesChgRC]==rcBin) continue; // row cross
1056     Int_t ily(coord[kBC]-1);
1057     h3[ily]->AddBinContent(h3[ily]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1058   }
1059
1060   TF1 fg("fg", "gaus", ay->GetXmin(), ay->GetXmax());
1061   if(!fProj){
1062     AliInfo("Building array of projections ...");
1063     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1064   }
1065   TObjArray *arr(NULL);
1066   fProj->AddAt(arr = new TObjArray(fgNproj[cidx]), cidx);
1067
1068   TH2F *h2(NULL);
1069   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1070     arr->AddAt(h2 = new TH2F(Form("h2CLs%d", ily),
1071             Form("Cl Resolution Ly[%d];%s;%s;#sigmay [#mum]", ily, aeta->GetTitle(), aphi->GetTitle()),
1072             neta, aeta->GetXmin(), aeta->GetXmax(),
1073             nphi, aphi->GetXmin(), aphi->GetXmax()), ily);
1074     TAxis *ax = h2->GetZaxis();
1075     ax->CenterTitle(); ax->SetTitleOffset(1.3);
1076     ax->SetRangeUser(250, 500);
1077     for(Int_t iphi(1); iphi<=nphi; iphi++){
1078       for(Int_t ieta(1); ieta<=neta; ieta++){
1079         TH1 *h = h3[ily]->ProjectionZ(Form("h1CLs%d", ily), ieta, ieta, iphi, iphi);
1080         Int_t ne(Int_t(h->GetEntries()));
1081         if(ne<100.) h2->SetBinContent(ieta, iphi, -999);
1082         else {
1083           fg.SetParameter(0, ne);fg.SetParameter(1, 0.);fg.SetParameter(2, 0.05);
1084           h->Fit(&fg, "QW0");
1085           Float_t val=TMath::Max(250., 1.e4*fg.GetParameter(2));
1086           h2->SetBinContent(ieta, iphi, val);
1087         }
1088       }
1089     }
1090   }
1091   for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++) delete h3[iproj];
1092   return kTRUE;
1093 }
1094
1095 //________________________________________________________
1096 Bool_t AliTRDresolution::MakeProjectionTracklet()
1097 {
1098 // Analyse tracklet
1099   Int_t cidx = kTracklet;
1100   if(fProj && fProj->At(cidx)) return kTRUE;
1101   if(!fContainer){
1102     AliError("Missing data container.");
1103     return kFALSE;
1104   }
1105   THnSparse *H(NULL);
1106   if(!(H = (THnSparse*)fContainer->At(cidx))){
1107 //    AliError(Form("Missing/Wrong data @ %d.", cidx));
1108     return kFALSE;
1109   }
1110   return kTRUE;
1111 }
1112
1113 //________________________________________________________
1114 Bool_t AliTRDresolution::MakeProjectionTrackIn()
1115 {
1116 // Analyse track in
1117
1118   Int_t cidx = kTrackIn;
1119   if(fProj && fProj->At(cidx)) return kTRUE;
1120   if(!fContainer){
1121     AliError("Missing data container.");
1122     return kFALSE;
1123   }
1124   THnSparse *H(NULL);
1125   if(!(H = (THnSparse*)fContainer->At(cidx))){
1126     AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1127     return kFALSE;
1128   }
1129
1130   Int_t  coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1131   TAxis //*abc(H->GetAxis(kBC)),
1132         *aphi(H->GetAxis(kPhi)),
1133         *aeta(H->GetAxis(kEta)),
1134         //*as(H->GetAxis(kSpeciesChgRC)),
1135         //*apt(H->GetAxis(kPt)),
1136         *ay(H->GetAxis(kYrez)),
1137         *az(H->GetAxis(kZrez)),
1138         *ap(H->GetAxis(kPrez));
1139   Int_t neta(aeta->GetNbins()), nphi(aphi->GetNbins());
1140   TH3I *h3[fgNproj[cidx]];
1141   h3[0] = new TH3I("h3TI0", Form("r-#phi residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
1142             neta, aeta->GetXmin(), aeta->GetXmax(),
1143             nphi, aphi->GetXmin(), aphi->GetXmax(),
1144             ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1145   h3[1] = (TH3I*)h3[0]->Clone("h3TI1"); h3[1]->SetTitle("r-#phi residuals for pos tracks");
1146   h3[2] = new TH3I("h3TI2", Form("z residuals for row cross;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), az->GetTitle()),
1147             neta, aeta->GetXmin(), aeta->GetXmax(),
1148             nphi, aphi->GetXmin(), aphi->GetXmax(),
1149             az->GetNbins(), az->GetXmin(), az->GetXmax());
1150   h3[3] = new TH3I("h3TI3", Form("angular residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ap->GetTitle()),
1151             neta, aeta->GetXmin(), aeta->GetXmax(),
1152             nphi, aphi->GetXmin(), aphi->GetXmax(),
1153             ap->GetNbins(), ap->GetXmin(), ap->GetXmax());
1154   h3[4] = (TH3I*)h3[3]->Clone("h3TI4"); h3[4]->SetTitle("angular residuals for pos tracks");
1155   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1156     v = H->GetBinContent(ib, coord);
1157     if(v<1.) continue;
1158     if(coord[kBC]>1) continue; // bunch cross cut
1159     // species selection
1160     if(coord[kSpeciesChgRC]<6){
1161       h3[0]->AddBinContent(
1162         h3[0]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1163       h3[3]->AddBinContent(
1164         h3[3]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1165     } else if(coord[kSpeciesChgRC]==6) {
1166       h3[2]->AddBinContent(
1167         h3[2]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1168     } else if(coord[kSpeciesChgRC]>6) {
1169       h3[1]->AddBinContent(
1170         h3[1]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1171       h3[4]->AddBinContent(
1172         h3[4]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1173     }
1174   }
1175   if(!fProj){
1176     AliInfo("Building array of projections ...");
1177     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1178   }
1179   TObjArray *arr(NULL);
1180   fProj->AddAt(arr = new TObjArray(fgNproj[cidx]), cidx);
1181
1182   TH2F *h2(NULL);
1183   for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++){
1184     TAxis *ax(h3[iproj]->GetZaxis());
1185     Float_t zm(ax->GetXmin()/3.), zM(ax->GetXmax()/3.), dz=(zM-zm)/50;
1186     arr->AddAt(h2 = new TH2F(Form("h2TI%d", iproj),
1187             Form("%s;%s;%s;%s", h3[iproj]->GetTitle(), aeta->GetTitle(), aphi->GetTitle(), ax->GetTitle()),
1188             neta, aeta->GetXmin(), aeta->GetXmax(),
1189             nphi, aphi->GetXmin(), aphi->GetXmax()), iproj);
1190     h2->SetContour(50);
1191     h2->GetZaxis()->CenterTitle();
1192     h2->GetZaxis()->SetRangeUser(zm, zM);
1193     zm+=dz; zM-=dz;
1194     for(Int_t iphi(1); iphi<=nphi; iphi++){
1195       for(Int_t ieta(1); ieta<=neta; ieta++){
1196         TH1 *h = h3[iproj]->ProjectionZ(Form("hy%d", iproj), ieta, ieta, iphi, iphi);
1197         h2->SetBinContent(ieta, iphi, GetMeanWithBoundary(h, zm, zM, dz));
1198       }
1199     }
1200   }
1201 //   h2[5] = (TH2F*)h2[0]->Clone("h25");
1202 //   h2[5]->SetTitle("Systematic shift between neg/pos tracks");
1203 //   h2[5]->SetZTitle("#Delta(#Delta^{-}y - #Delta^{+}y) [cm]"); h2[5]->Reset();
1204 //   h2[6] = (TH2F*)h2[1]->Clone("h26");
1205 //   h2[6]->SetTitle("Average shift of pos&neg tracks");
1206 //   h2[6]->SetZTitle("<#Delta^{-}y, #Delta^{+}y> [cm]"); h2[6]->Reset();
1207 //   for(Int_t iphi(1); iphi<=nphi; iphi++){
1208 //     for(Int_t ieta(1); ieta<=neta; ieta++){
1209 //       Float_t neg = h2[0]->GetBinContent(ieta, iphi),
1210 //               pos = h2[1]->GetBinContent(ieta, iphi);
1211 //       if(neg<-100 || pos<-100){
1212 //         h2[5]->SetBinContent(ieta, iphi, -999.);
1213 //         h2[6]->SetBinContent(ieta, iphi, -999.);
1214 //       } else {
1215 //         h2[5]->SetBinContent(ieta, iphi, neg-pos);
1216 //         h2[6]->SetBinContent(ieta, iphi, 0.5*(neg+pos));
1217 //       }
1218 //     }
1219 //   }
1220
1221   for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++) delete h3[iproj];
1222   return kTRUE;
1223 }
1224
1225
1226
1227 //________________________________________________________
1228 Bool_t AliTRDresolution::PostProcess()
1229 {
1230 // Fit, Project, Combine, Extract values from the containers filled during execution
1231
1232   if (!fContainer) {
1233     AliError("ERROR: list not available");
1234     return kFALSE;
1235   }
1236
1237   // DEFINE MODELS
1238   // simple gauss
1239   TF1 fg("fGauss", "gaus", -.5, .5);  
1240   // Landau for charge resolution
1241   TF1 fch("fClCh", "landau", 0., 1000.);  
1242   // Landau for e+- pt resolution
1243   TF1 fpt("fPt", "landau", -0.1, 0.2);  
1244
1245   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1246   // Clusters residuals
1247   if(!MakeProjectionCluster()) return kFALSE;
1248   fNRefFigures = 3;
1249   // Tracklet residual/pulls
1250   if(!MakeProjectionTracklet()) return kFALSE;
1251   fNRefFigures = 7;
1252   // TRDin residual/pulls
1253   if(!MakeProjectionTrackIn()) return kFALSE;
1254   fNRefFigures = 11;
1255   // TRDout residual/pulls
1256 //  if(!MakeProjectionTrackOut()) return kFALSE;
1257   fNRefFigures = 15;
1258
1259   if(!HasMCdata()) return kTRUE;
1260
1261
1262   //PROCESS MC RESIDUAL DISTRIBUTIONS
1263
1264   // CLUSTER Y RESOLUTION/PULLS
1265 //  if(!MakeProjectionClusterMC()) return kFALSE;
1266   fNRefFigures = 17;
1267
1268   // TRACKLET RESOLUTION/PULLS
1269 //  if(!MakeProjectionTrackletMC()) return kFALSE;
1270   fNRefFigures = 21;
1271
1272   // TRACK RESOLUTION/PULLS
1273 /*  Process3Darray(kMCtrack, 0, &fg, 1.e4);   // y
1274   Process3DlinkedArray(kMCtrack, 1, &fg);   // y PULL
1275   Process3Darray(kMCtrack, 2, &fg, 1.e4);   // z
1276   Process3Darray(kMCtrack, 3, &fg);         // z PULL
1277   Process2Darray(kMCtrack, 4, &fg, 1.e3);   // phi
1278   Process2Darray(kMCtrack, 5, &fg);         // snp PULL
1279   Process2Darray(kMCtrack, 6, &fg, 1.e3);   // theta
1280   Process2Darray(kMCtrack, 7, &fg);         // tgl PULL
1281   Process3Darray(kMCtrack, 8, &fg, 1.e2);   // pt resolution
1282   Process3Darray(kMCtrack, 9, &fg);         // 1/pt pulls
1283   Process3Darray(kMCtrack, 10, &fg, 1.e2);  // p resolution*/
1284 //  if(!MakeProjectionTrackMC(kMCtrack)) return kFALSE;
1285   fNRefFigures+=16;
1286
1287   // TRACK TRDin RESOLUTION/PULLS
1288 //  if(!MakeProjectionTrackMC(kMCtrackIn)) return kFALSE;
1289   fNRefFigures+=8;
1290
1291   // TRACK TRDout RESOLUTION/PULLS
1292 //  if(!MakeProjectionTrackMC(kMCtrackOut)) return kFALSE;
1293   fNRefFigures+=8;
1294
1295   return kTRUE;
1296 }
1297
1298
1299 //________________________________________________________
1300 void AliTRDresolution::Terminate(Option_t *opt)
1301 {
1302   AliTRDrecoTask::Terminate(opt);
1303   if(HasPostProcess()) PostProcess();
1304 }
1305
1306 //________________________________________________________
1307 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1308 {
1309 // Helper function to avoid duplication of code
1310 // Make first guesses on the fit parameters
1311
1312   // find the intial parameters of the fit !! (thanks George)
1313   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1314   Double_t sum = 0.;
1315   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1316   f->SetParLimits(0, 0., 3.*sum);
1317   f->SetParameter(0, .9*sum);
1318   Double_t rms = h->GetRMS();
1319   f->SetParLimits(1, -rms, rms);
1320   f->SetParameter(1, h->GetMean());
1321
1322   f->SetParLimits(2, 0., 2.*rms);
1323   f->SetParameter(2, rms);
1324   if(f->GetNpar() <= 4) return;
1325
1326   f->SetParLimits(3, 0., sum);
1327   f->SetParameter(3, .1*sum);
1328
1329   f->SetParLimits(4, -.3, .3);
1330   f->SetParameter(4, 0.);
1331
1332   f->SetParLimits(5, 0., 1.e2);
1333   f->SetParameter(5, 2.e-1);
1334 }
1335
1336 //________________________________________________________
1337 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
1338 {
1339 // Build performance histograms for AliTRDcluster.vs TRD track or MC
1340 //  - y reziduals/pulls
1341
1342   TObjArray *arr = new TObjArray(2);
1343   arr->SetName(name); arr->SetOwner();
1344   TH1 *h(NULL); char hname[100], htitle[300];
1345
1346   // tracklet resolution/pull in y direction
1347   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
1348   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
1349   Float_t rr = range<0.?fDyRange:range;
1350   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1351     Int_t nybins=50;
1352     if(expand) nybins*=2;
1353     h = new TH3S(hname, htitle, 
1354                  48, -.48, .48,            // phi
1355                  60, -rr, rr,              // dy
1356                  nybins, -0.5, nybins-0.5);// segment
1357   } else h->Reset();
1358   arr->AddAt(h, 0);
1359   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
1360   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
1361   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1362     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
1363   } else h->Reset();
1364   arr->AddAt(h, 1);
1365
1366   return arr;
1367 }
1368
1369 //________________________________________________________
1370 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
1371 {
1372 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1373 //  - y reziduals/pulls
1374 //  - z reziduals/pulls
1375 //  - phi reziduals
1376   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05); 
1377   arr->Expand(5);
1378   TH1 *h(NULL); char hname[100], htitle[300];
1379
1380   // tracklet resolution/pull in z direction
1381   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
1382   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1383   if(!(h = (TH2S*)gROOT->FindObject(hname))){
1384     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
1385   } else h->Reset();
1386   arr->AddAt(h, 2);
1387   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1388   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
1389   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1390     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
1391     h->GetZaxis()->SetBinLabel(1, "no RC");
1392     h->GetZaxis()->SetBinLabel(2, "RC");
1393   } else h->Reset();
1394   arr->AddAt(h, 3);
1395
1396   // tracklet to track phi resolution
1397   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
1398   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1399   Int_t nsgms=540;
1400   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1401     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
1402   } else h->Reset();
1403   arr->AddAt(h, 4);
1404
1405   return arr;
1406 }
1407
1408 //________________________________________________________
1409 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1410 {
1411 // Build performance histograms for AliExternalTrackParam.vs MC
1412 //  - y resolution/pulls
1413 //  - z resolution/pulls
1414 //  - phi resolution, snp pulls
1415 //  - theta resolution, tgl pulls
1416 //  - pt resolution, 1/pt pulls, p resolution
1417
1418   TObjArray *arr = BuildMonitorContainerTracklet(name); 
1419   arr->Expand(11);
1420   TH1 *h(NULL); char hname[100], htitle[300];
1421   TAxis *ax(NULL);
1422
1423   // snp pulls
1424   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1425   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
1426   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1427     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1428   } else h->Reset();
1429   arr->AddAt(h, 5);
1430
1431   // theta resolution
1432   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1433   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
1434   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1435     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1436   } else h->Reset();
1437   arr->AddAt(h, 6);
1438   // tgl pulls
1439   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1440   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
1441   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1442     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1443   } else h->Reset();
1444   arr->AddAt(h, 7);
1445   
1446   const Int_t kNdpt(150); 
1447   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1448   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
1449   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
1450   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1451   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1452   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
1453
1454   // Pt resolution
1455   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1456   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
1457   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1458     h = new TH3S(hname, htitle, 
1459                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1460     //ax = h->GetZaxis();
1461     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1462   } else h->Reset();
1463   arr->AddAt(h, 8);
1464   // 1/Pt pulls
1465   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1466   snprintf(htitle, 300, "#splitline{1/P_{t} pull for}{\"%s\" @ %s};1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
1467   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1468     h = new TH3S(hname, htitle, 
1469                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1470     ax = h->GetZaxis();
1471     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1472   } else h->Reset();
1473   arr->AddAt(h, 9);
1474   // P resolution
1475   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1476   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
1477   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1478     h = new TH3S(hname, htitle, 
1479                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1480     ax = h->GetZaxis();
1481     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1482   } else h->Reset();
1483   arr->AddAt(h, 10);
1484
1485   return arr;
1486 }
1487
1488
1489 //________________________________________________________
1490 TObjArray* AliTRDresolution::Histos()
1491 {
1492   //
1493   // Define histograms
1494   //
1495
1496   if(fContainer) return fContainer;
1497
1498   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1499   THnSparse *H(NULL);
1500   const Int_t nhn(100); Char_t hn[nhn]; TString st;
1501
1502   //++++++++++++++++++++++
1503   // cluster to track residuals/pulls
1504   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1505   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1506     const Char_t *clTitle[5/*kNdim*/] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], "chg*Q/vd/angle", fgkTitle[kYrez]/*, "#Deltax [cm]", "Q</Q", "#Phi^{*} - ExB [deg]"*/};
1507     const Int_t clNbins[5/*kNdim*/]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], 51, fgkNbins[kYrez]/*, 33, 10, 30*/};
1508     const Double_t clMin[5/*kNdim*/]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], -102., fgkMin[kYrez]/3./*, 0., 0.1, -30*/},
1509                    clMax[5/*kNdim*/]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], 102., fgkMax[kYrez]/3./*, 3.3, 2.1, 30*/};
1510     st = "cluster spatial&charge resolution;";
1511     for(Int_t idim(0); idim<5/*kNdim*/; idim++){ st += clTitle[idim]; st+=";";}
1512     H = new THnSparseI(hn, st.Data(), kNdim, clNbins, clMin, clMax);
1513   } else H->Reset();
1514   fContainer->AddAt(H, kCluster);
1515   //++++++++++++++++++++++
1516   // tracklet to TRD track
1517   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1518   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1519     const Char_t *trTitle[kNdim+1] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kSpeciesChgRC], fgkTitle[kYrez], "#Deltaz [cm]/#Phi^{*} - ExB [rad]", fgkTitle[kPrez], "dq/dl [a.u.]"/*, fgkTitle[kPt]*/};
1520     const Int_t trNbins[kNdim+1]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kSpeciesChgRC], fgkNbins[kYrez], fgkNbins[kZrez], fgkNbins[kPrez], 30/*, fgkNbins[kPt]*/};
1521     const Double_t trMin[kNdim+1]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kSpeciesChgRC], fgkMin[kYrez], fgkMin[kZrez], fgkMin[kPrez], 0./*, fgkMin[kPt]*/},
1522                    trMax[kNdim+1]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kSpeciesChgRC], fgkMax[kYrez], fgkMax[kZrez], fgkMax[kPrez], 3000./*, fgkMax[kPt]*/};
1523     st = "tracklet spatial&charge resolution;";
1524     for(Int_t idim(0); idim<kNdim+1; idim++){ st += trTitle[idim]; st+=";";}
1525     H = new THnSparseI(hn, st.Data(), kNdim+1, trNbins, trMin, trMax);
1526   } else H->Reset();
1527   fContainer->AddAt(H, kTracklet);
1528   //++++++++++++++++++++++
1529   // tracklet to TRDin
1530   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1531   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1532     st = "r-#phi/z/angular residuals @ TRD entry;";
1533     for(Int_t idim(0); idim<kNdim; idim++){ st += fgkTitle[idim]; st+=";";}
1534     H = new THnSparseI(hn, st.Data(), kNdim, fgkNbins, fgkMin, fgkMax);
1535   } else H->Reset();
1536   fContainer->AddAt(H, kTrackIn);
1537   // tracklet to TRDout
1538   fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1539
1540
1541   // Resolution histos
1542   if(!HasMCdata()) return fContainer;
1543
1544   // cluster resolution 
1545   fContainer->AddAt(BuildMonitorContainerCluster("MCcl"),  kMCcluster);
1546
1547   // tracklet resolution
1548   fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
1549
1550   // track resolution
1551   TObjArray *arr(NULL);
1552   fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
1553   arr->SetName("MCtrk");
1554   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
1555
1556   // TRDin TRACK RESOLUTION
1557   fContainer->AddAt(H, kMCtrackIn);
1558
1559   // TRDout TRACK RESOLUTION
1560   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1561
1562   return fContainer;
1563 }
1564
1565 //________________________________________________________
1566 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1567 {
1568 // Robust function to process sigma/mean for 2D plot dy(x)
1569 // For each x bin a gauss fit is performed on the y projection and the range
1570 // with the minimum chi2/ndf is choosen
1571
1572   if(!h2) {
1573     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
1574     return kFALSE;
1575   }
1576   if(!Int_t(h2->GetEntries())){
1577     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
1578     return kFALSE;
1579   }
1580   if(!g || !g[0]|| !g[1]) {
1581     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
1582     return kFALSE;
1583   }
1584   // prepare
1585   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1586   Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
1587   TF1 f("f", "gaus", ymin, ymax);
1588   Int_t n(0);
1589   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1590   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1591   TH1D *h(NULL);
1592   if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
1593   Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
1594   
1595
1596   // do actual loop
1597   Float_t chi2OverNdf(0.);
1598   for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
1599     x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
1600     ymin = ay->GetXmin(); ymax = ay->GetXmax();
1601
1602     h = h2->ProjectionY("py", ix, ix);
1603     if((n=(Int_t)h->GetEntries())<stat){
1604       if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
1605       continue;
1606     }
1607     // looking for a first order mean value
1608     f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
1609     h->Fit(&f, "QNW");
1610     chi2OverNdf = f.GetChisquare()/f.GetNDF();
1611     printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
1612     y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
1613     ey  = f.GetParError(1);
1614     sy = f.GetParameter(2); esy = f.GetParError(2);
1615 //     // looking for the best chi2/ndf value
1616 //     while(ymin<y0 && ymax>y1){
1617 //       f.SetParameter(1, y);
1618 //       f.SetParameter(2, sy);
1619 //       h->Fit(&f, "QNW", "", y0, y1);
1620 //       printf("   range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
1621 //       if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
1622 //         chi2OverNdf = f.GetChisquare()/f.GetNDF();
1623 //         y  = f.GetParameter(1); ey  = f.GetParError(1);
1624 //         sy = f.GetParameter(2); esy = f.GetParError(2);
1625 //         printf("    set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
1626 //       }
1627 //       y0-=dy; y1+=dy;
1628 //     }
1629     g[0]->SetPoint(np, x, y);
1630     g[0]->SetPointError(np, ex, ey);
1631     g[1]->SetPoint(np, x, sy);
1632     g[1]->SetPointError(np, ex, esy);
1633     np++;
1634   }
1635   return kTRUE;
1636 }
1637
1638
1639 //________________________________________________________
1640 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1641 {
1642   //
1643   // Do the processing
1644   //
1645
1646   Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
1647   Int_t n = 0;
1648   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1649   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1650   if(Int_t(h2->GetEntries())){ 
1651     AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
1652   } else {
1653     AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
1654     fIdxPlot++;
1655     return kTRUE;
1656   }
1657
1658   const Int_t kINTEGRAL=1;
1659   for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
1660     Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
1661     Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
1662     TH1D *h = h2->ProjectionY(pn, abin, bbin);
1663     if((n=(Int_t)h->GetEntries())<150){ 
1664       AliDebug(4, Form("  x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
1665       continue;
1666     }
1667     h->Fit(f, "QN");
1668     Int_t ip = g[0]->GetN();
1669     AliDebug(4, Form("  x_%d[%f] range[%d %d] stat[%d] M[%f] Sgm[%f]", ip, x, abin, bbin, n, f->GetParameter(1), f->GetParameter(2)));
1670     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1671     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1672     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1673     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1674 /*  
1675     g[0]->SetPoint(ip, x, k*h->GetMean());
1676     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1677     g[1]->SetPoint(ip, x, k*h->GetRMS());
1678     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1679   }
1680   fIdxPlot++;
1681   return kTRUE;
1682 }
1683
1684
1685 //____________________________________________________________________
1686 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
1687 {
1688 //
1689 // Fit track with a staight line using the "np" clusters stored in the array "points".
1690 // The following particularities are stored in the clusters from points:
1691 //   1. pad tilt as cluster charge
1692 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
1693 // The parameters of the straight line fit are stored in the array "param" in the following order :
1694 //     param[0] - x0 reference radial position
1695 //     param[1] - y0 reference r-phi position @ x0
1696 //     param[2] - z0 reference z position @ x0
1697 //     param[3] - slope dy/dx
1698 //     param[4] - slope dz/dx
1699 //
1700 // Attention :
1701 // Function should be used to refit tracks for B=0T
1702 //
1703
1704   if(np<40){
1705     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
1706     return kFALSE;
1707   }
1708   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
1709
1710   Double_t x0(0.);
1711   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
1712   x0/=Float_t(np);
1713
1714   Double_t x, y, z, dx, tilt(0.);
1715   for(Int_t ip(0); ip<np; ip++){
1716     x = points[ip].GetX(); z = points[ip].GetZ();
1717     dx = x - x0;
1718     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
1719   }
1720   if(zfitter.Eval() != 0) return kFALSE;
1721
1722   Double_t z0    = zfitter.GetParameter(0);
1723   Double_t dzdx  = zfitter.GetParameter(1);
1724   for(Int_t ip(0); ip<np; ip++){
1725     if(points[ip].GetClusterType()) continue;
1726     x    = points[ip].GetX();
1727     dx   = x - x0;
1728     y    = points[ip].GetY();
1729     z    = points[ip].GetZ();
1730     tilt = points[ip].GetCharge();
1731     y -= tilt*(-dzdx*dx + z - z0);
1732     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
1733     yfitter.AddPoint(&dx, y, 1.);
1734   }
1735   if(yfitter.Eval() != 0) return kFALSE;
1736   Double_t y0   = yfitter.GetParameter(0);
1737   Double_t dydx = yfitter.GetParameter(1);
1738
1739   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
1740   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
1741   return kTRUE;
1742 }
1743
1744 //____________________________________________________________________
1745 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
1746 {
1747 //
1748 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
1749 // See function FitTrack for the data stored in the "clusters" array
1750
1751 // The parameters of the straight line fit are stored in the array "param" in the following order :
1752 //     par[0] - x0 reference radial position
1753 //     par[1] - y0 reference r-phi position @ x0
1754 //     par[2] - slope dy/dx
1755 //
1756 // Attention :
1757 // Function should be used to refit tracks for B=0T
1758 //
1759
1760   TLinearFitter yfitter(2, "pol1");
1761
1762   // grep data for tracklet
1763   Double_t x0(0.), x[60], y[60], dy[60];
1764   Int_t nly(0);
1765   for(Int_t ip(0); ip<np; ip++){
1766     if(points[ip].GetClusterType()) continue;
1767     if(points[ip].GetVolumeID() != ly) continue;
1768     Float_t xt(points[ip].GetX())
1769            ,yt(param[1] + param[3] * (xt - param[0]));
1770     x[nly] = xt;
1771     y[nly] = points[ip].GetY();
1772     dy[nly]= y[nly]-yt;
1773     x0    += xt;
1774     nly++;
1775   }
1776   if(nly<10){
1777     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
1778     return kFALSE;
1779   }
1780   // set radial reference for fit
1781   x0 /= Float_t(nly);
1782
1783   // find tracklet core
1784   Double_t mean(0.), sig(1.e3);
1785   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
1786
1787   // simple cluster error parameterization
1788   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
1789
1790   // fit tracklet core
1791   for(Int_t jly(0); jly<nly; jly++){
1792     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
1793     Double_t dx(x[jly]-x0);
1794     yfitter.AddPoint(&dx, y[jly], 1.);
1795   }
1796   if(yfitter.Eval() != 0) return kFALSE;
1797   par[0] = x0;
1798   par[1] = yfitter.GetParameter(0);
1799   par[2] = yfitter.GetParameter(1);
1800   return kTRUE;
1801 }
1802
1803 //____________________________________________________________________
1804 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
1805 {
1806 //
1807 // Global selection mechanism of tracksbased on cluster to fit residuals
1808 // The parameters are the same as used ni function FitTrack().
1809
1810   const Float_t kS(0.6), kM(0.2);
1811   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
1812   Float_t dy, dz, s, m;
1813   for(Int_t ip(0); ip<np; ip++){
1814     if(points[ip].GetClusterType()) continue;
1815     Float_t x0(points[ip].GetX())
1816           ,y0(param[1] + param[3] * (x0 - param[0]))
1817           ,z0(param[2] + param[4] * (x0 - param[0]));
1818     dy=points[ip].GetY() - y0; h.Fill(dy);
1819     dz=points[ip].GetZ() - z0;
1820   }
1821   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
1822   fg.SetParameter(1, 0.);
1823   fg.SetParameter(2, 2.e-2);
1824   h.Fit(&fg, "QN");
1825   m=fg.GetParameter(1); s=fg.GetParameter(2);
1826   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
1827   return kTRUE;
1828 }
1829
1830 //________________________________________________________
1831 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
1832 {
1833   //
1834   // Get the most probable value and the full width half mean 
1835   // of a Landau distribution
1836   //
1837
1838   const Float_t dx = 1.;
1839   mpv = f->GetParameter(1);
1840   Float_t fx, max = f->Eval(mpv);
1841
1842   xm = mpv - dx;
1843   while((fx = f->Eval(xm))>.5*max){
1844     if(fx>max){ 
1845       max = fx;
1846       mpv = xm;
1847     }
1848     xm -= dx;
1849   }
1850
1851   xM += 2*(mpv - xm);
1852   while((fx = f->Eval(xM))>.5*max) xM += dx;
1853 }
1854
1855
1856 // #include "TFile.h"
1857 // //________________________________________________________
1858 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
1859 // {
1860 //   if(!file){
1861 //     AliWarning("Use cluster position as in reconstruction.");
1862 //     SetLoadCorrection();
1863 //     return kTRUE;
1864 //   }
1865 //   TDirectory *cwd(gDirectory);
1866 //   TString fileList;
1867 //   FILE *filePtr = fopen(file, "rt");
1868 //   if(!filePtr){
1869 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
1870 //     SetLoadCorrection();
1871 //     return kFALSE;
1872 //   }
1873 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
1874 //   while(fileList.Gets(filePtr)){
1875 //     if(!TFile::Open(fileList.Data())) {
1876 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
1877 //       continue;
1878 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
1879 // 
1880 //     TTree *tSys = (TTree*)gFile->Get("tSys");
1881 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
1882 //     tSys->Draw("det:t>>h2", "dx", "goff");
1883 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1884 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
1885 //     }
1886 //     h2->SetDirectory(cwd);
1887 //     gFile->Close();
1888 //   }
1889 //   cwd->cd();
1890 // 
1891 //   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
1892 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1893 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
1894 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1895 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
1896 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
1897 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1898 //       printf("%03d|", idet);
1899 //       fprintf(fout, "    {");
1900 //       for(Int_t it(0); it<30; it++){
1901 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
1902 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
1903 //       }
1904 //       printf("\n");
1905 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
1906 //     }
1907 //     fprintf(fout, "  };\n");
1908 //   }
1909 //   SetLoadCorrection();
1910 //   return kTRUE;
1911 // }