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