fix coverity
[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)
95
3ceb45ae 96Int_t const AliTRDresolution::fgkNbins[kNdim] = {
97 Int_t(kNbunchCross)/*bc*/,
98 180/*phi*/,
99 50/*eta*/,
100 Int_t(kNcharge)*AliPID::kSPECIES+1/*chg*species*/,
3ceb45ae 101 50/*dy*/,
102 50/*dz*/,
103 40/*dphi*/
0d69a485 104// kNpt/*pt*/,
3ceb45ae 105}; //! no of bins/projection
106Double_t const AliTRDresolution::fgkMin[kNdim] = {
107 -0.5,
108 -TMath::Pi(),
109 -1.,
110 -AliPID::kSPECIES-0.5,
3ceb45ae 111 -1.5,
112 -2.5,
113 -10.
0d69a485 114// -0.5,
3ceb45ae 115}; //! low limits for projections
116Double_t const AliTRDresolution::fgkMax[kNdim] = {
117 Int_t(kNbunchCross)-0.5,
118 TMath::Pi(),
119 1.,
120 AliPID::kSPECIES+0.5,
3ceb45ae 121 1.5,
122 2.5,
123 10.
0d69a485 124// kNpt-0.5,
3ceb45ae 125}; //! high limits for projections
126Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
127 "bunch cross",
128 "#phi [rad]",
129 "#eta",
130 "chg*spec*rc",
3ceb45ae 131 "#Deltay [cm]",
132 "#Deltaz [cm]",
133 "#Delta#phi [deg]"
0d69a485 134// "bin_p_{t}",
3ceb45ae 135}; //! title of projection
136
137UChar_t const AliTRDresolution::fgNproj[kNclasses] = {
138 6, 5, 5, 5,
92d6d80c 139 2, 5, 11, 11, 11
1ee39b3a 140};
3ceb45ae 141Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
142 "Cluster2Track"
1ee39b3a 143 ,"Tracklet2Track"
a310e49b 144 ,"Tracklet2TRDin"
145 ,"Tracklet2TRDout"
1ee39b3a 146 ,"Cluster2MC"
147 ,"Tracklet2MC"
92d6d80c 148 ,"TRDin2MC"
149 ,"TRDout2MC"
1ee39b3a 150 ,"TRD2MC"
151};
3ceb45ae 152Float_t AliTRDresolution::fgPtBin[kNpt+1] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
1ee39b3a 153
154//________________________________________________________
155AliTRDresolution::AliTRDresolution()
f2e89a4c 156 :AliTRDrecoTask()
1ee39b3a 157 ,fIdxPlot(0)
92c40e64 158 ,fIdxFrame(0)
3ba3e0a4 159 ,fPtThreshold(1.)
9dcc64cc 160 ,fDyRange(0.75)
3ceb45ae 161 ,fProj(NULL)
92c40e64 162 ,fDBPDG(NULL)
b91fdd71 163 ,fCl(NULL)
b91fdd71 164 ,fMCcl(NULL)
f8f46e4d 165{
166 //
167 // Default constructor
168 //
f2e89a4c 169 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
3ceb45ae 170 MakePtSegmentation();
f8f46e4d 171}
172
705f8b0a 173//________________________________________________________
f8f46e4d 174AliTRDresolution::AliTRDresolution(char* name)
f2e89a4c 175 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
f8f46e4d 176 ,fIdxPlot(0)
177 ,fIdxFrame(0)
3ba3e0a4 178 ,fPtThreshold(1.)
9dcc64cc 179 ,fDyRange(0.75)
3ceb45ae 180 ,fProj(NULL)
f8f46e4d 181 ,fDBPDG(NULL)
f8f46e4d 182 ,fCl(NULL)
f8f46e4d 183 ,fMCcl(NULL)
1ee39b3a 184{
185 //
186 // Default constructor
187 //
188
1ee39b3a 189 InitFunctorList();
3ceb45ae 190 MakePtSegmentation();
1ee39b3a 191
705f8b0a 192 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
705f8b0a 193 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
1ee39b3a 194}
195
196//________________________________________________________
197AliTRDresolution::~AliTRDresolution()
198{
199 //
200 // Destructor
201 //
202
3ceb45ae 203 if(fProj){fProj->Delete(); delete fProj;}
1ee39b3a 204 if(fCl){fCl->Delete(); delete fCl;}
1ee39b3a 205 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
1ee39b3a 206}
207
208
209//________________________________________________________
f8f46e4d 210void AliTRDresolution::UserCreateOutputObjects()
1ee39b3a 211{
212 // spatial resolution
5935a6da 213
068e2c00 214 AliTRDrecoTask::UserCreateOutputObjects();
83b44483 215 InitExchangeContainers();
216}
1ee39b3a 217
83b44483 218//________________________________________________________
219void AliTRDresolution::InitExchangeContainers()
220{
61f6b45e 221// Init containers for subsequent tasks (AliTRDclusterResolution)
222
3ceb45ae 223 fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
224 fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
e3147647 225 PostData(kClToTrk, fCl);
226 PostData(kClToMC, fMCcl);
1ee39b3a 227}
228
229//________________________________________________________
f8f46e4d 230void AliTRDresolution::UserExec(Option_t *opt)
1ee39b3a 231{
232 //
233 // Execution part
234 //
235
236 fCl->Delete();
1ee39b3a 237 fMCcl->Delete();
b4414720 238 AliTRDrecoTask::UserExec(opt);
1ee39b3a 239}
240
553528eb 241//________________________________________________________
6475ec36 242Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
553528eb 243{
244// Helper function to calculate pulls in the yz plane
245// using proper tilt rotation
246// Uses functionality defined by AliTRDseedV1.
247
6475ec36 248 return kTRUE;
249/*
553528eb 250 Double_t t2(tilt*tilt);
9dcc64cc 251 // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
553528eb 252
253 // rotate along pad
254 Double_t cc[3];
255 cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
256 cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
257 cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
258 // do sqrt
259 Double_t sqr[3]={0., 0., 0.};
260 if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
261 Double_t invsqr[3]={0., 0., 0.};
262 if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
263 Double_t tmp(dyz[0]);
264 dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
265 dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
266 return kTRUE;
6475ec36 267*/
553528eb 268}
269
1ee39b3a 270//________________________________________________________
3ceb45ae 271TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
1ee39b3a 272{
273 //
3ceb45ae 274 // Plot the cluster distributions
1ee39b3a 275 //
276
277 if(track) fkTrack = track;
278 if(!fkTrack){
3d2a3dff 279 AliDebug(4, "No Track defined.");
b91fdd71 280 return NULL;
1ee39b3a 281 }
3ceb45ae 282 if(fkESD->GetTOFbc() > 1){
283 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
b91fdd71 284 return NULL;
1ee39b3a 285 }
3ceb45ae 286 if(fPt<fPtThreshold){
287 AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
b91fdd71 288 return NULL;
1ee39b3a 289 }
3ceb45ae 290 THnSparse *H(NULL);
291 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
1ee39b3a 292 AliWarning("No output container defined.");
b91fdd71 293 return NULL;
1ee39b3a 294 }
3ceb45ae 295
296 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
297 Double_t val[kNdim]; Float_t exb, vd, t0, s2, dl, dt, corr;
298 TObjArray *clInfoArr(NULL);
299 AliTRDseedV1 *fTracklet(NULL);
0d69a485 300 AliTRDcluster *c(NULL)/*, *cc(NULL)*/;
1ee39b3a 301 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
302 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
303 if(!fTracklet->IsOK()) continue;
3ceb45ae 304 fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
305 val[kBC] = ily;
306 val[kPhi] = fPhi;
307 val[kEta] = fEta;
0d69a485 308 //val[kPt] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb))*TMath::RadToDeg();
3ceb45ae 309 corr = fkTrack->Charge()/TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1))/vd;
310 fTracklet->ResetClusterIter(kTRUE);
311 while((c = fTracklet->NextCluster())){
0d69a485 312 Float_t xc(c->GetX()); //Int_t tb(c->GetLocalTimeBin());
3ceb45ae 313
314 val[kYrez] = c->GetY()-fTracklet->GetYat(xc);
0d69a485 315 //val[kZrez] = fTracklet->GetX0()-xc;
316 //val[kPrez] = 0.; Int_t ic(0);
317 //if((cc = fTracklet->GetClusters(tb-1))) {val[kPrez] += cc->GetQ(); ic++;}
318 //if((cc = fTracklet->GetClusters(tb-2))) {val[kPrez] += cc->GetQ(); ic++;}
319 //if(ic) val[kPrez] /= (ic*c->GetQ());
3ceb45ae 320 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:(c->GetQ()*corr);
321 H->Fill(val);
322/* // tilt rotation of covariance for clusters
553528eb 323 Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
a47a87c5 324 cov[0] = (sy2+t2*sz2)*corr;
325 cov[1] = tilt*(sz2 - sy2)*corr;
326 cov[2] = (t2*sy2+sz2)*corr;
553528eb 327 // sum with track covariance
a47a87c5 328 cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
553528eb 329 Double_t dyz[2]= {dy[1], dz[1]};
3ceb45ae 330 Pulls(dyz, cov, tilt);*/
1ee39b3a 331
dfd7d48b 332 // Get z-position with respect to anode wire
3ceb45ae 333 Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
334 zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
fd7ffd88 335 Int_t istk = geo->GetStack(c->GetDetector());
336 AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
dfd7d48b 337 Float_t row0 = pp->GetRow0();
94b94be0 338 Float_t d = row0 - zt + pp->GetAnodeWireOffset();
dfd7d48b 339 d -= ((Int_t)(2 * d)) / 2.0;
340 if (d > 0.25) d = 0.5 - d;
341
8ee59659 342 AliTRDclusterInfo *clInfo(NULL);
343 clInfo = new AliTRDclusterInfo;
dfd7d48b 344 clInfo->SetCluster(c);
3ceb45ae 345 //Float_t covF[] = {cov[0], cov[1], cov[2]};
346 clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
347 clInfo->SetResolution(val[kYrez]);
dfd7d48b 348 clInfo->SetAnisochronity(d);
3ceb45ae 349 clInfo->SetDriftLength(val[kZrez]);
350 clInfo->SetTilt(fTracklet->GetTilt());
83b44483 351 if(fCl) fCl->Add(clInfo);
352 else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
353
3ba3e0a4 354 if(DebugLevel()>=1){
d14a8ac2 355 if(!clInfoArr){
356 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
357 clInfoArr->SetOwner(kFALSE);
358 }
dfd7d48b 359 clInfoArr->Add(clInfo);
df0514f6 360 }
1ee39b3a 361 }
3ba3e0a4 362 if(DebugLevel()>=1 && clInfoArr){
3ceb45ae 363 ULong_t status = fkESD->GetStatus();
dfd7d48b 364 (*DebugStream()) << "cluster"
365 <<"status=" << status
366 <<"clInfo.=" << clInfoArr
367 << "\n";
d14a8ac2 368 clInfoArr->Clear();
dfd7d48b 369 }
1ee39b3a 370 }
d14a8ac2 371 if(clInfoArr) delete clInfoArr;
00a56108 372
3ceb45ae 373 return NULL;//H->Projection(kEta, kPhi);
1ee39b3a 374}
375
376
377//________________________________________________________
378TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
379{
380// Plot normalized residuals for tracklets to track.
381//
382// We start from the result that if X=N(|m|, |Cov|)
383// BEGIN_LATEX
384// (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
385// END_LATEX
386// in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
387// reference position.
388 if(track) fkTrack = track;
389 if(!fkTrack){
3d2a3dff 390 AliDebug(4, "No Track defined.");
b91fdd71 391 return NULL;
1ee39b3a 392 }
3ceb45ae 393 if(fkESD->GetTOFbc()>1){
394 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
395 return NULL;
396 }
397 THnSparse *H(NULL);
398 if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
1ee39b3a 399 AliWarning("No output container defined.");
b91fdd71 400 return NULL;
1ee39b3a 401 }
dbb2e0a7 402// return NULL;
3ceb45ae 403 Double_t val[kNdim+1];
404 AliTRDseedV1 *fTracklet(NULL);
3ba3e0a4 405 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
1ee39b3a 406 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
407 if(!fTracklet->IsOK()) continue;
3ceb45ae 408 val [kBC] = il;
409 val[kPhi] = fPhi;
410 val[kEta] = fEta;
411 val[kSpeciesChgRC]= fSpecies;
0d69a485 412 //val[kPt] = GetPtBin(fTracklet->GetPt());
3ceb45ae 413 Double_t dyt(fTracklet->GetYref(0) - fTracklet->GetYfit(0)),
414 dzt(fTracklet->GetZref(0) - fTracklet->GetZfit(0)),
415 dydx(fTracklet->GetYfit(1)),
416 tilt(fTracklet->GetTilt());
417 // correct for tilt rotation
418 val[kYrez] = dyt - dzt*tilt;
419 val[kZrez] = dzt + dyt*tilt;
420 dydx+= tilt*fTracklet->GetZref(1);
421 val[kPrez] = TMath::ATan((fTracklet->GetYref(1) - dydx)/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
422 if(fTracklet->IsRowCross()){
423 val[kSpeciesChgRC]= 0.;
424 val[kPrez] = fkTrack->Charge(); // may be better defined
425 } else {
426 Float_t exb, vd, t0, s2, dl, dt;
427 fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
428 val[kZrez] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb));
429 }
430 val[kNdim] = fTracklet->GetdQdl();
dbb2e0a7 431 if(DebugLevel()>=1) H->Fill(val);
3ceb45ae 432
433// // compute covariance matrix
434// fTracklet->GetCovAt(x, cov);
435// fTracklet->GetCovRef(covR);
436// cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
437// Double_t dyz[2]= {dy[1], dz[1]};
438// Pulls(dyz, cov, tilt);
439// ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
440// ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
dfd7d48b 441
dbb2e0a7 442 if(DebugLevel()>=2){
3ceb45ae 443 Bool_t rc(fTracklet->IsRowCross());
dfd7d48b 444 UChar_t err(fTracklet->GetErrorMsg());
3ceb45ae 445 Double_t x(fTracklet->GetX()),
446 pt(fTracklet->GetPt()),
447 yt(fTracklet->GetYref(0)),
448 zt(fTracklet->GetZref(0)),
449 phi(fTracklet->GetYref(1)),
450 tht(fTracklet->GetZref(1));
451 Int_t ncl(fTracklet->GetN()),
452 det(fTracklet->GetDetector());
dfd7d48b 453 (*DebugStream()) << "tracklet"
3ba3e0a4 454 <<"pt=" << pt
3ceb45ae 455 <<"x=" << x
6475ec36 456 <<"yt=" << yt
3ceb45ae 457 <<"zt=" << zt
3ba3e0a4 458 <<"phi=" << phi
459 <<"tht=" << tht
3ceb45ae 460 <<"det=" << det
6475ec36 461 <<"n=" << ncl
3ceb45ae 462 <<"dy0=" << dyt
463 <<"dz0=" << dzt
464 <<"dy=" << val[kYrez]
465 <<"dz=" << val[kZrez]
466 <<"dphi="<< val[kPrez]
467 <<"dQ ="<< val[kNdim]
dfd7d48b 468 <<"rc=" << rc
469 <<"err=" << err
470 << "\n";
471 }
1ee39b3a 472 }
3ceb45ae 473 return NULL;//H->Projection(kEta, kPhi);
1ee39b3a 474}
475
476
477//________________________________________________________
a310e49b 478TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
1ee39b3a 479{
480// Store resolution/pulls of Kalman before updating with the TRD information
481// at the radial position of the first tracklet. The following points are used
482// for comparison
483// - the (y,z,snp) of the first TRD tracklet
484// - the (y, z, snp, tgl, pt) of the MC track reference
485//
486// Additionally the momentum resolution/pulls are calculated for usage in the
487// PID calculation.
488
489 if(track) fkTrack = track;
3ceb45ae 490 if(!fkTrack){
3d2a3dff 491 AliDebug(4, "No Track defined.");
b91fdd71 492 return NULL;
1ee39b3a 493 }
3ceb45ae 494 // check container
495 THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
496 if(!H){
497 AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
83b44483 498 return NULL;
499 }
3ceb45ae 500 // check input track status
501 AliExternalTrackParam *tin(NULL);
a310e49b 502 if(!(tin = fkTrack->GetTrackIn())){
3ceb45ae 503 AliError("Track did not entered TRD fiducial volume.");
b91fdd71 504 return NULL;
1ee39b3a 505 }
3ceb45ae 506 // check first tracklet
507 AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
508 if(!fTracklet){
509 AliDebug(3, "No Tracklet in ly[0]. Skip track.");
b91fdd71 510 return NULL;
1ee39b3a 511 }
3ceb45ae 512 // check radial position
513 Double_t x = tin->GetX();
514 if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
515 AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
516 return NULL;
1ee39b3a 517 }
518
3ceb45ae 519 Int_t bc(fkESD->GetTOFbc()%2);
520 const Double_t *parR(tin->GetParameter());
521 Double_t dyt(parR[0] - fTracklet->GetYfit(0)), dzt(parR[1] - fTracklet->GetZfit(0)),
522 phit(fTracklet->GetYfit(1)),
523 tilt(fTracklet->GetTilt());
524
525 // correct for tilt rotation
526 Double_t dy = dyt - dzt*tilt,
527 dz = dzt + dyt*tilt;
528 phit += tilt*parR[3];
529 Double_t dphi = TMath::ASin(parR[2])-TMath::ATan(phit);
530
531 Double_t val[kNdim];
532 val[kBC] = bc;
533 val[kPhi] = fPhi;
534 val[kEta] = fEta;
535 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fSpecies;
0d69a485 536// val[kPt] = GetPtBin(fPt);
3ceb45ae 537 val[kYrez] = dy;
538 val[kZrez] = dz;
539 val[kPrez] = dphi*TMath::RadToDeg();
540 H->Fill(val);
541
542 if(!HasMCdata()) return NULL; // H->Projection(kEta, kPhi);
543 if(!(H = (THnSparseI*)fContainer->At(kMCtrackIn))) {
544 AliError(Form("Missing container @ %d", Int_t(kMCtrackIn)));
545 return NULL;
546 }
1ee39b3a 547
3ceb45ae 548 // get MC info
1ee39b3a 549 UChar_t s;
3ceb45ae 550 Float_t pt0, eta, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
551 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, s)) return NULL;
552 dyt = y0 - fTracklet->GetYfit(0);
553 dzt = z0 - fTracklet->GetZfit(0);
554 phit= fTracklet->GetYfit(1) + tilt*dzdx0;
555 Float_t phi = TMath::ATan2(y0, x0);
556 dy = dyt - dzt*tilt;
557 dz = dzt + dyt*tilt;
558 dphi= TMath::ASin(dydx0)-TMath::ATan(phit);
559
92c40e64 560 Int_t pdg = fkMC->GetPDG(),
561 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
562 sign(0);
563 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
564 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
565 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
566
83b44483 567
3ceb45ae 568 val[kBC] = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
569 val[kPhi] = phi;
570 val[kEta] = eta;
571 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:sign*(sIdx+1);
0d69a485 572// val[kPt] = GetPtBin(pt0);
3ceb45ae 573 val[kYrez] = dy;
574 val[kZrez] = dz;
575 val[kPrez] = dphi*TMath::RadToDeg();
576 H->Fill(val);
577
578 return NULL; //H->Projection(kEta, kPhi);
1ee39b3a 579}
580
a310e49b 581//________________________________________________________
582TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
583{
584// Store resolution/pulls of Kalman after last update with the TRD information
585// at the radial position of the first tracklet. The following points are used
586// for comparison
587// - the (y,z,snp) of the first TRD tracklet
588// - the (y, z, snp, tgl, pt) of the MC track reference
589//
590// Additionally the momentum resolution/pulls are calculated for usage in the
591// PID calculation.
592
593 if(track) fkTrack = track;
596f82fb 594 return NULL;
a310e49b 595}
596
1ee39b3a 597//________________________________________________________
598TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
599{
600 //
601 // Plot MC distributions
602 //
603
604 if(!HasMCdata()){
e9d62bd2 605 AliDebug(2, "No MC defined. Results will not be available.");
b91fdd71 606 return NULL;
1ee39b3a 607 }
608 if(track) fkTrack = track;
609 if(!fkTrack){
3d2a3dff 610 AliDebug(4, "No Track defined.");
b91fdd71 611 return NULL;
1ee39b3a 612 }
83b44483 613 if(!fContainer){
614 AliWarning("No output container defined.");
615 return NULL;
616 }
92c40e64 617 // retriev track characteristics
618 Int_t pdg = fkMC->GetPDG(),
619 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
620 sign(0),
3ba3e0a4 621 sgm[3],
3ceb45ae 622 label(fkMC->GetLabel()),
623 fSegmentLevel(0);
92c40e64 624 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
625 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
626 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
92c40e64 627
628 TObjArray *arr(NULL);TH1 *h(NULL);
3ceb45ae 629 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
630 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
1ee39b3a 631 UChar_t s;
1ee39b3a 632 Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
633 Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
634 Double_t covR[7]/*, cov[3]*/;
3ceb45ae 635
3ba3e0a4 636 if(DebugLevel()>=3){
3ceb45ae 637 // get first detector
638 Int_t det = -1;
639 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
640 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
641 det = fTracklet->GetDetector();
642 break;
643 }
644 if(det>=0){
645 TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
646 Double_t m(-1.);
647 m = fkTrack->GetMass();
648 if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
649 (*DebugStream()) << "MCkalman"
650 << "pdg=" << pdg
651 << "det=" << det
652 << "x=" << &X
653 << "y=" << &Y
654 << "z=" << &Z
655 << "dx=" << &dX
656 << "dy=" << &dY
657 << "dz=" << &dZ
658 << "pt=" << &vPt
659 << "dpt=" << &dPt
660 << "bgt=" << &budget
661 << "cov=" << &cCOV
662 << "\n";
663 }
664 }
1ee39b3a 665 }
1ee39b3a 666 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
667 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
668 !fTracklet->IsOK())*/ continue;
669
3ba3e0a4 670 sgm[2] = fTracklet->GetDetector();
671 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
672 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
6558fd69 673 Double_t tilt(fTracklet->GetTilt())
674 ,t2(tilt*tilt)
675 ,corr(1./(1. + t2))
676 ,cost(TMath::Sqrt(corr));
1ee39b3a 677 x0 = fTracklet->GetX0();
678 //radial shift with respect to the MC reference (radial position of the pad plane)
679 x= fTracklet->GetX();
3ceb45ae 680 Bool_t rc(fTracklet->IsRowCross()); Float_t eta;
681 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, s)) continue;
1ee39b3a 682 xAnode = fTracklet->GetX0();
683
684 // MC track position at reference radial position
685 dx = x0 - x;
3ba3e0a4 686 if(DebugLevel()>=4){
1ee39b3a 687 (*DebugStream()) << "MC"
3ba3e0a4 688 << "det=" << sgm[2]
1ee39b3a 689 << "pdg=" << pdg
92c40e64 690 << "sgn=" << sign
1ee39b3a 691 << "pt=" << pt0
692 << "x=" << x0
693 << "y=" << y0
694 << "z=" << z0
695 << "dydx=" << dydx0
696 << "dzdx=" << dzdx0
697 << "\n";
698 }
699 Float_t ymc = y0 - dx*dydx0;
700 Float_t zmc = z0 - dx*dzdx0;
701 //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
702
703 // Kalman position at reference radial position
704 dx = xAnode - x;
705 dydx = fTracklet->GetYref(1);
706 dzdx = fTracklet->GetZref(1);
707 dzdl = fTracklet->GetTgl();
708 y = fTracklet->GetYref(0) - dx*dydx;
709 dy = y - ymc;
710 z = fTracklet->GetZref(0) - dx*dzdx;
711 dz = z - zmc;
712 pt = TMath::Abs(fTracklet->GetPt());
713 fTracklet->GetCovRef(covR);
714
92d6d80c 715 arr = (TObjArray*)((TObjArray*)fContainer->At(kMCtrack))->At(ily);
1ee39b3a 716 // y resolution/pulls
2589cf64 717 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
718 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(covR[0]), dz/TMath::Sqrt(covR[2]));
1ee39b3a 719 // z resolution/pulls
3ceb45ae 720 ((TH2S*)arr->At(2))->Fill(dzdx0, dz);
81979445 721 ((TH3S*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]), 0);
1ee39b3a 722 // phi resolution/ snp pulls
723 Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
724 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
725 Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
726 ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
727 // theta resolution/ tgl pulls
728 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
729 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
730 ((TH2I*)arr->At(6))->Fill(dzdl0,
731 TMath::ATan(dtgl));
732 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
733 // pt resolution \\ 1/pt pulls \\ p resolution for PID
92c40e64 734 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
735 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
afca20ef 736 ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
737 ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
738 ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);
1ee39b3a 739
740 // Fill Debug stream for Kalman track
3ba3e0a4 741 if(DebugLevel()>=4){
1ee39b3a 742 (*DebugStream()) << "MCtrack"
743 << "pt=" << pt
744 << "x=" << x
745 << "y=" << y
746 << "z=" << z
747 << "dydx=" << dydx
748 << "dzdx=" << dzdx
749 << "s2y=" << covR[0]
750 << "s2z=" << covR[2]
751 << "\n";
752 }
753
754 // recalculate tracklet based on the MC info
755 AliTRDseedV1 tt(*fTracklet);
756 tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
757 tt.SetZref(1, dzdx0);
fd7ffd88 758 tt.SetReconstructor(AliTRDinfoGen::Reconstructor());
26e4e234 759 tt.Fit(1);
1ee39b3a 760 x= tt.GetX();y= tt.GetY();z= tt.GetZ();
761 dydx = tt.GetYfit(1);
762 dx = x0 - x;
763 ymc = y0 - dx*dydx0;
764 zmc = z0 - dx*dzdx0;
6558fd69 765 dy = y-ymc;
766 dz = z-zmc;
81979445 767 Float_t dphi = (dydx - dydx0);
768 dphi /= (1.- dydx*dydx0);
6558fd69 769
1ee39b3a 770 // add tracklet residuals for y and dydx
771 arr = (TObjArray*)fContainer->At(kMCtracklet);
81979445 772
2589cf64 773 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
774 if(tt.GetS2Y()>0. && tt.GetS2Z()>0.) ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(tt.GetS2Y()), dz/TMath::Sqrt(tt.GetS2Z()));
81979445 775 ((TH3S*)arr->At(2))->Fill(dzdl0, dz, rc);
776 if(tt.GetS2Z()>0.) ((TH3S*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()), rc);
777 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
1ee39b3a 778
779 // Fill Debug stream for tracklet
3ba3e0a4 780 if(DebugLevel()>=4){
1ee39b3a 781 Float_t s2y = tt.GetS2Y();
782 Float_t s2z = tt.GetS2Z();
783 (*DebugStream()) << "MCtracklet"
784 << "rc=" << rc
785 << "x=" << x
786 << "y=" << y
787 << "z=" << z
788 << "dydx=" << dydx
789 << "s2y=" << s2y
790 << "s2z=" << s2z
791 << "\n";
792 }
793
fd7ffd88 794 AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
1ee39b3a 795 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
1ee39b3a 796 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
797
798 arr = (TObjArray*)fContainer->At(kMCcluster);
b91fdd71 799 AliTRDcluster *c = NULL;
83b44483 800 tt.ResetClusterIter(kFALSE);
801 while((c = tt.PrevCluster())){
1ee39b3a 802 Float_t q = TMath::Abs(c->GetQ());
3ceb45ae 803 x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()]; y = c->GetY();z = c->GetZ();
1ee39b3a 804 dx = x0 - x;
805 ymc= y0 - dx*dydx0;
806 zmc= z0 - dx*dzdx0;
6558fd69 807 dy = cost*(y - ymc - tilt*(z-zmc));
808 dz = cost*(z - zmc + tilt*(y-ymc));
809
1ee39b3a 810 // Fill Histograms
83b44483 811 if(q>20. && q<250. && pt0>fPtThreshold && c->IsInChamber()){
2589cf64 812 ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
813 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(c->GetSigmaY2()), dz/TMath::Sqrt(c->GetSigmaZ2()));
1ee39b3a 814 }
815
816 // Fill calibration container
817 Float_t d = zr0 - zmc;
818 d -= ((Int_t)(2 * d)) / 2.0;
819 if (d > 0.25) d = 0.5 - d;
820 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
1ee39b3a 821 clInfo->SetCluster(c);
822 clInfo->SetMC(pdg, label);
823 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
824 clInfo->SetResolution(dy);
825 clInfo->SetAnisochronity(d);
2589cf64 826 clInfo->SetDriftLength(dx);
1ee39b3a 827 clInfo->SetTilt(tilt);
83b44483 828 if(fMCcl) fMCcl->Add(clInfo);
829 else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
3ba3e0a4 830 if(DebugLevel()>=5){
d14a8ac2 831 if(!clInfoArr){
832 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
833 clInfoArr->SetOwner(kFALSE);
834 }
d43e2ad1 835 clInfoArr->Add(clInfo);
1ee39b3a 836 }
837 }
d43e2ad1 838 // Fill Debug Tree
3ba3e0a4 839 if(DebugLevel()>=5 && clInfoArr){
d43e2ad1 840 (*DebugStream()) << "MCcluster"
841 <<"clInfo.=" << clInfoArr
842 << "\n";
92c40e64 843 clInfoArr->Clear();
d43e2ad1 844 }
1ee39b3a 845 }
92c40e64 846 if(clInfoArr) delete clInfoArr;
1ee39b3a 847 return h;
848}
849
850
3ceb45ae 851//__________________________________________________________________________
852Int_t AliTRDresolution::GetPtBin(Float_t pt)
853{
854// Find pt bin according to local pt segmentation
855 Int_t ipt(0);
856 while(ipt<AliTRDresolution::kNpt){
857 if(pt<fgPtBin[ipt]) break;
858 ipt++;
859 }
860 return TMath::Max(0,ipt);
861}
862
863//________________________________________________________
864Float_t AliTRDresolution::GetMeanWithBoundary(TH1 *h, Float_t zm, Float_t zM, Float_t dz) const
865{
866// return mean of histogram "h"
867// if histo is empty returns -infinity
868// if few entries returns zM+epsilon
869// if mean less than zm returns zm-epsilon
870
0d69a485 871 Int_t ne(Int_t(h->GetEntries()));
3ceb45ae 872 if(ne==0) return -1.e+5;
873 else if(ne<20) return zM+0.5*dz;
874 else{
875 Float_t val(h->GetMean());
876 if(val<zm) return zm-0.5*dz;
877 else if(val>zM) return zM-0.5*dz;
878 else return val;
879 }
880}
881
1ee39b3a 882//________________________________________________________
883Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
884{
885 //
886 // Get the reference figures
887 //
888
1ee39b3a 889 if(!gPad){
890 AliWarning("Please provide a canvas to draw results.");
891 return kFALSE;
892 }
3ceb45ae 893/* Int_t selection[100], n(0), selStart(0); //
a310e49b 894 Int_t ly0(0), dly(5);
3ceb45ae 895 TList *l(NULL); TVirtualPad *pad(NULL); */
1ee39b3a 896 switch(ifig){
3ceb45ae 897 case 0:
1ee39b3a 898 break;
1ee39b3a 899 }
900 AliWarning(Form("Reference plot [%d] missing result", ifig));
901 return kFALSE;
902}
903
3ceb45ae 904
905//________________________________________________________
906void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
907{
908// Build pt segments
909 for(Int_t j(0); j<=kNpt; j++){
910 pt0+=(TMath::Exp(j*j*dpt)-1.);
911 fgPtBin[j]=pt0;
912 }
913}
914
00a3f7a4 915//________________________________________________________
916void AliTRDresolution::MakeSummary()
917{
918// Build summary plots
919
3ceb45ae 920 if(!fProj){
00a3f7a4 921 AliError("Missing results");
922 return;
923 }
3ceb45ae 924 Int_t iSumPlot(0);
925 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
926 TObjArray *arr(NULL);
068e2c00 927
3ceb45ae 928 // cluster resolution
929 // define palette
930 gStyle->SetPalette(1);
931 cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot++), "Cluster Resolution", 1024, 768);
932 cOut->Divide(3,2, 2.e-3, 2.e-3);
933 arr = (TObjArray*)fProj->At(kCluster);
934 for(Int_t iplot(0); iplot<fgNproj[kCluster]; iplot++){
935 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
936 ((TH2*)arr->At(iplot))->Draw("colz");
937 }
068e2c00 938 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
00a3f7a4 939
068e2c00 940
3ceb45ae 941 // trackIn systematic
942 // define palette
943 Int_t palette[50];
944 for (int i=1;i<49;i++) palette[i] = 50+i; palette[0]=kMagenta; palette[49]=kBlack;
945 gStyle->SetPalette(50, palette);
946 cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot++), "Track IN Resolution", 1024, 768);
947 cOut->Divide(3,2, 2.e-3, 2.e-3);
948 arr = (TObjArray*)fProj->At(kTrackIn);
949 for(Int_t iplot(0); iplot<fgNproj[kTrackIn]; iplot++){
950 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
951 ((TH2*)arr->At(iplot))->Draw("colz");
952 }
068e2c00 953 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
3ceb45ae 954
068e2c00 955 delete cOut;
3ceb45ae 956 gStyle->SetPalette(1);
068e2c00 957}
958
959//________________________________________________________
960void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
961{
962// Returns the range of the bulk of data in histogram h2. Removes outliers.
963// The "range" vector should be initialized with 2 elements
964// Option "mod" can be any of
965// - 0 : gaussian like distribution
966// - 1 : tailed distribution
967
968 Int_t nx(h2->GetNbinsX())
969 , ny(h2->GetNbinsY())
970 , n(nx*ny);
971 Double_t *data=new Double_t[n];
972 for(Int_t ix(1), in(0); ix<=nx; ix++){
973 for(Int_t iy(1); iy<=ny; iy++)
974 data[in++] = h2->GetBinContent(ix, iy);
00a3f7a4 975 }
068e2c00 976 Double_t mean, sigm;
977 AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
978
979 range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
980 if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
981 AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
982 TH1S h1("h1SF0", "", 100, range[0], range[1]);
983 h1.FillN(n,data,0);
984 delete [] data;
985
986 switch(mod){
987 case 0:// gaussian distribution
988 {
989 TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
990 h1.Fit(&fg, "QN");
991 mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
992 range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
993 AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1]));
994 break;
00a3f7a4 995 }
068e2c00 996 case 1:// tailed distribution
997 {
998 Int_t bmax(h1.GetMaximumBin());
999 Int_t jBinMin(1), jBinMax(100);
1000 for(Int_t ibin(bmax); ibin--;){
61f6b45e 1001 if(h1.GetBinContent(ibin)<1.){
068e2c00 1002 jBinMin=ibin; break;
1003 }
1004 }
1005 for(Int_t ibin(bmax); ibin++;){
61f6b45e 1006 if(h1.GetBinContent(ibin)<1.){
068e2c00 1007 jBinMax=ibin; break;
1008 }
00a3f7a4 1009 }
068e2c00 1010 range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1011 AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1]));
1012 break;
1013 }
00a3f7a4 1014 }
00a3f7a4 1015
068e2c00 1016 return;
1017}
1018
3ceb45ae 1019
068e2c00 1020//________________________________________________________
3ceb45ae 1021Bool_t AliTRDresolution::MakeProjectionCluster()
068e2c00 1022{
3ceb45ae 1023// Analyse cluster
1024 Int_t cidx = kCluster;
1025 if(fProj && fProj->At(cidx)) return kTRUE;
1026 if(!fContainer){
1027 AliError("Missing data container.");
1028 return kFALSE;
1029 }
1030 THnSparse *H(NULL);
1031 if(!(H = (THnSparse*)fContainer->At(cidx))){
1032 AliError(Form("Missing/Wrong data @ %d.", cidx));
1033 return kFALSE;
1034 }
1035
1036 TAxis *aphi(H->GetAxis(kPhi)),
1037 *aeta(H->GetAxis(kEta)),
1038 *as(H->GetAxis(kSpeciesChgRC)),
1039 //*apt(H->GetAxis(kPt)),
1040 *ay(H->GetAxis(kYrez));
1041 //*az(H->GetAxis(kZrez)),
1042 //*ap(H->GetAxis(kPrez));
1043 Int_t neta(aeta->GetNbins()), nphi(aphi->GetNbins()), rcBin(as->GetNbins()/2 + 1);
1044 TH3I *h3[fgNproj[cidx]];
1045 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1046 h3[ily] = new TH3I(Form("h3CL%d", ily), Form("r-#phi residuals @ ly[%d];%s;%s;%s", ily, aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
1047 neta, aeta->GetXmin(), aeta->GetXmax(),
1048 nphi, aphi->GetXmin(), aphi->GetXmax(),
1049 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1050 }
1051 Int_t coord[AliTRDresolution::kNdim]; memset(coord, 0, sizeof(Int_t) * AliTRDresolution::kNdim); Double_t v = 0.;
1052 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1053 v = H->GetBinContent(ib, coord);
1054 if(v<1.) continue;
1055 if(coord[kSpeciesChgRC]==rcBin) continue; // row cross
1056 Int_t ily(coord[kBC]-1);
1057 h3[ily]->AddBinContent(h3[ily]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1058 }
1059
1060 TF1 fg("fg", "gaus", ay->GetXmin(), ay->GetXmax());
1061 if(!fProj){
1062 AliInfo("Building array of projections ...");
1063 fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1064 }
1065 TObjArray *arr(NULL);
1066 fProj->AddAt(arr = new TObjArray(fgNproj[cidx]), cidx);
1067
1068 TH2F *h2(NULL);
1069 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1070 arr->AddAt(h2 = new TH2F(Form("h2CLs%d", ily),
1071 Form("Cl Resolution Ly[%d];%s;%s;#sigmay [#mum]", ily, aeta->GetTitle(), aphi->GetTitle()),
1072 neta, aeta->GetXmin(), aeta->GetXmax(),
1073 nphi, aphi->GetXmin(), aphi->GetXmax()), ily);
1074 TAxis *ax = h2->GetZaxis();
1075 ax->CenterTitle(); ax->SetTitleOffset(1.3);
1076 ax->SetRangeUser(250, 500);
1077 for(Int_t iphi(1); iphi<=nphi; iphi++){
1078 for(Int_t ieta(1); ieta<=neta; ieta++){
1079 TH1 *h = h3[ily]->ProjectionZ(Form("h1CLs%d", ily), ieta, ieta, iphi, iphi);
0d69a485 1080 Int_t ne(Int_t(h->GetEntries()));
3ceb45ae 1081 if(ne<100.) h2->SetBinContent(ieta, iphi, -999);
1082 else {
1083 fg.SetParameter(0, ne);fg.SetParameter(1, 0.);fg.SetParameter(2, 0.05);
1084 h->Fit(&fg, "QW0");
1085 Float_t val=TMath::Max(250., 1.e4*fg.GetParameter(2));
1086 h2->SetBinContent(ieta, iphi, val);
1087 }
1088 }
00a3f7a4 1089 }
1090 }
3ceb45ae 1091 for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++) delete h3[iproj];
1092 return kTRUE;
00a3f7a4 1093}
1094
3ceb45ae 1095//________________________________________________________
1096Bool_t AliTRDresolution::MakeProjectionTracklet()
1097{
1098// Analyse tracklet
3ceb45ae 1099 Int_t cidx = kTracklet;
1100 if(fProj && fProj->At(cidx)) return kTRUE;
1101 if(!fContainer){
1102 AliError("Missing data container.");
1103 return kFALSE;
1104 }
1105 THnSparse *H(NULL);
1106 if(!(H = (THnSparse*)fContainer->At(cidx))){
dbb2e0a7 1107// AliError(Form("Missing/Wrong data @ %d.", cidx));
3ceb45ae 1108 return kFALSE;
1109 }
395d3507 1110 AliDebug(2, Form("%s[%d]", H->GetName(), H->GetNdimensions()));
3ceb45ae 1111 return kTRUE;
1112}
068e2c00 1113
00a3f7a4 1114//________________________________________________________
3ceb45ae 1115Bool_t AliTRDresolution::MakeProjectionTrackIn()
1ee39b3a 1116{
3ceb45ae 1117// Analyse track in
61f6b45e 1118
3ceb45ae 1119 Int_t cidx = kTrackIn;
1120 if(fProj && fProj->At(cidx)) return kTRUE;
1121 if(!fContainer){
1122 AliError("Missing data container.");
1123 return kFALSE;
1124 }
1125 THnSparse *H(NULL);
1126 if(!(H = (THnSparse*)fContainer->At(cidx))){
1127 AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1ee39b3a 1128 return kFALSE;
1129 }
61f6b45e 1130
3ceb45ae 1131 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1132 TAxis //*abc(H->GetAxis(kBC)),
1133 *aphi(H->GetAxis(kPhi)),
1134 *aeta(H->GetAxis(kEta)),
1135 //*as(H->GetAxis(kSpeciesChgRC)),
1136 //*apt(H->GetAxis(kPt)),
1137 *ay(H->GetAxis(kYrez)),
1138 *az(H->GetAxis(kZrez)),
1139 *ap(H->GetAxis(kPrez));
1140 Int_t neta(aeta->GetNbins()), nphi(aphi->GetNbins());
1141 TH3I *h3[fgNproj[cidx]];
1142 h3[0] = new TH3I("h3TI0", Form("r-#phi residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
1143 neta, aeta->GetXmin(), aeta->GetXmax(),
1144 nphi, aphi->GetXmin(), aphi->GetXmax(),
1145 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1146 h3[1] = (TH3I*)h3[0]->Clone("h3TI1"); h3[1]->SetTitle("r-#phi residuals for pos tracks");
1147 h3[2] = new TH3I("h3TI2", Form("z residuals for row cross;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), az->GetTitle()),
1148 neta, aeta->GetXmin(), aeta->GetXmax(),
1149 nphi, aphi->GetXmin(), aphi->GetXmax(),
1150 az->GetNbins(), az->GetXmin(), az->GetXmax());
1151 h3[3] = new TH3I("h3TI3", Form("angular residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ap->GetTitle()),
1152 neta, aeta->GetXmin(), aeta->GetXmax(),
1153 nphi, aphi->GetXmin(), aphi->GetXmax(),
1154 ap->GetNbins(), ap->GetXmin(), ap->GetXmax());
1155 h3[4] = (TH3I*)h3[3]->Clone("h3TI4"); h3[4]->SetTitle("angular residuals for pos tracks");
1156 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1157 v = H->GetBinContent(ib, coord);
1158 if(v<1.) continue;
1159 if(coord[kBC]>1) continue; // bunch cross cut
1160 // species selection
1161 if(coord[kSpeciesChgRC]<6){
1162 h3[0]->AddBinContent(
1163 h3[0]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1164 h3[3]->AddBinContent(
1165 h3[3]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1166 } else if(coord[kSpeciesChgRC]==6) {
1167 h3[2]->AddBinContent(
0d69a485 1168 h3[2]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
3ceb45ae 1169 } else if(coord[kSpeciesChgRC]>6) {
1170 h3[1]->AddBinContent(
1171 h3[1]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1172 h3[4]->AddBinContent(
1173 h3[4]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1174 }
1175 }
1176 if(!fProj){
1177 AliInfo("Building array of projections ...");
1178 fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1179 }
1180 TObjArray *arr(NULL);
1181 fProj->AddAt(arr = new TObjArray(fgNproj[cidx]), cidx);
1182
1183 TH2F *h2(NULL);
1184 for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++){
1185 TAxis *ax(h3[iproj]->GetZaxis());
1186 Float_t zm(ax->GetXmin()/3.), zM(ax->GetXmax()/3.), dz=(zM-zm)/50;
1187 arr->AddAt(h2 = new TH2F(Form("h2TI%d", iproj),
1188 Form("%s;%s;%s;%s", h3[iproj]->GetTitle(), aeta->GetTitle(), aphi->GetTitle(), ax->GetTitle()),
1189 neta, aeta->GetXmin(), aeta->GetXmax(),
1190 nphi, aphi->GetXmin(), aphi->GetXmax()), iproj);
1191 h2->SetContour(50);
1192 h2->GetZaxis()->CenterTitle();
1193 h2->GetZaxis()->SetRangeUser(zm, zM);
1194 zm+=dz; zM-=dz;
1195 for(Int_t iphi(1); iphi<=nphi; iphi++){
1196 for(Int_t ieta(1); ieta<=neta; ieta++){
1197 TH1 *h = h3[iproj]->ProjectionZ(Form("hy%d", iproj), ieta, ieta, iphi, iphi);
1198 h2->SetBinContent(ieta, iphi, GetMeanWithBoundary(h, zm, zM, dz));
1ee39b3a 1199 }
1200 }
1201 }
3ceb45ae 1202// h2[5] = (TH2F*)h2[0]->Clone("h25");
1203// h2[5]->SetTitle("Systematic shift between neg/pos tracks");
1204// h2[5]->SetZTitle("#Delta(#Delta^{-}y - #Delta^{+}y) [cm]"); h2[5]->Reset();
1205// h2[6] = (TH2F*)h2[1]->Clone("h26");
1206// h2[6]->SetTitle("Average shift of pos&neg tracks");
1207// h2[6]->SetZTitle("<#Delta^{-}y, #Delta^{+}y> [cm]"); h2[6]->Reset();
1208// for(Int_t iphi(1); iphi<=nphi; iphi++){
1209// for(Int_t ieta(1); ieta<=neta; ieta++){
1210// Float_t neg = h2[0]->GetBinContent(ieta, iphi),
1211// pos = h2[1]->GetBinContent(ieta, iphi);
1212// if(neg<-100 || pos<-100){
1213// h2[5]->SetBinContent(ieta, iphi, -999.);
1214// h2[6]->SetBinContent(ieta, iphi, -999.);
1215// } else {
1216// h2[5]->SetBinContent(ieta, iphi, neg-pos);
1217// h2[6]->SetBinContent(ieta, iphi, 0.5*(neg+pos));
1218// }
1219// }
1220// }
1221
1222 for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++) delete h3[iproj];
1223 return kTRUE;
1224}
1ee39b3a 1225
3ceb45ae 1226
1227
1228//________________________________________________________
1229Bool_t AliTRDresolution::PostProcess()
1230{
1231// Fit, Project, Combine, Extract values from the containers filled during execution
1232
1233 if (!fContainer) {
1234 AliError("ERROR: list not available");
1235 return kFALSE;
1236 }
1ee39b3a 1237
1238 // DEFINE MODELS
1239 // simple gauss
1240 TF1 fg("fGauss", "gaus", -.5, .5);
1241 // Landau for charge resolution
92c40e64 1242 TF1 fch("fClCh", "landau", 0., 1000.);
1243 // Landau for e+- pt resolution
1244 TF1 fpt("fPt", "landau", -0.1, 0.2);
1ee39b3a 1245
1246 //PROCESS EXPERIMENTAL DISTRIBUTIONS
1ee39b3a 1247 // Clusters residuals
3ceb45ae 1248 if(!MakeProjectionCluster()) return kFALSE;
6558fd69 1249 fNRefFigures = 3;
1ee39b3a 1250 // Tracklet residual/pulls
3ceb45ae 1251 if(!MakeProjectionTracklet()) return kFALSE;
6558fd69 1252 fNRefFigures = 7;
a310e49b 1253 // TRDin residual/pulls
3ceb45ae 1254 if(!MakeProjectionTrackIn()) return kFALSE;
6558fd69 1255 fNRefFigures = 11;
a310e49b 1256 // TRDout residual/pulls
3ceb45ae 1257// if(!MakeProjectionTrackOut()) return kFALSE;
6558fd69 1258 fNRefFigures = 15;
1ee39b3a 1259
1260 if(!HasMCdata()) return kTRUE;
1261
1262
1263 //PROCESS MC RESIDUAL DISTRIBUTIONS
1264
1265 // CLUSTER Y RESOLUTION/PULLS
3ceb45ae 1266// if(!MakeProjectionClusterMC()) return kFALSE;
6558fd69 1267 fNRefFigures = 17;
1ee39b3a 1268
1269 // TRACKLET RESOLUTION/PULLS
3ceb45ae 1270// if(!MakeProjectionTrackletMC()) return kFALSE;
d25116b6 1271 fNRefFigures = 21;
1ee39b3a 1272
1273 // TRACK RESOLUTION/PULLS
3ceb45ae 1274/* Process3Darray(kMCtrack, 0, &fg, 1.e4); // y
81979445 1275 Process3DlinkedArray(kMCtrack, 1, &fg); // y PULL
1276 Process3Darray(kMCtrack, 2, &fg, 1.e4); // z
1277 Process3Darray(kMCtrack, 3, &fg); // z PULL
92d6d80c 1278 Process2Darray(kMCtrack, 4, &fg, 1.e3); // phi
1279 Process2Darray(kMCtrack, 5, &fg); // snp PULL
1280 Process2Darray(kMCtrack, 6, &fg, 1.e3); // theta
1281 Process2Darray(kMCtrack, 7, &fg); // tgl PULL
1282 Process3Darray(kMCtrack, 8, &fg, 1.e2); // pt resolution
1283 Process3Darray(kMCtrack, 9, &fg); // 1/pt pulls
3ceb45ae 1284 Process3Darray(kMCtrack, 10, &fg, 1.e2); // p resolution*/
1285// if(!MakeProjectionTrackMC(kMCtrack)) return kFALSE;
d25116b6 1286 fNRefFigures+=16;
92d6d80c 1287
1288 // TRACK TRDin RESOLUTION/PULLS
3ceb45ae 1289// if(!MakeProjectionTrackMC(kMCtrackIn)) return kFALSE;
d25116b6 1290 fNRefFigures+=8;
92d6d80c 1291
1292 // TRACK TRDout RESOLUTION/PULLS
3ceb45ae 1293// if(!MakeProjectionTrackMC(kMCtrackOut)) return kFALSE;
d25116b6 1294 fNRefFigures+=8;
1ee39b3a 1295
1296 return kTRUE;
1297}
1298
1299
1300//________________________________________________________
1301void AliTRDresolution::Terminate(Option_t *opt)
1302{
1303 AliTRDrecoTask::Terminate(opt);
1304 if(HasPostProcess()) PostProcess();
1305}
1306
1307//________________________________________________________
1308void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1309{
1310// Helper function to avoid duplication of code
1311// Make first guesses on the fit parameters
1312
1313 // find the intial parameters of the fit !! (thanks George)
1314 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1315 Double_t sum = 0.;
1316 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1317 f->SetParLimits(0, 0., 3.*sum);
1318 f->SetParameter(0, .9*sum);
1319 Double_t rms = h->GetRMS();
1320 f->SetParLimits(1, -rms, rms);
1321 f->SetParameter(1, h->GetMean());
1322
1323 f->SetParLimits(2, 0., 2.*rms);
1324 f->SetParameter(2, rms);
1325 if(f->GetNpar() <= 4) return;
1326
1327 f->SetParLimits(3, 0., sum);
1328 f->SetParameter(3, .1*sum);
1329
1330 f->SetParLimits(4, -.3, .3);
1331 f->SetParameter(4, 0.);
1332
1333 f->SetParLimits(5, 0., 1.e2);
1334 f->SetParameter(5, 2.e-1);
1335}
1336
afca20ef 1337//________________________________________________________
9dcc64cc 1338TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
afca20ef 1339{
3d2a3dff 1340// Build performance histograms for AliTRDcluster.vs TRD track or MC
afca20ef 1341// - y reziduals/pulls
3d2a3dff 1342
1343 TObjArray *arr = new TObjArray(2);
afca20ef 1344 arr->SetName(name); arr->SetOwner();
1345 TH1 *h(NULL); char hname[100], htitle[300];
1346
1347 // tracklet resolution/pull in y direction
7fe4e88b 1348 snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
3ceb45ae 1349 snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
9dcc64cc 1350 Float_t rr = range<0.?fDyRange:range;
3d2a3dff 1351 if(!(h = (TH3S*)gROOT->FindObject(hname))){
3ceb45ae 1352 Int_t nybins=50;
2589cf64 1353 if(expand) nybins*=2;
3d2a3dff 1354 h = new TH3S(hname, htitle,
6859821f 1355 48, -.48, .48, // phi
9dcc64cc 1356 60, -rr, rr, // dy
6859821f 1357 nybins, -0.5, nybins-0.5);// segment
afca20ef 1358 } else h->Reset();
1359 arr->AddAt(h, 0);
7fe4e88b 1360 snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
3ceb45ae 1361 snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, "Detector");
81979445 1362 if(!(h = (TH3S*)gROOT->FindObject(hname))){
3ceb45ae 1363 h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
afca20ef 1364 } else h->Reset();
1365 arr->AddAt(h, 1);
1366
3d2a3dff 1367 return arr;
1368}
1369
1370//________________________________________________________
2589cf64 1371TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
3d2a3dff 1372{
1373// Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1374// - y reziduals/pulls
1375// - z reziduals/pulls
1376// - phi reziduals
9dcc64cc 1377 TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
3d2a3dff 1378 arr->Expand(5);
1379 TH1 *h(NULL); char hname[100], htitle[300];
1380
afca20ef 1381 // tracklet resolution/pull in z direction
7fe4e88b 1382 snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
9dcc64cc 1383 snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1384 if(!(h = (TH2S*)gROOT->FindObject(hname))){
1385 h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
afca20ef 1386 } else h->Reset();
1387 arr->AddAt(h, 2);
7fe4e88b 1388 snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1389 snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
81979445 1390 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1391 h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
dfd7d48b 1392 h->GetZaxis()->SetBinLabel(1, "no RC");
1393 h->GetZaxis()->SetBinLabel(2, "RC");
afca20ef 1394 } else h->Reset();
1395 arr->AddAt(h, 3);
1396
1397 // tracklet to track phi resolution
7fe4e88b 1398 snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
3ceb45ae 1399 snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1400 Int_t nsgms=540;
9dcc64cc 1401 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1402 h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
afca20ef 1403 } else h->Reset();
1404 arr->AddAt(h, 4);
1405
1406 return arr;
1407}
1408
1409//________________________________________________________
1410TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1411{
1412// Build performance histograms for AliExternalTrackParam.vs MC
1413// - y resolution/pulls
1414// - z resolution/pulls
1415// - phi resolution, snp pulls
1416// - theta resolution, tgl pulls
1417// - pt resolution, 1/pt pulls, p resolution
1418
afca20ef 1419 TObjArray *arr = BuildMonitorContainerTracklet(name);
1420 arr->Expand(11);
3d2a3dff 1421 TH1 *h(NULL); char hname[100], htitle[300];
395d3507 1422 //TAxis *ax(NULL);
3d2a3dff 1423
afca20ef 1424 // snp pulls
7fe4e88b 1425 snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1426 snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
afca20ef 1427 if(!(h = (TH2I*)gROOT->FindObject(hname))){
1428 h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1429 } else h->Reset();
1430 arr->AddAt(h, 5);
1431
1432 // theta resolution
7fe4e88b 1433 snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1434 snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
afca20ef 1435 if(!(h = (TH2I*)gROOT->FindObject(hname))){
1436 h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1437 } else h->Reset();
1438 arr->AddAt(h, 6);
1439 // tgl pulls
7fe4e88b 1440 snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1441 snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
afca20ef 1442 if(!(h = (TH2I*)gROOT->FindObject(hname))){
1443 h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1444 } else h->Reset();
1445 arr->AddAt(h, 7);
1446
afca20ef 1447 const Int_t kNdpt(150);
1448 const Int_t kNspc = 2*AliPID::kSPECIES+1;
61f6b45e 1449 Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
afca20ef 1450 Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
61f6b45e 1451 for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1452 for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1453 for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
afca20ef 1454
1455 // Pt resolution
7fe4e88b 1456 snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1457 snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
afca20ef 1458 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1459 h = new TH3S(hname, htitle,
1460 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
dbb2e0a7 1461 //ax = h->GetZaxis();
3ceb45ae 1462 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
afca20ef 1463 } else h->Reset();
1464 arr->AddAt(h, 8);
1465 // 1/Pt pulls
7fe4e88b 1466 snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1467 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 1468 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1469 h = new TH3S(hname, htitle,
1470 kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
395d3507 1471 //ax = h->GetZaxis();
3ceb45ae 1472 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
afca20ef 1473 } else h->Reset();
1474 arr->AddAt(h, 9);
1475 // P resolution
7fe4e88b 1476 snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1477 snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
afca20ef 1478 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1479 h = new TH3S(hname, htitle,
1480 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
395d3507 1481 //ax = h->GetZaxis();
3ceb45ae 1482 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
afca20ef 1483 } else h->Reset();
1484 arr->AddAt(h, 10);
1485
1486 return arr;
1487}
1488
1489
1ee39b3a 1490//________________________________________________________
1491TObjArray* AliTRDresolution::Histos()
1492{
1493 //
1494 // Define histograms
1495 //
1496
1497 if(fContainer) return fContainer;
1498
3ceb45ae 1499 fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1500 THnSparse *H(NULL);
1501 const Int_t nhn(100); Char_t hn[nhn]; TString st;
1ee39b3a 1502
3ceb45ae 1503 //++++++++++++++++++++++
1ee39b3a 1504 // cluster to track residuals/pulls
3ceb45ae 1505 snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1506 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
0d69a485 1507 const Char_t *clTitle[5/*kNdim*/] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], "chg*Q/vd/angle", fgkTitle[kYrez]/*, "#Deltax [cm]", "Q</Q", "#Phi^{*} - ExB [deg]"*/};
1508 const Int_t clNbins[5/*kNdim*/] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], 51, fgkNbins[kYrez]/*, 33, 10, 30*/};
1509 const Double_t clMin[5/*kNdim*/] = {-0.5, fgkMin[kPhi], fgkMin[kEta], -102., fgkMin[kYrez]/3./*, 0., 0.1, -30*/},
1510 clMax[5/*kNdim*/] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], 102., fgkMax[kYrez]/3./*, 3.3, 2.1, 30*/};
3ceb45ae 1511 st = "cluster spatial&charge resolution;";
0d69a485 1512 for(Int_t idim(0); idim<5/*kNdim*/; idim++){ st += clTitle[idim]; st+=";";}
3ceb45ae 1513 H = new THnSparseI(hn, st.Data(), kNdim, clNbins, clMin, clMax);
1514 } else H->Reset();
1515 fContainer->AddAt(H, kCluster);
1516 //++++++++++++++++++++++
afca20ef 1517 // tracklet to TRD track
3ceb45ae 1518 snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1519 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
0d69a485 1520 const Char_t *trTitle[kNdim+1] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kSpeciesChgRC], fgkTitle[kYrez], "#Deltaz [cm]/#Phi^{*} - ExB [rad]", fgkTitle[kPrez], "dq/dl [a.u.]"/*, fgkTitle[kPt]*/};
1521 const Int_t trNbins[kNdim+1] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kSpeciesChgRC], fgkNbins[kYrez], fgkNbins[kZrez], fgkNbins[kPrez], 30/*, fgkNbins[kPt]*/};
1522 const Double_t trMin[kNdim+1] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kSpeciesChgRC], fgkMin[kYrez], fgkMin[kZrez], fgkMin[kPrez], 0./*, fgkMin[kPt]*/},
1523 trMax[kNdim+1] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kSpeciesChgRC], fgkMax[kYrez], fgkMax[kZrez], fgkMax[kPrez], 3000./*, fgkMax[kPt]*/};
3ceb45ae 1524 st = "tracklet spatial&charge resolution;";
1525 for(Int_t idim(0); idim<kNdim+1; idim++){ st += trTitle[idim]; st+=";";}
1526 H = new THnSparseI(hn, st.Data(), kNdim+1, trNbins, trMin, trMax);
1527 } else H->Reset();
1528 fContainer->AddAt(H, kTracklet);
1529 //++++++++++++++++++++++
afca20ef 1530 // tracklet to TRDin
3ceb45ae 1531 snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1532 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1533 st = "r-#phi/z/angular residuals @ TRD entry;";
1534 for(Int_t idim(0); idim<kNdim; idim++){ st += fgkTitle[idim]; st+=";";}
1535 H = new THnSparseI(hn, st.Data(), kNdim, fgkNbins, fgkMin, fgkMax);
1536 } else H->Reset();
1537 fContainer->AddAt(H, kTrackIn);
a310e49b 1538 // tracklet to TRDout
1539 fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1ee39b3a 1540
1541
1542 // Resolution histos
1543 if(!HasMCdata()) return fContainer;
1544
3d2a3dff 1545 // cluster resolution
1546 fContainer->AddAt(BuildMonitorContainerCluster("MCcl"), kMCcluster);
1ee39b3a 1547
3d2a3dff 1548 // tracklet resolution
1549 fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
1ee39b3a 1550
3d2a3dff 1551 // track resolution
3ceb45ae 1552 TObjArray *arr(NULL);
92d6d80c 1553 fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
3d2a3dff 1554 arr->SetName("MCtrk");
1555 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
31c8fa8a 1556
afca20ef 1557 // TRDin TRACK RESOLUTION
3ceb45ae 1558 fContainer->AddAt(H, kMCtrackIn);
1ee39b3a 1559
afca20ef 1560 // TRDout TRACK RESOLUTION
92d6d80c 1561 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1ee39b3a 1562
1563 return fContainer;
1564}
1565
5468a262 1566//________________________________________________________
1567Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1568{
b37d601d 1569// Robust function to process sigma/mean for 2D plot dy(x)
1570// For each x bin a gauss fit is performed on the y projection and the range
1571// with the minimum chi2/ndf is choosen
5468a262 1572
1573 if(!h2) {
1574 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
1575 return kFALSE;
1576 }
1577 if(!Int_t(h2->GetEntries())){
1578 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
1579 return kFALSE;
1580 }
1581 if(!g || !g[0]|| !g[1]) {
1582 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
1583 return kFALSE;
1584 }
1585 // prepare
1586 TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
b37d601d 1587 Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
1588 TF1 f("f", "gaus", ymin, ymax);
5468a262 1589 Int_t n(0);
1590 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1591 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1592 TH1D *h(NULL);
1593 if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
b37d601d 1594 Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
1595
5468a262 1596
1597 // do actual loop
b37d601d 1598 Float_t chi2OverNdf(0.);
5468a262 1599 for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
b37d601d 1600 x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
1601 ymin = ay->GetXmin(); ymax = ay->GetXmax();
1602
5468a262 1603 h = h2->ProjectionY("py", ix, ix);
1604 if((n=(Int_t)h->GetEntries())<stat){
1605 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
1606 continue;
1607 }
b37d601d 1608 // looking for a first order mean value
1609 f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
1610 h->Fit(&f, "QNW");
1611 chi2OverNdf = f.GetChisquare()/f.GetNDF();
1612 printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
1613 y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
1614 ey = f.GetParError(1);
1615 sy = f.GetParameter(2); esy = f.GetParError(2);
1616// // looking for the best chi2/ndf value
1617// while(ymin<y0 && ymax>y1){
1618// f.SetParameter(1, y);
1619// f.SetParameter(2, sy);
1620// h->Fit(&f, "QNW", "", y0, y1);
1621// printf(" range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
1622// if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
1623// chi2OverNdf = f.GetChisquare()/f.GetNDF();
1624// y = f.GetParameter(1); ey = f.GetParError(1);
1625// sy = f.GetParameter(2); esy = f.GetParError(2);
1626// printf(" set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
1627// }
1628// y0-=dy; y1+=dy;
1629// }
1630 g[0]->SetPoint(np, x, y);
1631 g[0]->SetPointError(np, ex, ey);
1632 g[1]->SetPoint(np, x, sy);
1633 g[1]->SetPointError(np, ex, esy);
5468a262 1634 np++;
1635 }
1636 return kTRUE;
1637}
1638
2589cf64 1639
1ee39b3a 1640//________________________________________________________
1641Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1642{
1643 //
1644 // Do the processing
1645 //
1646
7fe4e88b 1647 Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
1ee39b3a 1648 Int_t n = 0;
1649 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1650 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2589cf64 1651 if(Int_t(h2->GetEntries())){
1652 AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
1653 } else {
1654 AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
1655 fIdxPlot++;
1656 return kTRUE;
1657 }
1ee39b3a 1658
dfd7d48b 1659 const Int_t kINTEGRAL=1;
1660 for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
1661 Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
1662 Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
1663 TH1D *h = h2->ProjectionY(pn, abin, bbin);
068e2c00 1664 if((n=(Int_t)h->GetEntries())<150){
2589cf64 1665 AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
1666 continue;
1667 }
1ee39b3a 1668 h->Fit(f, "QN");
1ee39b3a 1669 Int_t ip = g[0]->GetN();
2589cf64 1670 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 1671 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1672 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1673 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1674 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1ee39b3a 1675/*
1676 g[0]->SetPoint(ip, x, k*h->GetMean());
1677 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1678 g[1]->SetPoint(ip, x, k*h->GetRMS());
1679 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1680 }
1681 fIdxPlot++;
1682 return kTRUE;
1683}
1684
1ee39b3a 1685
08bdd9d1 1686//____________________________________________________________________
1687Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
1688{
1689//
1690// Fit track with a staight line using the "np" clusters stored in the array "points".
1691// The following particularities are stored in the clusters from points:
1692// 1. pad tilt as cluster charge
1693// 2. pad row cross or vertex constrain as fake cluster with cluster type 1
1694// The parameters of the straight line fit are stored in the array "param" in the following order :
1695// param[0] - x0 reference radial position
1696// param[1] - y0 reference r-phi position @ x0
1697// param[2] - z0 reference z position @ x0
1698// param[3] - slope dy/dx
1699// param[4] - slope dz/dx
1700//
1701// Attention :
1702// Function should be used to refit tracks for B=0T
1703//
1704
1705 if(np<40){
b37d601d 1706 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
08bdd9d1 1707 return kFALSE;
1708 }
1709 TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
1710
1711 Double_t x0(0.);
1712 for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
1713 x0/=Float_t(np);
1714
00a56108 1715 Double_t x, y, z, dx, tilt(0.);
08bdd9d1 1716 for(Int_t ip(0); ip<np; ip++){
1717 x = points[ip].GetX(); z = points[ip].GetZ();
1718 dx = x - x0;
1719 zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
1720 }
1721 if(zfitter.Eval() != 0) return kFALSE;
1722
1723 Double_t z0 = zfitter.GetParameter(0);
1724 Double_t dzdx = zfitter.GetParameter(1);
1725 for(Int_t ip(0); ip<np; ip++){
1726 if(points[ip].GetClusterType()) continue;
1727 x = points[ip].GetX();
1728 dx = x - x0;
1729 y = points[ip].GetY();
1730 z = points[ip].GetZ();
1731 tilt = points[ip].GetCharge();
1732 y -= tilt*(-dzdx*dx + z - z0);
1733 Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
1734 yfitter.AddPoint(&dx, y, 1.);
1735 }
1736 if(yfitter.Eval() != 0) return kFALSE;
1737 Double_t y0 = yfitter.GetParameter(0);
1738 Double_t dydx = yfitter.GetParameter(1);
1739
1740 param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
b37d601d 1741 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 1742 return kTRUE;
1743}
1744
5f7f1c07 1745//____________________________________________________________________
b37d601d 1746Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
5f7f1c07 1747{
1748//
1749// Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
1750// See function FitTrack for the data stored in the "clusters" array
1751
1752// The parameters of the straight line fit are stored in the array "param" in the following order :
1753// par[0] - x0 reference radial position
1754// par[1] - y0 reference r-phi position @ x0
1755// par[2] - slope dy/dx
1756//
1757// Attention :
1758// Function should be used to refit tracks for B=0T
1759//
1760
1761 TLinearFitter yfitter(2, "pol1");
1762
1763 // grep data for tracklet
1764 Double_t x0(0.), x[60], y[60], dy[60];
1765 Int_t nly(0);
1766 for(Int_t ip(0); ip<np; ip++){
1767 if(points[ip].GetClusterType()) continue;
1768 if(points[ip].GetVolumeID() != ly) continue;
1769 Float_t xt(points[ip].GetX())
1770 ,yt(param[1] + param[3] * (xt - param[0]));
1771 x[nly] = xt;
1772 y[nly] = points[ip].GetY();
1773 dy[nly]= y[nly]-yt;
1774 x0 += xt;
1775 nly++;
1776 }
1777 if(nly<10){
b37d601d 1778 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
5f7f1c07 1779 return kFALSE;
1780 }
1781 // set radial reference for fit
1782 x0 /= Float_t(nly);
1783
1784 // find tracklet core
1785 Double_t mean(0.), sig(1.e3);
1786 AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
1787
1788 // simple cluster error parameterization
1789 Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
1790
1791 // fit tracklet core
1792 for(Int_t jly(0); jly<nly; jly++){
1793 if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
1794 Double_t dx(x[jly]-x0);
1795 yfitter.AddPoint(&dx, y[jly], 1.);
1796 }
1797 if(yfitter.Eval() != 0) return kFALSE;
1798 par[0] = x0;
1799 par[1] = yfitter.GetParameter(0);
1800 par[2] = yfitter.GetParameter(1);
1801 return kTRUE;
1802}
1803
00a56108 1804//____________________________________________________________________
b37d601d 1805Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
00a56108 1806{
1807//
1808// Global selection mechanism of tracksbased on cluster to fit residuals
1809// The parameters are the same as used ni function FitTrack().
1810
1811 const Float_t kS(0.6), kM(0.2);
1812 TH1S h("h1", "", 100, -5.*kS, 5.*kS);
1813 Float_t dy, dz, s, m;
1814 for(Int_t ip(0); ip<np; ip++){
1815 if(points[ip].GetClusterType()) continue;
1816 Float_t x0(points[ip].GetX())
1817 ,y0(param[1] + param[3] * (x0 - param[0]))
1818 ,z0(param[2] + param[4] * (x0 - param[0]));
1819 dy=points[ip].GetY() - y0; h.Fill(dy);
1820 dz=points[ip].GetZ() - z0;
1821 }
1822 TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
1823 fg.SetParameter(1, 0.);
1824 fg.SetParameter(2, 2.e-2);
1825 h.Fit(&fg, "QN");
1826 m=fg.GetParameter(1); s=fg.GetParameter(2);
1827 if(s>kS || TMath::Abs(m)>kM) return kFALSE;
1828 return kTRUE;
1829}
1830
1ee39b3a 1831//________________________________________________________
1832void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
1833{
1834 //
1835 // Get the most probable value and the full width half mean
1836 // of a Landau distribution
1837 //
1838
1839 const Float_t dx = 1.;
1840 mpv = f->GetParameter(1);
1841 Float_t fx, max = f->Eval(mpv);
1842
1843 xm = mpv - dx;
1844 while((fx = f->Eval(xm))>.5*max){
1845 if(fx>max){
1846 max = fx;
1847 mpv = xm;
1848 }
1849 xm -= dx;
1850 }
1851
1852 xM += 2*(mpv - xm);
1853 while((fx = f->Eval(xM))>.5*max) xM += dx;
1854}
1855
1856
3ceb45ae 1857// #include "TFile.h"
1858// //________________________________________________________
1859// Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
1860// {
1861// if(!file){
1862// AliWarning("Use cluster position as in reconstruction.");
1863// SetLoadCorrection();
1864// return kTRUE;
1865// }
1866// TDirectory *cwd(gDirectory);
1867// TString fileList;
1868// FILE *filePtr = fopen(file, "rt");
1869// if(!filePtr){
1870// AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
1871// SetLoadCorrection();
1872// return kFALSE;
1873// }
1874// TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
1875// while(fileList.Gets(filePtr)){
1876// if(!TFile::Open(fileList.Data())) {
1877// AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
1878// continue;
1879// } else AliInfo(Form("\"%s\"", fileList.Data()));
1880//
1881// TTree *tSys = (TTree*)gFile->Get("tSys");
1882// h2->SetDirectory(gDirectory); h2->Reset("ICE");
1883// tSys->Draw("det:t>>h2", "dx", "goff");
1884// for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1885// for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
1886// }
1887// h2->SetDirectory(cwd);
1888// gFile->Close();
1889// }
1890// cwd->cd();
1891//
1892// if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
1893// for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1894// printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
1895// for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1896// FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
1897// fprintf(fout, " static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
1898// for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1899// printf("%03d|", idet);
1900// fprintf(fout, " {");
1901// for(Int_t it(0); it<30; it++){
1902// printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
1903// fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
1904// }
1905// printf("\n");
1906// fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
1907// }
1908// fprintf(fout, " };\n");
1909// }
1910// SetLoadCorrection();
1911// return kTRUE;
1912// }