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