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