updates in the PWGPP/TRD
[u/mrichter/AliRoot.git] / PWGPP / 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 <TFile.h>
52 #include <TSystem.h>
53 #include <TStyle.h>
54 #include <TROOT.h>
55 #include <TObjArray.h>
56 #include <TH3.h>
57 #include <TH2.h>
58 #include <TH1.h>
59 #include <THnSparse.h>
60 #include <TF1.h>
61 #include <TCanvas.h>
62 #include <TGaxis.h>
63 #include <TBox.h>
64 #include <TLegend.h>
65 #include <TGraphErrors.h>
66 #include <TGraphAsymmErrors.h>
67 #include <TLinearFitter.h>
68 #include <TMath.h>
69 #include <TMatrixT.h>
70 #include <TVectorT.h>
71 #include <TTreeStream.h>
72 #include <TGeoManager.h>
73 #include <TDatabasePDG.h>
74
75 #include "AliPID.h"
76 #include "AliLog.h"
77 #include "AliESDtrack.h"
78 #include "AliMathBase.h"
79 #include "AliTrackPointArray.h"
80
81 #include "AliTRDresolution.h"
82 #include "AliTRDgeometry.h"
83 #include "AliTRDtransform.h"
84 #include "AliTRDpadPlane.h"
85 #include "AliTRDcluster.h"
86 #include "AliTRDseedV1.h"
87 #include "AliTRDtrackV1.h"
88 #include "AliTRDReconstructor.h"
89 #include "AliTRDrecoParam.h"
90 #include "AliTRDpidUtil.h"
91 #include "AliTRDinfoGen.h"
92 #include "AliAnalysisManager.h"
93 #include "info/AliTRDclusterInfo.h"
94 #include "info/AliTRDeventInfo.h"
95
96 ClassImp(AliTRDresolution)
97 ClassImp(AliTRDresolution::AliTRDrecoProjection)
98
99 Int_t const   AliTRDresolution::fgkNbins[kNdim] = {
100   Int_t(kNbunchCross)/*bc*/,
101   144/*phi*/,
102   45/*eta*/,
103   50/*dy*/,
104   40/*dphi*/,
105   50/*dz*/,
106   5/*chg*species*/,
107   kNpt/*pt*/
108 };  //! no of bins/projection
109 Double_t const AliTRDresolution::fgkMin[kNdim] = {
110   -1.5,
111   -TMath::Pi(),
112   -0.9,
113   -1.5,
114   -5.,
115   -2.5,
116   -2.5,
117   -0.5
118 };    //! low limits for projections
119 Double_t const AliTRDresolution::fgkMax[kNdim] = {
120   1.5,
121   TMath::Pi(),
122   0.9,
123   1.5,
124   5.,
125   2.5,
126   2.5,
127   kNpt-0.5
128 };    //! high limits for projections
129 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
130   "bunch cross",
131   "#phi [rad]",
132   "#eta",
133   "#Deltay [cm]",
134   "#Delta#phi [deg]",
135   "#Deltaz [cm]",
136   "chg*spec*rc",
137   "bin_p_{t}"
138 };  //! title of projection
139
140 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
141     "Det2Cluster"
142     ,"Cluster2Track"
143     ,"Tracklet2Track"
144     ,"Tracklet2TRDin"
145     ,"Cluster2MC"
146     ,"Tracklet2MC"
147     ,"TRDin2MC"
148     ,"TRD2MC"
149 //    ,"Tracklet2TRDout"
150 //    ,"TRDout2MC"
151 };
152
153 //________________________________________________________
154 AliTRDresolution::AliTRDresolution()
155   :AliTRDrecoTask()
156   ,fSteer(0)
157   ,fPtThreshold(.3)
158   ,fBCbinTOF(0)
159   ,fBCbinFill(0)
160   ,fBsign(kFALSE)
161   ,fProj(NULL)
162   ,fDBPDG(NULL)
163   ,fCl(NULL)
164   ,fMCcl(NULL)
165 {
166   //
167   // Default constructor
168   //
169   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
170   SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
171 }
172
173 //________________________________________________________
174 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
175   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
176   ,fSteer(0)
177   ,fPtThreshold(.3)
178   ,fBCbinTOF(0)
179   ,fBCbinFill(0)
180   ,fBsign(kFALSE)
181   ,fProj(NULL)
182   ,fDBPDG(NULL)
183   ,fCl(NULL)
184   ,fMCcl(NULL)
185 {
186   //
187   // Default constructor
188   //
189
190   InitFunctorList();
191   SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
192   if(xchange){
193     SetUseExchangeContainers();
194     DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
195     DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
196   }
197 }
198
199 //________________________________________________________
200 AliTRDresolution::~AliTRDresolution()
201 {
202   //
203   // Destructor
204   //
205   if (AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
206   if(fProj){fProj->Delete(); delete fProj;}
207   if(fCl){fCl->Delete(); delete fCl;}
208   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
209 }
210
211
212 //________________________________________________________
213 void AliTRDresolution::UserCreateOutputObjects()
214 {
215   // spatial resolution
216
217   AliTRDrecoTask::UserCreateOutputObjects();
218   if(UseExchangeContainers()) InitExchangeContainers();
219 }
220
221 //________________________________________________________
222 void AliTRDresolution::InitExchangeContainers()
223 {
224 // Init containers for subsequent tasks (AliTRDclusterResolution)
225
226   fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
227   fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
228   PostData(kClToTrk, fCl);
229   PostData(kClToMC, fMCcl);
230 }
231
232 //________________________________________________________
233 void AliTRDresolution::UserExec(Option_t *opt)
234 {
235   //
236   // Execution part
237   //
238
239   if(fCl) fCl->Delete();
240   if(fMCcl) fMCcl->Delete();
241   AliTRDrecoTask::UserExec(opt);
242 }
243
244 //________________________________________________________
245 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
246 {
247 // Helper function to calculate pulls in the yz plane
248 // using proper tilt rotation
249 // Uses functionality defined by AliTRDseedV1.
250
251   return kTRUE;
252 /*
253   Double_t t2(tilt*tilt);
254   // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
255
256   // rotate along pad
257   Double_t cc[3];
258   cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
259   cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
260   cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
261   // do sqrt
262   Double_t sqr[3]={0., 0., 0.};
263   if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
264   Double_t invsqr[3]={0., 0., 0.};
265   if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
266   Double_t tmp(dyz[0]);
267   dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
268   dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
269   return kTRUE;
270 */
271 }
272
273 //________________________________________________________
274 TH1* AliTRDresolution::DetCluster(const TObjArray *cls)
275 {
276   //
277   // Plot the cluster distributions
278   //
279
280   if(cls) fkClusters = cls;
281   if(!fkClusters){
282     AliDebug(4, "No Clusters defined.");
283     return NULL;
284   }
285   Int_t ncl(0);
286   if(!(ncl = fkClusters->GetEntriesFast())){
287     AliDebug(1, "No RecPoints defined.");
288     return NULL;
289   }
290   Int_t det(-1);
291   AliTRDcluster *cl(NULL);
292   for(Int_t icl(0); icl<ncl; icl++){
293     if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])) continue;
294     det = cl->GetDetector(); break;
295   }
296   if(det<0){
297     AliDebug(1, "No useful clusters defined.");
298     return NULL;
299   }
300   THnSparse *H(NULL);
301   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kDetector)))){
302     AliWarning("No output container defined.");
303     return NULL;
304   }
305   Int_t ly(AliTRDgeometry::GetLayer(det)),
306         stk(AliTRDgeometry::GetStack(det));
307   Double_t val[kNdim],
308            alpha((0.5+AliTRDgeometry::GetSector(det))*AliTRDgeometry::GetAlpha()),
309            cs(TMath::Cos(alpha)),
310            sn(TMath::Sin(alpha));
311   for(Int_t icl(0); icl<ncl; icl++){
312     if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])){
313       AliDebug(2, Form("Missing clusters @ %d", icl));
314       continue;
315     }
316     if(cl->IsMasked()) printf("Mask %03d[%02d_%d_%d] %d[%d]", det, AliTRDgeometry::GetSector(det), stk, ly, cl->GetPadMaskedPosition(), cl->GetPadMaskedStatus());
317     val[kBC]  = ly;
318     val[kPhi] = TMath::ATan2(cl->GetX()*sn + cl->GetY()*cs, cl->GetX()*cs - cl->GetY()*sn);
319     val[kEta] = (5-stk)*16-cl->GetPadRow()-1-(stk<3?4:0);
320     val[kYrez]= fEvent->GetMultiplicity();
321     val[4]    = TMath::Abs(cl->GetQ());
322     val[5]    = fTriggerList?fTriggerSlot:(cl->IsFivePad()?1:cl->GetNPads());
323     H->Fill(val);
324   }
325   return NULL;
326 }
327
328 //________________________________________________________
329 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
330 {
331   //
332   // Plot the cluster distributions
333   //
334
335   if(track) fkTrack = track;
336   if(!fkTrack){
337     AliDebug(4, "No Track defined.");
338     return NULL;
339   }
340   if(fkESD && TMath::Abs(fkESD->GetTOFbc()) > 1){
341     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
342     return NULL;
343   }
344   if(fkESD && fPt<fPtThreshold){
345     AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
346     return NULL;
347   }
348   THnSparse *H(NULL);
349   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
350     AliWarning("No output container defined.");
351     return NULL;
352   }
353
354   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
355   Double_t val[kNdim+2],
356            alpha(0.), cs(-2.), sn(0.); //Float_t exb, vd, t0, s2, dl, dt;
357   TObjArray     *clInfoArr(NULL);
358   AliTRDseedV1  *fTracklet(NULL);
359   AliTRDcluster *c(NULL), *cc(NULL);
360   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
361     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
362     if(!fTracklet->IsOK()) continue;
363     //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
364     val[kBC]  = ily;
365     if(cs<-1.){
366       alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
367       cs    = TMath::Cos(alpha);
368       sn    = TMath::Sin(alpha);
369     }
370     val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
371     Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
372     val[kEta] = -TMath::Log(TMath::Tan(0.5 *  (0.5*TMath::Pi() - TMath::ATan(tgl))));
373     val[kPt]  = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
374     Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
375     Int_t row0(-1);
376     Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
377     fTracklet->ResetClusterIter(kTRUE);
378     while((c = fTracklet->NextCluster())){
379       Float_t xc(c->GetX()),
380               q(TMath::Abs(c->GetQ()));
381       if(row0<0) row0 = c->GetPadRow();
382
383       val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
384       val[kPrez] = fTracklet->GetX0()-xc;
385       val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
386       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
387       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
388       if(ic) val[kZrez] /= (ic*q);
389       val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
390       val[kNdim]   = fEvent?fEvent->GetMultiplicity():0.;
391       val[kNdim+1] = fTriggerList?fTriggerSlot:(c->IsFivePad()?1:c->GetNPads());
392       H->Fill(val);
393 /*      // tilt rotation of covariance for clusters
394       Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
395       cov[0] = (sy2+t2*sz2)*corr;
396       cov[1] = tilt*(sz2 - sy2)*corr;
397       cov[2] = (t2*sy2+sz2)*corr;
398       // sum with track covariance
399       cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
400       Double_t dyz[2]= {dy[1], dz[1]};
401       Pulls(dyz, cov, tilt);*/
402
403       // Get z-position with respect to anode wire
404       Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
405               zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
406       Int_t istk = geo->GetStack(c->GetDetector());
407       AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
408       Float_t rowZ = pp->GetRow0();
409       Float_t d  = rowZ - zt + pp->GetAnodeWireOffset();
410       d -= ((Int_t)(2 * d)) / 2.0;
411       if (d > 0.25) d  = 0.5 - d;
412
413       AliTRDclusterInfo *clInfo(NULL);
414       clInfo = new AliTRDclusterInfo;
415       clInfo->SetCluster(c);
416       //Float_t covF[] = {cov[0], cov[1], cov[2]};
417       clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
418       clInfo->SetResolution(val[kYrez]);
419       clInfo->SetAnisochronity(d);
420       clInfo->SetDriftLength(val[kZrez]);
421       clInfo->SetTilt(fTracklet->GetTilt());
422       if(fCl) fCl->Add(clInfo);
423       //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
424
425       if(DebugLevel()>=2){
426         if(!clInfoArr){
427           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
428           clInfoArr->SetOwner(kFALSE);
429         }
430         clInfoArr->Add(clInfo);
431       }
432     }
433     if(DebugLevel()>=2 && clInfoArr){
434       ULong_t status = fkESD->GetStatus();
435       (*DebugStream()) << "cluster"
436         <<"status="  << status
437         <<"clInfo.=" << clInfoArr
438         << "\n";
439       clInfoArr->Clear();
440     }
441   }
442   if(clInfoArr) delete clInfoArr;
443
444   if(!track) return NULL;
445   // special care for EVE usage
446   TH1 *h(NULL);
447   if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
448   return H->Projection(kYrez);
449 }
450
451
452 //________________________________________________________
453 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
454 {
455 // Plot normalized residuals for tracklets to track.
456 //
457 // We start from the result that if X=N(|m|, |Cov|)
458 // BEGIN_LATEX
459 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
460 // END_LATEX
461 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
462 // reference position.
463   if(track) fkTrack = track;
464   if(!fkTrack){
465     AliDebug(4, "No Track defined.");
466     return NULL;
467   }
468   if(fkESD && TMath::Abs(fkESD->GetTOFbc())>1){
469     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
470     //return NULL;
471   }
472   THnSparse *H(NULL);
473   if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
474     AliWarning("No output container defined.");
475     return NULL;
476   }
477
478   // down scale PID resolution
479   Int_t spc(fSpecies); if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
480   const Int_t ndim(kNdim+8);
481   Double_t val[ndim],
482            alpha(0.), cs(-2.), sn(0.);
483 //  Float_t sz[AliTRDseedV1::kNtb], pos[AliTRDseedV1::kNtb]; Int_t ntbGap[AliTRDseedV1::kNtb];
484   AliTRDseedV1 *fTracklet(NULL);
485   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
486     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
487     if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()) continue;
488     val [kBC] = il;
489     if(cs<-1.){
490       alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
491       cs    = TMath::Cos(alpha);
492       sn    = TMath::Sin(alpha);
493     }
494     val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
495     Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());//fTracklet->GetTgl();
496     val[kEta] = -TMath::Log(TMath::Tan(0.5 *  (0.5*TMath::Pi() - TMath::ATan(tgl))));
497
498     val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:spc;// fSpecies;
499     val[kPt]  = GetPtBin(fPt); //fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fTracklet->GetMomentum());
500     Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)),
501              dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)),
502              dydx(fTracklet->GetYfit(1)),
503              tilt(fTracklet->GetTilt());
504     // correct for tilt rotation
505     val[kYrez] = dyt - dzt*tilt;
506     val[kZrez] = fTracklet->IsRowCross()?(dzt + dyt*tilt):(fTracklet->GetdQdl()*3.e-4-1.5);
507     dydx+= tilt*fTracklet->GetZref(1);
508     val[kPrez] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
509     val[kNdim] = fEvent?fEvent->GetMultiplicity():0;
510     val[kNdim+1] = fTriggerList?fTriggerSlot:(1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb);
511 /*    Int_t n = fTracklet->GetChargeGaps(sz, pos, ntbGap);
512     val[kNdim+2] = 0.; for(Int_t igap(0); igap<n; igap++) val[kNdim+2] += sz[igap];
513     for(Int_t ifill(0); ifill<3; ifill++){ val[kNdim+3+ifill]=0.;val[kNdim+4+ifill]=0.;}
514     for(Int_t igap(0), ifill(0); igap<n&&ifill<2; igap++){
515       if(ntbGap[igap]<2) continue;
516       val[kNdim+3+ifill] = sz[igap];
517       val[kNdim+4+ifill] = pos[igap];
518       ifill++;
519     }*/
520     H->Fill(val);
521
522 //     // compute covariance matrix
523 //     fTracklet->GetCovAt(x, cov);
524 //     fTracklet->GetCovRef(covR);
525 //     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
526 //     Double_t dyz[2]= {dy[1], dz[1]};
527 //     Pulls(dyz, cov, tilt);
528 //     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
529 //     ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
530
531     if(DebugLevel()>=3){
532       Bool_t rc(fTracklet->IsRowCross());
533       UChar_t err(fTracklet->GetErrorMsg());
534       Double_t x(fTracklet->GetX()),
535                pt(fTracklet->GetPt()),
536                yt(fTracklet->GetYref(0)),
537                zt(fTracklet->GetZref(0)),
538                phi(fTracklet->GetYref(1)),
539                tht(fTracklet->GetZref(1));
540       Int_t ncl(fTracklet->GetN()), det(fTracklet->GetDetector());
541       (*DebugStream()) << "tracklet"
542         <<"pt="  << pt
543         <<"x="   << x
544         <<"yt="  << yt
545         <<"zt="  << zt
546         <<"phi=" << phi
547         <<"tht=" << tht
548         <<"det=" << det
549         <<"n="   << ncl
550         <<"dy0=" << dyt
551         <<"dz0=" << dzt
552         <<"dy="  << val[kYrez]
553         <<"dz="  << val[kZrez]
554         <<"dphi="<< val[kPrez]
555         <<"dQ  ="<< val[kNdim]
556         <<"rc="  << rc
557         <<"err=" << err
558         << "\n";
559     }
560   }
561   if(!track) return NULL;
562   // special care for EVE usage
563   TH1 *h(NULL);
564   if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
565   return H->Projection(kYrez);
566 }
567
568
569 //________________________________________________________
570 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
571 {
572 // Store resolution/pulls of Kalman before updating with the TRD information
573 // at the radial position of the first tracklet. The following points are used
574 // for comparison
575 //  - the (y,z,snp) of the first TRD tracklet
576 //  - the (y, z, snp, tgl, pt) of the MC track reference
577 //
578 // Additionally the momentum resolution/pulls are calculated for usage in the
579 // PID calculation.
580  
581   if(track) fkTrack = track;
582   if(!fkTrack){
583     AliDebug(4, "No Track defined.");
584     return NULL;
585   }
586   //fkTrack->Print();
587   TH1 *h(NULL); // EVE projection
588   // check container
589   THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
590   if(!H){
591     AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
592     return NULL;
593   }
594   // check input track status
595   AliExternalTrackParam *tin(NULL);
596   if(!(tin = fkTrack->GetTrackIn())){
597     AliError("Track did not entered TRD fiducial volume.");
598     return NULL;
599   }
600   // check first tracklet
601   AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
602   if(!fTracklet){
603     AliDebug(3, "No Tracklet in ly[0]. Skip track.");
604     if(!track) return NULL;
605     // special care for EVE usage
606     if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
607     return H->Projection(kYrez);
608   }
609   if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()){
610     AliDebug(3, "Tracklet or Chamber not OK. Skip track.");
611     if(!track) return NULL;
612     // special care for EVE usage
613     if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
614     return H->Projection(kYrez);
615   }
616   // check radial position
617   Double_t x = tin->GetX();
618   if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
619     AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
620     if(!track) return NULL;
621     // special care for EVE usage
622     if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
623     return H->Projection(kYrez);
624   }
625
626   // down scale PID resolution
627   Int_t spc(fSpecies); if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
628
629   if(DebugLevel()>=3){
630     Float_t tpc(fkESD->GetTPCdedx());
631     Float_t tof(fkESD->GetTOFbeta());
632     (*DebugStream()) << "trackIn"
633       <<"spc="        << spc
634       <<"tpc="        << tpc
635       <<"tof="        << tof
636       <<"tracklet.="  << fTracklet
637       <<"trackIn.="   << tin
638       << "\n";
639   }
640
641   //Int_t bc(fkESD?fkESD->GetTOFbc()/2:0);
642   const Double_t *parR(tin->GetParameter());
643   Double_t dyt(fTracklet->GetYfit(0)-parR[0]), dzt(fTracklet->GetZfit(0)-parR[1]),
644             phit(fTracklet->GetYfit(1)),
645             tilt(fTracklet->GetTilt()),
646             norm(1./TMath::Sqrt((1.-parR[2])*(1.+parR[2])));
647
648   // correct for tilt rotation
649   Double_t dy  = dyt - dzt*tilt,
650            dz  = dzt + dyt*tilt,
651            dx  = dy/(parR[2]*norm-parR[3]*norm*tilt);
652   phit       += tilt*parR[3];
653   Double_t dphi = TMath::ATan(phit) - TMath::ASin(parR[2]);
654
655   Double_t val[kNdim+2];
656   val[kBC]          = fTriggerSlot;//bc==0?0:(bc<0?-1.:1.);
657   Double_t alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha(),
658            cs    = TMath::Cos(alpha),
659            sn    = TMath::Sin(alpha);
660   val[kPhi]         = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
661   Float_t tgl       = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
662   val[kEta]         = -TMath::Log(TMath::Tan(0.5 *  (0.5*TMath::Pi() - TMath::ATan(tgl))));
663   val[kYrez]        = dy;
664   val[kZrez]        = fTracklet->IsRowCross()?dz:(fTracklet->GetdQdl()*5.e-4 - 2.5);
665   val[kPrez]        = dphi*TMath::RadToDeg();
666   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:spc;
667   val[kPt]          = GetPtBin(fPt);
668   val[kNdim]        = GetPtBin(fTracklet->GetMomentum());
669   val[kNdim+1]      = dx;
670   //val[kNdim+2]      = fEvent?fEvent->GetBunchFill():0;
671   H->Fill(val);
672
673   if(!track) return NULL;
674   // special care for EVE usage
675   if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
676   return H->Projection(kYrez);
677 }
678
679 /*
680 //________________________________________________________
681 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
682 {
683 // Store resolution/pulls of Kalman after last update with the TRD information
684 // at the radial position of the first tracklet. The following points are used
685 // for comparison
686 //  - the (y,z,snp) of the first TRD tracklet
687 //  - the (y, z, snp, tgl, pt) of the MC track reference
688 //
689 // Additionally the momentum resolution/pulls are calculated for usage in the
690 // PID calculation.
691
692   if(track) fkTrack = track;
693   return NULL;
694 }
695 */
696 //________________________________________________________
697 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
698 {
699   //
700   // Plot MC distributions
701   //
702
703   if(!HasMCdata()) return NULL;
704   if(track) fkTrack = track;
705   if(!fkTrack){
706     AliDebug(4, "No Track defined.");
707     return NULL;
708   }
709   Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
710
711   THnSparse *H(NULL);
712   if(!fContainer){
713     AliWarning("No output container defined.");
714     return NULL;
715   }
716   // retriev track characteristics
717   Int_t pdg = fkMC->GetPDG(),
718         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
719         sign(0),
720 //        sgm[3],
721         label(fkMC->GetLabel());
722 //        fSegmentLevel(0);
723   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
724   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
725   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
726   if(sIdx==3) sIdx=2; if(sIdx>3) sIdx=3; // downscale PID
727
728   TH1 *h(NULL);
729   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
730   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
731   UChar_t s;
732   Double_t x, y, z, pt, dydx, dzdx, dzdl;
733   Float_t pt0, p0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
734   Double_t covR[7]/*, cov[3]*/;
735
736 /*  if(DebugLevel()>=3){
737     // get first detector
738     Int_t det = -1;
739     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
740       if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
741       det = fTracklet->GetDetector();
742       break;
743     }
744     if(det>=0){
745       TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
746       Double_t m(-1.);
747       m = fkTrack->GetMass();
748       if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
749         (*DebugStream()) << "MCkalman"
750           << "pdg=" << pdg
751           << "det=" << det
752           << "x="   << &X
753           << "y="   << &Y
754           << "z="   << &Z
755           << "dx="  << &dX
756           << "dy="  << &dY
757           << "dz="  << &dZ
758           << "pt="  << &vPt
759           << "dpt=" << &dPt
760           << "bgt=" << &budget
761           << "cov=" << &cCOV
762           << "\n";
763       }
764     }
765   }*/
766   AliTRDcluster *c(NULL);
767   Double_t val[kNdim+2];
768   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
769     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
770        !fTracklet->IsOK())*/ continue;
771
772     x= x0 = fTracklet->GetX();
773     Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
774     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, p0, eta, phi, s)) continue;
775
776     // MC track position at reference radial position
777     dx  = x0 - x;
778     Float_t ymc = y0 - dx*dydx0;
779     Float_t zmc = z0 - dx*dzdx0;
780     //phi -= TMath::Pi();
781
782     val[kBC]  = ily;
783     val[kPhi] = phi;
784     val[kEta] = eta;
785     val[kSpeciesChgRC]= rc?0.:sign*sIdx;
786     val[kPt]  = GetPtBin(pt0);
787     val[kNdim]= GetPtBin(p0);
788     Double_t tilt(fTracklet->GetTilt());
789 //             ,t2(tilt*tilt)
790 //             ,corr(1./(1. + t2))
791 //             ,cost(TMath::Sqrt(corr));
792
793     AliExternalTrackParam *tin(fkTrack->GetTrackIn());
794     if(ily==0 && tin){ // trackIn residuals
795       // check radial position
796       if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
797       else{
798         val[kBC]          = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
799         val[kYrez]        = tin->GetY()-ymc;
800         val[kZrez]        = rc?(tin->GetZ()-zmc):(fTracklet->GetdQdl()*5e-4 - 2.5);
801         val[kPrez]        = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
802         val[kNdim+1]      = 0.;//dx;
803         if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
804       }
805     }
806     //if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
807
808     // track residuals
809     dydx = fTracklet->GetYref(1);
810     dzdx = fTracklet->GetZref(1);
811     dzdl = fTracklet->GetTgl();
812     y  = fTracklet->GetYref(0);
813     dy = y - ymc;
814     z  = fTracklet->GetZref(0);
815     dz = z - zmc;
816     pt = TMath::Abs(fTracklet->GetPt());
817     fTracklet->GetCovRef(covR);
818
819     val[kYrez] = dy;
820     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
821     val[kZrez] = dz;
822     val[kPt]   = GetPtBin(pt0);
823     val[kNdim] = 1.e2*(pt/pt0-1.);
824     if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
825 /*      // theta resolution/ tgl pulls
826       Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
827                 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
828       ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
829       ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
830       // pt resolution  \\ 1/pt pulls \\ p resolution for PID
831       Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
832               p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
833       ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
834       ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
835       ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
836
837     // Fill Debug stream for MC track
838     if(DebugLevel()>=4){
839       Int_t det(fTracklet->GetDetector());
840       (*DebugStream()) << "MC"
841         << "det="     << det
842         << "pdg="     << pdg
843         << "sgn="     << sign
844         << "pt="      << pt0
845         << "x="       << x0
846         << "y="       << y0
847         << "z="       << z0
848         << "dydx="    << dydx0
849         << "dzdx="    << dzdx0
850         << "\n";
851
852       // Fill Debug stream for Kalman track
853       (*DebugStream()) << "MCtrack"
854         << "pt="      << pt
855         << "x="       << x
856         << "y="       << y
857         << "z="       << z
858         << "dydx="    << dydx
859         << "dzdx="    << dzdx
860         << "s2y="     << covR[0]
861         << "s2z="     << covR[2]
862         << "\n";
863     }
864
865     // tracklet residuals
866     dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
867     dzdx = fTracklet->GetZfit(1);
868     y  = fTracklet->GetYfit(0);
869     dy = y - ymc;
870     z  = fTracklet->GetZfit(0);
871     dz = z - zmc;
872     val[kYrez] = dy - dz*tilt;
873     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
874     val[kZrez] = rc?(dz + dy*tilt):(fTracklet->GetdQdl()*3.e-4 - 1.5);
875     val[kNdim] = fEvent?fEvent->GetMultiplicity():0;
876     val[kNdim+1] = 1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb;
877 //      val[kNdim] = pt/pt0-1.;
878     if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
879
880
881     // Fill Debug stream for tracklet
882     if(DebugLevel()>=4){
883       Float_t s2y = fTracklet->GetS2Y();
884       Float_t s2z = fTracklet->GetS2Z();
885       (*DebugStream()) << "MCtracklet"
886         << "rc="    << rc
887         << "x="     << x
888         << "y="     << y
889         << "z="     << z
890         << "dydx="  << dydx
891         << "s2y="   << s2y
892         << "s2z="   << s2z
893         << "\n";
894     }
895
896     AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
897     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
898     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
899
900     H = (THnSparse*)fContainer->At(kMCcluster);
901     val[kPt]  = TMath::ATan(dydx0)*TMath::RadToDeg();
902     //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
903     Int_t row0(-1);
904     Float_t padCorr(tilt*fTracklet->GetPadLength()),
905             corr(1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0));
906     fTracklet->ResetClusterIter(kTRUE);
907     while((c = fTracklet->NextCluster())){
908       if(row0<0) row0 = c->GetPadRow();
909       x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
910       y = c->GetY()  + padCorr*(c->GetPadRow() - row0);
911       z = c->GetZ();
912       dx = x0 - x;
913       ymc= y0 - dx*dydx0;
914       zmc= z0 - dx*dzdx0;
915       dy = y - ymc;
916       dz = z - zmc;
917       val[kYrez] = dy - dz*tilt;
918       val[kPrez] = dx;
919       val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t  q(TMath::Abs(c->GetQ()));
920       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
921       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
922       if(ic) val[kZrez] /= (ic*q);
923       val[kSpeciesChgRC]= rc?0.:(TMath::Max(q*corr, Float_t(3.)));
924       val[kNdim+1] = c->IsFivePad()?1:c->GetNPads();
925       if(H) H->Fill(val);
926
927
928       // Fill calibration container
929       Float_t d = zr0 - zmc;
930       d -= ((Int_t)(2 * d)) / 2.0;
931       if (d > 0.25) d  = 0.5 - d;
932       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
933       clInfo->SetCluster(c);
934       clInfo->SetMC(pdg, label);
935       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
936       clInfo->SetResolution(dy);
937       clInfo->SetAnisochronity(d);
938       clInfo->SetDriftLength(dx);
939       clInfo->SetTilt(tilt);
940       if(fMCcl) fMCcl->Add(clInfo);
941       else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
942       if(DebugLevel()>=5){
943         if(!clInfoArr){
944           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
945           clInfoArr->SetOwner(kFALSE);
946         }
947         clInfoArr->Add(clInfo);
948       }
949     }
950     // Fill Debug Tree
951     if(DebugLevel()>=5 && clInfoArr){
952       (*DebugStream()) << "MCcluster"
953         <<"clInfo.=" << clInfoArr
954         << "\n";
955       clInfoArr->Clear();
956     }
957   }
958   if(clInfoArr) delete clInfoArr;
959   if(!track) return NULL;
960   // special care for EVE usage
961   if(H && (h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
962   return H?H->Projection(kYrez):NULL;
963 }
964
965 //________________________________________________________
966 void AliTRDresolution::GetRangeZ(TH2 *h2, Float_t &min, Float_t &max)
967 {
968 // Get range on Z axis such to avoid outliers
969
970   Double_t cnt[20000], c, m, s;
971   Int_t nx(h2->GetXaxis()->GetNbins()), ny(h2->GetYaxis()->GetNbins()), nc(0);
972   for(Int_t ix(1); ix<=nx; ix++){
973     for(Int_t iy(1); iy<=ny; iy++){
974       if((c = h2->GetBinContent(ix, iy))<10) continue;
975       cnt[nc++] = c;
976       if(nc==20000) break;
977     }
978     if(nc==20000) break;
979   }
980   AliMathBase::EvaluateUni(nc, cnt, m, s, 0);
981   min = m-s; max = m+2.*s;
982 }
983
984 //________________________________________________________
985 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
986 {
987   //
988   // Get the reference figures
989   //
990
991   if(!gPad){
992     AliWarning("Please provide a canvas to draw results.");
993     return kFALSE;
994   }
995 /*  Int_t selection[100], n(0), selStart(0); //
996   Int_t ly0(0), dly(5);
997   TList *l(NULL); TVirtualPad *pad(NULL); */
998   switch(ifig){
999   case 0:
1000     break;
1001   }
1002   AliWarning(Form("Reference plot [%d] missing result", ifig));
1003   return kFALSE;
1004 }
1005
1006 //________________________________________________________
1007 void AliTRDresolution::MakeSummary()
1008 {
1009 // Build summary plots
1010
1011   if(!fProj){
1012     AliError("Missing results");
1013     return;
1014   }
1015   TVirtualPad *p(NULL); TCanvas *cOut(NULL);
1016   TObjArray *arr(NULL); TH2 *h2(NULL);
1017
1018   // cluster resolution
1019   // define palette
1020   gStyle->SetPalette(1);
1021   const Int_t nClViews(9);
1022   const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh"};
1023   const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1};
1024   const Float_t vClSignMin[2] = {2.6e2, 4.4e2},
1025                 vClSignMax[2] = {4.4e2, 6.2e2},
1026                 vClMin[nClViews] = {3.2e2, vClSignMin[Int_t(fBsign)], vClSignMin[Int_t(!fBsign)], 0., 0., 0., 0., 0., 3.2e2},
1027                 vClMax[nClViews] = {5.e2, vClSignMax[Int_t(fBsign)], vClSignMax[Int_t(!fBsign)], 0., 0., 0., 0., 0., 5.e2};
1028   const Int_t nTrkltViews(20);
1029   const Char_t *vTrkltName[nTrkltViews] = {
1030     "TrkltY", "TrkltYn", "TrkltYp", "TrkltY", "TrkltYn", "TrkltYp",
1031     "TrkltRCZ",
1032     "TrkltPh", "TrkltPhn", "TrkltPhp",
1033     "TrkltQ", "TrkltQn", "TrkltQp",
1034     "TrkltQS", "TrkltQSn", "TrkltQSp",
1035     "TrkltTBn", "TrkltTBp", "TrkltTBn", "TrkltTBp"
1036 //    "TrkltYnl0", "TrkltYpl0", "TrkltPhnl0", "TrkltPhpl0", "TrkltQnl0", "TrkltQpl0",  // electrons low pt
1037 /*    "TrkltYnl1", "TrkltYpl1", "TrkltPhnl1", "TrkltPhpl1", "TrkltQnl1", "TrkltQpl1",  // muons low pt
1038     "TrkltYnl2", "TrkltYpl2", "TrkltPhnl2", "TrkltPhpl2", "TrkltQnl2", "TrkltQpl2"  // pions low pt*/
1039   };
1040   const UChar_t vTrkltOpt[nTrkltViews] = {
1041     0, 0, 0, 1, 1, 1,
1042     0,
1043     0, 0, 0,
1044     0, 0, 0,
1045     0, 0, 0,
1046     0, 0, 2, 2
1047   };
1048   const Int_t nTrkInViews(5);
1049   const Char_t *vTrkInName[nTrkInViews][6] = {
1050     {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1051     {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1052     {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"},
1053     {"TrkInXnl", "TrkInXpl", "TrkInXl", "TrkInYnh", "TrkInYph", "TrkInYh"},
1054     {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"},
1055     //{"TrkInRCX", "TrkInRCY", "TrkInRCPh", "TrkInRCZl", "TrkInRCZi", "TrkInRCZh"}
1056   };
1057   const UChar_t vTrkInOpt[nTrkInViews] = {0, 1, 0, 0, 0};
1058   const Float_t min[6] = {0.15, 0.15, 0.15, 0.15, 0.5, 0.5};
1059   const Float_t max[6] = {0.6, 0.6, 0.6, 0.6, 2.3, 2.3};
1060   const Char_t *ttt[6] = {"#sigma(#Deltay) [cm]", "#sigma(#Deltay) [cm]", "#sigma(#Deltay) [cm]", "#sigma(#Deltaz) [cm]", "#sigma(#Delta#phi) [deg]", "#sigma(#Delta#phi) [deg]"};
1061
1062   const Int_t nTrkViews(27);
1063   const Char_t *vTrkName[nTrkViews] = {
1064     "TrkY", "TrkYn", "TrkYp",
1065     "TrkPh", "TrkPhn", "TrkPhp",
1066     "TrkDPt", "TrkDPtn", "TrkDPtp",
1067     "TrkYnl0", "TrkYpl0", "TrkPhnl0", "TrkPhpl0", "TrkDPtnl0", "TrkDPtpl0",  // electrons low pt
1068     "TrkYnl1", "TrkYpl1", "TrkPhnl1", "TrkPhpl1", "TrkDPtnl1", "TrkDPtpl1",  // muons low pt
1069     "TrkYnl2", "TrkYpl2", "TrkPhnl2", "TrkPhpl2", "TrkDPtnl2", "TrkDPtpl2"  // pions low pt
1070   };
1071   const Char_t *typName[] = {"", "MC"};
1072   const Int_t nx(2048), ny(1536);
1073
1074   if((arr = (TObjArray*)fProj->FindObject("hDet2Cluster"))){
1075     cOut = new TCanvas(Form("%s_DetOccupancy", GetName()), "Detector performance", 2*nx, 2*ny);
1076     cOut->Divide(AliTRDgeometry::kNlayer,AliTRDeventInfo::kCentralityClasses, 1.e-5, 1.e-5);
1077     Int_t n=0;
1078     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1079       for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1080         p=cOut->cd(icen*AliTRDgeometry::kNlayer+ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1081         if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d%dEn", ily, icen)))) continue;
1082         if(SetNormZ(h2, 1, -1, 1, -1, 10.) < 1.e3) continue;  // cut all bins with < 10 entries
1083         SetRangeZ(h2, -90., 90, -200.);
1084         //h2->GetXaxis()->SetNdivisions(1010);
1085         h2->GetZaxis()->SetTitle("Rel. Det. Occup. [%]");
1086         h2->GetZaxis()->CenterTitle();
1087         h2->SetContour(9); h2->Draw("colz"); n++;
1088         MakeDetectorPlot(ily, "pad");
1089       }
1090     }
1091     if(n>=AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1092
1093     cOut = new TCanvas(Form("%s_DetCharge", GetName()), "Detector performance", nx, ny);
1094     cOut->Divide(3,2, 1.e-5, 1.e-5);
1095     for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1096       p=cOut->cd(ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1097       if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d_2D", ily)))) continue;
1098       SetNormZ(h2, 1, -1, 1, -1, 10.); // cut all <q> < 10
1099       SetRangeZ(h2, -30., 30., -200.);
1100       //h2->GetXaxis()->SetNdivisions(1010);
1101       h2->GetZaxis()->SetTitle("Rel. Mean(q) [%]");
1102       h2->GetZaxis()->CenterTitle();
1103       h2->Draw("colz");
1104       MakeDetectorPlot(ily, "pad");
1105     }
1106     cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1107   }
1108   for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
1109     if((arr = (TObjArray*)fProj->FindObject(ityp?"hCluster2MC":"hCluster2Track"))){
1110       for(Int_t iview(0); iview<nClViews; iview++){
1111         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vClName[iview], vClOpt[iview]), "Cluster Resolution", nx, ny);
1112         cOut->Divide(3,2, 1.e-5, 1.e-5);
1113         Int_t nplot(0);
1114         for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1115           p=cOut->cd(ily+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1116           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], ily)))) continue;
1117           nplot++;
1118           if(vClOpt[iview]==0) h2->Draw("colz");
1119           else if(vClOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [#mum]", vClMin[iview], vClMax[iview], 1.e4);
1120           if(iview<5) MakeDetectorPlot(ily);
1121         }
1122         if(nplot==AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1123         else delete cOut;
1124       }
1125     }
1126     // tracklet systematic
1127     if((arr = (TObjArray*)fProj->FindObject(ityp?"hTracklet2MC":"hTracklet2Track"))){
1128       for(Int_t iview(0); iview<nTrkltViews; iview++){
1129         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkltName[iview], vTrkltOpt[iview]), "Tracklet Resolution", nx, ny);
1130         cOut->Divide(3,2, 1.e-5, 1.e-5);
1131         Int_t nplot(0);
1132         for(Int_t iplot(0); iplot<6; iplot++){
1133           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1134           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
1135           nplot++;
1136           if(vTrkltOpt[iview]==0) h2->Draw("colz");
1137           else if (vTrkltOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [cm]", .15, .4);
1138           else if (vTrkltOpt[iview]==2) DrawSigma(h2, "#sigma(occupancy) [%]", 10.5, 15.);
1139           MakeDetectorPlot(iplot);
1140         }
1141         if(nplot==6){
1142           cOut->Modified();cOut->Update();
1143           cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1144         }
1145         delete cOut;
1146       }
1147     }
1148     // trackIn systematic
1149     if((arr = (TObjArray*)fProj->FindObject(ityp?"hTRDin2MC":"hTracklet2TRDin"))){
1150       for(Int_t iview(0); iview<nTrkInViews; iview++){
1151         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkInName[iview][0], vTrkInOpt[iview]), "Track IN Resolution", nx, ny);
1152         cOut->Divide(3,2, 1.e-5, 1.e-5);
1153         Int_t nplot(0);
1154         for(Int_t iplot(0); iplot<6; iplot++){
1155           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1156           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])))){
1157             AliInfo(Form("Missing H%s%s_2D", typName[ityp], vTrkInName[iview][iplot]));
1158             continue;
1159           }
1160           nplot++;
1161           if(vTrkInOpt[iview]==0) h2->Draw("colz");
1162           else DrawSigma(h2, ttt[iplot], min[iplot], max[iplot]);
1163           MakeDetectorPlot(0);
1164         }
1165         if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1166         else delete cOut;
1167       }
1168       // species
1169       const Float_t zmin[] = {1., 61., 15.},
1170                     zmax[] = {10., 79., 33.};
1171       cOut = new TCanvas(Form("%s_%sTrkInSpc", GetName(), typName[ityp]), "Track IN PID", Int_t(1.5*ny), Int_t(1.5*ny));
1172       cOut->Divide(3,3, 1.e-5, 1.e-5);
1173       Int_t nplot(0); const Char_t *chName[] = {"p", "n", ""};
1174       for(Int_t ich(0), ipad(1); ich<3; ich++){
1175         TH2 *h2s(NULL);
1176         if(!(h2s = (TH2*)arr->FindObject(Form("H%sTrkInY%sEn", typName[ityp], chName[ich])))) {
1177           AliInfo(Form("Missing H%sTrkIn%sEn", typName[ityp], chName[ich]));
1178           continue;
1179         }
1180         Int_t irebin(0), dxBin(1), dyBin(1);
1181         const Int_t nrebin(5); Int_t rebinX[nrebin] = {1, 2, 1, 2, 1}, rebinY[nrebin] = {2, 1, 2, 1, 2};
1182         if(h2s->GetNbinsY()==180){ // old binning
1183           rebinX[0] = 1; rebinY[0] = 2;
1184           rebinX[1] = 2; rebinY[1] = 1;
1185           rebinX[2] = 2; rebinY[2] = 1;
1186           rebinX[3] = 1; rebinY[3] = 5;
1187           rebinX[4] = 1; rebinY[4] = 1; // dummy
1188         }
1189         while(irebin<nrebin && (AliTRDrecoTask::GetMeanStat(h2s, .5, ">")<200)){
1190           h2s->Rebin2D(rebinX[irebin], rebinY[irebin]);
1191           dxBin*=rebinX[irebin];dyBin*=rebinY[irebin];irebin++;
1192         }
1193         AliDebug(2, Form("Rebin level[%d] @ chg[%d]. Binning[%2dx%2d]", irebin, ich, h2s->GetNbinsX(), h2s->GetNbinsY()));
1194         for(Int_t ispec(0); ispec<kNspc; ispec++){
1195           p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1196           if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkInY%s%dEn", typName[ityp], chName[ich], ispec)))) {
1197             AliInfo(Form("Missing H%sTrkIn%s%dEn", typName[ityp], chName[ich], ispec));
1198             continue;
1199           }
1200           nplot++;
1201           h2->Rebin2D(dxBin,dyBin);
1202           h2->Divide(h2, h2s, 1.e2);
1203           h2->SetContour(9);
1204           h2->GetZaxis()->SetRangeUser(zmin[ispec], zmax[ispec]);
1205           h2->GetZaxis()->SetTitle("Rel. Abundancy [%]");h2->GetZaxis()->CenterTitle();
1206           TString tit(h2->GetTitle()); TObjArray *atit = tit.Tokenize("::");
1207           h2->SetTitle(Form("%s :: Relative Abundancy", ((*atit)[0])->GetName()));
1208           atit->Delete(); delete atit;
1209           h2->Draw("colz");
1210           MakeDetectorPlot(0);
1211         }
1212       }
1213       if(nplot==9) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1214       else delete cOut;
1215
1216       const char *chQ[] = {"Q", "QS"};
1217       for(Int_t iq(0); iq<2; iq++){
1218         cOut = new TCanvas(Form("%s_%sTrkIn%s", GetName(), typName[ityp], chQ[iq]), "Track IN PID", Int_t(1.5*ny), Int_t(1.5*ny));
1219         cOut->Divide(3,3, 1.e-5, 1.e-5);
1220         nplot=0;
1221         for(Int_t ich(0), ipad(1); ich<3; ich++){
1222           for(Int_t ispec(0); ispec<kNspc; ispec++){
1223             p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1224             if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkIn%s%s%d_2D", typName[ityp], chQ[iq], chName[ich], ispec)))) {
1225               AliInfo(Form("Missing H%sTrkIn%s%s%d_2D", typName[ityp], chQ[iq], chName[ich], ispec));
1226               continue;
1227             }
1228             nplot++;
1229             h2->Draw("colz");
1230             MakeDetectorPlot(0);
1231           }
1232         }
1233         if(nplot==9) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1234         else delete cOut;
1235       }
1236     }
1237   }
1238   // track MC systematic
1239   if((arr = (TObjArray*)fProj->FindObject("hTRD2MC"))) {
1240     for(Int_t iview(0); iview<nTrkViews; iview++){
1241       cOut = new TCanvas(Form("%s_MC%s", GetName(), vTrkName[iview]), "Track Resolution", nx, ny);
1242       cOut->Divide(3,2, 1.e-5, 1.e-5);
1243       Int_t nplot(0);
1244       for(Int_t iplot(0); iplot<6; iplot++){
1245         p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1246         if(!(h2 = (TH2*)arr->FindObject(Form("HMC%s%d_2D", vTrkName[iview], iplot)))) continue;
1247         h2->Draw("colz"); nplot++;
1248       }
1249       if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1250       else delete cOut;
1251     }
1252   }
1253
1254
1255   gStyle->SetPalette(1);
1256 }
1257
1258 //________________________________________________________
1259 void AliTRDresolution::DrawSigma(TH2 *h2, const Char_t *title, Float_t m, Float_t M, Float_t scale)
1260 {
1261   // Draw error bars scaled with "scale" instead of content values
1262   //use range [m,M] if limits are specified
1263
1264   if(!h2) return;
1265   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1266   TH2F *h2e = new TH2F(Form("%s_E", h2->GetName()),
1267                 Form("%s;%s;%s;%s", h2->GetTitle(), ax->GetTitle(), ay->GetTitle(), title),
1268                 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
1269                 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1270   h2e->SetContour(9);
1271   TAxis *az(h2e->GetZaxis());
1272   if(M>m) az->SetRangeUser(m, M);
1273   az->CenterTitle();
1274   az->SetTitleOffset(1.5);
1275   for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1276     for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1277       if(h2->GetBinContent(ix, iy)<-100.) continue;
1278       Float_t v(scale*h2->GetBinError(ix, iy));
1279       if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1280       h2e->SetBinContent(ix, iy, v);
1281     }
1282   }
1283   h2e->Draw("colz");
1284 }
1285
1286 //________________________________________________________
1287 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1288 {
1289 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1290 // The "range" vector should be initialized with 2 elements
1291 // Option "mod" can be any of
1292 //   - 0 : gaussian like distribution
1293 //   - 1 : tailed distribution
1294
1295   Int_t nx(h2->GetNbinsX())
1296       , ny(h2->GetNbinsY())
1297       , n(nx*ny);
1298   Double_t *data=new Double_t[n];
1299   for(Int_t ix(1), in(0); ix<=nx; ix++){
1300     for(Int_t iy(1); iy<=ny; iy++)
1301       data[in++] = h2->GetBinContent(ix, iy);
1302   }
1303   Double_t mean, sigm;
1304   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1305
1306   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1307   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1308   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1309   TH1S h1("h1SF0", "", 100, range[0], range[1]);
1310   h1.FillN(n,data,0);
1311   delete [] data;
1312
1313   switch(mod){
1314   case 0:// gaussian distribution
1315   {
1316     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1317     h1.Fit(&fg, "QN");
1318     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1319     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1320     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
1321     break;
1322   }
1323   case 1:// tailed distribution
1324   {
1325     Int_t bmax(h1.GetMaximumBin());
1326     Int_t jBinMin(1), jBinMax(100);
1327     for(Int_t ibin(bmax); ibin--;){
1328       if(h1.GetBinContent(ibin)<1.){
1329         jBinMin=ibin; break;
1330       }
1331     }
1332     for(Int_t ibin(bmax); ibin++;){
1333       if(h1.GetBinContent(ibin)<1.){
1334         jBinMax=ibin; break;
1335       }
1336     }
1337     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1338     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
1339     break;
1340   }
1341   }
1342
1343   return;
1344 }
1345
1346 //________________________________________________________
1347 Bool_t AliTRDresolution::MakeProjectionDetector()
1348 {
1349 // Analyse cluster
1350   const Int_t kNcontours(9);
1351   const Int_t kNstat(100);
1352   if(fProj && fProj->At(kDetector)) return kTRUE;
1353   if(!fContainer){
1354     AliError("Missing data container.");
1355     return kFALSE;
1356   }
1357   THnSparse *H(NULL);
1358   if(!(H = (THnSparse*)fContainer->FindObject("hDet2Cluster"))){
1359     AliInfo(Form("Missing/Wrong data @ hDet2Cluster."));
1360     return kTRUE;
1361   }
1362   Int_t ndim(H->GetNdimensions());
1363   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1364   TAxis *aa[kNdim], *an(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1365   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1366   if(ndim < 5) aa[4] = new TAxis(1, -0.5, 0.5);
1367   Int_t nPad(1);
1368   if(ndim > 5){
1369     an = H->GetAxis(5);
1370     nPad = kNpads+1;
1371   }
1372   // build list of projections
1373   const Int_t nsel(8*AliTRDgeometry::kNlayer*AliTRDeventInfo::kCentralityClasses);
1374   // define rebinning strategy
1375   const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 1, 2, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1376   if(aa[1]->GetNbins()==180){ // old binning
1377     rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1378     rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1379     rebinEtaPhiX[2] = 2; rebinEtaPhiY[2] = 1;
1380     rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1381     rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1382   }
1383   //const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1384   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1385   AliTRDrecoProjection hp[kDetNproj];  TObjArray php(kDetNproj);
1386   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1387   for(Int_t ipad(0); ipad<nPad; ipad++){
1388     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1389       for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1390         isel++; // new selection
1391         hp[ih].Build(Form("HDet%d%d%d", ily, icen, ipad),
1392                     Form("Detectors :: Ly[%d] Cen[%s] Pads[%s]", ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1393                     kEta, kPhi, 4, aa);
1394         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1395         hp[ih].SetShowRange(10., 55.);
1396         php.AddLast(&hp[ih++]); np[isel]++;
1397       }
1398     }
1399   }
1400   AliInfo(Form("Build %3d 3D projections.", ih));
1401
1402   Int_t ly(0), cen(0), npad(0);
1403   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1404     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1405     ly = coord[kBC]-1;    // layer selection
1406     cen = coord[kYrez]-1; // centrality selection
1407     npad = 0;             // no. of pads selection
1408     if(an) npad = TMath::Min(coord[5]-1, Int_t(kNpads));
1409     isel = npad*AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer+cen*AliTRDgeometry::kNlayer+ly;
1410     ((AliTRDrecoProjection*)php.At(isel))->Increment(coord, v);
1411     //Int_t ioff=isel;for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1412   }
1413   TObjArray *arr(NULL);
1414   fProj->AddAt(arr = new TObjArray(kDetNproj), kDetector);
1415   arr->SetName("hDet2Cluster");
1416
1417   TH2 *h2(NULL);  Int_t jh(0);
1418   for(; ih--; ){
1419     if(!hp[ih].H()) continue;
1420     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1421     arr->AddAt(h2, jh++);
1422     if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].H()->GetName())))) continue;
1423     arr->AddAt(h2, jh++);
1424   }
1425   AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1426   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1427     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1428       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, 0)))){
1429         for(Int_t ipad(1); ipad<nPad; ipad++){
1430           if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, ipad)))){
1431             (*pr0)+=(*pr1);
1432           }
1433         }
1434         pr0->H()->SetNameTitle(Form("HDet%d%d", ily, icen), Form("Detectors :: Ly[%d] Cen[%s]", ily, cenName[icen]));
1435         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0, kFALSE))) arr->AddAt(h2, jh++);
1436         if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
1437         if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))) (*pr1)+=(*pr0);
1438       }
1439     }
1440     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))){
1441       pr0->H()->SetNameTitle(Form("HDet%d", ily), Form("Detectors :: Ly[%d]", ily));
1442       if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1443     }
1444   }
1445   AliInfo(Form("Done %3d 2D projections.", jh));
1446   return kTRUE;
1447 }
1448
1449 //________________________________________________________
1450 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1451 {
1452 // Analyse cluster
1453   const Int_t kNcontours(9);
1454   const Int_t kNstat(100);
1455   Int_t cidx=mc?kMCcluster:kCluster;
1456   if(fProj && fProj->At(cidx)) return kTRUE;
1457   if(!fContainer){
1458     AliError("Missing data container.");
1459     return kFALSE;
1460   }
1461   const Char_t *projName[] = {"hCluster2Track", "hCluster2MC"};
1462   THnSparse *H(NULL);
1463   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1464     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1465     return kTRUE;
1466   }
1467   Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimCl));
1468   Int_t coord[10]; memset(coord, 0, sizeof(Int_t) * 10); Double_t v = 0.;
1469   TAxis *aa[kNdim], *as(NULL), *apt(NULL), *acen(NULL), *anp(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1470   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1471   if(ndim > Int_t(kPt)) apt = H->GetAxis(kPt);
1472   if(ndim > Int_t(kSpeciesChgRC)) as  = H->GetAxis(kSpeciesChgRC);
1473   if(ndim > Int_t(kNdim)) acen  = H->GetAxis(kNdim);
1474   if(ndim > Int_t(kNdim)+1) anp  = H->GetAxis(kNdim+1);
1475   // calculate size depending on debug level
1476   const Int_t nCh(apt?2:1);
1477   const Int_t nCen(acen?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1478   const Int_t nNpad(anp?(Int_t(kNpads)+1):1);
1479   
1480   // build list of projections
1481   const Int_t nsel(AliTRDgeometry::kNlayer*kNcharge*(kNpads+1)*AliTRDeventInfo::kCentralityClasses);
1482   // define rebinning strategy
1483   const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1484   if(aa[1]->GetNbins()==180){ // old binning
1485     rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1486     rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1487     rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
1488     rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1489     rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1490   }
1491   AliTRDrecoProjection hp[kClNproj];  TObjArray php(kClNproj);
1492   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1493   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1494   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1495   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1496     for(Int_t ich(0); ich<nCh; ich++){
1497       for(Int_t icen(0); icen<nCen; icen++){
1498         for(Int_t ipad(0); ipad<nNpad; ipad++){
1499           isel++; // new selection
1500           hp[ih].Build(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1501                       Form("Clusters[%c] :: #Deltay Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1502                       kEta, kPhi, kYrez, aa);
1503           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1504           php.AddLast(&hp[ih++]); np[isel]++;
1505           if(!debug) break;
1506           hp[ih].Build(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1507                       Form("Clusters[%c] :: Q Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1508                       kEta, kPhi, kSpeciesChgRC, aa);
1509           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1510           hp[ih].SetShowRange(24., 33.);
1511           php.AddLast(&hp[ih++]); np[isel]++;
1512           hp[ih].Build(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1513                       Form("Clusters[%c] :: #Deltay(x,TC) Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1514                       kPrez, kZrez, kYrez, aa);
1515           php.AddLast(&hp[ih++]); np[isel]++;
1516           hp[ih].Build(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1517                         Form("Clusters[%c] :: #Deltay(x,#Phi) Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1518                         kPrez, kPt, kYrez, aa);
1519           php.AddLast(&hp[ih++]); np[isel]++;
1520         }
1521       }
1522     }
1523   }
1524   AliInfo(Form("Build %3d 3D %s projections.", ih, mc?"MC":""));
1525
1526   AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1527   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1), ioff(0), cen(0), npad(0);
1528   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1529     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1530     ly = coord[kBC]-1;
1531     // RC selection
1532     if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1533
1534     // charge selection
1535     ch = 0; // [-] track
1536     if(chBin>0 && coord[kPt] > chBin) ch = 1;  // [+] track
1537     cen = 0; // centrality selection
1538     if(acen) cen = coord[kNdim]-1;
1539     npad = 0;             // no. of pads selection
1540     if(anp) npad = TMath::Min(coord[kNdim+1]-1, Int_t(kNpads));
1541
1542     if(debug){
1543       isel = ly*nCh*nCen*nNpad
1544             +ch*nCen*nNpad
1545             +cen*nNpad
1546             +npad;
1547       ioff=isel*4;
1548     } else {
1549       isel = ly; ioff = isel;
1550     }
1551     if(ioff>=ih){
1552       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1553       return kFALSE;
1554     }
1555     if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
1556       AliError(Form("Missing projection %d", ioff));
1557       return kFALSE;
1558     }
1559     if(strcmp(pr0->H()->GetName(), Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ch], ly, cen, npad))!=0){
1560       AliError(Form("Projection mismatch :: request[H%sClY%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ly, cen, npad, pr0->H()->GetName()));
1561       return kFALSE;
1562     }
1563     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1564   }
1565   if(HasDump3DFor(kCluster)){
1566     TDirectory *cwd = gDirectory;
1567     TFile::Open(Form("ResDump_%s.root", H->GetName()), "RECREATE");
1568     for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
1569       if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
1570       if(!pr0->H()) continue;
1571       TH3 *h3=(TH3*)pr0->H()->Clone();
1572       h3->Write();
1573     }
1574     gFile->Close();
1575     cwd->cd();
1576   }
1577
1578   TObjArray *arr(NULL);
1579   fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1580   arr->SetName(projName[Int_t(mc)]);
1581
1582   TH2 *h2(NULL);  Int_t jh(0);
1583   for(; ih--; ){
1584     if(!hp[ih].H()) continue;
1585     if(strchr(hp[ih].H()->GetName(), 'Q')){
1586       if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1587       arr->AddAt(h2, jh++);
1588       if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].H()->GetName())))) continue;
1589       arr->AddAt(h2, jh++);
1590     } else {
1591       if((h2 = hp[ih].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1592     }
1593   }
1594   Double_t m(0.), s(0.), trend(0.);
1595   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1596     for(Int_t ich(0); ich<nCh; ich++){
1597       for(Int_t icen(0); icen<nCen; icen++){
1598         /*!dy*/
1599         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1600           for(Int_t ipad(1); ipad<nNpad; ipad++){
1601             if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1602             (*pr0)+=(*pr1);
1603           }
1604           pr0->H()->SetNameTitle(Form("H%sClY%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1605           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1606           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1607         }
1608         /*!Q*/
1609         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1610           for(Int_t ipad(1); ipad<nNpad; ipad++){
1611             if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1612             (*pr0)+=(*pr1);
1613           }
1614           pr0->H()->SetNameTitle(Form("H%sClQ%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: Q Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1615           if((h2 = pr0->Projection2D(kNstat, kNcontours, 2, kFALSE))) arr->AddAt(h2, jh++);
1616           if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
1617           pr0->H()->SetName(Form("H%sClQS%c%d%d", mc?"MC":"", chName[ich], ily, icen));
1618           if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1619           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1620         }
1621         /*!YXTC*/
1622         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1623           for(Int_t ipad(1); ipad<nNpad; ipad++){
1624             if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1625             (*pr0)+=(*pr1);
1626           }
1627           pr0->H()->SetNameTitle(Form("H%sClYXTC%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay(x,TC) Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1628           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1629           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1630         }
1631         /*!YXPh*/
1632         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1633           for(Int_t ipad(1); ipad<nNpad; ipad++){
1634             if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1635             (*pr0)+=(*pr1);
1636           }
1637           pr0->H()->SetNameTitle(Form("H%sClYXPh%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay(x,#Phi) Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1638           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1639           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1640         }
1641       } // end centrality integration
1642       /*!dy*/
1643       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1644         pr0->H()->SetNameTitle(Form("H%sClY%c%d", mc?"MC":"", chName[ich], ily),
1645                               Form("Clusters[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
1646         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1647         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1648       }
1649       /*!Q*/
1650       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1651         pr0->H()->SetNameTitle(Form("H%sClQ%c%d", mc?"MC":"", chName[ich], ily),
1652                               Form("Clusters[%c]:: Q Ly[%d]", chSgn[ich], ily));
1653         if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
1654         pr0->H()->SetName(Form("H%sClQS%c%d", mc?"MC":"", chName[ich], ily));
1655         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1656         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1657       }
1658       /*!YXTC*/
1659       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1660         pr0->H()->SetNameTitle(Form("H%sClYXTC%c%d", mc?"MC":"", chName[ich], ily),
1661                               Form("Clusters[%c]:: #Deltay(x,TC) Ly[%d]", chSgn[ich], ily));
1662         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1663         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1664       }
1665       /*!YXPh*/
1666       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1667         pr0->H()->SetNameTitle(Form("H%sClYXPh%c%d", mc?"MC":"", chName[ich], ily),
1668                               Form("Clusters[%c]:: #Deltay(x,#Phi) Ly[%d]", chSgn[ich], ily));
1669         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1670         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1671       }
1672     } // end charge integration
1673     /*!dy*/
1674     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1675       pr0->H()->SetNameTitle(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1676       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1677       if((trend=pr0->GetTrendValue(1,&m,&s))>-100.) PutTrendValue(Form("%sClS%d", mc?"MC":"", ily), s);
1678     }
1679     /*!YXPh*/
1680     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1681       pr0->H()->SetNameTitle(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1682       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1683     }
1684
1685   }
1686   AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
1687   return kTRUE;
1688 }
1689
1690 //________________________________________________________
1691 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1692 {
1693 // Analyse tracklet
1694   const Int_t kNcontours(9);
1695   const Int_t kNstat(30);
1696   const Int_t kNstatQ(30);
1697   Int_t cidx=mc?kMCtracklet:kTracklet;
1698   if(fProj && fProj->At(cidx)) return kTRUE;
1699   if(!fContainer){
1700     AliError("Missing data container.");
1701     return kFALSE;
1702   }
1703   const Char_t *projName[] = {"hTracklet2Track", "hTracklet2MC"};
1704   THnSparse *H(NULL);
1705   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1706     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1707     return kTRUE;
1708   }
1709   const Int_t mdim(kNdim+8);
1710   Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimTrklt));
1711   Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
1712   TAxis *aa[mdim], *as(NULL), *ap(NULL), *ac(NULL); memset(aa, 0, sizeof(TAxis*) * mdim);
1713   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1714   if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC); // init species/charge selection
1715   if(ndim > Int_t(kPt))           ap = H->GetAxis(kPt);           // init pt selection
1716   if(ndim > Int_t(kNdim))         ac = H->GetAxis(kNdim);         // init centrality selection
1717   // calculate size depending on debug level
1718   const Int_t nCen(debug?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1719   const Int_t nPt(debug?Int_t(kNpt+2):1);
1720   const Int_t nSpc(1);//ndim>kNdimTrklt?fgkNbins[kSpeciesChgRC]:1);
1721   const Int_t nCh(debug?Int_t(kNcharge):1);
1722
1723   // build list of projections
1724   const Int_t nsel(AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer*(kNpt+2)*(AliPID::kSPECIES*kNcharge + 1));
1725   // define rebinning strategy
1726   const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1727   if(aa[1]->GetNbins()==180){ // old binning
1728     rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1729     rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1730     rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
1731     rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1732     rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1733   }
1734   AliTRDrecoProjection hp[kTrkltNproj]; TObjArray php(kTrkltNproj);
1735   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1736   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1737   const Char_t ptName[kNpt+2] = {'L', 'l', 'i', 'h', 'H'};
1738   const Char_t *ptCut[kNpt+2] = {"p_{t}[GeV/c]<0.5", "0.5<=p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "1.5<=p_{t}[GeV/c]<5.0", "p_{t}[GeV/c]>=5.0"};
1739 //  const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1740   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1741   for(Int_t icen(0); icen<nCen; icen++){
1742     for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1743       for(Int_t ipt(0); ipt<nPt; ipt++){
1744         for(Int_t isp(0); isp<nSpc; isp++){
1745           for(Int_t ich(0); ich<nCh; ich++){
1746             isel++; // new selection
1747             hp[ih].Build(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1748                         Form("Tracklets[%s%c]:: #Deltay{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1749                         kEta, kPhi, kYrez, aa);
1750             //hp[ih].SetShowRange(-0.1,0.1);
1751             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1752             php.AddLast(&hp[ih++]); np[isel]++;
1753             hp[ih].Build(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1754                         Form("Tracklets[%s%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1755                         kEta, kPhi, kPrez, aa);
1756             //hp[ih].SetShowRange(-0.5,0.5);
1757             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1758             php.AddLast(&hp[ih++]); np[isel]++;
1759             hp[ih].Build(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1760                         Form("Tracklets[%s%c]:: dQdl{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1761                         kEta, kPhi, kZrez, aa);
1762             hp[ih].SetShowRange(1.,2.3);
1763             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1764             php.AddLast(&hp[ih++]); np[isel]++;
1765             hp[ih].Build(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1766                         Form("Tracklets[%s%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1767                         kEta, kPhi, kNdim+1, aa);
1768             hp[ih].SetShowRange(30., 70.);
1769             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1770             php.AddLast(&hp[ih++]); np[isel]++;
1771           }
1772         }
1773         if(ndim==kNdimTrklt) continue;
1774
1775         isel++; // new selection
1776         hp[ih].Build(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1777                     Form("Tracklets[RC]:: #Deltaz{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1778                     kEta, kPhi, kZrez, aa);
1779   //      hp[ih].SetShowRange(-0.1,0.1);
1780         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1781         php.AddLast(&hp[ih++]); np[isel]++;
1782 /*        hp[ih].Build(Form("H%sTrkltRCY%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1783                     Form("Tracklets[RC]:: #Deltay{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1784                     kEta, kPhi, kYrez, aa);
1785         //hp[ih].SetShowRange(-0.1,0.1);
1786         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1787         php.AddLast(&hp[ih++]); np[isel]++;
1788         hp[ih].Build(Form("H%sTrkltRCPh%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1789                     Form("Tracklets[RC]:: #Delta#phi{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1790                     kEta, kPhi, kPrez, aa);
1791         //hp[ih].SetShowRange(-0.1,0.1);
1792         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1793         php.AddLast(&hp[ih++]); np[isel]++;
1794         hp[ih].Build(Form("H%sTrkltRCQ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1795                     Form("Tracklets[RC]:: dQdl{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1796                     kEta, kPhi, kNdim, aa);
1797         //hp[ih].SetShowRange(-0.1,0.1);
1798         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1799         php.AddLast(&hp[ih++]); np[isel]++;*/
1800       }
1801     }
1802   }
1803   AliInfo(Form("Build %3d 3D %s projections.", ih, mc?"MC":""));
1804
1805   AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1806   Int_t ly(0), ch(0), sp(2), rcBin(as?as->FindBin(0.):-1), pt(0), cen(0), ioff(0), jspc(nSpc*nCh+1), kspc(nSpc*nCh*4+1);
1807   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1808     v = H->GetBinContent(ib, coord);
1809     if(v<1.) continue;
1810     ly = coord[kBC]-1; // layer selection
1811     // charge selection
1812     ch = 0; sp=0;// [e-] track [dafault]
1813     if(rcBin>0){ // debug mode in which species/charge are also saved
1814       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1815       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1816       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1817     }
1818     // pt selection
1819     pt = 0; // low pt
1820     if(ap) pt = TMath::Min(coord[kPt], Int_t(kNpt)+1);
1821     // centrality selection
1822     cen = 0; // default
1823     if(ac) cen = coord[kNdim]-1;
1824     // global selection
1825     if(ndim==kNdimTrklt){
1826       ioff = ly*4;
1827       isel = ly;
1828     } else {
1829       isel = cen*AliTRDgeometry::kNlayer*nPt*jspc+ly*nPt*jspc+pt*jspc; isel+=sp<0?(nSpc*nCh):ch;
1830       ioff = cen*AliTRDgeometry::kNlayer*nPt*kspc+ly*nPt*kspc+pt*kspc; ioff+=sp<0?((nSpc*nCh)*4):(ch*4);
1831     }
1832     if(ioff>=ih){
1833       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1834       return kFALSE;
1835     }
1836     if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
1837       AliError(Form("Missing projection %d", ioff));
1838       return kFALSE;
1839     }
1840     if(sp>=0){
1841       if(strcmp(pr0->H()->GetName(), Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ch], ptName[pt], sp, ly, cen))!=0){
1842         AliError(Form("Projection mismatch :: request[H%sTrkltY%c%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, ly, cen, pr0->H()->GetName()));
1843         return kFALSE;
1844       }
1845     } else {
1846       if(strcmp(pr0->H()->GetName(), Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[pt], ly, cen))!=0){
1847         AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c%d%d] found[%s]", mc?"MC":"", ptName[pt], ly, cen, pr0->H()->GetName()));
1848         return kFALSE;
1849       }
1850     }
1851     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1852   }
1853   if(HasDump3DFor(kTracklet)){
1854     TDirectory *cwd = gDirectory;
1855     TFile::Open(Form("ResDump_%s.root", H->GetName()), "RECREATE");
1856     for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
1857       if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
1858       if(!pr0->H()) continue;
1859       TH3 *h3=(TH3*)pr0->H()->Clone();
1860       h3->Write();
1861     }
1862     gFile->Close();
1863     cwd->cd();
1864   }
1865
1866   TObjArray *arr(NULL);
1867   fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1868   arr->SetName(projName[Int_t(mc)]);
1869
1870   TH2 *h2(NULL); Int_t jh(0);
1871   for(; ih--; ){
1872     if(!hp[ih].H()) continue;
1873     Int_t mid(0), nstat(kNstat);
1874     if(strchr(hp[ih].H()->GetName(), 'Q')){ mid=2; nstat=kNstatQ;}
1875     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1876     arr->AddAt(h2, jh++);
1877   }
1878   // build combined performance plots
1879   Double_t m(0.), s(0.), trend(0.);
1880   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1881     for(Int_t ich(0); ich<nCh; ich++){
1882       for(Int_t icen(0); icen<nCen; icen++){
1883         for(Int_t ipt(0); ipt<nPt; ipt++){
1884           /*!dy*/
1885           if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1886             for(Int_t isp(1); isp<nSpc; isp++){
1887               if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1888               (*pr0)+=(*pr1);
1889             }
1890             pr0->H()->SetNameTitle(Form("H%sTrkltY%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1891                                       Form("Tracklets[%c]:: #Deltay{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1892             if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1893             if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1894           }
1895           /*!dphi*/
1896           if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1897             for(Int_t isp(1); isp<nSpc; isp++){
1898               if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1899               (*pr0)+=(*pr1);
1900             }
1901             pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1902                                       Form("Tracklets[%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1903             if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1904             if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1905           }
1906           /*!dQ/dl*/
1907           if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1908             for(Int_t isp(1); isp<nSpc; isp++){
1909               if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1910               (*pr0)+=(*pr1);
1911             }
1912             pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1913                                       Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1914             if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1915             pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1916                                       Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1917             pr0->SetShowRange(2.4, 5.1);
1918             if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1919             if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1920           }
1921           /*!TB occupancy*/
1922           if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1923             for(Int_t isp(1); isp<nSpc; isp++){
1924               if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1925               (*pr0)+=(*pr1);
1926             }
1927             pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1928                                       Form("Tracklets[%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1929             if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1930             if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1931           }
1932         } // end pt integration
1933         /*!dy*/
1934         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1935           pr0->H()->SetNameTitle(Form("H%sTrkltY%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1936                                 Form("Tracklets[%c]:: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1937           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1938           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1939         }
1940         /*!dphi*/
1941         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1942           pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1943                                 Form("Tracklets[%c]:: #Delta#phi Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1944           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1945           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1946         }
1947         /*!dQ/dl*/
1948         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1949           pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1950                                 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1951           pr0->SetShowRange(1.,2.3);
1952           if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1953           pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1954                                 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1955           pr0->SetShowRange(2.4,5.1);
1956           if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1957           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1958         }
1959         /*!TB occupancy*/
1960         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1961           pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1962                                 Form("Tracklets[%c]:: OccupancyTB Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1963           if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1964           if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1965         }
1966       } // end centrality integration
1967       /*!dy*/
1968       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1969         pr0->H()->SetNameTitle(Form("H%sTrkltY%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Deltay Ly[%d]", chSgn[ich], ily));
1970         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1971         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1972       }
1973       /*!dphi*/
1974       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1975         pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Delta#phi Ly[%d]", chSgn[ich], ily));
1976         pr0->SetShowRange(-.9,.9);
1977         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1978         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1979       }
1980       /*!dQ/dl*/
1981       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1982         pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1983         pr0->SetShowRange(1.,2.3);
1984         if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1985         pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1986         pr0->SetShowRange(2.4,5.1);
1987         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1988         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1989       }
1990       /*!TB occupancy*/
1991       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1992         pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: OccupancyTB Ly[%d]", chSgn[ich], ily));
1993         if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1994         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1995       }
1996     } // end charge integration
1997     /*!dy*/
1998     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1999       pr0->H()->SetNameTitle(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets :: #Deltay Ly[%d]", ily));
2000       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2001       if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2002         PutTrendValue(Form("%sTrkltY%d", mc?"MC":"", ily), trend);
2003         PutTrendValue(Form("%sTrkltYS%d", mc?"MC":"", ily), s);
2004       }
2005     }
2006     /*!dphi*/
2007     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
2008       pr0->H()->SetNameTitle(Form("H%sTrkltPh%d", mc?"MC":"", ily), Form("Tracklets :: #Delta#phi Ly[%d]", ily));
2009       pr0->SetShowRange(-.45,.45);
2010       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2011     }
2012     /*!dQdl*/
2013     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
2014       pr0->H()->SetNameTitle(Form("H%sTrkltQ%d", mc?"MC":"", ily), Form("Tracklets :: dQdl Ly[%d]", ily));
2015       pr0->SetShowRange(1.,2.3);
2016       if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2017       pr0->H()->SetName(Form("H%sTrkltQS%d", mc?"MC":"", ily));
2018       pr0->SetShowRange(2.4,5.1);
2019       if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2020     }
2021     /*!TB occupancy*/
2022     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
2023       pr0->H()->SetNameTitle(Form("H%sTrkltTB%d", mc?"MC":"", ily), Form("Tracklets :: OccupancyTB Ly[%d]", ily));
2024       if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
2025     }
2026
2027     /*! Row Cross processing*/
2028     for(Int_t icen(0); icen<nCen; icen++){
2029       /*!RC dz*/
2030       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, icen)))){
2031         for(Int_t ipt(0); ipt<kNpt; ipt++){
2032           if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen)))) continue;
2033           (*pr0)+=(*pr1);
2034         }
2035         pr0->H()->SetNameTitle(Form("H%sTrkltRCZ%d%d", mc?"MC":"", ily, icen), Form("Tracklets[RC]:: #Deltaz Ly[%d] Cen[%s]", ily, cenName[icen]));
2036         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2037         if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))) (*pr1)+=(*pr0);
2038       }
2039     } // end centrality integration for row cross
2040     /*!RC dz*/
2041     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))){
2042       pr0->H()->SetNameTitle(Form("H%sTrkltRCZ%d", mc?"MC":"", ily), Form("Tracklets[RC] :: #Deltaz Ly[%d]", ily));
2043       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2044     }
2045   } // end layer loop
2046   AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
2047   return kTRUE;
2048 }
2049
2050 //________________________________________________________
2051 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
2052 {
2053 // Analyse track in
2054   const Int_t kNcontours(9);
2055   const Int_t kNstat(30);
2056   Int_t cidx=mc?kMCtrackIn:kTrackIn;
2057   if(fProj && fProj->At(cidx)) return kTRUE;
2058   if(!fContainer){
2059     AliError("Missing data container.");
2060     return kFALSE;
2061   }
2062   const Char_t *projName[] = {"hTracklet2TRDin", "hTRDin2MC"};
2063   THnSparse *H(NULL);
2064   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
2065     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
2066     return kTRUE;
2067   }
2068
2069   const Int_t mdim(kNdim+3);
2070   Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
2071   Int_t ndim(H->GetNdimensions());
2072   TAxis *aa[mdim], *as(NULL), *ap(NULL), *apt(NULL), *ax(NULL), *abf(NULL); memset(aa, 0, sizeof(TAxis*) * (mdim));
2073   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2074   if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC);
2075   if(ndim > Int_t(kPt))          apt = H->GetAxis(kPt);
2076   if(ndim > Int_t(kNdim))         ap = H->GetAxis(kNdim);
2077   if(ndim > Int_t(kNdim)+1)       ax = H->GetAxis(kNdim+1);
2078   if(ndim > Int_t(kNdim)+2)      abf = H->GetAxis(kNdim+2);
2079   //AliInfo(Form("Using : Species[%c] Pt[%c] BunchFill[%c]", as?'y':'n', ap?'y':'n', abf?'y':'n'));
2080   const Int_t nPt(ndim>Int_t(kNdimTrkIn)?Int_t(kNpt+2):1);
2081
2082   // build list of projections
2083   const Int_t nsel((kNpt+2)*(kNspc*kNcharge + 1));
2084   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2085   const Char_t *spcName[kNspc] = {"e", "#mu#pi", "Kp"};
2086   const Char_t ptName[kNpt+2] = {'L', 'l', 'i', 'h', 'H'};
2087   const Char_t *ptCut[kNpt+2] = {"p_{t}[GeV/c]<0.5", "0.5<=p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "1.5<=p_{t}[GeV/c]<5.0", "p_{t}[GeV/c]>=5.0"};
2088   const Char_t *pCut[kNpt+2] = {"p[GeV/c]<0.5", "0.5<=p[GeV/c]<0.8", "0.8<=p[GeV/c]<1.5", "1.5<=p[GeV/c]<5.0", "p[GeV/c]>=5.0"};
2089   // define rebinning strategy
2090   const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
2091   if(aa[1]->GetNbins()==180){ // old binning
2092     rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
2093     rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
2094     rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
2095     rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
2096     rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
2097   }
2098   AliTRDrecoProjection hp[kMCTrkInNproj]; TObjArray php(kMCTrkInNproj+2);
2099   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2100   // define list of projections
2101   for(Int_t ipt(0); ipt<nPt; ipt++){
2102     for(Int_t isp(0); isp<kNspc; isp++){
2103       for(Int_t ich(0); ich<kNcharge; ich++){
2104         isel++; // new selection
2105         hp[ih].Build(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2106                      Form("TrackIn[%s%c]:: #Deltay{%s}", spcName[isp], chSgn[ich], ptCut[ipt]),
2107                      kEta, kPhi, kYrez, aa);
2108         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2109         php.AddLast(&hp[ih++]); np[isel]++;
2110         hp[ih].Build(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2111                      Form("TrackIn[%s%c]:: #Delta#phi{%s}", spcName[isp], chSgn[ich], ptCut[ipt]),
2112                      kEta, kPhi, kPrez, aa);
2113         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2114         php.AddLast(&hp[ih++]); np[isel]++;
2115         hp[ih].Build(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2116                      Form("TrackIn[%s%c]:: dQdl {%s}", spcName[isp], chSgn[ich], pCut[ipt]),
2117                      kEta, kPhi, kZrez, aa);
2118         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2119         php.AddLast(&hp[ih++]); np[isel]++;
2120         if(!ax) continue;
2121         hp[ih].Build(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2122                      Form("TrackIn[%s%c]:: #Deltax{%s}", spcName[isp], chSgn[ich], ptCut[ipt]),
2123                      kEta, kPhi, kNdim+1, aa);
2124         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2125         php.AddLast(&hp[ih++]); np[isel]++;
2126       }
2127     }
2128     isel++; // RC tracks
2129     hp[ih].Build(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt]),
2130                   Form("TrackIn[RC]:: #Deltaz{%s}", ptCut[ipt]),
2131                   kEta, kPhi, kZrez, aa);
2132     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2133     php.AddLast(&hp[ih++]); np[isel]++;
2134     hp[ih].Build(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt]),
2135                   Form("TrackIn[RC]:: #Deltay{%s}", ptCut[ipt]),
2136                   kEta, kPhi, kYrez, aa);
2137     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2138     php.AddLast(&hp[ih++]); np[isel]++;
2139     hp[ih].Build(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt]),
2140                   Form("TrackIn[RC]:: #Delta#phi{%s}", ptCut[ipt]),
2141                   kEta, kPhi, kPrez, aa);
2142     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2143     php.AddLast(&hp[ih++]); np[isel]++;
2144     if(!ax) continue;
2145     hp[ih].Build(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt]),
2146                   Form("TrackIn[RC]:: #Deltax{%s}", ptCut[ipt]),
2147                   kEta, kPhi, kNdim+1, aa);
2148     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2149     php.AddLast(&hp[ih++]); np[isel]++;
2150   }
2151   AliInfo(Form("Build %3d 3D %s projections.", ih, mc?"MC":""));
2152
2153   // fill projections
2154   Int_t ch(0), pt(0), p(0), sp(1), rcBin(as?as->FindBin(0.):-1), ioff(0), joff(0), jsel(0);
2155   AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
2156   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2157     v = H->GetBinContent(ib, coord);
2158     if(v<1.) continue;
2159     if(fBCbinTOF>0 && coord[kBC]!=fBCbinTOF) continue; // TOF bunch cross cut
2160     if(fBCbinFill>0 && abf && coord[kNdim+2]!=fBCbinTOF) continue; // Fill bunch cut
2161     if(coord[kBC]==3) continue;
2162     // charge selection
2163     ch = 0; sp=1;// [pi-] track
2164     if(rcBin>0){ // debug mode in which species are also saved
2165       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2166       // take care of old data format (2*AliPID::kSPECIES+1)
2167       if(as->GetNbins() == kNcharge*AliPID::kSPECIES+1){
2168         if(sp>2) sp=2;
2169         else if(sp==2) sp=1;
2170       }
2171       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
2172       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
2173     }
2174     // pt selection
2175     pt = 0; p = 0; // low pt
2176     if(apt) pt = TMath::Min(coord[kPt], Int_t(kNpt)+1);
2177     if(ap ) p  = TMath::Min(coord[kNdim], Int_t(kNpt)+1);
2178     // global selection
2179     isel = pt*(kNspc*kNcharge+1); isel+=(ch==2?(kNspc*kNcharge):(sp*kNcharge+ch)); ioff = isel*(ax?4:3);
2180     jsel = p *(kNspc*kNcharge+1); jsel+=(ch==2?(kNspc*kNcharge):(sp*kNcharge+ch)); joff = jsel*(ax?4:3);
2181     if(ioff>=ih){
2182       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
2183       return kFALSE;
2184     }
2185     if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
2186       AliError(Form("Missing projection %d", ioff));
2187       return kFALSE;
2188     }
2189     if(ch<2){
2190       if(strcmp(pr0->H()->GetName(), Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ch], ptName[pt], sp))!=0){
2191         AliError(Form("Projection mismatch :: request[H%sTrkInY%c%c%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, pr0->H()->GetName()));
2192         return kFALSE;
2193       }
2194     } else {
2195       if(strcmp(pr0->H()->GetName(), Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[pt]))!=0){
2196         AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c] found[%s]", mc?"MC":"", ptName[pt], pr0->H()->GetName()));
2197         return kFALSE;
2198       }
2199     }
2200     AliDebug(2, Form("Found %s for selection sp[%d] ch[%d] pt[%d]", pr0->H()->GetName(), sp, ch, pt));
2201     for(Int_t jh(0); jh<np[isel]; jh++){
2202       if(ch<2 && jh==2) ((AliTRDrecoProjection*)php.At(joff+jh))->Increment(coord, v); // special care for dQdl=f(p)
2203       else ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v); // all = f(p_t)
2204     }
2205   }
2206   if(HasDump3DFor(kTrackIn)){
2207     TDirectory *cwd = gDirectory;
2208     TFile::Open(Form("ResDump_%s.root", H->GetName()), "RECREATE");
2209     for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
2210       if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
2211       if(!pr0->H()) continue;
2212       TH3 *h3=(TH3*)pr0->H()->Clone();
2213       h3->Write();
2214     }
2215     gFile->Close();
2216     cwd->cd();
2217   }
2218
2219   TObjArray *arr(NULL);
2220   fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx);
2221   arr->SetName(projName[Int_t(mc)]);
2222
2223   TH2 *h2(NULL); Int_t jh(0);
2224   for(; ih--; ){
2225     if(!hp[ih].H()) continue;
2226     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2227     arr->AddAt(h2, jh++);
2228   }
2229   // build combined performance plots
2230   // combine up the tree of projections
2231   Double_t m(0.), s(0.), trend(0.);
2232   AliTRDrecoProjection xlow[2], specY[kNcharge*kNspc], specPh[kNcharge*kNspc], specQ[kNcharge*kNspc];
2233   for(Int_t ich(0); ich<kNcharge; ich++){
2234     // PID dependency - summation over pt
2235     for(Int_t isp(0); isp<kNspc; isp++){
2236       /*!dy*/
2237       Int_t idx(ich*kNspc+isp);
2238       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2239         specY[idx] = (*pr0);
2240         specY[idx].SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2241         specY[idx].H()->SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp),
2242                               Form("TrackIn[%s%c]:: #Deltay", spcName[isp], chSgn[ich]));
2243         if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2244           PutTrendValue(Form("%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), trend);
2245           PutTrendValue(Form("%sTrkInYS%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), s);
2246         }
2247         for(Int_t ipt(1); ipt<nPt; ipt++){
2248           if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2249           if((trend=pr1->GetTrendValue(1,&m,&s))>-100.){
2250             PutTrendValue(Form("%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), trend);
2251             PutTrendValue(Form("%sTrkInYS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), s);
2252           }
2253           specY[idx]+=(*pr1);
2254         }
2255         php.AddLast(&specY[idx]);
2256         if((h2 = specY[idx].Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2257         if((trend=specY[idx].GetTrendValue(1,&m,&s))>-100.){
2258           PutTrendValue(Form("%sTrkInY%c%d", mc?"MC":"", chName[ich], isp), trend);
2259           PutTrendValue(Form("%sTrkInYS%c%d", mc?"MC":"", chName[ich], isp), s);
2260         }
2261         if((h2 = (TH2*)gDirectory->Get(Form("%sEn", specY[idx].H()->GetName())))) arr->AddAt(h2, jh++);
2262         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specY[idx];
2263       }
2264       /*!dphi*/
2265       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2266         specPh[idx] = (*pr0);
2267         specPh[idx].SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2268         specPh[idx].H()->SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp),
2269                               Form("TrackIn[%s%c]:: #Delta#phi", spcName[isp], chSgn[ich]));
2270         specPh[idx].SetShowRange(-1.5, 1.5);
2271         if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2272           PutTrendValue(Form("%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), trend);
2273           PutTrendValue(Form("%sTrkInPhS%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), s);
2274         }
2275         for(Int_t ipt(1); ipt<nPt; ipt++){
2276           if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2277           if((trend=pr1->GetTrendValue(1,&m,&s))>-100.){
2278             PutTrendValue(Form("%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), trend);
2279             PutTrendValue(Form("%sTrkInPhS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), s);
2280           }
2281           specPh[idx]+=(*pr1);
2282         }
2283         php.AddLast(&specPh[idx]);
2284         if((h2 = specPh[idx].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2285         if((trend=specPh[idx].GetTrendValue(1,&m,&s))>-100.){
2286           PutTrendValue(Form("%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp), trend);
2287           PutTrendValue(Form("%sTrkInPhS%c%d", mc?"MC":"", chName[ich], isp), s);
2288         }
2289         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specPh[idx];
2290       }
2291       /*!dQdl*/
2292       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2293         specQ[idx] = (*pr0);
2294         specQ[idx].SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp), "Sum over p");
2295         specQ[idx].H()->SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp),
2296                               Form("TrackIn[%s%c]:: dQdl", spcName[isp], chSgn[ich]));
2297         specQ[idx].SetShowRange(-2.2, -1.75);
2298         specQ[idx].H()->GetZaxis()->SetTitle("dQdl [a.u.]");
2299         if((trend = pr0->GetTrendValue(2, &m))>-100.){
2300           PutTrendValue(Form("%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), trend);
2301           PutTrendValue(Form("%sTrkInQS%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), m);
2302         }
2303         for(Int_t ipt(1); ipt<nPt; ipt++){
2304           if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2305           if((trend=pr1->GetTrendValue(2, &m))>-100.){
2306             PutTrendValue(Form("%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), trend);
2307             PutTrendValue(Form("%sTrkInQS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), m);
2308           }
2309           specQ[idx]+=(*pr1);
2310         }
2311         php.AddLast(&specQ[idx]);
2312         if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2313         specQ[idx].H()->SetName(Form("H%sTrkInQS%c%d", mc?"MC":"", chName[ich], isp));
2314         specQ[idx].SetShowRange(-1.85, -1.4);
2315         if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2316         if((trend=specQ[idx].GetTrendValue(2, &m))>-100.){
2317           PutTrendValue(Form("%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp), trend);
2318           PutTrendValue(Form("%sTrkInQS%c%d", mc?"MC":"", chName[ich], isp), m);
2319         }
2320         if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specQ[idx];
2321       }
2322     } // end PID loop for pt integration
2323
2324     // pt dependency - summation over PID
2325     for(Int_t ipt(0); ipt<nPt; ipt++){
2326       /*!dy*/
2327       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2328         for(Int_t isp(1); isp<kNspc; isp++){
2329           if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2330           (*pr0)+=(*pr1);
2331         }
2332         pr0->H()->SetNameTitle(Form("H%sTrkInY%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2333                                   Form("TrackIn[%c]:: #Deltay{%s}", chSgn[ich], ptCut[ipt]));
2334         pr0->SetShowRange(-0.3, 0.3);
2335         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2336         if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2337           PutTrendValue(Form("%sTrkInY%c%c", mc?"MC":"", chName[ich], ptName[ipt]), trend);
2338           PutTrendValue(Form("%sTrkInYS%c%c", mc?"MC":"", chName[ich], ptName[ipt]), s);
2339         }
2340         if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2341       }
2342       /*!dphi*/
2343       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2344         for(Int_t isp(1); isp<kNspc; isp++){
2345           if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2346           (*pr0)+=(*pr1);
2347         }
2348         pr0->H()->SetNameTitle(Form("H%sTrkInPh%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2349                                   Form("TrackIn[%c]:: #Delta#phi{%s}", chSgn[ich], ptCut[ipt]));
2350         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2351         if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2352           PutTrendValue(Form("%sTrkInPh%c%c", mc?"MC":"", chName[ich], ptName[ipt]), trend);
2353           PutTrendValue(Form("%sTrkInPhS%c%c", mc?"MC":"", chName[ich], ptName[ipt]), s);
2354         }
2355         if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2356       }
2357       /*!dx*/
2358       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2359         for(Int_t isp(1); isp<kNspc; isp++){
2360           if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2361           (*pr0)+=(*pr1);
2362         }
2363         pr0->H()->SetNameTitle(Form("H%sTrkInX%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2364                                   Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[ipt]));
2365         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2366         PutTrendValue(Form("%sTrkInX%c%c", mc?"MC":"", chName[ich], ptName[ipt]), pr0->GetTrendValue(1));
2367         if(!ipt){
2368           xlow[ich] = (*pr0);
2369           xlow[ich].SetNameTitle(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5),
2370                                  Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[0]));
2371           php.AddLast(&xlow[ich]);
2372         }
2373         if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2374       }
2375     } // end pt loop for PID integration
2376
2377     /*!dy*/
2378     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2379       pr0->H()->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", chName[ich]),
2380                             Form("TrackIn[%c]:: #Deltay", chSgn[ich]));
2381       pr0->SetShowRange(-0.3, 0.3);
2382       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2383       if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2384       if((trend=pr0->GetTrendValue(1, &m, &s))>-100.){
2385         PutTrendValue(Form("%sTrkInY%c", mc?"MC":"", chName[ich]), trend);
2386         PutTrendValue(Form("%sTrkInYS%c", mc?"MC":"", chName[ich]), s);
2387       }
2388       if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2389     }
2390     /*!dy high pt*/
2391     if(ich && (pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[3], 0)))){
2392       if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[3], 0)))){
2393         (*pr0)+=(*pr1);
2394         pr0->H()->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", ptName[3]), Form("TrackIn :: #Deltay{%s}", ptCut[3]));
2395         pr0->SetShowRange(-0.3, 0.3);
2396         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2397       }
2398     }
2399     /*!dphi*/
2400     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2401       pr0->H()->SetNameTitle(Form("H%sTrkInPh%c", mc?"MC":"", chName[ich]),
2402                             Form("TrackIn[%c]:: #Delta#phi", chSgn[ich]));
2403       pr0->SetShowRange(-1., 1.);
2404       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2405       if((trend=pr0->GetTrendValue(1, &m, &s))>-100.){
2406         PutTrendValue(Form("%sTrkInPh%c", mc?"MC":"", chName[ich]), trend);
2407         PutTrendValue(Form("%sTrkInPhS%c", mc?"MC":"", chName[ich]), s);
2408       }
2409       if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2410     }
2411     /*!dx*/
2412     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2413       pr0->H()->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", chName[ich]),
2414                             Form("TrackIn[%c]:: #Deltax", chSgn[ich]));
2415       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2416       PutTrendValue(Form("%sTrkInX%c", mc?"MC":"", chName[ich]), pr0->GetTrendValue(1));
2417       if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2418     }
2419     /*!dx low pt*/
2420     if(ich && (pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[1], 5)))){
2421       if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[1], 5)))){
2422         (*pr0)+=(*pr1);
2423         pr0->H()->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", ptName[1]), Form("TrackIn :: #Deltax{%s}", ptCut[1]));
2424         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2425       }
2426     }
2427   } // end charge loop
2428
2429   for(Int_t isp(0); isp<kNspc; isp++){
2430     /*!dy*/
2431     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))){
2432       pr0->H()->SetNameTitle(Form("H%sTrkInY%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Deltay", spcName[isp]));
2433       pr0->SetShowRange(-0.3, 0.3);
2434       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2435       PutTrendValue(Form("%sTrkInY%d", mc?"MC":"", isp), pr0->GetTrendValue(1));
2436       if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2437     }
2438     /*!dphi*/
2439     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))){
2440       pr0->H()->SetNameTitle(Form("H%sTrkInPh%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Delta#phi", spcName[isp]));
2441       pr0->SetShowRange(-1., 1.);
2442       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2443       PutTrendValue(Form("%sTrkInPh%d", mc?"MC":"", isp), pr0->GetTrendValue(1));
2444     }
2445     /*!dQdl*/
2446     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))){
2447       pr0->H()->SetNameTitle(Form("H%sTrkInQ%d", mc?"MC":"", isp), Form("TrackIn[%s] :: dQdl", spcName[isp]));
2448       pr0->SetShowRange(-2.2, -1.75);
2449       if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2450       pr0->H()->SetName(Form("H%sTrkInQS%d", mc?"MC":"", isp));
2451       pr0->SetShowRange(-1.85, -1.4);
2452       if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2453       if((trend=pr0->GetTrendValue(2, &m))>-100.){
2454         PutTrendValue(Form("%sTrkInQ%d", mc?"MC":"", isp), trend);
2455         PutTrendValue(Form("%sTrkInQS%d", mc?"MC":"", isp), m);
2456       }
2457     }
2458   } // end PID processing
2459
2460   /*!dy*/
2461   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2462     pr0->H()->SetNameTitle(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: #Deltay");
2463     pr0->SetShowRange(-0.3, 0.3);
2464     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2465     if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2466   }
2467   /*!dphi*/
2468   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2469     pr0->H()->SetNameTitle(Form("H%sTrkInPh", mc?"MC":""), "TrackIn :: #Delta#phi");
2470     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2471   }
2472   /*!dx*/
2473   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2474     pr0->H()->SetNameTitle(Form("H%sTrkInX", mc?"MC":""), "TrackIn :: #Deltax");
2475     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2476   }
2477
2478   // Row Cross processing
2479   /*!RC dz*/
2480   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[0])))){
2481     for(Int_t ipt(0); ipt<kNpt; ipt++){
2482       if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt])))) continue;
2483       (*pr0)+=(*pr1);
2484     }
2485     pr0->H()->SetNameTitle(Form("H%sTrkInRCZ", mc?"MC":""), "TrackIn[RC]:: #Deltaz");
2486     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2487     if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2488       PutTrendValue(Form("%sTrkInRCZ", mc?"MC":""), trend);
2489       PutTrendValue(Form("%sTrkInRCZS", mc?"MC":""), s);
2490     }
2491   }
2492   /*!RC dy*/
2493   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[0])))){
2494     for(Int_t ipt(0); ipt<kNpt; ipt++){
2495       if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt])))) continue;
2496       (*pr0)+=(*pr1);
2497     }
2498     pr0->H()->SetNameTitle(Form("H%sTrkInRCY", mc?"MC":""), "TrackIn[RC]:: #Deltay");
2499     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2500   }
2501   /*!RC dphi*/
2502   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[0])))){
2503     for(Int_t ipt(0); ipt<kNpt; ipt++){
2504       if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt])))) continue;
2505       (*pr0)+=(*pr1);
2506     }
2507     pr0->H()->SetNameTitle(Form("H%sTrkInRCPh", mc?"MC":""), "TrackIn[RC]:: #Delta#phi");
2508     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2509   }
2510   /*!RC dx*/
2511   if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[0])))){
2512     for(Int_t ipt(0); ipt<kNpt; ipt++){
2513       if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt])))) continue;
2514       (*pr0)+=(*pr1);
2515     }
2516     pr0->H()->SetNameTitle(Form("H%sTrkInRCX", mc?"MC":""), "TrackIn[RC]:: #Deltax");
2517     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2518   }
2519   AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
2520   return kTRUE;
2521 }
2522
2523
2524 //________________________________________________________
2525 Bool_t AliTRDresolution::MakeProjectionTrack()
2526 {
2527 // Analyse tracklet
2528   const Int_t kNcontours(9);
2529   const Int_t kNstat(30);
2530   Int_t cidx(kMCtrack);
2531   if(fProj && fProj->At(cidx)) return kTRUE;
2532   if(!fContainer){
2533     AliError("Missing data container.");
2534     return kFALSE;
2535   }
2536   THnSparse *H(NULL);
2537   if(!(H = (THnSparse*)fContainer->FindObject("hTRD2MC"))){
2538     AliError("Missing/Wrong data @ hTRD2MC.");
2539     return kTRUE;
2540   }
2541   Int_t ndim(H->GetNdimensions());
2542   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
2543   TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
2544   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2545   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
2546   if(ndim > kPt) ap = H->GetAxis(kPt);
2547
2548   // build list of projections
2549   const Int_t nsel(AliTRDgeometry::kNlayer*kNpt*AliPID::kSPECIES*7);//, npsel(3);
2550   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2551   const Char_t ptName[kNpt] = {'l', 'i', 'h'};
2552   const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"};
2553   // define rebinning strategy
2554   const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
2555   if(aa[1]->GetNbins()==180){ // old binning
2556     rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
2557     rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
2558     rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
2559     rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
2560     rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
2561   }
2562   AliTRDrecoProjection hp[kTrkNproj]; TObjArray php(kTrkNproj);
2563   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2564   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2565     for(Int_t ipt(0); ipt<kNpt; ipt++){
2566       for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2567         for(Int_t ich(0); ich<kNcharge; ich++){
2568           isel++; // new selection
2569           hp[ih].Build(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2570                        Form("Tracks[%s%c]:: #Deltay{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2571                        kEta, kPhi, kYrez, aa);
2572           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2573           php.AddLast(&hp[ih++]); np[isel]++;
2574           hp[ih].Build(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2575                        Form("Tracks[%s%c]:: #Delta#phi{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2576                        kEta, kPhi, kPrez, aa);
2577           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2578           php.AddLast(&hp[ih++]); np[isel]++;
2579           hp[ih].Build(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2580                        Form("Tracks[%s%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2581                        kEta, kPhi, kNdim, aa);
2582           hp[ih].SetShowRange(0.,10.);
2583           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2584           php.AddLast(&hp[ih++]); np[isel]++;
2585         }
2586       }
2587       isel++; // new selection
2588       hp[ih].Build(Form("HMCTrkZ%c%d", ptName[ipt], ily),
2589                     Form("Tracks[RC]:: #Deltaz{%s} Ly[%d]", ptCut[ipt], ily),
2590                     kEta, kPhi, kZrez, aa);
2591       hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2592       php.AddLast(&hp[ih++]); np[isel]++;
2593     }
2594   }
2595   AliInfo(Form("Build %3d 3D projections.", ih));
2596
2597   Int_t ly(0), ch(0), pt(0), sp(2), rcBin(as?as->FindBin(0.):-1);
2598   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2599     v = H->GetBinContent(ib, coord);
2600     if(v<1.) continue;
2601     ly = coord[kBC]-1; // layer selection
2602     // charge selection
2603     ch=0; sp=2;// [pi-] track [dafault]
2604     if(rcBin>0){ // debug mode in which species are also saved
2605       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2606       if(coord[kSpeciesChgRC] > rcBin) ch = 1;        // [+] track
2607       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
2608     }
2609     // pt selection
2610     pt = 0; // low pt [default]
2611     if(ap) pt = coord[kPt]-1;
2612     // global selection
2613     Int_t ioff = ly*kNpt*31+pt*31; ioff+=3*(sp<0?10:(sp*kNcharge+ch));
2614     isel = ly*kNpt*11+pt*11; isel+=sp<0?10:(sp*kNcharge+ch);
2615     AliDebug(4, Form("SELECTION[%d] :: ch[%c] pt[%c] sp[%d] ly[%d]\n", np[isel], ch==2?'Z':chName[ch], ptName[pt], sp, ly));
2616     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
2617   }
2618   TObjArray *arr(NULL);
2619   fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
2620   arr->SetName("hTRD2MC");
2621
2622   TH2 *h2(NULL); Int_t jh(0);
2623   for(; ih--; ){
2624     if(!hp[ih].H()) continue;
2625     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2626     arr->AddAt(h2, jh++);
2627   }
2628
2629   // combine up the tree of projections
2630   AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
2631   //Int_t iproj(0), jproj(0);
2632   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2633     for(Int_t ich(0); ich<kNcharge; ich++){
2634       for(Int_t ipt(0); ipt<kNpt; ipt++){
2635         /*!dy*/
2636         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2637           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2638             if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2639             (*pr0)+=(*pr1);
2640           }
2641           AliDebug(2, Form("Rename %s to HMCTrkY%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], ily));
2642           pr0->H()->SetNameTitle(Form("HMCTrkY%c%c%d", chName[ich], ptName[ipt], ily),
2643                                     Form("Tracks[%c]:: #Deltay{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2644           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2645           if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2646         }
2647         /*!dphi*/
2648         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2649           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2650             if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2651             (*pr0)+=(*pr1);
2652           }
2653           AliDebug(2, Form("Rename %s to HMCTrkPh%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], ily));
2654           pr0->H()->SetNameTitle(Form("HMCTrkPh%c%c%d", chName[ich], ptName[ipt], ily),
2655                                     Form("Tracks[%c]:: #Delta#phi{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2656           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2657           if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2658         }
2659
2660         /*!dpt/pt*/
2661         if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2662           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2663             if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2664             (*pr0)+=(*pr1);
2665           }
2666           AliDebug(2, Form("Rename %s to HMCTrkDPt%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], ily));
2667           pr0->H()->SetNameTitle(Form("HMCTrkDPt%c%c%d", chName[ich], ptName[ipt], ily),
2668                                     Form("Tracks[%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2669           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2670           if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2671         }
2672       }
2673       /*!dy*/
2674       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2675         pr0->H()->SetNameTitle(Form("HMCTrkY%c%d", chName[ich], ily),
2676                               Form("Tracks[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
2677         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2678         if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2679       }
2680       /*!dphi*/
2681       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2682         pr0->H()->SetNameTitle(Form("HMCTrkPh%c%d", chName[ich], ily),
2683                               Form("Tracks[%c]:: #Delta#phi Ly[%d]", chSgn[ich], ily));
2684         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2685         if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2686       }
2687       /*!dpt/pt*/
2688       if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2689         pr0->H()->SetNameTitle(Form("HMCTrkDPt%c%d", chName[ich], ily),
2690                               Form("Tracks[%c]:: #Deltap_{t}/p_{t} Ly[%d]", chSgn[ich], ily));
2691         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2692         if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2693       }
2694     }
2695     /*!dy*/
2696     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2697       pr0->H()->SetNameTitle(Form("HMCTrkY%d", ily), Form("Tracks :: #Deltay Ly[%d]", ily));
2698       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2699     }
2700     /*!dphi*/
2701     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2702       pr0->H()->SetNameTitle(Form("HMCTrkPh%d", ily), Form("Tracks :: #Delta#phi Ly[%d]", ily));
2703       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2704     }
2705     /*!dpt/pt*/
2706     if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2707       pr0->H()->SetNameTitle(Form("HMCTrkDPt%d", ily), Form("Tracks :: #Deltap_{t}/p_{t} Ly[%d]", ily));
2708       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2709     }
2710   }
2711   AliInfo(Form("Done %3d 2D projections.", jh));
2712   return kTRUE;
2713 }
2714
2715 //________________________________________________________
2716 Bool_t AliTRDresolution::PostProcess()
2717 {
2718 // Fit, Project, Combine, Extract values from the containers filled during execution
2719
2720   if (!fContainer) {
2721     AliError("ERROR: list not available");
2722     return kFALSE;
2723   }
2724   if(!fProj){
2725     AliInfo("Building array of projections ...");
2726     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
2727   }
2728   // set pt/p segmentation. guess from data
2729   THnSparse *H(NULL);
2730   if(!(H = (THnSparse*)fContainer->FindObject("hTracklet2TRDin"))){
2731     AliError("Missing/Wrong data @ hTracklet2TRDin.");
2732     return kFALSE;
2733   }
2734   fNpt=H->GetAxis(kPt)->GetNbins()+1;
2735   if(!MakeMomSegmentation()) return kFALSE;
2736
2737   //PROCESS EXPERIMENTAL DISTRIBUTIONS
2738   // Clusters detector
2739   if(HasProcess(kDetector)) if(!MakeProjectionDetector()) return kFALSE;
2740   // Clusters residuals
2741   if(HasProcess(kCluster)) if(!MakeProjectionCluster()) return kFALSE;
2742   fNRefFigures = 3;
2743   // Tracklet residual/pulls
2744   if(HasProcess(kTracklet)) if(!MakeProjectionTracklet()) return kFALSE;
2745   fNRefFigures = 7;
2746   // TRDin residual/pulls
2747   if(HasProcess(kTrackIn)) if(!MakeProjectionTrackIn()) return kFALSE;
2748   fNRefFigures = 11;
2749
2750   if(!HasMCdata()) return kTRUE;
2751   //PROCESS MC RESIDUAL DISTRIBUTIONS
2752
2753   // CLUSTER Y RESOLUTION/PULLS
2754   if(HasProcess(kCluster)) if(!MakeProjectionCluster(kTRUE)) return kFALSE;
2755   fNRefFigures = 17;
2756
2757   // TRACKLET RESOLUTION/PULLS
2758   if(HasProcess(kTracklet)) if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
2759   fNRefFigures = 21;
2760
2761   // TRACK RESOLUTION/PULLS
2762   if(HasProcess(kTracklet)) if(!MakeProjectionTrack()) return kFALSE;
2763   fNRefFigures+=16;
2764
2765   // TRACK TRDin RESOLUTION/PULLS
2766   if(HasProcess(kTrackIn)) if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
2767   fNRefFigures+=8;
2768
2769   return kTRUE;
2770 }
2771
2772
2773 //________________________________________________________
2774 void AliTRDresolution::Terminate(Option_t *opt)
2775 {
2776   AliTRDrecoTask::Terminate(opt);
2777   if(HasPostProcess()) PostProcess();
2778 }
2779
2780 //________________________________________________________
2781 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2782 {
2783 // Helper function to avoid duplication of code
2784 // Make first guesses on the fit parameters
2785
2786   // find the intial parameters of the fit !! (thanks George)
2787   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2788   Double_t sum = 0.;
2789   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2790   f->SetParLimits(0, 0., 3.*sum);
2791   f->SetParameter(0, .9*sum);
2792   Double_t rms = h->GetRMS();
2793   f->SetParLimits(1, -rms, rms);
2794   f->SetParameter(1, h->GetMean());
2795
2796   f->SetParLimits(2, 0., 2.*rms);
2797   f->SetParameter(2, rms);
2798   if(f->GetNpar() <= 4) return;
2799
2800   f->SetParLimits(3, 0., sum);
2801   f->SetParameter(3, .1*sum);
2802
2803   f->SetParLimits(4, -.3, .3);
2804   f->SetParameter(4, 0.);
2805
2806   f->SetParLimits(5, 0., 1.e2);
2807   f->SetParameter(5, 2.e-1);
2808 }
2809
2810 //________________________________________________________
2811 TObjArray* AliTRDresolution::Histos()
2812 {
2813   //
2814   // Define histograms
2815   //
2816
2817   if(fContainer) return fContainer;
2818
2819   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
2820   THnSparse *H(NULL);
2821   const Int_t nhn(100); Char_t hn[nhn]; TString st;
2822   Int_t triggerBins(fTriggerList?(TMath::Power(2.,fTriggerList->GetEntriesFast())-1):0);
2823
2824   //++++++++++++++++++++++
2825   // cluster to detector
2826   snprintf(hn, nhn, "h%s", fgPerformanceName[kDetector]);
2827   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2828     const Int_t mdim(6);
2829     const Char_t *cldTitle[mdim] = {"layer", fgkTitle[kPhi], "pad row", "centrality", "q [a.u.]", triggerBins?"trigger":"no. of pads"/*, "tb [10 Hz]"*/};
2830     const Int_t cldNbins[mdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], 76, AliTRDeventInfo::kCentralityClasses, 50, triggerBins?triggerBins:kNpads};
2831     const Double_t cldMin[mdim]  = {-0.5, fgkMin[kPhi], -0.5, -0.5, 0., 0.5},
2832                    cldMax[mdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], 75.5, AliTRDeventInfo::kCentralityClasses - 0.5, 1200., triggerBins?(triggerBins+.5):(kNpads+.5)};
2833     st = "cluster proprieties;";
2834     // define minimum info to be saved in non debug mode
2835     Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimDet);
2836     for(Int_t idim(0); idim<ndim; idim++){ st += cldTitle[idim]; st+=";";}
2837     H = new THnSparseI(hn, st.Data(), ndim, cldNbins, cldMin, cldMax);
2838   } else H->Reset();
2839   fContainer->AddAt(H, kDetector);
2840   //++++++++++++++++++++++
2841   // cluster to tracklet residuals/pulls
2842   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
2843   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2844     const Int_t mdim(10);
2845     const Char_t *clTitle[mdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]", "centrality", triggerBins?"trigger":"no. of pads"};
2846     const Int_t clNbins[mdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 30, 30, 15, AliTRDeventInfo::kCentralityClasses, triggerBins?triggerBins:kNpads};
2847     const Double_t clMin[mdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45, -0.5, 0.5},
2848                    clMax[mdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45, AliTRDeventInfo::kCentralityClasses - 0.5, triggerBins?(triggerBins+.5):(kNpads+.5)};
2849     st = "cluster spatial&charge resolution;";
2850     // define minimum info to be saved in non debug mode
2851     Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimCl);
2852     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2853     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2854   } else H->Reset();
2855   fContainer->AddAt(H, kCluster);
2856   //++++++++++++++++++++++
2857   // tracklet to TRD track
2858   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
2859   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2860     const Int_t mdim(kNdim+8);
2861     Char_t *trTitle[mdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2862     Int_t trNbins[mdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2863     Double_t trMin[mdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2864     Double_t trMax[mdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2865     // set specific fields
2866     trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
2867     trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
2868     trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
2869     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2870     Int_t jdim(kNdim);
2871     trTitle[jdim]=StrDup("centrality"); trNbins[jdim] = AliTRDeventInfo::kCentralityClasses; trMin[jdim] = -.5; trMax[jdim] = AliTRDeventInfo::kCentralityClasses - 0.5;
2872     jdim++; 
2873     if(triggerBins){
2874       trTitle[jdim]=StrDup("trigger"); trNbins[jdim] = triggerBins;
2875       trMin[jdim] = 0.5; trMax[jdim] = triggerBins+.5;
2876     } else {
2877       trTitle[jdim]=StrDup("occupancy [%]"); trNbins[jdim] = 12;
2878       trMin[jdim] = 25.; trMax[jdim] = 85.;
2879     }
2880
2881     st = "tracklet spatial&charge resolution;";
2882     // define minimum info to be saved in non debug mode
2883     Int_t ndim=DebugLevel()>=1?(jdim+1):kNdimTrklt;
2884     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2885     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2886   } else H->Reset();
2887   fContainer->AddAt(H, kTracklet);
2888   //++++++++++++++++++++++
2889   // tracklet to TRDin
2890   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
2891   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2892     // set specific fields
2893     const Int_t mdim(kNdim+2);
2894     Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*));
2895     Int_t trinNbins[mdim];   memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t));
2896     Double_t trinMin[mdim];  memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t));
2897     Double_t trinMax[mdim];  memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t));
2898     if(triggerBins){
2899       trinTitle[kBC]=StrDup("trigger"); trinNbins[kBC] = triggerBins;
2900       trinMin[kBC] = 0.5; trinMax[kBC] = triggerBins+.5;
2901     }
2902 //    trinNbins[kSpeciesChgRC] = Int_t(kNcharge)*(kNspc-1)+1; trinMin[kSpeciesChgRC] = -kNspc+0.5; trinMax[kSpeciesChgRC] = kNspc-0.5;
2903     trinNbins[kPt]=fNpt-1; trinMax[kPt] = trinNbins[kPt]-.5;
2904     trinTitle[kNdim]=StrDup("bin_p"); trinNbins[kNdim] = fNpt-1; trinMin[kNdim] = -0.5; trinMax[kNdim] = trinNbins[kNdim]-.5;
2905     trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4;
2906     //trinTitle[kNdim+2]=StrDup("Fill Bunch"); trinNbins[kNdim+2]=3500; trinMin[kNdim+2]=-0.5; trinMax[kNdim+2]=3499.5;
2907     st = "r-#phi/z/angular residuals @ TRD entry;";
2908     // define minimum info to be saved in non debug mode
2909     Int_t ndim=(DebugLevel()>=1?mdim:(kNdim+1));//kNdimTrkIn;
2910     for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";}
2911     H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax);
2912   } else H->Reset();
2913   fContainer->AddAt(H, kTrackIn);
2914   // tracklet to TRDout
2915 //  fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2916
2917
2918   // Resolution histos
2919   if(!HasMCdata()) return fContainer;
2920
2921   //++++++++++++++++++++++
2922   // cluster to TrackRef residuals/pulls
2923   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
2924   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2925     H = (THnSparseI*)((THnSparseI*)fContainer->At(kCluster))->Clone(hn);
2926     H->SetTitle("MC cluster spatial resolution");
2927   } else H->Reset();
2928   fContainer->AddAt(H, kMCcluster);
2929   //++++++++++++++++++++++
2930   // tracklet to TrackRef
2931   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
2932   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2933     H = (THnSparseI*)((THnSparseI*)fContainer->At(kTracklet))->Clone(hn);
2934     H->SetTitle("MC tracklet spatial resolution");
2935   } else H->Reset();
2936   fContainer->AddAt(H, kMCtracklet);
2937   //++++++++++++++++++++++
2938   // TRDin to TrackRef
2939   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
2940   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2941     H = (THnSparseI*)((THnSparseI*)fContainer->At(kTrackIn))->Clone(hn);
2942     H->SetTitle("MC r-#phi/z/angular residuals @ TRD entry");
2943   } else H->Reset();
2944   fContainer->AddAt(H, kMCtrackIn);
2945   //++++++++++++++++++++++
2946   // track to TrackRef
2947   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
2948   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2949     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2950     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2951     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2952     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2953     // set specific fields
2954     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2955     trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
2956     trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
2957     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2958     trNbins[kPt] =fNpt-1;trMax[kPt]   = trNbins[kPt]-.5;
2959     //trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2960     trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5;
2961
2962     st = "MC track spatial&p_{t} resolution;";
2963     // define minimum info to be saved in non debug mode
2964     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
2965     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2966     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2967   } else H->Reset();
2968   fContainer->AddAt(H, kMCtrack);
2969
2970   return fContainer;
2971 }
2972
2973 //____________________________________________________________________
2974 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2975 {
2976 //
2977 // Fit track with a staight line using the "np" clusters stored in the array "points".
2978 // The following particularities are stored in the clusters from points:
2979 //   1. pad tilt as cluster charge
2980 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
2981 // The parameters of the straight line fit are stored in the array "param" in the following order :
2982 //     param[0] - x0 reference radial position
2983 //     param[1] - y0 reference r-phi position @ x0
2984 //     param[2] - z0 reference z position @ x0
2985 //     param[3] - slope dy/dx
2986 //     param[4] - slope dz/dx
2987 //
2988 // Attention :
2989 // Function should be used to refit tracks for B=0T
2990 //
2991
2992   if(np<40){
2993     if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2994     return kFALSE;
2995   }
2996   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2997
2998   Double_t x0(0.);
2999   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
3000   x0/=Float_t(np);
3001
3002   Double_t x, y, z, dx, tilt(0.);
3003   for(Int_t ip(0); ip<np; ip++){
3004     x = points[ip].GetX(); z = points[ip].GetZ();
3005     dx = x - x0;
3006     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
3007   }
3008   if(zfitter.Eval() != 0) return kFALSE;
3009
3010   Double_t z0    = zfitter.GetParameter(0);
3011   Double_t dzdx  = zfitter.GetParameter(1);
3012   for(Int_t ip(0); ip<np; ip++){
3013     if(points[ip].GetClusterType()) continue;
3014     x    = points[ip].GetX();
3015     dx   = x - x0;
3016     y    = points[ip].GetY();
3017     z    = points[ip].GetZ();
3018     tilt = points[ip].GetCharge();
3019     y -= tilt*(-dzdx*dx + z - z0);
3020     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
3021     yfitter.AddPoint(&dx, y, 1.);
3022   }
3023   if(yfitter.Eval() != 0) return kFALSE;
3024   Double_t y0   = yfitter.GetParameter(0);
3025   Double_t dydx = yfitter.GetParameter(1);
3026
3027   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
3028   if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
3029   return kTRUE;
3030 }
3031
3032 //____________________________________________________________________
3033 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
3034 {
3035 //
3036 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
3037 // See function FitTrack for the data stored in the "clusters" array
3038
3039 // The parameters of the straight line fit are stored in the array "param" in the following order :
3040 //     par[0] - x0 reference radial position
3041 //     par[1] - y0 reference r-phi position @ x0
3042 //     par[2] - slope dy/dx
3043 //
3044 // Attention :
3045 // Function should be used to refit tracks for B=0T
3046 //
3047
3048   TLinearFitter yfitter(2, "pol1");
3049
3050   // grep data for tracklet
3051   Double_t x0(0.), x[60], y[60], dy[60];
3052   Int_t nly(0);
3053   for(Int_t ip(0); ip<np; ip++){
3054     if(points[ip].GetClusterType()) continue;
3055     if(points[ip].GetVolumeID() != ly) continue;
3056     Float_t xt(points[ip].GetX())
3057            ,yt(param[1] + param[3] * (xt - param[0]));
3058     x[nly] = xt;
3059     y[nly] = points[ip].GetY();
3060     dy[nly]= y[nly]-yt;
3061     x0    += xt;
3062     nly++;
3063   }
3064   if(nly<10){
3065     if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
3066     return kFALSE;
3067   }
3068   // set radial reference for fit
3069   x0 /= Float_t(nly);
3070
3071   // find tracklet core
3072   Double_t mean(0.), sig(1.e3);
3073   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
3074
3075   // simple cluster error parameterization
3076   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
3077
3078   // fit tracklet core
3079   for(Int_t jly(0); jly<nly; jly++){
3080     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
3081     Double_t dx(x[jly]-x0);
3082     yfitter.AddPoint(&dx, y[jly], 1.);
3083   }
3084   if(yfitter.Eval() != 0) return kFALSE;
3085   par[0] = x0;
3086   par[1] = yfitter.GetParameter(0);
3087   par[2] = yfitter.GetParameter(1);
3088   return kTRUE;
3089 }
3090
3091 //____________________________________________________________________
3092 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
3093 {
3094 //
3095 // Global selection mechanism of tracksbased on cluster to fit residuals
3096 // The parameters are the same as used ni function FitTrack().
3097
3098   const Float_t kS(0.6), kM(0.2);
3099   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
3100   Float_t dy, dz, s, m;
3101   for(Int_t ip(0); ip<np; ip++){
3102     if(points[ip].GetClusterType()) continue;
3103     Float_t x0(points[ip].GetX())
3104           ,y0(param[1] + param[3] * (x0 - param[0]))
3105           ,z0(param[2] + param[4] * (x0 - param[0]));
3106     dy=points[ip].GetY() - y0; h.Fill(dy);
3107     dz=points[ip].GetZ() - z0;
3108   }
3109   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
3110   fg.SetParameter(1, 0.);
3111   fg.SetParameter(2, 2.e-2);
3112   h.Fit(&fg, "QN");
3113   m=fg.GetParameter(1); s=fg.GetParameter(2);
3114   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
3115   return kTRUE;
3116 }
3117
3118 //________________________________________________________
3119 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
3120 {
3121   //
3122   // Get the most probable value and the full width half mean
3123   // of a Landau distribution
3124   //
3125
3126   const Float_t dx = 1.;
3127   mpv = f->GetParameter(1);
3128   Float_t fx, max = f->Eval(mpv);
3129
3130   xm = mpv - dx;
3131   while((fx = f->Eval(xm))>.5*max){
3132     if(fx>max){
3133       max = fx;
3134       mpv = xm;
3135     }
3136     xm -= dx;
3137   }
3138
3139   xM += 2*(mpv - xm);
3140   while((fx = f->Eval(xM))>.5*max) xM += dx;
3141 }
3142
3143
3144 // #include "TFile.h"
3145 // //________________________________________________________
3146 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
3147 // {
3148 //   if(!file){
3149 //     AliWarning("Use cluster position as in reconstruction.");
3150 //     SetLoadCorrection();
3151 //     return kTRUE;
3152 //   }
3153 //   TDirectory *cwd(gDirectory);
3154 //   TString fileList;
3155 //   FILE *filePtr = fopen(file, "rt");
3156 //   if(!filePtr){
3157 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
3158 //     SetLoadCorrection();
3159 //     return kFALSE;
3160 //   }
3161 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
3162 //   while(fileList.Gets(filePtr)){
3163 //     if(!TFile::Open(fileList.Data())) {
3164 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
3165 //       continue;
3166 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
3167 //
3168 //     TTree *tSys = (TTree*)gFile->Get("tSys");
3169 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
3170 //     tSys->Draw("det:t>>h2", "dx", "goff");
3171 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3172 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
3173 //     }
3174 //     h2->SetDirectory(cwd);
3175 //     gFile->Close();
3176 //   }
3177 //   cwd->cd();
3178 //
3179 //   if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>=2){
3180 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3181 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
3182 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3183 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
3184 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
3185 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3186 //       printf("%03d|", idet);
3187 //       fprintf(fout, "    {");
3188 //       for(Int_t it(0); it<30; it++){
3189 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
3190 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
3191 //       }
3192 //       printf("\n");
3193 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
3194 //     }
3195 //     fprintf(fout, "  };\n");
3196 //   }
3197 //   SetLoadCorrection();
3198 //   return kTRUE;
3199 // }
3200
3201 //________________________________________________________
3202 void AliTRDresolution::SetProcesses(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3203 {
3204 // steer processes
3205   if(det) SETBIT(fSteer, kDetector); else CLRBIT(fSteer, kDetector);
3206   if(cl)  SETBIT(fSteer, kCluster); else CLRBIT(fSteer, kCluster);
3207   if(trklt) SETBIT(fSteer, kTracklet); else CLRBIT(fSteer, kTracklet);
3208   if(trkin) SETBIT(fSteer, kTrackIn); else CLRBIT(fSteer, kTrackIn);
3209 }
3210
3211
3212 //________________________________________________________
3213 void AliTRDresolution::SetDump3D(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3214 {
3215 // steer processes
3216   if(det) SETBIT(fSteer, 4+kDetector); else CLRBIT(fSteer, 4+kDetector);
3217   if(cl)  SETBIT(fSteer, 4+kCluster); else CLRBIT(fSteer, 4+kCluster);
3218   if(trklt) SETBIT(fSteer, 4+kTracklet); else CLRBIT(fSteer, 4+kTracklet);
3219   if(trkin) SETBIT(fSteer, 4+kTrackIn); else CLRBIT(fSteer, 4+kTrackIn);
3220 }
3221