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