fix coverity
[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   AliDebug(2, Form("%s[%d]", H->GetName(), H->GetNdimensions()));
1111   return kTRUE;
1112 }
1113
1114 //________________________________________________________
1115 Bool_t AliTRDresolution::MakeProjectionTrackIn()
1116 {
1117 // Analyse track in
1118
1119   Int_t cidx = kTrackIn;
1120   if(fProj && fProj->At(cidx)) return kTRUE;
1121   if(!fContainer){
1122     AliError("Missing data container.");
1123     return kFALSE;
1124   }
1125   THnSparse *H(NULL);
1126   if(!(H = (THnSparse*)fContainer->At(cidx))){
1127     AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1128     return kFALSE;
1129   }
1130
1131   Int_t  coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1132   TAxis //*abc(H->GetAxis(kBC)),
1133         *aphi(H->GetAxis(kPhi)),
1134         *aeta(H->GetAxis(kEta)),
1135         //*as(H->GetAxis(kSpeciesChgRC)),
1136         //*apt(H->GetAxis(kPt)),
1137         *ay(H->GetAxis(kYrez)),
1138         *az(H->GetAxis(kZrez)),
1139         *ap(H->GetAxis(kPrez));
1140   Int_t neta(aeta->GetNbins()), nphi(aphi->GetNbins());
1141   TH3I *h3[fgNproj[cidx]];
1142   h3[0] = new TH3I("h3TI0", Form("r-#phi residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
1143             neta, aeta->GetXmin(), aeta->GetXmax(),
1144             nphi, aphi->GetXmin(), aphi->GetXmax(),
1145             ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1146   h3[1] = (TH3I*)h3[0]->Clone("h3TI1"); h3[1]->SetTitle("r-#phi residuals for pos tracks");
1147   h3[2] = new TH3I("h3TI2", Form("z residuals for row cross;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), az->GetTitle()),
1148             neta, aeta->GetXmin(), aeta->GetXmax(),
1149             nphi, aphi->GetXmin(), aphi->GetXmax(),
1150             az->GetNbins(), az->GetXmin(), az->GetXmax());
1151   h3[3] = new TH3I("h3TI3", Form("angular residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ap->GetTitle()),
1152             neta, aeta->GetXmin(), aeta->GetXmax(),
1153             nphi, aphi->GetXmin(), aphi->GetXmax(),
1154             ap->GetNbins(), ap->GetXmin(), ap->GetXmax());
1155   h3[4] = (TH3I*)h3[3]->Clone("h3TI4"); h3[4]->SetTitle("angular residuals for pos tracks");
1156   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1157     v = H->GetBinContent(ib, coord);
1158     if(v<1.) continue;
1159     if(coord[kBC]>1) continue; // bunch cross cut
1160     // species selection
1161     if(coord[kSpeciesChgRC]<6){
1162       h3[0]->AddBinContent(
1163         h3[0]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1164       h3[3]->AddBinContent(
1165         h3[3]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1166     } else if(coord[kSpeciesChgRC]==6) {
1167       h3[2]->AddBinContent(
1168         h3[2]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1169     } else if(coord[kSpeciesChgRC]>6) {
1170       h3[1]->AddBinContent(
1171         h3[1]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1172       h3[4]->AddBinContent(
1173         h3[4]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1174     }
1175   }
1176   if(!fProj){
1177     AliInfo("Building array of projections ...");
1178     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1179   }
1180   TObjArray *arr(NULL);
1181   fProj->AddAt(arr = new TObjArray(fgNproj[cidx]), cidx);
1182
1183   TH2F *h2(NULL);
1184   for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++){
1185     TAxis *ax(h3[iproj]->GetZaxis());
1186     Float_t zm(ax->GetXmin()/3.), zM(ax->GetXmax()/3.), dz=(zM-zm)/50;
1187     arr->AddAt(h2 = new TH2F(Form("h2TI%d", iproj),
1188             Form("%s;%s;%s;%s", h3[iproj]->GetTitle(), aeta->GetTitle(), aphi->GetTitle(), ax->GetTitle()),
1189             neta, aeta->GetXmin(), aeta->GetXmax(),
1190             nphi, aphi->GetXmin(), aphi->GetXmax()), iproj);
1191     h2->SetContour(50);
1192     h2->GetZaxis()->CenterTitle();
1193     h2->GetZaxis()->SetRangeUser(zm, zM);
1194     zm+=dz; zM-=dz;
1195     for(Int_t iphi(1); iphi<=nphi; iphi++){
1196       for(Int_t ieta(1); ieta<=neta; ieta++){
1197         TH1 *h = h3[iproj]->ProjectionZ(Form("hy%d", iproj), ieta, ieta, iphi, iphi);
1198         h2->SetBinContent(ieta, iphi, GetMeanWithBoundary(h, zm, zM, dz));
1199       }
1200     }
1201   }
1202 //   h2[5] = (TH2F*)h2[0]->Clone("h25");
1203 //   h2[5]->SetTitle("Systematic shift between neg/pos tracks");
1204 //   h2[5]->SetZTitle("#Delta(#Delta^{-}y - #Delta^{+}y) [cm]"); h2[5]->Reset();
1205 //   h2[6] = (TH2F*)h2[1]->Clone("h26");
1206 //   h2[6]->SetTitle("Average shift of pos&neg tracks");
1207 //   h2[6]->SetZTitle("<#Delta^{-}y, #Delta^{+}y> [cm]"); h2[6]->Reset();
1208 //   for(Int_t iphi(1); iphi<=nphi; iphi++){
1209 //     for(Int_t ieta(1); ieta<=neta; ieta++){
1210 //       Float_t neg = h2[0]->GetBinContent(ieta, iphi),
1211 //               pos = h2[1]->GetBinContent(ieta, iphi);
1212 //       if(neg<-100 || pos<-100){
1213 //         h2[5]->SetBinContent(ieta, iphi, -999.);
1214 //         h2[6]->SetBinContent(ieta, iphi, -999.);
1215 //       } else {
1216 //         h2[5]->SetBinContent(ieta, iphi, neg-pos);
1217 //         h2[6]->SetBinContent(ieta, iphi, 0.5*(neg+pos));
1218 //       }
1219 //     }
1220 //   }
1221
1222   for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++) delete h3[iproj];
1223   return kTRUE;
1224 }
1225
1226
1227
1228 //________________________________________________________
1229 Bool_t AliTRDresolution::PostProcess()
1230 {
1231 // Fit, Project, Combine, Extract values from the containers filled during execution
1232
1233   if (!fContainer) {
1234     AliError("ERROR: list not available");
1235     return kFALSE;
1236   }
1237
1238   // DEFINE MODELS
1239   // simple gauss
1240   TF1 fg("fGauss", "gaus", -.5, .5);  
1241   // Landau for charge resolution
1242   TF1 fch("fClCh", "landau", 0., 1000.);  
1243   // Landau for e+- pt resolution
1244   TF1 fpt("fPt", "landau", -0.1, 0.2);  
1245
1246   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1247   // Clusters residuals
1248   if(!MakeProjectionCluster()) return kFALSE;
1249   fNRefFigures = 3;
1250   // Tracklet residual/pulls
1251   if(!MakeProjectionTracklet()) return kFALSE;
1252   fNRefFigures = 7;
1253   // TRDin residual/pulls
1254   if(!MakeProjectionTrackIn()) return kFALSE;
1255   fNRefFigures = 11;
1256   // TRDout residual/pulls
1257 //  if(!MakeProjectionTrackOut()) return kFALSE;
1258   fNRefFigures = 15;
1259
1260   if(!HasMCdata()) return kTRUE;
1261
1262
1263   //PROCESS MC RESIDUAL DISTRIBUTIONS
1264
1265   // CLUSTER Y RESOLUTION/PULLS
1266 //  if(!MakeProjectionClusterMC()) return kFALSE;
1267   fNRefFigures = 17;
1268
1269   // TRACKLET RESOLUTION/PULLS
1270 //  if(!MakeProjectionTrackletMC()) return kFALSE;
1271   fNRefFigures = 21;
1272
1273   // TRACK RESOLUTION/PULLS
1274 /*  Process3Darray(kMCtrack, 0, &fg, 1.e4);   // y
1275   Process3DlinkedArray(kMCtrack, 1, &fg);   // y PULL
1276   Process3Darray(kMCtrack, 2, &fg, 1.e4);   // z
1277   Process3Darray(kMCtrack, 3, &fg);         // z PULL
1278   Process2Darray(kMCtrack, 4, &fg, 1.e3);   // phi
1279   Process2Darray(kMCtrack, 5, &fg);         // snp PULL
1280   Process2Darray(kMCtrack, 6, &fg, 1.e3);   // theta
1281   Process2Darray(kMCtrack, 7, &fg);         // tgl PULL
1282   Process3Darray(kMCtrack, 8, &fg, 1.e2);   // pt resolution
1283   Process3Darray(kMCtrack, 9, &fg);         // 1/pt pulls
1284   Process3Darray(kMCtrack, 10, &fg, 1.e2);  // p resolution*/
1285 //  if(!MakeProjectionTrackMC(kMCtrack)) return kFALSE;
1286   fNRefFigures+=16;
1287
1288   // TRACK TRDin RESOLUTION/PULLS
1289 //  if(!MakeProjectionTrackMC(kMCtrackIn)) return kFALSE;
1290   fNRefFigures+=8;
1291
1292   // TRACK TRDout RESOLUTION/PULLS
1293 //  if(!MakeProjectionTrackMC(kMCtrackOut)) return kFALSE;
1294   fNRefFigures+=8;
1295
1296   return kTRUE;
1297 }
1298
1299
1300 //________________________________________________________
1301 void AliTRDresolution::Terminate(Option_t *opt)
1302 {
1303   AliTRDrecoTask::Terminate(opt);
1304   if(HasPostProcess()) PostProcess();
1305 }
1306
1307 //________________________________________________________
1308 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1309 {
1310 // Helper function to avoid duplication of code
1311 // Make first guesses on the fit parameters
1312
1313   // find the intial parameters of the fit !! (thanks George)
1314   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1315   Double_t sum = 0.;
1316   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1317   f->SetParLimits(0, 0., 3.*sum);
1318   f->SetParameter(0, .9*sum);
1319   Double_t rms = h->GetRMS();
1320   f->SetParLimits(1, -rms, rms);
1321   f->SetParameter(1, h->GetMean());
1322
1323   f->SetParLimits(2, 0., 2.*rms);
1324   f->SetParameter(2, rms);
1325   if(f->GetNpar() <= 4) return;
1326
1327   f->SetParLimits(3, 0., sum);
1328   f->SetParameter(3, .1*sum);
1329
1330   f->SetParLimits(4, -.3, .3);
1331   f->SetParameter(4, 0.);
1332
1333   f->SetParLimits(5, 0., 1.e2);
1334   f->SetParameter(5, 2.e-1);
1335 }
1336
1337 //________________________________________________________
1338 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
1339 {
1340 // Build performance histograms for AliTRDcluster.vs TRD track or MC
1341 //  - y reziduals/pulls
1342
1343   TObjArray *arr = new TObjArray(2);
1344   arr->SetName(name); arr->SetOwner();
1345   TH1 *h(NULL); char hname[100], htitle[300];
1346
1347   // tracklet resolution/pull in y direction
1348   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
1349   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
1350   Float_t rr = range<0.?fDyRange:range;
1351   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1352     Int_t nybins=50;
1353     if(expand) nybins*=2;
1354     h = new TH3S(hname, htitle, 
1355                  48, -.48, .48,            // phi
1356                  60, -rr, rr,              // dy
1357                  nybins, -0.5, nybins-0.5);// segment
1358   } else h->Reset();
1359   arr->AddAt(h, 0);
1360   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
1361   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
1362   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1363     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
1364   } else h->Reset();
1365   arr->AddAt(h, 1);
1366
1367   return arr;
1368 }
1369
1370 //________________________________________________________
1371 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
1372 {
1373 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1374 //  - y reziduals/pulls
1375 //  - z reziduals/pulls
1376 //  - phi reziduals
1377   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05); 
1378   arr->Expand(5);
1379   TH1 *h(NULL); char hname[100], htitle[300];
1380
1381   // tracklet resolution/pull in z direction
1382   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
1383   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1384   if(!(h = (TH2S*)gROOT->FindObject(hname))){
1385     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
1386   } else h->Reset();
1387   arr->AddAt(h, 2);
1388   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1389   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
1390   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1391     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
1392     h->GetZaxis()->SetBinLabel(1, "no RC");
1393     h->GetZaxis()->SetBinLabel(2, "RC");
1394   } else h->Reset();
1395   arr->AddAt(h, 3);
1396
1397   // tracklet to track phi resolution
1398   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
1399   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1400   Int_t nsgms=540;
1401   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1402     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
1403   } else h->Reset();
1404   arr->AddAt(h, 4);
1405
1406   return arr;
1407 }
1408
1409 //________________________________________________________
1410 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1411 {
1412 // Build performance histograms for AliExternalTrackParam.vs MC
1413 //  - y resolution/pulls
1414 //  - z resolution/pulls
1415 //  - phi resolution, snp pulls
1416 //  - theta resolution, tgl pulls
1417 //  - pt resolution, 1/pt pulls, p resolution
1418
1419   TObjArray *arr = BuildMonitorContainerTracklet(name); 
1420   arr->Expand(11);
1421   TH1 *h(NULL); char hname[100], htitle[300];
1422   //TAxis *ax(NULL);
1423
1424   // snp pulls
1425   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1426   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
1427   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1428     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1429   } else h->Reset();
1430   arr->AddAt(h, 5);
1431
1432   // theta resolution
1433   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1434   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
1435   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1436     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1437   } else h->Reset();
1438   arr->AddAt(h, 6);
1439   // tgl pulls
1440   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1441   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
1442   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1443     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1444   } else h->Reset();
1445   arr->AddAt(h, 7);
1446   
1447   const Int_t kNdpt(150); 
1448   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1449   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
1450   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
1451   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1452   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1453   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
1454
1455   // Pt resolution
1456   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1457   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
1458   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1459     h = new TH3S(hname, htitle, 
1460                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1461     //ax = h->GetZaxis();
1462     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1463   } else h->Reset();
1464   arr->AddAt(h, 8);
1465   // 1/Pt pulls
1466   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1467   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);
1468   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1469     h = new TH3S(hname, htitle, 
1470                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1471     //ax = h->GetZaxis();
1472     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1473   } else h->Reset();
1474   arr->AddAt(h, 9);
1475   // P resolution
1476   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1477   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
1478   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1479     h = new TH3S(hname, htitle, 
1480                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1481     //ax = h->GetZaxis();
1482     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1483   } else h->Reset();
1484   arr->AddAt(h, 10);
1485
1486   return arr;
1487 }
1488
1489
1490 //________________________________________________________
1491 TObjArray* AliTRDresolution::Histos()
1492 {
1493   //
1494   // Define histograms
1495   //
1496
1497   if(fContainer) return fContainer;
1498
1499   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1500   THnSparse *H(NULL);
1501   const Int_t nhn(100); Char_t hn[nhn]; TString st;
1502
1503   //++++++++++++++++++++++
1504   // cluster to track residuals/pulls
1505   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1506   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1507     const Char_t *clTitle[5/*kNdim*/] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], "chg*Q/vd/angle", fgkTitle[kYrez]/*, "#Deltax [cm]", "Q</Q", "#Phi^{*} - ExB [deg]"*/};
1508     const Int_t clNbins[5/*kNdim*/]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], 51, fgkNbins[kYrez]/*, 33, 10, 30*/};
1509     const Double_t clMin[5/*kNdim*/]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], -102., fgkMin[kYrez]/3./*, 0., 0.1, -30*/},
1510                    clMax[5/*kNdim*/]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], 102., fgkMax[kYrez]/3./*, 3.3, 2.1, 30*/};
1511     st = "cluster spatial&charge resolution;";
1512     for(Int_t idim(0); idim<5/*kNdim*/; idim++){ st += clTitle[idim]; st+=";";}
1513     H = new THnSparseI(hn, st.Data(), kNdim, clNbins, clMin, clMax);
1514   } else H->Reset();
1515   fContainer->AddAt(H, kCluster);
1516   //++++++++++++++++++++++
1517   // tracklet to TRD track
1518   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1519   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1520     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]*/};
1521     const Int_t trNbins[kNdim+1]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kSpeciesChgRC], fgkNbins[kYrez], fgkNbins[kZrez], fgkNbins[kPrez], 30/*, fgkNbins[kPt]*/};
1522     const Double_t trMin[kNdim+1]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kSpeciesChgRC], fgkMin[kYrez], fgkMin[kZrez], fgkMin[kPrez], 0./*, fgkMin[kPt]*/},
1523                    trMax[kNdim+1]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kSpeciesChgRC], fgkMax[kYrez], fgkMax[kZrez], fgkMax[kPrez], 3000./*, fgkMax[kPt]*/};
1524     st = "tracklet spatial&charge resolution;";
1525     for(Int_t idim(0); idim<kNdim+1; idim++){ st += trTitle[idim]; st+=";";}
1526     H = new THnSparseI(hn, st.Data(), kNdim+1, trNbins, trMin, trMax);
1527   } else H->Reset();
1528   fContainer->AddAt(H, kTracklet);
1529   //++++++++++++++++++++++
1530   // tracklet to TRDin
1531   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1532   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1533     st = "r-#phi/z/angular residuals @ TRD entry;";
1534     for(Int_t idim(0); idim<kNdim; idim++){ st += fgkTitle[idim]; st+=";";}
1535     H = new THnSparseI(hn, st.Data(), kNdim, fgkNbins, fgkMin, fgkMax);
1536   } else H->Reset();
1537   fContainer->AddAt(H, kTrackIn);
1538   // tracklet to TRDout
1539   fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1540
1541
1542   // Resolution histos
1543   if(!HasMCdata()) return fContainer;
1544
1545   // cluster resolution 
1546   fContainer->AddAt(BuildMonitorContainerCluster("MCcl"),  kMCcluster);
1547
1548   // tracklet resolution
1549   fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
1550
1551   // track resolution
1552   TObjArray *arr(NULL);
1553   fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
1554   arr->SetName("MCtrk");
1555   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
1556
1557   // TRDin TRACK RESOLUTION
1558   fContainer->AddAt(H, kMCtrackIn);
1559
1560   // TRDout TRACK RESOLUTION
1561   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1562
1563   return fContainer;
1564 }
1565
1566 //________________________________________________________
1567 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1568 {
1569 // Robust function to process sigma/mean for 2D plot dy(x)
1570 // For each x bin a gauss fit is performed on the y projection and the range
1571 // with the minimum chi2/ndf is choosen
1572
1573   if(!h2) {
1574     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
1575     return kFALSE;
1576   }
1577   if(!Int_t(h2->GetEntries())){
1578     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
1579     return kFALSE;
1580   }
1581   if(!g || !g[0]|| !g[1]) {
1582     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
1583     return kFALSE;
1584   }
1585   // prepare
1586   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1587   Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
1588   TF1 f("f", "gaus", ymin, ymax);
1589   Int_t n(0);
1590   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1591   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1592   TH1D *h(NULL);
1593   if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
1594   Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
1595   
1596
1597   // do actual loop
1598   Float_t chi2OverNdf(0.);
1599   for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
1600     x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
1601     ymin = ay->GetXmin(); ymax = ay->GetXmax();
1602
1603     h = h2->ProjectionY("py", ix, ix);
1604     if((n=(Int_t)h->GetEntries())<stat){
1605       if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
1606       continue;
1607     }
1608     // looking for a first order mean value
1609     f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
1610     h->Fit(&f, "QNW");
1611     chi2OverNdf = f.GetChisquare()/f.GetNDF();
1612     printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
1613     y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
1614     ey  = f.GetParError(1);
1615     sy = f.GetParameter(2); esy = f.GetParError(2);
1616 //     // looking for the best chi2/ndf value
1617 //     while(ymin<y0 && ymax>y1){
1618 //       f.SetParameter(1, y);
1619 //       f.SetParameter(2, sy);
1620 //       h->Fit(&f, "QNW", "", y0, y1);
1621 //       printf("   range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
1622 //       if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
1623 //         chi2OverNdf = f.GetChisquare()/f.GetNDF();
1624 //         y  = f.GetParameter(1); ey  = f.GetParError(1);
1625 //         sy = f.GetParameter(2); esy = f.GetParError(2);
1626 //         printf("    set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
1627 //       }
1628 //       y0-=dy; y1+=dy;
1629 //     }
1630     g[0]->SetPoint(np, x, y);
1631     g[0]->SetPointError(np, ex, ey);
1632     g[1]->SetPoint(np, x, sy);
1633     g[1]->SetPointError(np, ex, esy);
1634     np++;
1635   }
1636   return kTRUE;
1637 }
1638
1639
1640 //________________________________________________________
1641 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1642 {
1643   //
1644   // Do the processing
1645   //
1646
1647   Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
1648   Int_t n = 0;
1649   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1650   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1651   if(Int_t(h2->GetEntries())){ 
1652     AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
1653   } else {
1654     AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
1655     fIdxPlot++;
1656     return kTRUE;
1657   }
1658
1659   const Int_t kINTEGRAL=1;
1660   for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
1661     Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
1662     Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
1663     TH1D *h = h2->ProjectionY(pn, abin, bbin);
1664     if((n=(Int_t)h->GetEntries())<150){ 
1665       AliDebug(4, Form("  x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
1666       continue;
1667     }
1668     h->Fit(f, "QN");
1669     Int_t ip = g[0]->GetN();
1670     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)));
1671     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1672     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1673     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1674     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1675 /*  
1676     g[0]->SetPoint(ip, x, k*h->GetMean());
1677     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1678     g[1]->SetPoint(ip, x, k*h->GetRMS());
1679     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1680   }
1681   fIdxPlot++;
1682   return kTRUE;
1683 }
1684
1685
1686 //____________________________________________________________________
1687 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
1688 {
1689 //
1690 // Fit track with a staight line using the "np" clusters stored in the array "points".
1691 // The following particularities are stored in the clusters from points:
1692 //   1. pad tilt as cluster charge
1693 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
1694 // The parameters of the straight line fit are stored in the array "param" in the following order :
1695 //     param[0] - x0 reference radial position
1696 //     param[1] - y0 reference r-phi position @ x0
1697 //     param[2] - z0 reference z position @ x0
1698 //     param[3] - slope dy/dx
1699 //     param[4] - slope dz/dx
1700 //
1701 // Attention :
1702 // Function should be used to refit tracks for B=0T
1703 //
1704
1705   if(np<40){
1706     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
1707     return kFALSE;
1708   }
1709   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
1710
1711   Double_t x0(0.);
1712   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
1713   x0/=Float_t(np);
1714
1715   Double_t x, y, z, dx, tilt(0.);
1716   for(Int_t ip(0); ip<np; ip++){
1717     x = points[ip].GetX(); z = points[ip].GetZ();
1718     dx = x - x0;
1719     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
1720   }
1721   if(zfitter.Eval() != 0) return kFALSE;
1722
1723   Double_t z0    = zfitter.GetParameter(0);
1724   Double_t dzdx  = zfitter.GetParameter(1);
1725   for(Int_t ip(0); ip<np; ip++){
1726     if(points[ip].GetClusterType()) continue;
1727     x    = points[ip].GetX();
1728     dx   = x - x0;
1729     y    = points[ip].GetY();
1730     z    = points[ip].GetZ();
1731     tilt = points[ip].GetCharge();
1732     y -= tilt*(-dzdx*dx + z - z0);
1733     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
1734     yfitter.AddPoint(&dx, y, 1.);
1735   }
1736   if(yfitter.Eval() != 0) return kFALSE;
1737   Double_t y0   = yfitter.GetParameter(0);
1738   Double_t dydx = yfitter.GetParameter(1);
1739
1740   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
1741   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);
1742   return kTRUE;
1743 }
1744
1745 //____________________________________________________________________
1746 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
1747 {
1748 //
1749 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
1750 // See function FitTrack for the data stored in the "clusters" array
1751
1752 // The parameters of the straight line fit are stored in the array "param" in the following order :
1753 //     par[0] - x0 reference radial position
1754 //     par[1] - y0 reference r-phi position @ x0
1755 //     par[2] - slope dy/dx
1756 //
1757 // Attention :
1758 // Function should be used to refit tracks for B=0T
1759 //
1760
1761   TLinearFitter yfitter(2, "pol1");
1762
1763   // grep data for tracklet
1764   Double_t x0(0.), x[60], y[60], dy[60];
1765   Int_t nly(0);
1766   for(Int_t ip(0); ip<np; ip++){
1767     if(points[ip].GetClusterType()) continue;
1768     if(points[ip].GetVolumeID() != ly) continue;
1769     Float_t xt(points[ip].GetX())
1770            ,yt(param[1] + param[3] * (xt - param[0]));
1771     x[nly] = xt;
1772     y[nly] = points[ip].GetY();
1773     dy[nly]= y[nly]-yt;
1774     x0    += xt;
1775     nly++;
1776   }
1777   if(nly<10){
1778     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
1779     return kFALSE;
1780   }
1781   // set radial reference for fit
1782   x0 /= Float_t(nly);
1783
1784   // find tracklet core
1785   Double_t mean(0.), sig(1.e3);
1786   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
1787
1788   // simple cluster error parameterization
1789   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
1790
1791   // fit tracklet core
1792   for(Int_t jly(0); jly<nly; jly++){
1793     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
1794     Double_t dx(x[jly]-x0);
1795     yfitter.AddPoint(&dx, y[jly], 1.);
1796   }
1797   if(yfitter.Eval() != 0) return kFALSE;
1798   par[0] = x0;
1799   par[1] = yfitter.GetParameter(0);
1800   par[2] = yfitter.GetParameter(1);
1801   return kTRUE;
1802 }
1803
1804 //____________________________________________________________________
1805 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
1806 {
1807 //
1808 // Global selection mechanism of tracksbased on cluster to fit residuals
1809 // The parameters are the same as used ni function FitTrack().
1810
1811   const Float_t kS(0.6), kM(0.2);
1812   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
1813   Float_t dy, dz, s, m;
1814   for(Int_t ip(0); ip<np; ip++){
1815     if(points[ip].GetClusterType()) continue;
1816     Float_t x0(points[ip].GetX())
1817           ,y0(param[1] + param[3] * (x0 - param[0]))
1818           ,z0(param[2] + param[4] * (x0 - param[0]));
1819     dy=points[ip].GetY() - y0; h.Fill(dy);
1820     dz=points[ip].GetZ() - z0;
1821   }
1822   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
1823   fg.SetParameter(1, 0.);
1824   fg.SetParameter(2, 2.e-2);
1825   h.Fit(&fg, "QN");
1826   m=fg.GetParameter(1); s=fg.GetParameter(2);
1827   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
1828   return kTRUE;
1829 }
1830
1831 //________________________________________________________
1832 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
1833 {
1834   //
1835   // Get the most probable value and the full width half mean 
1836   // of a Landau distribution
1837   //
1838
1839   const Float_t dx = 1.;
1840   mpv = f->GetParameter(1);
1841   Float_t fx, max = f->Eval(mpv);
1842
1843   xm = mpv - dx;
1844   while((fx = f->Eval(xm))>.5*max){
1845     if(fx>max){ 
1846       max = fx;
1847       mpv = xm;
1848     }
1849     xm -= dx;
1850   }
1851
1852   xM += 2*(mpv - xm);
1853   while((fx = f->Eval(xM))>.5*max) xM += dx;
1854 }
1855
1856
1857 // #include "TFile.h"
1858 // //________________________________________________________
1859 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
1860 // {
1861 //   if(!file){
1862 //     AliWarning("Use cluster position as in reconstruction.");
1863 //     SetLoadCorrection();
1864 //     return kTRUE;
1865 //   }
1866 //   TDirectory *cwd(gDirectory);
1867 //   TString fileList;
1868 //   FILE *filePtr = fopen(file, "rt");
1869 //   if(!filePtr){
1870 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
1871 //     SetLoadCorrection();
1872 //     return kFALSE;
1873 //   }
1874 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
1875 //   while(fileList.Gets(filePtr)){
1876 //     if(!TFile::Open(fileList.Data())) {
1877 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
1878 //       continue;
1879 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
1880 // 
1881 //     TTree *tSys = (TTree*)gFile->Get("tSys");
1882 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
1883 //     tSys->Draw("det:t>>h2", "dx", "goff");
1884 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1885 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
1886 //     }
1887 //     h2->SetDirectory(cwd);
1888 //     gFile->Close();
1889 //   }
1890 //   cwd->cd();
1891 // 
1892 //   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
1893 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1894 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
1895 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1896 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
1897 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
1898 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1899 //       printf("%03d|", idet);
1900 //       fprintf(fout, "    {");
1901 //       for(Int_t it(0); it<30; it++){
1902 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
1903 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
1904 //       }
1905 //       printf("\n");
1906 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
1907 //     }
1908 //     fprintf(fout, "  };\n");
1909 //   }
1910 //   SetLoadCorrection();
1911 //   return kTRUE;
1912 // }