bbdd9472eca709a63516fc83df0b76ad3dc14301
[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 //ClassImp(AliTRDresolution::AliTRDresolutionProjection)
96
97 Int_t const   AliTRDresolution::fgkNbins[kNdim] = {
98   Int_t(kNbunchCross)/*bc*/,
99   180/*phi*/,
100   50/*eta*/,
101   50/*dy*/,
102   40/*dphi*/,
103   50/*dz*/,
104   Int_t(kNcharge)*AliPID::kSPECIES+1/*chg*species*/,
105   kNpt/*pt*/
106 };  //! no of bins/projection
107 Double_t const AliTRDresolution::fgkMin[kNdim] = {
108   -0.5,
109   -TMath::Pi(),
110   -1.,
111   -1.5,
112   -10.,
113   -2.5,
114   -AliPID::kSPECIES-0.5,
115   -0.5
116 };    //! low limits for projections
117 Double_t const AliTRDresolution::fgkMax[kNdim] = {
118   Int_t(kNbunchCross)-0.5,
119   TMath::Pi(),
120   1.,
121   1.5,
122   10.,
123   2.5,
124   AliPID::kSPECIES+0.5,
125   kNpt-0.5
126 };    //! high limits for projections
127 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
128   "bunch cross",
129   "#phi [rad]",
130   "#eta",
131   "#Deltay [cm]",
132   "#Delta#phi [deg]",
133   "#Deltaz [cm]",
134   "chg*spec*rc",
135   "bin_p_{t}"
136 };  //! title of projection
137
138 Int_t const AliTRDresolution::fgkNproj[kNclasses] = {
139   48, 72, 8    // clusters/tracklet/trackIn
140  ,48, 72, 8, 72 // MCclusters/MCtracklet/MCtrackIn/MCtrack
141 // ,5, 11 // trackOut/MCtrackOut
142 };
143 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
144     "Cluster2Track"
145     ,"Tracklet2Track"
146     ,"Tracklet2TRDin" 
147     ,"Cluster2MC"
148     ,"Tracklet2MC"
149     ,"TRDin2MC"
150     ,"TRD2MC"
151 //    ,"Tracklet2TRDout"
152 //    ,"TRDout2MC"
153 };
154 Float_t AliTRDresolution::fgPtBin[kNpt+1];
155
156 //________________________________________________________
157 AliTRDresolution::AliTRDresolution()
158   :AliTRDrecoTask()
159   ,fIdxPlot(0)
160   ,fIdxFrame(0)
161   ,fPtThreshold(1.)
162   ,fDyRange(0.75)
163   ,fProj(NULL)
164   ,fDBPDG(NULL)
165   ,fCl(NULL)
166   ,fMCcl(NULL)
167 {
168   //
169   // Default constructor
170   //
171   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
172   MakePtSegmentation();
173 }
174
175 //________________________________________________________
176 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
177   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
178   ,fIdxPlot(0)
179   ,fIdxFrame(0)
180   ,fPtThreshold(1.)
181   ,fDyRange(0.75)
182   ,fProj(NULL)
183   ,fDBPDG(NULL)
184   ,fCl(NULL)
185   ,fMCcl(NULL)
186 {
187   //
188   // Default constructor
189   //
190
191   InitFunctorList();
192   MakePtSegmentation();
193   if(xchange){
194     SetUseExchangeContainers();
195     DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
196     DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
197   }
198 }
199
200 //________________________________________________________
201 AliTRDresolution::~AliTRDresolution()
202 {
203   //
204   // Destructor
205   //
206
207   if(fProj){fProj->Delete(); delete fProj;}
208   if(fCl){fCl->Delete(); delete fCl;}
209   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
210 }
211
212
213 //________________________________________________________
214 void AliTRDresolution::UserCreateOutputObjects()
215 {
216   // spatial resolution
217
218   AliTRDrecoTask::UserCreateOutputObjects();
219   if(UseExchangeContainers()) InitExchangeContainers();
220 }
221
222 //________________________________________________________
223 void AliTRDresolution::InitExchangeContainers()
224 {
225 // Init containers for subsequent tasks (AliTRDclusterResolution)
226
227   fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
228   fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
229   PostData(kClToTrk, fCl);
230   PostData(kClToMC, fMCcl);
231 }
232
233 //________________________________________________________
234 void AliTRDresolution::UserExec(Option_t *opt)
235 {
236   //
237   // Execution part
238   //
239
240   if(fCl) fCl->Delete();
241   if(fMCcl) fMCcl->Delete();
242   AliTRDrecoTask::UserExec(opt);
243 }
244
245 //________________________________________________________
246 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
247 {
248 // Helper function to calculate pulls in the yz plane 
249 // using proper tilt rotation
250 // Uses functionality defined by AliTRDseedV1.
251
252   return kTRUE;
253 /*
254   Double_t t2(tilt*tilt);
255   // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
256
257   // rotate along pad
258   Double_t cc[3];
259   cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2]; 
260   cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
261   cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
262   // do sqrt
263   Double_t sqr[3]={0., 0., 0.}; 
264   if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
265   Double_t invsqr[3]={0., 0., 0.}; 
266   if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
267   Double_t tmp(dyz[0]);
268   dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
269   dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
270   return kTRUE;
271 */
272 }
273
274 //________________________________________________________
275 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
276 {
277   //
278   // Plot the cluster distributions
279   //
280
281   if(track) fkTrack = track;
282   if(!fkTrack){
283     AliDebug(4, "No Track defined.");
284     return NULL;
285   }
286   if(TMath::Abs(fkESD->GetTOFbc()) > 1){
287     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
288     return NULL;
289   }
290   if(fPt<fPtThreshold){
291     AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
292     return NULL;
293   }
294   THnSparse *H(NULL);
295   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
296     AliWarning("No output container defined.");
297     return NULL;
298   }
299
300   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
301   Double_t val[kNdim]; //Float_t exb, vd, t0, s2, dl, dt;
302   TObjArray     *clInfoArr(NULL);
303   AliTRDseedV1  *fTracklet(NULL);
304   AliTRDcluster *c(NULL), *cc(NULL);
305   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
306     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
307     if(!fTracklet->IsOK()) continue;
308     //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
309     val[kBC]  = ily;
310     val[kPhi] = fPhi;
311     val[kEta] = fEta;
312     val[kPt]  = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
313     Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
314     Int_t row0(-1);
315     Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
316     fTracklet->ResetClusterIter(kTRUE);
317     while((c = fTracklet->NextCluster())){
318       Float_t xc(c->GetX()),
319               q(TMath::Abs(c->GetQ()));
320       if(row0<0) row0 = c->GetPadRow();
321
322       val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
323       val[kPrez] = fTracklet->GetX0()-xc;
324       val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
325       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
326       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
327       if(ic) val[kZrez] /= (ic*q);
328       val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
329       H->Fill(val);
330 /*      // tilt rotation of covariance for clusters
331       Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
332       cov[0] = (sy2+t2*sz2)*corr;
333       cov[1] = tilt*(sz2 - sy2)*corr;
334       cov[2] = (t2*sy2+sz2)*corr;
335       // sum with track covariance
336       cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
337       Double_t dyz[2]= {dy[1], dz[1]};
338       Pulls(dyz, cov, tilt);*/
339   
340       // Get z-position with respect to anode wire
341       Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
342               zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
343       Int_t istk = geo->GetStack(c->GetDetector());
344       AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
345       Float_t rowZ = pp->GetRow0();
346       Float_t d  = rowZ - zt + pp->GetAnodeWireOffset();
347       d -= ((Int_t)(2 * d)) / 2.0;
348       if (d > 0.25) d  = 0.5 - d;
349
350       AliTRDclusterInfo *clInfo(NULL);
351       clInfo = new AliTRDclusterInfo;
352       clInfo->SetCluster(c);
353       //Float_t covF[] = {cov[0], cov[1], cov[2]};
354       clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
355       clInfo->SetResolution(val[kYrez]);
356       clInfo->SetAnisochronity(d);
357       clInfo->SetDriftLength(val[kZrez]);
358       clInfo->SetTilt(fTracklet->GetTilt());
359       if(fCl) fCl->Add(clInfo);
360       //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
361
362       if(DebugLevel()>=2){
363         if(!clInfoArr){ 
364           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
365           clInfoArr->SetOwner(kFALSE);
366         }
367         clInfoArr->Add(clInfo);
368       }
369     }
370     if(DebugLevel()>=2 && clInfoArr){
371       ULong_t status = fkESD->GetStatus();
372       (*DebugStream()) << "cluster"
373         <<"status="  << status
374         <<"clInfo.=" << clInfoArr
375         << "\n";
376       clInfoArr->Clear();
377     }
378   }
379   if(clInfoArr) delete clInfoArr;
380
381   return NULL;//H->Projection(kEta, kPhi);
382 }
383
384
385 //________________________________________________________
386 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
387 {
388 // Plot normalized residuals for tracklets to track. 
389 // 
390 // We start from the result that if X=N(|m|, |Cov|)
391 // BEGIN_LATEX
392 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
393 // END_LATEX
394 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial 
395 // reference position. 
396   if(track) fkTrack = track;
397   if(!fkTrack){
398     AliDebug(4, "No Track defined.");
399     return NULL;
400   }
401   if(TMath::Abs(fkESD->GetTOFbc())>1){
402     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
403     return NULL;
404   }
405   THnSparse *H(NULL);
406   if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
407     AliWarning("No output container defined.");
408     return NULL;
409   }
410 //  return NULL;
411   Double_t val[kNdim+1];
412   AliTRDseedV1 *fTracklet(NULL);
413   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
414     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
415     if(!fTracklet->IsOK()) continue;
416     val [kBC] = il; 
417     val[kPhi] = fPhi;
418     val[kEta] = fEta;
419     val[kSpeciesChgRC]= fSpecies;
420     val[kPt]  = GetPtBin(fTracklet->GetMomentum());
421     Double_t dyt(fTracklet->GetYref(0) - fTracklet->GetYfit(0)),
422              dzt(fTracklet->GetZref(0) - fTracklet->GetZfit(0)),
423              dydx(fTracklet->GetYfit(1)),
424              tilt(fTracklet->GetTilt());
425     // correct for tilt rotation
426     val[kYrez] = dyt - dzt*tilt;
427     val[kZrez] = dzt + dyt*tilt;
428     dydx+= tilt*fTracklet->GetZref(1);
429     val[kPrez] = TMath::ATan((fTracklet->GetYref(1) - dydx)/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
430     if(fTracklet->IsRowCross()){
431       val[kSpeciesChgRC]= 0.;
432 //      val[kPrez] = fkTrack->Charge(); // may be better defined
433     }/* else {
434       Float_t exb, vd, t0, s2, dl, dt;
435       fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
436       val[kZrez] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb));
437     }*/
438     val[kNdim] = fTracklet->GetdQdl();
439     if(DebugLevel()>=1) H->Fill(val);
440
441 //     // compute covariance matrix
442 //     fTracklet->GetCovAt(x, cov);
443 //     fTracklet->GetCovRef(covR);
444 //     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2]; 
445 //     Double_t dyz[2]= {dy[1], dz[1]};
446 //     Pulls(dyz, cov, tilt);
447 //     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
448 //     ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
449
450     if(DebugLevel()>=3){
451       Bool_t rc(fTracklet->IsRowCross());
452       UChar_t err(fTracklet->GetErrorMsg());
453       Double_t x(fTracklet->GetX()),
454                pt(fTracklet->GetPt()),
455                yt(fTracklet->GetYref(0)),
456                zt(fTracklet->GetZref(0)),
457                phi(fTracklet->GetYref(1)),
458                tht(fTracklet->GetZref(1));
459       Int_t ncl(fTracklet->GetN()),
460             det(fTracklet->GetDetector());
461       (*DebugStream()) << "tracklet"
462         <<"pt="  << pt
463         <<"x="   << x
464         <<"yt="  << yt
465         <<"zt="  << zt
466         <<"phi=" << phi
467         <<"tht=" << tht
468         <<"det=" << det
469         <<"n="   << ncl
470         <<"dy0=" << dyt
471         <<"dz0=" << dzt
472         <<"dy="  << val[kYrez]
473         <<"dz="  << val[kZrez]
474         <<"dphi="<< val[kPrez]
475         <<"dQ  ="<< val[kNdim]
476         <<"rc="  << rc
477         <<"err=" << err
478         << "\n";
479     }
480   }
481   return NULL;//H->Projection(kEta, kPhi);
482 }
483
484
485 //________________________________________________________
486 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
487 {
488 // Store resolution/pulls of Kalman before updating with the TRD information 
489 // at the radial position of the first tracklet. The following points are used 
490 // for comparison  
491 //  - the (y,z,snp) of the first TRD tracklet
492 //  - the (y, z, snp, tgl, pt) of the MC track reference
493 // 
494 // Additionally the momentum resolution/pulls are calculated for usage in the 
495 // PID calculation. 
496   //printf("AliTRDresolution::PlotTrackIn() :: track[%p]\n", (void*)track);
497   
498   if(track) fkTrack = track;
499   if(!fkTrack){
500     AliDebug(4, "No Track defined.");
501     return NULL;
502   }
503   //fkTrack->Print();
504   // check container
505   THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
506   if(!H){
507     AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
508     return NULL;
509   }
510   // check input track status
511   AliExternalTrackParam *tin(NULL);
512   if(!(tin = fkTrack->GetTrackIn())){
513     AliError("Track did not entered TRD fiducial volume.");
514     return NULL;
515   }
516   // check first tracklet
517   AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
518   if(!fTracklet){
519     AliDebug(3, "No Tracklet in ly[0]. Skip track.");
520     return NULL;
521   }
522   // check radial position
523   Double_t x = tin->GetX();
524   if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
525     AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
526     return NULL;
527   }
528   //printf("USE y[%+f] dydx[%+f]\n", fTracklet->GetYfit(0), fTracklet->GetYfit(1));
529
530   Int_t bc(TMath::Abs(fkESD->GetTOFbc())%2);
531   const Double_t *parR(tin->GetParameter());
532   Double_t dyt(parR[0] - fTracklet->GetYfit(0)), dzt(parR[1] - fTracklet->GetZfit(0)),
533             phit(fTracklet->GetYfit(1)),
534             tilt(fTracklet->GetTilt());
535
536   // correct for tilt rotation
537   Double_t dy  = dyt - dzt*tilt,
538            dz  = dzt + dyt*tilt;
539   phit       += tilt*parR[3];
540   Double_t dphi = TMath::ASin(parR[2])-TMath::ATan(phit);
541
542   Double_t val[kNdim];
543   val[kBC]          = bc;
544   val[kPhi]         = fPhi;
545   val[kEta]         = fEta;
546   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fSpecies;
547   val[kPt]          = GetPtBin(fPt);
548   val[kYrez]        = dy;
549   val[kZrez]        = dz;
550   val[kPrez]        = dphi*TMath::RadToDeg();
551   H->Fill(val);
552   if(DebugLevel()>=3){
553     (*DebugStream()) << "trackIn"
554       <<"tracklet.="  << fTracklet
555       <<"trackIn.="   << tin
556       << "\n";
557   }
558
559   return NULL; // H->Projection(kEta, kPhi);
560 }
561
562 /*
563 //________________________________________________________
564 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
565 {
566 // Store resolution/pulls of Kalman after last update with the TRD information 
567 // at the radial position of the first tracklet. The following points are used 
568 // for comparison  
569 //  - the (y,z,snp) of the first TRD tracklet
570 //  - the (y, z, snp, tgl, pt) of the MC track reference
571 // 
572 // Additionally the momentum resolution/pulls are calculated for usage in the 
573 // PID calculation. 
574
575   if(track) fkTrack = track;
576   return NULL;
577 }
578 */
579 //________________________________________________________
580 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
581 {
582   //
583   // Plot MC distributions
584   //
585
586   if(!HasMCdata()){ 
587     AliDebug(2, "No MC defined. Results will not be available.");
588     return NULL;
589   }
590   if(track) fkTrack = track;
591   if(!fkTrack){
592     AliDebug(4, "No Track defined.");
593     return NULL;
594   }
595   Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
596
597   THnSparse *H(NULL);
598   if(!fContainer){
599     AliWarning("No output container defined.");
600     return NULL;
601   }
602   // retriev track characteristics
603   Int_t pdg = fkMC->GetPDG(),
604         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
605         sign(0),
606 //        sgm[3],
607         label(fkMC->GetLabel());
608 //        fSegmentLevel(0);
609   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
610   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
611   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
612
613   TH1 *h(NULL);
614   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
615   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
616   UChar_t s;
617   Double_t x, y, z, pt, dydx, dzdx, dzdl;
618   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
619   Double_t covR[7]/*, cov[3]*/;
620   
621 /*  if(DebugLevel()>=3){
622     // get first detector
623     Int_t det = -1;
624     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
625       if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
626       det = fTracklet->GetDetector();
627       break;
628     }
629     if(det>=0){
630       TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
631       Double_t m(-1.);
632       m = fkTrack->GetMass();
633       if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
634         (*DebugStream()) << "MCkalman"
635           << "pdg=" << pdg
636           << "det=" << det
637           << "x="   << &X
638           << "y="   << &Y
639           << "z="   << &Z
640           << "dx="  << &dX
641           << "dy="  << &dY
642           << "dz="  << &dZ
643           << "pt="  << &vPt
644           << "dpt=" << &dPt
645           << "bgt=" << &budget
646           << "cov=" << &cCOV
647           << "\n";
648       }
649     }
650   }*/
651   AliTRDcluster *c(NULL);
652   Double_t val[kNdim+1];
653   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
654     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
655        !fTracklet->IsOK())*/ continue;
656
657     x= x0 = fTracklet->GetX();
658     Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
659     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, phi, s)) continue;
660
661     // MC track position at reference radial position
662     dx  = x0 - x;
663     Float_t ymc = y0 - dx*dydx0;
664     Float_t zmc = z0 - dx*dzdx0;
665     phi -= TMath::Pi();
666
667     val[kBC]  = ily;
668     val[kPhi] = phi;
669     val[kEta] = eta;
670     val[kSpeciesChgRC]= rc?0.:sign*sIdx;
671     val[kPt]  = GetPtBin(pt0);
672     Double_t tilt(fTracklet->GetTilt());
673 //             ,t2(tilt*tilt)
674 //             ,corr(1./(1. + t2))
675 //             ,cost(TMath::Sqrt(corr));
676
677     AliExternalTrackParam *tin(fkTrack->GetTrackIn());
678     if(ily==0 && tin){ // trackIn residuals
679       // check radial position
680       if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
681       else{
682       //  Float_t phi = TMath::ATan2(y0, x0);
683         val[kBC]          = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
684         val[kYrez]        = tin->GetY()-ymc;
685         val[kZrez]        = tin->GetZ()-zmc;
686         val[kPrez]        = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
687         if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
688       }
689     }
690     if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
691
692     // track residuals
693     dydx = fTracklet->GetYref(1);
694     dzdx = fTracklet->GetZref(1);
695     dzdl = fTracklet->GetTgl();
696     y  = fTracklet->GetYref(0);
697     dy = y - ymc;
698     z  = fTracklet->GetZref(0);
699     dz = z - zmc;
700     pt = TMath::Abs(fTracklet->GetPt());
701     fTracklet->GetCovRef(covR);
702
703     val[kYrez] = dy;
704     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
705     val[kZrez] = dz;
706     val[kNdim] = pt/pt0-1.;
707     if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
708 /*      // theta resolution/ tgl pulls
709       Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
710                 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
711       ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
712       ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
713       // pt resolution  \\ 1/pt pulls \\ p resolution for PID
714       Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
715               p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
716       ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
717       ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
718       ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
719
720     // Fill Debug stream for MC track
721     if(DebugLevel()>=4){
722       Int_t det(fTracklet->GetDetector());
723       (*DebugStream()) << "MC"
724         << "det="     << det
725         << "pdg="     << pdg
726         << "sgn="     << sign
727         << "pt="      << pt0
728         << "x="       << x0
729         << "y="       << y0
730         << "z="       << z0
731         << "dydx="    << dydx0
732         << "dzdx="    << dzdx0
733         << "\n";
734     
735       // Fill Debug stream for Kalman track
736       (*DebugStream()) << "MCtrack"
737         << "pt="      << pt
738         << "x="       << x
739         << "y="       << y
740         << "z="       << z
741         << "dydx="    << dydx
742         << "dzdx="    << dzdx
743         << "s2y="     << covR[0]
744         << "s2z="     << covR[2]
745         << "\n";
746     }
747
748     // tracklet residuals
749     dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
750     dzdx = fTracklet->GetZfit(1);
751     y  = fTracklet->GetYfit(0);
752     dy = y - ymc;
753     z  = fTracklet->GetZfit(0);
754     dz = z - zmc;
755     val[kYrez] = dy - dz*tilt;
756     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
757     val[kZrez] = dz + dy*tilt;
758 //      val[kNdim] = pt/pt0-1.;
759     if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
760     
761
762     // Fill Debug stream for tracklet
763     if(DebugLevel()>=4){
764       Float_t s2y = fTracklet->GetS2Y();
765       Float_t s2z = fTracklet->GetS2Z();
766       (*DebugStream()) << "MCtracklet"
767         << "rc="    << rc
768         << "x="     << x
769         << "y="     << y
770         << "z="     << z
771         << "dydx="  << dydx
772         << "s2y="   << s2y
773         << "s2z="   << s2z
774         << "\n";
775     }
776
777     AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
778     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
779     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
780
781     H = (THnSparse*)fContainer->At(kMCcluster);
782     val[kPt]  = TMath::ATan(dydx0)*TMath::RadToDeg();
783     //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
784     Int_t row0(-1);
785     Float_t padCorr(tilt*fTracklet->GetPadLength());
786     fTracklet->ResetClusterIter(kTRUE);
787     while((c = fTracklet->NextCluster())){
788       if(row0<0) row0 = c->GetPadRow();
789       x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
790       y = c->GetY()  + padCorr*(c->GetPadRow() - row0);
791       z = c->GetZ();
792       dx = x0 - x; 
793       ymc= y0 - dx*dydx0;
794       zmc= z0 - dx*dzdx0;
795       dy = y - ymc;
796       dz = z - zmc;
797       val[kYrez] = dy - dz*tilt;
798       val[kPrez] = dx;
799       val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t  q(TMath::Abs(c->GetQ()));
800       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
801       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
802       if(ic) val[kZrez] /= (ic*q);
803       if(H) H->Fill(val);
804
805
806       // Fill calibration container
807       Float_t d = zr0 - zmc;
808       d -= ((Int_t)(2 * d)) / 2.0;
809       if (d > 0.25) d  = 0.5 - d;
810       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
811       clInfo->SetCluster(c);
812       clInfo->SetMC(pdg, label);
813       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
814       clInfo->SetResolution(dy);
815       clInfo->SetAnisochronity(d);
816       clInfo->SetDriftLength(dx);
817       clInfo->SetTilt(tilt);
818       if(fMCcl) fMCcl->Add(clInfo);
819       else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
820       if(DebugLevel()>=5){ 
821         if(!clInfoArr){ 
822           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
823           clInfoArr->SetOwner(kFALSE);
824         }
825         clInfoArr->Add(clInfo);
826       }
827     }
828     // Fill Debug Tree
829     if(DebugLevel()>=5 && clInfoArr){
830       (*DebugStream()) << "MCcluster"
831         <<"clInfo.=" << clInfoArr
832         << "\n";
833       clInfoArr->Clear();
834     }
835   }
836   if(clInfoArr) delete clInfoArr;
837   return h;
838 }
839
840
841 //__________________________________________________________________________
842 Int_t AliTRDresolution::GetPtBin(Float_t pt)
843 {
844 // Find pt bin according to local pt segmentation
845   Int_t ipt(-1);
846   while(ipt<AliTRDresolution::kNpt){
847     if(pt<fgPtBin[ipt+1]) break;
848     ipt++;
849   }
850   return ipt;
851 }
852
853 //________________________________________________________
854 Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)
855 {
856 // return mean number of entries/bin of histogram "h"
857 // if option "opt" is given the following values are accepted:
858 //   "<" : consider only entries less than "cut"
859 //   ">" : consider only entries greater than "cut"
860
861   //Int_t dim(h->GetDimension());
862   Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());
863   Double_t sum(0.); Int_t n(0);
864   for(Int_t ix(1); ix<=nbx; ix++)
865     for(Int_t iy(1); iy<=nby; iy++)
866       for(Int_t iz(1); iz<=nbz; iz++){
867         if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}
868         else{
869           if(strcmp(opt, "<")==0) {
870             if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
871           } else if(strcmp(opt, ">")==0){
872             if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
873           } else {sum += h->GetBinContent(ix, iy, iz); n++;}
874         }
875       }
876   return n>0?sum/n:0.;
877 }
878
879 //________________________________________________________
880 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
881 {
882   //
883   // Get the reference figures
884   //
885
886   if(!gPad){
887     AliWarning("Please provide a canvas to draw results.");
888     return kFALSE;
889   }
890 /*  Int_t selection[100], n(0), selStart(0); // 
891   Int_t ly0(0), dly(5);
892   TList *l(NULL); TVirtualPad *pad(NULL); */
893   switch(ifig){
894   case 0:
895     break;
896   }
897   AliWarning(Form("Reference plot [%d] missing result", ifig));
898   return kFALSE;
899 }
900
901
902 //________________________________________________________
903 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
904 {
905 // Build pt segments
906   for(Int_t j(0); j<=kNpt; j++){
907     pt0+=(TMath::Exp(j*j*dpt)-1.);
908     fgPtBin[j]=pt0;
909   }
910 }
911
912 //________________________________________________________
913 void AliTRDresolution::MakeSummary()
914 {
915 // Build summary plots
916
917   if(!fProj){
918     AliError("Missing results");
919     return;
920   }  
921   TVirtualPad *p(NULL); TCanvas *cOut(NULL);
922   TObjArray *arr(NULL); TH2 *h2(NULL);
923
924   // cluster resolution
925   // define palette
926   gStyle->SetPalette(1);
927   const Int_t nClViews(11);
928   const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh", "ClY", "ClYn"};
929   const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0};
930   const Int_t nTrkltViews(10);
931   const Char_t *vTrkltName[nTrkltViews] = {"TrkltY", "TrkltYn", "TrkltYp", "TrkltPhn", "TrkltPhp", "TrkltZ", "TrkltQn", "TrkltQp", "TrkltPn", "TrkltPp"};
932   const Int_t nTrkInViews(6);
933   const Char_t *vTrkInName[nTrkInViews] = {"TrkInY", "TrkInYn", "TrkInYp", "TrkInZ", "TrkInPhn", "TrkInPhp"};
934   const Int_t nTrkViews(10);
935   const Char_t *vTrkName[nTrkViews] = {"TrkY", "TrkYn", "TrkYp", "TrkPhn", "TrkPhp", "TrkZ", "TrkQn", "TrkQp", "TrkPn", "TrkPp"};
936   const Char_t *typName[] = {"", "MC"};
937
938   for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
939     if((arr = (TObjArray*)fProj->At(ityp?kMCcluster:kCluster))){
940       for(Int_t iview(0); iview<nClViews; iview++){
941         cOut = new TCanvas(Form("TRDsummary%s_%sCl%02d", GetName(), typName[ityp], iview), "Cluster Resolution", 1024, 768);
942         cOut->Divide(3,2, 1.e-5, 1.e-5);
943         Int_t nplot(0);
944         for(Int_t iplot(0); iplot<6; iplot++){
945           p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
946           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], iplot)))) continue;
947           nplot++;
948           if(vClOpt[iview]==0) h2->Draw("colz");
949           else if(vClOpt[iview]==1) DrawSigma(h2, 1.e4, 2.e2, 6.5e2, "#sigma(#Deltay) [#mum]");
950         }
951         if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
952         else delete cOut;
953       }
954     }
955     // tracklet systematic
956     if((arr = (TObjArray*)fProj->At(ityp?kMCtracklet:kTracklet))){
957       for(Int_t iview(0); iview<nTrkltViews; iview++){
958         cOut = new TCanvas(Form("TRDsummary%s_%sTrklt%02d", GetName(), typName[ityp], iview), "Tracklet Resolution", 1024, 768);
959         cOut->Divide(3,2, 1.e-5, 1.e-5);
960         Int_t nplot(0);
961         for(Int_t iplot(0); iplot<6; iplot++){
962           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
963           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
964           h2->Draw("colz"); nplot++;
965         }
966         if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
967         else delete cOut;
968       }
969     }
970     // trackIn systematic
971     if((arr = (TObjArray*)fProj->At(ityp?kMCtrackIn:kTrackIn))){
972       cOut = new TCanvas(Form("TRDsummary%s_%sTrkIn", GetName(), typName[ityp]), "Track IN Resolution", 1024, 768);
973       cOut->Divide(3,2, 1.e-5, 1.e-5);
974       Int_t nplot(0);
975       for(Int_t iplot(0); iplot<nTrkInViews; iplot++){
976         p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
977           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iplot])))) continue;
978           h2->Draw("colz"); nplot++;
979       }
980       if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
981       else delete cOut;
982     }
983   }
984   // track MC systematic
985   if((arr = (TObjArray*)fProj->At(kMCtrack))) {
986     for(Int_t iview(0); iview<nTrkViews; iview++){
987       cOut = new TCanvas(Form("TRDsummary%s_MCTrk%02d", GetName(), iview), "Track Resolution", 1024, 768);
988       cOut->Divide(3,2, 1.e-5, 1.e-5);
989       Int_t nplot(0);
990       for(Int_t iplot(0); iplot<6; iplot++){
991         p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
992         if(!(h2 = (TH2*)arr->FindObject(Form("H%s%d_2D", vTrkName[iview], iplot)))) continue;
993         h2->Draw("colz"); nplot++;
994       }
995       if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
996       else delete cOut;
997     }
998   }
999
1000
1001   gStyle->SetPalette(1);
1002 }
1003
1004 //________________________________________________________
1005 void AliTRDresolution::DrawSigma(TH2 *h2, Float_t scale, Float_t m, Float_t M, const Char_t *title)
1006 {
1007   // Draw error bars scaled with "scale" instead of content values
1008   //use range [m,M] if limits are specified
1009
1010   if(!h2) return;
1011   TH2 *h2e = (TH2F*)h2->Clone(Form("%s_E", h2->GetName()));
1012   h2e->SetContour(10);
1013   if(M>m) h2e->GetZaxis()->SetRangeUser(m, M);
1014   if(title) h2e->GetZaxis()->SetTitle(title);
1015   
1016   for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1017     for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1018       if(h2->GetBinContent(ix, iy)<-100.) continue;
1019       Float_t v(scale*h2->GetBinError(ix, iy));
1020       if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1021       h2e->SetBinContent(ix, iy, v);
1022     }
1023   }
1024   h2e->Draw("colz");
1025 }
1026
1027 //________________________________________________________
1028 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1029 {
1030 // Returns the range of the bulk of data in histogram h2. Removes outliers. 
1031 // The "range" vector should be initialized with 2 elements
1032 // Option "mod" can be any of
1033 //   - 0 : gaussian like distribution 
1034 //   - 1 : tailed distribution 
1035
1036   Int_t nx(h2->GetNbinsX())
1037       , ny(h2->GetNbinsY())
1038       , n(nx*ny);
1039   Double_t *data=new Double_t[n];
1040   for(Int_t ix(1), in(0); ix<=nx; ix++){
1041     for(Int_t iy(1); iy<=ny; iy++)
1042       data[in++] = h2->GetBinContent(ix, iy);
1043   }
1044   Double_t mean, sigm;
1045   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1046
1047   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1048   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]); 
1049   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1050   TH1S h1("h1SF0", "", 100, range[0], range[1]);
1051   h1.FillN(n,data,0);
1052   delete [] data;
1053  
1054   switch(mod){
1055   case 0:// gaussian distribution  
1056   {
1057     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1058     h1.Fit(&fg, "QN");
1059     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1060     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1061     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
1062     break;
1063   }
1064   case 1:// tailed distribution  
1065   {  
1066     Int_t bmax(h1.GetMaximumBin());
1067     Int_t jBinMin(1), jBinMax(100);
1068     for(Int_t ibin(bmax); ibin--;){
1069       if(h1.GetBinContent(ibin)<1.){
1070         jBinMin=ibin; break;
1071       }
1072     }
1073     for(Int_t ibin(bmax); ibin++;){
1074       if(h1.GetBinContent(ibin)<1.){
1075         jBinMax=ibin; break;
1076       }
1077     }
1078     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1079     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
1080     break;
1081   }
1082   }
1083
1084   return;
1085 }
1086
1087
1088 //________________________________________________________
1089 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1090 {
1091 // Analyse cluster
1092   const Int_t kNcontours(9);
1093   const Int_t kNstat(300);
1094   Int_t cidx=mc?kMCcluster:kCluster;
1095   if(fProj && fProj->At(cidx)) return kTRUE;
1096   if(!fContainer){
1097     AliError("Missing data container.");
1098     return kFALSE;
1099   }
1100   THnSparse *H(NULL);
1101   if(!(H = (THnSparse*)fContainer->At(cidx))){
1102     AliError(Form("Missing/Wrong data @ %d.", cidx));
1103     return kFALSE;
1104   }
1105   Int_t ndim(H->GetNdimensions());
1106   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1107   TAxis *aa[kNdim], *as(NULL), *apt(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1108   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1109   if(ndim > kPt) apt = H->GetAxis(kPt);
1110   if(ndim > kSpeciesChgRC) as  = H->GetAxis(kSpeciesChgRC);
1111   // build list of projections
1112   const Int_t nsel(12), npsel(5);
1113   // define rebinning strategy
1114   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1115   AliTRDresolutionProjection hp[fgkNproj[cidx]], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1116   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1117   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1118     isel++; // new selection
1119     hp[ih].Build(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1120     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1121       php[isel][np[isel]++] = &hp[ih++];
1122     hp[ih].Build(Form("H%sClYn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1123     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1124       php[isel][np[isel]++] = &hp[ih++];
1125     hp[ih].Build(Form("H%sClQn%d", mc?"MC":"", ily), Form("Clusters[-]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1126     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1127     hp[ih].SetShowRange(20., 40.);
1128       php[isel][np[isel]++] = &hp[ih++];
1129     hp[ih].Build(Form("H%sClYXTCn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1130 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1131       php[isel][np[isel]++] = &hp[ih++];
1132     hp[ih].Build(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: r-#phi(x,#Phi) residuals ly%d", ily), kPrez, kPt, kYrez, aa);
1133 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1134       php[isel][np[isel]++] = &hp[ih++];
1135     isel++; // new selection
1136       php[isel][np[isel]++] = &hp[ih-5]; // relink HClY
1137     hp[ih].Build(Form("H%sClYp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1138     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1139       php[isel][np[isel]++] = &hp[ih++];
1140     hp[ih].Build(Form("H%sClQp%d", mc?"MC":"", ily), Form("Clusters[+]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1141     hp[ih].SetShowRange(20., 40.);
1142     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1143       php[isel][np[isel]++] = &hp[ih++];
1144     hp[ih].Build(Form("H%sClYXTCp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1145 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1146       php[isel][np[isel]++] = &hp[ih++];
1147       php[isel][np[isel]++] = &hp[ih-4]; // relink HClYXPh
1148   }
1149
1150   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1);
1151   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1152     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1153     ly = coord[kBC]-1;
1154     // RC selection
1155     if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1156
1157     // charge selection
1158     ch = 0; // [-] track
1159     if(chBin>0 && coord[kPt] > chBin) ch = 1;  // [+] track
1160
1161     isel = ly*2+ch;
1162     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1163   }
1164
1165   if(!fProj){
1166     AliInfo("Building array of projections ...");
1167     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1168   }
1169   TObjArray *arr(NULL);
1170   fProj->AddAt(arr = new TObjArray(fgkNproj[cidx]), cidx);
1171
1172   TH2 *h2(NULL);
1173   for(; ih--; ){
1174     if(!hp[ih].fH) continue;
1175     Int_t mid(1), nstat(kNstat);
1176     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; /*nstat=300;*/}
1177     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1178     arr->AddAt(h2, ih);
1179   }
1180
1181   return kTRUE;
1182 }
1183
1184 //________________________________________________________
1185 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1186 {
1187 // Analyse tracklet
1188   const Int_t kNcontours(9);
1189   const Int_t kNstat(100);
1190   Int_t cidx=mc?kMCtracklet:kTracklet;
1191   if(fProj && fProj->At(cidx)) return kTRUE;
1192   if(!fContainer){
1193     AliError("Missing data container.");
1194     return kFALSE;
1195   }
1196   THnSparse *H(NULL);
1197   if(!(H = (THnSparse*)fContainer->At(cidx))){
1198     AliError(Form("Missing/Wrong data @ %d.", cidx));
1199     return kFALSE;
1200   }
1201   Int_t ndim(H->GetNdimensions());
1202   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
1203   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1204   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1205   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1206   // build list of projections
1207   const Int_t nsel(18), npsel(6);
1208   // define rebinning strategy
1209   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1210   AliTRDresolutionProjection hp[fgkNproj[cidx]], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1211   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1212   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1213     isel++; // new selection
1214     hp[ih].Build(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets   :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1215     hp[ih].SetShowRange(-0.03,0.03);
1216     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1217       php[isel][np[isel]++] = &hp[ih++];
1218     hp[ih].Build(Form("H%sTrkltYn%d", mc?"MC":"", ily), Form("Tracklets[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1219     hp[ih].SetShowRange(-0.03,0.03);
1220     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1221       php[isel][np[isel]++] = &hp[ih++];
1222     hp[ih].Build(Form("H%sTrkltPhn%d", mc?"MC":"", ily), Form("Tracklets[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1223     hp[ih].SetShowRange(-0.5,0.5);
1224     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1225       php[isel][np[isel]++] = &hp[ih++];
1226     hp[ih].Build(Form("H%sTrkltPn%d", mc?"MC":"", ily), Form("Tracklets[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1227     hp[ih].SetShowRange(6.,12.);
1228     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1229       php[isel][np[isel]++] = &hp[ih++];
1230     hp[ih].Build(Form("H%sTrkltYPn%d", mc?"MC":"", ily), Form("Tracklets[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1231       php[isel][np[isel]++] = &hp[ih++];
1232     hp[ih].Build(Form("H%sTrkltQn%d", mc?"MC":"", ily), Form("Tracklets[-]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1233     hp[ih].SetShowRange(700.,1100.);
1234     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1235       php[isel][np[isel]++] = &hp[ih++];
1236     isel++; // new selection
1237     php[isel][np[isel]++] = &hp[ih-6]; // relink first histo
1238     hp[ih].Build(Form("H%sTrkltYp%d", mc?"MC":"", ily), Form("Tracklets[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1239     hp[ih].SetShowRange(-0.03,0.03);
1240     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1241       php[isel][np[isel]++] = &hp[ih++];
1242     hp[ih].Build(Form("H%sTrkltPhp%d", mc?"MC":"", ily), Form("Tracklets[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1243     hp[ih].SetShowRange(-0.5,0.5);
1244     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1245       php[isel][np[isel]++] = &hp[ih++];
1246     hp[ih].Build(Form("H%sTrkltPp%d", mc?"MC":"", ily), Form("Tracklets[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1247     hp[ih].SetShowRange(6.,12.);
1248     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1249       php[isel][np[isel]++] = &hp[ih++];
1250     hp[ih].Build(Form("H%sTrkltYPp%d", mc?"MC":"", ily), Form("Tracklets[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1251       php[isel][np[isel]++] = &hp[ih++];
1252     hp[ih].Build(Form("H%sTrkltQp%d", mc?"MC":"", ily), Form("Tracklets[+]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1253     hp[ih].SetShowRange(700.,1100.);
1254     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1255       php[isel][np[isel]++] = &hp[ih++];
1256     isel++; // new selection
1257     hp[ih].Build(Form("H%sTrkltZ%d", mc?"MC":"", ily), Form("Tracklets[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa);
1258     hp[ih].SetShowRange(-0.1,0.1);
1259     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1260       php[isel][np[isel]++] = &hp[ih++];
1261   }
1262
1263   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1);
1264   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1265     v = H->GetBinContent(ib, coord);
1266     if(v<1.) continue;
1267     ly = coord[kBC]-1; // layer selection
1268     // charge selection
1269     ch = 0; // [-] track
1270     if(rcBin>0){ // debug mode in which species are also saved
1271       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1272       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1273     }
1274     isel = ly*3+ch;
1275     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1276   }
1277
1278   if(!fProj){
1279     AliInfo("Building array of projections ...");
1280     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1281   }
1282   TObjArray *arr(NULL);
1283   fProj->AddAt(arr = new TObjArray(fgkNproj[cidx]), cidx);
1284
1285   TH2 *h2(NULL);
1286   for(; ih--; ){
1287     if(!hp[ih].fH) continue;
1288     Int_t mid(0), nstat(kNstat);
1289     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; /*nstat=300;*/}
1290     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1291     arr->AddAt(h2, ih);
1292   }
1293   return kTRUE;
1294 }
1295
1296 //________________________________________________________
1297 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
1298 {
1299 // Analyse track in
1300
1301   const Int_t kNcontours(9);
1302   const Int_t kNstat(30);
1303   Int_t cidx=mc?kMCtrackIn:kTrackIn;
1304   if(fProj && fProj->At(cidx)) return kTRUE;
1305   if(!fContainer){
1306     AliError("Missing data container.");
1307     return kFALSE;
1308   }
1309   THnSparse *H(NULL);
1310   if(!(H = (THnSparse*)fContainer->At(cidx))){
1311     AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1312     return kFALSE;
1313   }
1314
1315   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1316   Int_t ndim(H->GetNdimensions());
1317   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1318   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1319   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1320   // build list of projections
1321   const Int_t nsel(3), npsel(4);
1322   // define rebinning strategy
1323   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1324   const Int_t nPtPhi(2); Int_t rebinPtPhiX[nEtaPhi] = {1, 1}, rebinPtPhiY[nEtaPhi] = {2, 5};
1325   AliTRDresolutionProjection hp[fgkNproj[cidx]], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1326   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1327   // define list of projections
1328   isel++;  // negative tracks
1329   hp[ih].Build(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: r-#phi residuals", kEta, kPhi, kYrez, aa);
1330   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1331     php[isel][np[isel]++] = &hp[ih++];
1332   hp[ih].Build(Form("H%sTrkInYn", mc?"MC":""), "TrackIn[-]:: r-#phi residuals", kEta, kPhi, kYrez, aa);
1333   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1334     php[isel][np[isel]++] = &hp[ih++];
1335   hp[ih].Build(Form("H%sTrkInPhn", mc?"MC":""), "TrackIn[-]:: #Delta#phi residuals", kEta, kPhi, kPrez, aa);
1336   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1337     php[isel][np[isel]++] = &hp[ih++];
1338   hp[ih].Build(Form("H%sTrkInYPn", mc?"MC":""), "TrackIn[-]:: r-#phi/p_{t} residuals", kPt, kPhi, kYrez, aa);
1339   hp[ih].SetRebinStrategy(nPtPhi, rebinPtPhiX, rebinPtPhiY);
1340     php[isel][np[isel]++] = &hp[ih++];
1341   isel++; // positive tracks
1342   php[isel][np[isel]++] = &hp[ih-4]; // relink first histo
1343   hp[ih].Build(Form("H%sTrkInYp", mc?"MC":""), "TrackIn[+]:: r-#phi residuals", kEta, kPhi, kYrez, aa);
1344   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1345     php[isel][np[isel]++] = &hp[ih++];
1346   hp[ih].Build(Form("H%sTrkInPhp", mc?"MC":""), "TrackIn[+]:: #Delta#phi residuals", kEta, kPhi, kPrez, aa);
1347   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1348     php[isel][np[isel]++] = &hp[ih++];
1349   hp[ih].Build(Form("H%sTrkInYPp", mc?"MC":""), "TrackIn[+]:: r-#phi/p_{t} residuals", kPt, kPhi, kYrez, aa);
1350   hp[ih].SetRebinStrategy(nPtPhi, rebinPtPhiX, rebinPtPhiY);
1351     php[isel][np[isel]++] = &hp[ih++];
1352   isel++; // RC tracks
1353   hp[ih].Build(Form("H%sTrkInZ", mc?"MC":""), "TrackIn[RC]:: z residuals", kEta, kPhi, kZrez, aa);
1354   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1355     php[isel][np[isel]++] = &hp[ih++];
1356
1357   // fill projections
1358   Int_t ch(0), rcBin(as?as->FindBin(0.):-1);
1359   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1360     v = H->GetBinContent(ib, coord);
1361     if(v<1.) continue;
1362     if(coord[kBC]>1) continue; // bunch cross cut
1363     // charge selection
1364     ch = 0; // [-] track
1365     if(rcBin>0){ // debug mode in which species are also saved
1366       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1367       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1368     }
1369     for(Int_t jh(0); jh<np[ch]; jh++) php[ch][jh]->Increment(coord, v);
1370   }
1371   if(!fProj){
1372     AliInfo("Building array of projections ...");
1373     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1374   }
1375   TObjArray *arr(NULL);
1376   fProj->AddAt(arr = new TObjArray(fgkNproj[cidx]), cidx);
1377
1378   TH2 *h2(NULL);
1379   for(; ih--; ){
1380     if(!hp[ih].fH) continue;
1381     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours/*, mid*/))) continue;
1382     arr->AddAt(h2, ih);
1383   }
1384   return kTRUE;
1385 }
1386
1387
1388 //________________________________________________________
1389 Bool_t AliTRDresolution::MakeProjectionTrack()
1390 {
1391 // Analyse tracklet
1392   const Int_t kNcontours(9);
1393   const Int_t kNstat(100);
1394   Int_t cidx(kMCtrack);
1395   if(fProj && fProj->At(cidx)) return kTRUE;
1396   if(!fContainer){
1397     AliError("Missing data container.");
1398     return kFALSE;
1399   }
1400   THnSparse *H(NULL);
1401   if(!(H = (THnSparse*)fContainer->At(cidx))){
1402     AliError(Form("Missing/Wrong data @ %d.", cidx));
1403     return kFALSE;
1404   }
1405   Int_t ndim(H->GetNdimensions());
1406   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
1407   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1408   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1409   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1410   // build list of projections
1411   const Int_t nsel(18), npsel(6);
1412   // define rebinning strategy
1413   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1414   AliTRDresolutionProjection hp[fgkNproj[cidx]], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1415   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1416   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1417     isel++; // new selection
1418     hp[ih].Build(Form("HTrkY%d", ily), Form("Tracks   :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1419     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1420       php[isel][np[isel]++] = &hp[ih++];
1421     hp[ih].Build(Form("HTrkYn%d", ily), Form("Tracks[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1422     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1423       php[isel][np[isel]++] = &hp[ih++];
1424     hp[ih].Build(Form("HTrkPhn%d", ily), Form("Tracks[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1425     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1426       php[isel][np[isel]++] = &hp[ih++];
1427     hp[ih].Build(Form("HTrkPn%d", ily), Form("Tracks[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1428     hp[ih].SetShowRange(6.,12.);
1429     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1430       php[isel][np[isel]++] = &hp[ih++];
1431     hp[ih].Build(Form("HTrkYPn%d", ily), Form("Tracks[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1432       php[isel][np[isel]++] = &hp[ih++];
1433     hp[ih].Build(Form("HTrkQn%d", ily), Form("Tracks[-]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa);
1434     //hp[ih].SetShowRange(700.,1100.);
1435     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1436       php[isel][np[isel]++] = &hp[ih++];
1437     isel++; // new selection
1438     php[isel][np[isel]++] = &hp[ih-6]; // relink first histo
1439     hp[ih].Build(Form("HTrkYp%d", ily), Form("Tracks[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1440     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1441       php[isel][np[isel]++] = &hp[ih++];
1442     hp[ih].Build(Form("HTrkPhp%d", ily), Form("Tracks[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1443     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1444       php[isel][np[isel]++] = &hp[ih++];
1445     hp[ih].Build(Form("HTrkPp%d", ily), Form("Tracks[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1446     hp[ih].SetShowRange(6.,12.);
1447     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1448       php[isel][np[isel]++] = &hp[ih++];
1449     hp[ih].Build(Form("HTrkYPp%d", ily), Form("Tracks[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1450       php[isel][np[isel]++] = &hp[ih++];
1451     hp[ih].Build(Form("HTrkQp%d", ily), Form("Tracks[+]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa);
1452     //hp[ih].SetShowRange(700.,1100.);
1453     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1454       php[isel][np[isel]++] = &hp[ih++];
1455     isel++; // new selection
1456     hp[ih].Build(Form("HTrkZ%d", ily), Form("Tracks[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa);
1457     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1458       php[isel][np[isel]++] = &hp[ih++];
1459   }
1460
1461   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1);
1462   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1463     v = H->GetBinContent(ib, coord);
1464     if(v<1.) continue;
1465     ly = coord[kBC]-1; // layer selection
1466     // charge selection
1467     ch = 0; // [-] track
1468     if(rcBin>0){ // debug mode in which species are also saved
1469       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1470       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1471     }
1472     isel = ly*3+ch;
1473     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1474   }
1475
1476   if(!fProj){
1477     AliInfo("Building array of projections ...");
1478     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1479   }
1480   TObjArray *arr(NULL);
1481   fProj->AddAt(arr = new TObjArray(fgkNproj[cidx]), cidx);
1482
1483   TH2 *h2(NULL);
1484   for(; ih--; ){
1485     if(!hp[ih].fH) continue;
1486     Int_t mid(0), nstat(kNstat);
1487     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; /*nstat=300;*/}
1488     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1489     arr->AddAt(h2, ih);
1490   }
1491   return kTRUE;
1492 }
1493
1494 //________________________________________________________
1495 Bool_t AliTRDresolution::PostProcess()
1496 {
1497 // Fit, Project, Combine, Extract values from the containers filled during execution
1498
1499   if (!fContainer) {
1500     AliError("ERROR: list not available");
1501     return kFALSE;
1502   }
1503
1504   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1505   // Clusters residuals
1506   if(!MakeProjectionCluster()) return kFALSE;
1507   fNRefFigures = 3;
1508   // Tracklet residual/pulls
1509   if(!MakeProjectionTracklet()) return kFALSE;
1510   fNRefFigures = 7;
1511   // TRDin residual/pulls
1512   if(!MakeProjectionTrackIn()) return kFALSE;
1513   fNRefFigures = 11;
1514
1515   if(!HasMCdata()) return kTRUE;
1516   //PROCESS MC RESIDUAL DISTRIBUTIONS
1517
1518   // CLUSTER Y RESOLUTION/PULLS
1519   if(!MakeProjectionCluster(kTRUE)) return kFALSE;
1520   fNRefFigures = 17;
1521
1522   // TRACKLET RESOLUTION/PULLS
1523   if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
1524   fNRefFigures = 21;
1525
1526   // TRACK RESOLUTION/PULLS
1527   if(!MakeProjectionTrack()) return kFALSE;
1528   fNRefFigures+=16;
1529
1530   // TRACK TRDin RESOLUTION/PULLS
1531   if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
1532   fNRefFigures+=8;
1533
1534   return kTRUE;
1535 }
1536
1537
1538 //________________________________________________________
1539 void AliTRDresolution::Terminate(Option_t *opt)
1540 {
1541   AliTRDrecoTask::Terminate(opt);
1542   if(HasPostProcess()) PostProcess();
1543 }
1544
1545 //________________________________________________________
1546 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1547 {
1548 // Helper function to avoid duplication of code
1549 // Make first guesses on the fit parameters
1550
1551   // find the intial parameters of the fit !! (thanks George)
1552   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1553   Double_t sum = 0.;
1554   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1555   f->SetParLimits(0, 0., 3.*sum);
1556   f->SetParameter(0, .9*sum);
1557   Double_t rms = h->GetRMS();
1558   f->SetParLimits(1, -rms, rms);
1559   f->SetParameter(1, h->GetMean());
1560
1561   f->SetParLimits(2, 0., 2.*rms);
1562   f->SetParameter(2, rms);
1563   if(f->GetNpar() <= 4) return;
1564
1565   f->SetParLimits(3, 0., sum);
1566   f->SetParameter(3, .1*sum);
1567
1568   f->SetParLimits(4, -.3, .3);
1569   f->SetParameter(4, 0.);
1570
1571   f->SetParLimits(5, 0., 1.e2);
1572   f->SetParameter(5, 2.e-1);
1573 }
1574
1575 //________________________________________________________
1576 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
1577 {
1578 // Build performance histograms for AliTRDcluster.vs TRD track or MC
1579 //  - y reziduals/pulls
1580
1581   TObjArray *arr = new TObjArray(2);
1582   arr->SetName(name); arr->SetOwner();
1583   TH1 *h(NULL); char hname[100], htitle[300];
1584
1585   // tracklet resolution/pull in y direction
1586   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
1587   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
1588   Float_t rr = range<0.?fDyRange:range;
1589   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1590     Int_t nybins=50;
1591     if(expand) nybins*=2;
1592     h = new TH3S(hname, htitle, 
1593                  48, -.48, .48,            // phi
1594                  60, -rr, rr,              // dy
1595                  nybins, -0.5, nybins-0.5);// segment
1596   } else h->Reset();
1597   arr->AddAt(h, 0);
1598   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
1599   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
1600   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1601     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
1602   } else h->Reset();
1603   arr->AddAt(h, 1);
1604
1605   return arr;
1606 }
1607
1608 //________________________________________________________
1609 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
1610 {
1611 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1612 //  - y reziduals/pulls
1613 //  - z reziduals/pulls
1614 //  - phi reziduals
1615   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05); 
1616   arr->Expand(5);
1617   TH1 *h(NULL); char hname[100], htitle[300];
1618
1619   // tracklet resolution/pull in z direction
1620   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
1621   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1622   if(!(h = (TH2S*)gROOT->FindObject(hname))){
1623     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
1624   } else h->Reset();
1625   arr->AddAt(h, 2);
1626   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1627   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
1628   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1629     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
1630     h->GetZaxis()->SetBinLabel(1, "no RC");
1631     h->GetZaxis()->SetBinLabel(2, "RC");
1632   } else h->Reset();
1633   arr->AddAt(h, 3);
1634
1635   // tracklet to track phi resolution
1636   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
1637   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1638   Int_t nsgms=540;
1639   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1640     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
1641   } else h->Reset();
1642   arr->AddAt(h, 4);
1643
1644   return arr;
1645 }
1646
1647 //________________________________________________________
1648 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1649 {
1650 // Build performance histograms for AliExternalTrackParam.vs MC
1651 //  - y resolution/pulls
1652 //  - z resolution/pulls
1653 //  - phi resolution, snp pulls
1654 //  - theta resolution, tgl pulls
1655 //  - pt resolution, 1/pt pulls, p resolution
1656
1657   TObjArray *arr = BuildMonitorContainerTracklet(name); 
1658   arr->Expand(11);
1659   TH1 *h(NULL); char hname[100], htitle[300];
1660   //TAxis *ax(NULL);
1661
1662   // snp pulls
1663   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1664   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
1665   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1666     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1667   } else h->Reset();
1668   arr->AddAt(h, 5);
1669
1670   // theta resolution
1671   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1672   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
1673   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1674     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1675   } else h->Reset();
1676   arr->AddAt(h, 6);
1677   // tgl pulls
1678   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1679   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
1680   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1681     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1682   } else h->Reset();
1683   arr->AddAt(h, 7);
1684   
1685   const Int_t kNdpt(150); 
1686   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1687   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
1688   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
1689   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1690   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1691   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
1692
1693   // Pt resolution
1694   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1695   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
1696   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1697     h = new TH3S(hname, htitle, 
1698                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1699     //ax = h->GetZaxis();
1700     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1701   } else h->Reset();
1702   arr->AddAt(h, 8);
1703   // 1/Pt pulls
1704   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1705   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);
1706   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1707     h = new TH3S(hname, htitle, 
1708                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1709     //ax = h->GetZaxis();
1710     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1711   } else h->Reset();
1712   arr->AddAt(h, 9);
1713   // P resolution
1714   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1715   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
1716   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1717     h = new TH3S(hname, htitle, 
1718                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1719     //ax = h->GetZaxis();
1720     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1721   } else h->Reset();
1722   arr->AddAt(h, 10);
1723
1724   return arr;
1725 }
1726
1727
1728 //________________________________________________________
1729 TObjArray* AliTRDresolution::Histos()
1730 {
1731   //
1732   // Define histograms
1733   //
1734
1735   if(fContainer) return fContainer;
1736
1737   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1738   THnSparse *H(NULL);
1739   const Int_t nhn(100); Char_t hn[nhn]; TString st;
1740
1741   //++++++++++++++++++++++
1742   // cluster to tracklet residuals/pulls
1743   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1744   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1745     const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]"};
1746     const Int_t clNbins[kNdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 10, 30, 15};
1747     const Double_t clMin[kNdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45},
1748                    clMax[kNdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45};
1749     st = "cluster spatial&charge resolution;";
1750     // define minimum info to be saved in non debug mode
1751     Int_t ndim=DebugLevel()>=1?kNdim:4;
1752     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
1753     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
1754   } else H->Reset();
1755   fContainer->AddAt(H, kCluster);
1756   //++++++++++++++++++++++
1757   // tracklet to TRD track
1758   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1759   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1760     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1761     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1762     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1763     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1764     // set specific fields
1765     trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
1766     trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
1767     trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
1768     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1769     trTitle[kNdim]=StrDup("dq/dl [a.u.]"); trNbins[kNdim] = 30; trMin[kNdim] = 100.; trMax[kNdim] = 3100;
1770
1771     st = "tracklet spatial&charge resolution;";
1772     // define minimum info to be saved in non debug mode
1773     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
1774     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1775     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1776   } else H->Reset();
1777   fContainer->AddAt(H, kTracklet);
1778   //++++++++++++++++++++++
1779   // tracklet to TRDin
1780   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1781   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1782     st = "r-#phi/z/angular residuals @ TRD entry;";
1783     // define minimum info to be saved in non debug mode
1784     Int_t ndim=DebugLevel()>=1?kNdim:7;
1785     for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
1786     H = new THnSparseI(hn, st.Data(), ndim, fgkNbins, fgkMin, fgkMax);
1787   } else H->Reset();
1788   fContainer->AddAt(H, kTrackIn);
1789   // tracklet to TRDout
1790 //  fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1791
1792
1793   // Resolution histos
1794   if(!HasMCdata()) return fContainer;
1795
1796   //++++++++++++++++++++++
1797   // cluster to TrackRef residuals/pulls
1798   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
1799   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1800     const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", fgkTitle[kSpeciesChgRC], "#Phi [deg]"};
1801     const Int_t clNbins[kNdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 20, 10, fgkNbins[kSpeciesChgRC], 15};
1802     const Double_t clMin[kNdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, fgkMin[kSpeciesChgRC], -45},
1803                    clMax[kNdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, fgkMax[kSpeciesChgRC], 45};
1804     st = "MC cluster spatial resolution;";
1805     // define minimum info to be saved in non debug mode
1806     Int_t ndim=DebugLevel()>=1?kNdim:4;
1807     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
1808     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
1809   } else H->Reset();
1810   fContainer->AddAt(H, kMCcluster);
1811   //++++++++++++++++++++++
1812   // tracklet to TrackRef
1813   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
1814   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1815     Char_t *trTitle[kNdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1816     Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1817     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1818     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1819     // set specific fields
1820     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1821     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
1822     trMin[kPrez] = -4.5; trMax[kPrez] = -trMin[kPrez];
1823     trMin[kZrez] = -1.5; trMax[kZrez] = -trMin[kZrez];
1824
1825     st = "MC tracklet spatial resolution;";
1826     // define minimum info to be saved in non debug mode
1827     Int_t ndim=DebugLevel()>=1?kNdim:4;
1828     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1829     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1830   } else H->Reset();
1831   fContainer->AddAt(H, kMCtracklet);
1832   //++++++++++++++++++++++
1833   // TRDin to TrackRef
1834   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
1835   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1836     st = "MC r-#phi/z/angular residuals @ TRD entry;";
1837     // set specific fields
1838     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1839     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1840     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
1841     trMin[kPrez] = -2.4; trMax[kPrez] = -trMin[kPrez];
1842     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
1843     // define minimum info to be saved in non debug mode
1844     Int_t ndim=DebugLevel()>=1?kNdim:7;
1845     for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
1846     H = new THnSparseI(hn, st.Data(), ndim, fgkNbins, trMin, trMax);
1847   } else H->Reset();
1848   fContainer->AddAt(H, kMCtrackIn);
1849   //++++++++++++++++++++++
1850   // track to TrackRef
1851   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
1852   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1853     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1854     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1855     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1856     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1857     // set specific fields
1858     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1859     trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 30; trMin[kNdim] = -15.; trMax[kNdim] = 15.;
1860     trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
1861     trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
1862     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
1863
1864     st = "MC track spatial&p_{t} resolution;";
1865     // define minimum info to be saved in non debug mode
1866     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
1867     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1868     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1869   } else H->Reset();
1870   fContainer->AddAt(H, kMCtrack);
1871
1872 //   // cluster resolution
1873 //   fContainer->AddAt(BuildMonitorContainerCluster("MCcl"),  kMCcluster);
1874 //   // track resolution
1875 //   TObjArray *arr(NULL);
1876 //   fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
1877 //   arr->SetName("MCtrk");
1878 //   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
1879 //   // TRDin TRACK RESOLUTION
1880 //   fContainer->AddAt(H, kMCtrackIn);
1881 //   // TRDout TRACK RESOLUTION
1882 //   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1883
1884   return fContainer;
1885 }
1886
1887 //________________________________________________________
1888 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1889 {
1890 // Robust function to process sigma/mean for 2D plot dy(x)
1891 // For each x bin a gauss fit is performed on the y projection and the range
1892 // with the minimum chi2/ndf is choosen
1893
1894   if(!h2) {
1895     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
1896     return kFALSE;
1897   }
1898   if(!Int_t(h2->GetEntries())){
1899     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
1900     return kFALSE;
1901   }
1902   if(!g || !g[0]|| !g[1]) {
1903     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
1904     return kFALSE;
1905   }
1906   // prepare
1907   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1908   Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
1909   TF1 f("f", "gaus", ymin, ymax);
1910   Int_t n(0);
1911   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1912   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1913   TH1D *h(NULL);
1914   if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
1915   Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
1916   
1917
1918   // do actual loop
1919   Float_t chi2OverNdf(0.);
1920   for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
1921     x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
1922     ymin = ay->GetXmin(); ymax = ay->GetXmax();
1923
1924     h = h2->ProjectionY("py", ix, ix);
1925     if((n=(Int_t)h->GetEntries())<stat){
1926       if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
1927       continue;
1928     }
1929     // looking for a first order mean value
1930     f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
1931     h->Fit(&f, "QNW");
1932     chi2OverNdf = f.GetChisquare()/f.GetNDF();
1933     printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
1934     y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
1935     ey  = f.GetParError(1);
1936     sy = f.GetParameter(2); esy = f.GetParError(2);
1937 //     // looking for the best chi2/ndf value
1938 //     while(ymin<y0 && ymax>y1){
1939 //       f.SetParameter(1, y);
1940 //       f.SetParameter(2, sy);
1941 //       h->Fit(&f, "QNW", "", y0, y1);
1942 //       printf("   range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
1943 //       if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
1944 //         chi2OverNdf = f.GetChisquare()/f.GetNDF();
1945 //         y  = f.GetParameter(1); ey  = f.GetParError(1);
1946 //         sy = f.GetParameter(2); esy = f.GetParError(2);
1947 //         printf("    set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
1948 //       }
1949 //       y0-=dy; y1+=dy;
1950 //     }
1951     g[0]->SetPoint(np, x, y);
1952     g[0]->SetPointError(np, ex, ey);
1953     g[1]->SetPoint(np, x, sy);
1954     g[1]->SetPointError(np, ex, esy);
1955     np++;
1956   }
1957   return kTRUE;
1958 }
1959
1960
1961 //________________________________________________________
1962 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1963 {
1964   //
1965   // Do the processing
1966   //
1967
1968   Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
1969   Int_t n = 0;
1970   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1971   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1972   if(Int_t(h2->GetEntries())){ 
1973     AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
1974   } else {
1975     AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
1976     fIdxPlot++;
1977     return kTRUE;
1978   }
1979
1980   const Int_t kINTEGRAL=1;
1981   for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
1982     Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
1983     Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
1984     TH1D *h = h2->ProjectionY(pn, abin, bbin);
1985     if((n=(Int_t)h->GetEntries())<150){ 
1986       AliDebug(4, Form("  x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
1987       continue;
1988     }
1989     h->Fit(f, "QN");
1990     Int_t ip = g[0]->GetN();
1991     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)));
1992     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1993     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1994     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1995     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1996 /*  
1997     g[0]->SetPoint(ip, x, k*h->GetMean());
1998     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1999     g[1]->SetPoint(ip, x, k*h->GetRMS());
2000     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2001   }
2002   fIdxPlot++;
2003   return kTRUE;
2004 }
2005
2006
2007 //____________________________________________________________________
2008 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2009 {
2010 //
2011 // Fit track with a staight line using the "np" clusters stored in the array "points".
2012 // The following particularities are stored in the clusters from points:
2013 //   1. pad tilt as cluster charge
2014 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
2015 // The parameters of the straight line fit are stored in the array "param" in the following order :
2016 //     param[0] - x0 reference radial position
2017 //     param[1] - y0 reference r-phi position @ x0
2018 //     param[2] - z0 reference z position @ x0
2019 //     param[3] - slope dy/dx
2020 //     param[4] - slope dz/dx
2021 //
2022 // Attention :
2023 // Function should be used to refit tracks for B=0T
2024 //
2025
2026   if(np<40){
2027     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2028     return kFALSE;
2029   }
2030   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2031
2032   Double_t x0(0.);
2033   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
2034   x0/=Float_t(np);
2035
2036   Double_t x, y, z, dx, tilt(0.);
2037   for(Int_t ip(0); ip<np; ip++){
2038     x = points[ip].GetX(); z = points[ip].GetZ();
2039     dx = x - x0;
2040     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
2041   }
2042   if(zfitter.Eval() != 0) return kFALSE;
2043
2044   Double_t z0    = zfitter.GetParameter(0);
2045   Double_t dzdx  = zfitter.GetParameter(1);
2046   for(Int_t ip(0); ip<np; ip++){
2047     if(points[ip].GetClusterType()) continue;
2048     x    = points[ip].GetX();
2049     dx   = x - x0;
2050     y    = points[ip].GetY();
2051     z    = points[ip].GetZ();
2052     tilt = points[ip].GetCharge();
2053     y -= tilt*(-dzdx*dx + z - z0);
2054     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
2055     yfitter.AddPoint(&dx, y, 1.);
2056   }
2057   if(yfitter.Eval() != 0) return kFALSE;
2058   Double_t y0   = yfitter.GetParameter(0);
2059   Double_t dydx = yfitter.GetParameter(1);
2060
2061   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
2062   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);
2063   return kTRUE;
2064 }
2065
2066 //____________________________________________________________________
2067 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
2068 {
2069 //
2070 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
2071 // See function FitTrack for the data stored in the "clusters" array
2072
2073 // The parameters of the straight line fit are stored in the array "param" in the following order :
2074 //     par[0] - x0 reference radial position
2075 //     par[1] - y0 reference r-phi position @ x0
2076 //     par[2] - slope dy/dx
2077 //
2078 // Attention :
2079 // Function should be used to refit tracks for B=0T
2080 //
2081
2082   TLinearFitter yfitter(2, "pol1");
2083
2084   // grep data for tracklet
2085   Double_t x0(0.), x[60], y[60], dy[60];
2086   Int_t nly(0);
2087   for(Int_t ip(0); ip<np; ip++){
2088     if(points[ip].GetClusterType()) continue;
2089     if(points[ip].GetVolumeID() != ly) continue;
2090     Float_t xt(points[ip].GetX())
2091            ,yt(param[1] + param[3] * (xt - param[0]));
2092     x[nly] = xt;
2093     y[nly] = points[ip].GetY();
2094     dy[nly]= y[nly]-yt;
2095     x0    += xt;
2096     nly++;
2097   }
2098   if(nly<10){
2099     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
2100     return kFALSE;
2101   }
2102   // set radial reference for fit
2103   x0 /= Float_t(nly);
2104
2105   // find tracklet core
2106   Double_t mean(0.), sig(1.e3);
2107   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
2108
2109   // simple cluster error parameterization
2110   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
2111
2112   // fit tracklet core
2113   for(Int_t jly(0); jly<nly; jly++){
2114     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
2115     Double_t dx(x[jly]-x0);
2116     yfitter.AddPoint(&dx, y[jly], 1.);
2117   }
2118   if(yfitter.Eval() != 0) return kFALSE;
2119   par[0] = x0;
2120   par[1] = yfitter.GetParameter(0);
2121   par[2] = yfitter.GetParameter(1);
2122   return kTRUE;
2123 }
2124
2125 //____________________________________________________________________
2126 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
2127 {
2128 //
2129 // Global selection mechanism of tracksbased on cluster to fit residuals
2130 // The parameters are the same as used ni function FitTrack().
2131
2132   const Float_t kS(0.6), kM(0.2);
2133   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
2134   Float_t dy, dz, s, m;
2135   for(Int_t ip(0); ip<np; ip++){
2136     if(points[ip].GetClusterType()) continue;
2137     Float_t x0(points[ip].GetX())
2138           ,y0(param[1] + param[3] * (x0 - param[0]))
2139           ,z0(param[2] + param[4] * (x0 - param[0]));
2140     dy=points[ip].GetY() - y0; h.Fill(dy);
2141     dz=points[ip].GetZ() - z0;
2142   }
2143   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
2144   fg.SetParameter(1, 0.);
2145   fg.SetParameter(2, 2.e-2);
2146   h.Fit(&fg, "QN");
2147   m=fg.GetParameter(1); s=fg.GetParameter(2);
2148   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
2149   return kTRUE;
2150 }
2151
2152 //________________________________________________________
2153 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2154 {
2155   //
2156   // Get the most probable value and the full width half mean 
2157   // of a Landau distribution
2158   //
2159
2160   const Float_t dx = 1.;
2161   mpv = f->GetParameter(1);
2162   Float_t fx, max = f->Eval(mpv);
2163
2164   xm = mpv - dx;
2165   while((fx = f->Eval(xm))>.5*max){
2166     if(fx>max){ 
2167       max = fx;
2168       mpv = xm;
2169     }
2170     xm -= dx;
2171   }
2172
2173   xM += 2*(mpv - xm);
2174   while((fx = f->Eval(xM))>.5*max) xM += dx;
2175 }
2176
2177
2178 // #include "TFile.h"
2179 // //________________________________________________________
2180 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
2181 // {
2182 //   if(!file){
2183 //     AliWarning("Use cluster position as in reconstruction.");
2184 //     SetLoadCorrection();
2185 //     return kTRUE;
2186 //   }
2187 //   TDirectory *cwd(gDirectory);
2188 //   TString fileList;
2189 //   FILE *filePtr = fopen(file, "rt");
2190 //   if(!filePtr){
2191 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
2192 //     SetLoadCorrection();
2193 //     return kFALSE;
2194 //   }
2195 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
2196 //   while(fileList.Gets(filePtr)){
2197 //     if(!TFile::Open(fileList.Data())) {
2198 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
2199 //       continue;
2200 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
2201 // 
2202 //     TTree *tSys = (TTree*)gFile->Get("tSys");
2203 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
2204 //     tSys->Draw("det:t>>h2", "dx", "goff");
2205 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2206 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
2207 //     }
2208 //     h2->SetDirectory(cwd);
2209 //     gFile->Close();
2210 //   }
2211 //   cwd->cd();
2212 // 
2213 //   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
2214 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2215 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
2216 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2217 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
2218 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
2219 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2220 //       printf("%03d|", idet);
2221 //       fprintf(fout, "    {");
2222 //       for(Int_t it(0); it<30; it++){
2223 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
2224 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
2225 //       }
2226 //       printf("\n");
2227 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
2228 //     }
2229 //     fprintf(fout, "  };\n");
2230 //   }
2231 //   SetLoadCorrection();
2232 //   return kTRUE;
2233 // }
2234
2235 //________________________________________________________
2236 AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection()
2237   :fH(NULL)
2238   ,fNrebin(0)
2239   ,fRebinX(NULL)
2240   ,fRebinY(NULL)
2241 {
2242   // constructor
2243   memset(fAx, 0, 3*sizeof(Int_t));
2244   memset(fRange, 0, 4*sizeof(Float_t));
2245 }
2246
2247 //________________________________________________________
2248 AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection()
2249 {
2250   // destructor
2251   if(fH) delete fH;
2252 }
2253
2254 //________________________________________________________
2255 void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])
2256 {
2257 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"
2258   if(!aa[ix] || !aa[iy] || !aa[iz]) return;
2259   TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);
2260   fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),
2261     ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
2262     ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
2263     az->GetNbins(), az->GetXmin(), az->GetXmax());
2264   fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;
2265   fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;
2266 }
2267
2268 //________________________________________________________
2269 void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v)
2270 {
2271 // increment bin with value "v" pointed by general coord in "bin"
2272   if(!fH) return;
2273   fH->AddBinContent(
2274         fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), v);
2275 }
2276
2277 //________________________________________________________
2278 TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid)
2279 {
2280 // build the 2D projection and adjust binning
2281
2282   const Char_t *title[] = {"Mean", "#mu", "MPV"};
2283   if(!fH) return NULL;
2284   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());
2285   TH2 *h2s = (TH2*)fH->Project3D("yx");
2286   Int_t irebin(0), dxBin(1), dyBin(1);
2287   while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){
2288     h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);
2289     dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];
2290     irebin++;
2291   }
2292   Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
2293   if(h2s) delete h2s;
2294
2295   // start projection
2296   TH1 *h(NULL);
2297   Float_t dz=(fRange[1]-fRange[1])/ncol;
2298   TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));
2299   TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),
2300             Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], (*tokenTitle)[0]->GetName(), tokenTitle->GetEntriesFast()>1?(*tokenTitle)[1]->GetName():""),
2301             nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());
2302   h2->SetContour(ncol);
2303   h2->GetZaxis()->CenterTitle();
2304   h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);
2305   //printf("%s[%s] nx[%d] ny[%d]\n", h2->GetName(), h2->GetTitle(), nx, ny);
2306   for(Int_t iy(0); iy<ny; iy++){
2307     for(Int_t ix(0); ix<nx; ix++){
2308       h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);
2309       Int_t ne((Int_t)h->Integral());
2310       if(ne<nstat/2){
2311         h2->SetBinContent(ix+1, iy+1, -999);
2312         h2->SetBinError(ix+1, iy+1, 1.);
2313       }else{
2314         Float_t v(h->GetMean()), ve(h->GetRMS());
2315         if(mid==1){
2316           TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());
2317           fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);
2318           h->Fit(&fg, "WQ");
2319           v = fg.GetParameter(1); ve = fg.GetParameter(2);
2320         } else if (mid==2) {
2321           TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());
2322           fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);
2323           h->Fit(&fl, "WQ");
2324           v = fl.GetMaximumX(); ve = fl.GetParameter(2);
2325 /*          TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());
2326           fgle.SetParameter(0, fl.GetParameter(0));
2327           fgle.SetParameter(1, fl.GetParameter(1));
2328           fgle.SetParameter(2, fl.GetParameter(2));
2329           fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);
2330           h->Fit(&fgle, "WQ");
2331           v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/
2332         }
2333         if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);
2334         else h2->SetBinContent(ix+1, iy+1, v);
2335         h2->SetBinError(ix+1, iy+1, ve);
2336       }
2337     }
2338   }
2339   if(h) delete h;
2340   return h2;
2341 }
2342
2343 void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])
2344 {
2345 // define rebinning strategy for this projection
2346   fNrebin = n;
2347   fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));
2348   fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));
2349 }
2350
2351