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