1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // TRD tracking resolution //
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
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
32 // For stand alone usage please refer to the following example:
34 // gSystem->Load("libANALYSIS.so");
35 // gSystem->Load("libTRDqaRec.so");
36 // AliTRDresolution *res = new AliTRDresolution();
37 // //res->SetMCdata();
38 // //res->SetVerbose();
39 // //res->SetVisual();
41 // if(!res->PostProcess()) return;
42 // res->GetRefFigure(0);
46 // Alexandru Bercuci <A.Bercuci@gsi.de> //
47 // Markus Fasel <M.Fasel@gsi.de> //
49 ////////////////////////////////////////////////////////////////////////////
54 #include <TObjArray.h>
58 #include <THnSparse.h>
64 #include <TGraphErrors.h>
65 #include <TGraphAsymmErrors.h>
66 #include <TLinearFitter.h>
70 #include <TTreeStream.h>
71 #include <TGeoManager.h>
72 #include <TDatabasePDG.h>
76 #include "AliESDtrack.h"
77 #include "AliMathBase.h"
78 #include "AliTrackPointArray.h"
80 #include "AliTRDresolution.h"
81 #include "AliTRDgeometry.h"
82 #include "AliTRDtransform.h"
83 #include "AliTRDpadPlane.h"
84 #include "AliTRDcluster.h"
85 #include "AliTRDseedV1.h"
86 #include "AliTRDtrackV1.h"
87 #include "AliTRDReconstructor.h"
88 #include "AliTRDrecoParam.h"
89 #include "AliTRDpidUtil.h"
90 #include "AliTRDinfoGen.h"
91 #include "AliAnalysisManager.h"
92 #include "info/AliTRDclusterInfo.h"
93 #include "info/AliTRDeventInfo.h"
95 ClassImp(AliTRDresolution)
96 ClassImp(AliTRDresolution::AliTRDresolutionProjection)
98 Int_t const AliTRDresolution::fgkNbins[kNdim] = {
99 Int_t(kNbunchCross)/*bc*/,
107 }; //! no of bins/projection
108 Double_t const AliTRDresolution::fgkMin[kNdim] = {
117 }; //! low limits for projections
118 Double_t const AliTRDresolution::fgkMax[kNdim] = {
127 }; //! high limits for projections
128 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
137 }; //! title of projection
139 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
148 // ,"Tracklet2TRDout"
151 Float_t AliTRDresolution::fgPtBin[kNpt+1];
153 //________________________________________________________
154 AliTRDresolution::AliTRDresolution()
168 // Default constructor
170 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
171 MakePtSegmentation();
174 //________________________________________________________
175 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
176 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
189 // Default constructor
193 MakePtSegmentation();
195 SetUseExchangeContainers();
196 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
197 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
201 //________________________________________________________
202 AliTRDresolution::~AliTRDresolution()
207 if (AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
208 if(fProj){fProj->Delete(); delete fProj;}
209 if(fCl){fCl->Delete(); delete fCl;}
210 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
214 //________________________________________________________
215 void AliTRDresolution::UserCreateOutputObjects()
217 // spatial resolution
219 AliTRDrecoTask::UserCreateOutputObjects();
220 if(UseExchangeContainers()) InitExchangeContainers();
223 //________________________________________________________
224 void AliTRDresolution::InitExchangeContainers()
226 // Init containers for subsequent tasks (AliTRDclusterResolution)
228 fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
229 fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
230 PostData(kClToTrk, fCl);
231 PostData(kClToMC, fMCcl);
234 //________________________________________________________
235 void AliTRDresolution::UserExec(Option_t *opt)
241 if(fCl) fCl->Delete();
242 if(fMCcl) fMCcl->Delete();
243 AliTRDrecoTask::UserExec(opt);
246 //________________________________________________________
247 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
249 // Helper function to calculate pulls in the yz plane
250 // using proper tilt rotation
251 // Uses functionality defined by AliTRDseedV1.
255 Double_t t2(tilt*tilt);
256 // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
260 cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
261 cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
262 cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
264 Double_t sqr[3]={0., 0., 0.};
265 if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
266 Double_t invsqr[3]={0., 0., 0.};
267 if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
268 Double_t tmp(dyz[0]);
269 dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
270 dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
275 //________________________________________________________
276 TH1* AliTRDresolution::DetCluster(const TObjArray *cls)
279 // Plot the cluster distributions
282 if(cls) fkClusters = cls;
284 AliDebug(4, "No Clusters defined.");
288 if(!(ncl = fkClusters->GetEntriesFast())){
289 AliDebug(1, "No RecPoints defined.");
293 AliTRDcluster *cl(NULL);
294 for(Int_t icl(0); icl<ncl; icl++){
295 if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])) continue;
296 det = cl->GetDetector(); break;
299 AliDebug(1, "No useful clusters defined.");
303 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kDetector)))){
304 AliWarning("No output container defined.");
307 Int_t ly(AliTRDgeometry::GetLayer(det)),
308 stk(AliTRDgeometry::GetStack(det));
310 alpha((0.5+AliTRDgeometry::GetSector(det))*AliTRDgeometry::GetAlpha()),
311 cs(TMath::Cos(alpha)),
312 sn(TMath::Sin(alpha));
313 for(Int_t icl(0); icl<ncl; icl++){
314 if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])) continue;
316 val[kPhi] = TMath::ATan2(cl->GetX()*sn + cl->GetY()*cs, cl->GetX()*cs - cl->GetY()*sn);
317 val[kEta] = (5-stk)*16-cl->GetPadRow()-1-(stk<3?4:0);
318 val[kYrez]= TMath::Abs(cl->GetQ());
319 val[4] = fEvent->GetMultiplicity();
325 //________________________________________________________
326 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
329 // Plot the cluster distributions
332 if(track) fkTrack = track;
334 AliDebug(4, "No Track defined.");
337 if(fkESD && TMath::Abs(fkESD->GetTOFbc()) > 1){
338 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
341 if(fkESD && fPt<fPtThreshold){
342 AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
346 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
347 AliWarning("No output container defined.");
351 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
353 alpha(0.), cs(-2.), sn(0.); //Float_t exb, vd, t0, s2, dl, dt;
354 TObjArray *clInfoArr(NULL);
355 AliTRDseedV1 *fTracklet(NULL);
356 AliTRDcluster *c(NULL), *cc(NULL);
357 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
358 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
359 if(!fTracklet->IsOK()) continue;
360 //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
363 alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
364 cs = TMath::Cos(alpha);
365 sn = TMath::Sin(alpha);
367 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
368 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
369 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
370 val[kPt] = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
371 Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
373 Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
374 fTracklet->ResetClusterIter(kTRUE);
375 while((c = fTracklet->NextCluster())){
376 Float_t xc(c->GetX()),
377 q(TMath::Abs(c->GetQ()));
378 if(row0<0) row0 = c->GetPadRow();
380 val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
381 val[kPrez] = fTracklet->GetX0()-xc;
382 val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
383 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
384 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
385 if(ic) val[kZrez] /= (ic*q);
386 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
388 /* // tilt rotation of covariance for clusters
389 Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
390 cov[0] = (sy2+t2*sz2)*corr;
391 cov[1] = tilt*(sz2 - sy2)*corr;
392 cov[2] = (t2*sy2+sz2)*corr;
393 // sum with track covariance
394 cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
395 Double_t dyz[2]= {dy[1], dz[1]};
396 Pulls(dyz, cov, tilt);*/
398 // Get z-position with respect to anode wire
399 Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
400 zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
401 Int_t istk = geo->GetStack(c->GetDetector());
402 AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
403 Float_t rowZ = pp->GetRow0();
404 Float_t d = rowZ - zt + pp->GetAnodeWireOffset();
405 d -= ((Int_t)(2 * d)) / 2.0;
406 if (d > 0.25) d = 0.5 - d;
408 AliTRDclusterInfo *clInfo(NULL);
409 clInfo = new AliTRDclusterInfo;
410 clInfo->SetCluster(c);
411 //Float_t covF[] = {cov[0], cov[1], cov[2]};
412 clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
413 clInfo->SetResolution(val[kYrez]);
414 clInfo->SetAnisochronity(d);
415 clInfo->SetDriftLength(val[kZrez]);
416 clInfo->SetTilt(fTracklet->GetTilt());
417 if(fCl) fCl->Add(clInfo);
418 //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
422 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
423 clInfoArr->SetOwner(kFALSE);
425 clInfoArr->Add(clInfo);
428 if(DebugLevel()>=2 && clInfoArr){
429 ULong_t status = fkESD->GetStatus();
430 (*DebugStream()) << "cluster"
431 <<"status=" << status
432 <<"clInfo.=" << clInfoArr
437 if(clInfoArr) delete clInfoArr;
439 if(!track) return NULL;
440 // special care for EVE usage
442 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
443 return H->Projection(kYrez);
447 //________________________________________________________
448 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
450 // Plot normalized residuals for tracklets to track.
452 // We start from the result that if X=N(|m|, |Cov|)
454 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
456 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
457 // reference position.
458 if(track) fkTrack = track;
460 AliDebug(4, "No Track defined.");
463 if(fkESD && TMath::Abs(fkESD->GetTOFbc())>1){
464 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
468 if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
469 AliWarning("No output container defined.");
473 const Int_t ndim(kNdim+7);
475 alpha(0.), cs(-2.), sn(0.);
476 Float_t sz[AliTRDseedV1::kNtb], pos[AliTRDseedV1::kNtb];
477 Int_t ntbGap[AliTRDseedV1::kNtb];
478 AliTRDseedV1 *fTracklet(NULL);
479 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
480 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
481 if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()) continue;
484 alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
485 cs = TMath::Cos(alpha);
486 sn = TMath::Sin(alpha);
488 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
489 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());//fTracklet->GetTgl();
490 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
492 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fkTrack->Charge();// fSpecies;
493 val[kPt] = fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fTracklet->GetMomentum());
494 Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)),
495 dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)),
496 dydx(fTracklet->GetYfit(1)),
497 tilt(fTracklet->GetTilt());
498 // correct for tilt rotation
499 val[kYrez] = dyt - dzt*tilt;
500 val[kZrez] = dzt + dyt*tilt;
501 dydx+= tilt*fTracklet->GetZref(1);
502 val[kPrez] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
503 if(fTracklet->IsRowCross()) val[kSpeciesChgRC]= 0.;
504 val[kNdim] = fTracklet->GetdQdl()*2.e-3;
505 val[kNdim+1] = fEvent?fEvent->GetMultiplicity():0;
506 val[kNdim+2] = 1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb;
507 Int_t n = fTracklet->GetChargeGaps(sz, pos, ntbGap);
508 val[kNdim+3] = 0.; for(Int_t igap(0); igap<n; igap++) val[kNdim+3] += sz[igap];
509 for(Int_t ifill(0); ifill<3; ifill++){ val[kNdim+4+ifill]=0.;val[kNdim+5+ifill]=0.;}
510 for(Int_t igap(0), ifill(0); igap<n&&ifill<2; igap++){
511 if(ntbGap[igap]<2) continue;
512 val[kNdim+4+ifill] = sz[igap];
513 val[kNdim+5+ifill] = pos[igap];
518 // // compute covariance matrix
519 // fTracklet->GetCovAt(x, cov);
520 // fTracklet->GetCovRef(covR);
521 // cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
522 // Double_t dyz[2]= {dy[1], dz[1]};
523 // Pulls(dyz, cov, tilt);
524 // ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
525 // ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
528 Bool_t rc(fTracklet->IsRowCross());
529 UChar_t err(fTracklet->GetErrorMsg());
530 Double_t x(fTracklet->GetX()),
531 pt(fTracklet->GetPt()),
532 yt(fTracklet->GetYref(0)),
533 zt(fTracklet->GetZref(0)),
534 phi(fTracklet->GetYref(1)),
535 tht(fTracklet->GetZref(1));
536 Int_t ncl(fTracklet->GetN()), det(fTracklet->GetDetector());
537 (*DebugStream()) << "tracklet"
548 <<"dy=" << val[kYrez]
549 <<"dz=" << val[kZrez]
550 <<"dphi="<< val[kPrez]
551 <<"dQ ="<< val[kNdim]
557 if(!track) return NULL;
558 // special care for EVE usage
560 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
561 return H->Projection(kYrez);
565 //________________________________________________________
566 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
568 // Store resolution/pulls of Kalman before updating with the TRD information
569 // at the radial position of the first tracklet. The following points are used
571 // - the (y,z,snp) of the first TRD tracklet
572 // - the (y, z, snp, tgl, pt) of the MC track reference
574 // Additionally the momentum resolution/pulls are calculated for usage in the
576 //printf("AliTRDresolution::PlotTrackIn() :: track[%p]\n", (void*)track);
578 if(track) fkTrack = track;
580 AliDebug(4, "No Track defined.");
584 TH1 *h(NULL); // EVE projection
586 THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
588 AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
591 // check input track status
592 AliExternalTrackParam *tin(NULL);
593 if(!(tin = fkTrack->GetTrackIn())){
594 AliError("Track did not entered TRD fiducial volume.");
597 // check first tracklet
598 AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
600 AliDebug(3, "No Tracklet in ly[0]. Skip track.");
601 if(!track) return NULL;
602 // special care for EVE usage
603 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
604 return H->Projection(kYrez);
606 if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()){
607 AliDebug(3, "Tracklet or Chamber not OK. Skip track.");
608 if(!track) return NULL;
609 // special care for EVE usage
610 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
611 return H->Projection(kYrez);
613 // check radial position
614 Double_t x = tin->GetX();
615 if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
616 AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
617 if(!track) return NULL;
618 // special care for EVE usage
619 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
620 return H->Projection(kYrez);
622 //printf("USE y[%+f] dydx[%+f]\n", fTracklet->GetYfit(0), fTracklet->GetYfit(1));
624 Int_t bc(fkESD?fkESD->GetTOFbc()/2:0);
625 const Double_t *parR(tin->GetParameter());
626 Double_t dyt(fTracklet->GetYfit(0)-parR[0]), dzt(fTracklet->GetZfit(0)-parR[1]),
627 phit(fTracklet->GetYfit(1)),
628 tilt(fTracklet->GetTilt()),
629 norm(1./TMath::Sqrt((1.-parR[2])*(1.+parR[2])));
631 // correct for tilt rotation
632 Double_t dy = dyt - dzt*tilt,
634 dx = dy/(parR[2]*norm-parR[3]*norm*tilt);
635 phit += tilt*parR[3];
636 Double_t dphi = TMath::ATan(phit) - TMath::ASin(parR[2]);
638 Double_t val[kNdim+3];
639 val[kBC] = bc==0?0:(bc<0?-1.:1.);
640 Double_t alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha(),
641 cs = TMath::Cos(alpha),
642 sn = TMath::Sin(alpha);
643 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
644 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
645 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
647 val[kZrez] = fTracklet->IsRowCross()?dz:(fTracklet->GetdQdl()*5.e-4 - 2.5);
648 val[kPrez] = dphi*TMath::RadToDeg();
649 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fSpecies;
650 val[kPt] = fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fPt);
651 val[kNdim] = fTracklet->GetDetector();
653 val[kNdim+2] = fEvent?fEvent->GetBunchFill():0;
656 (*DebugStream()) << "trackIn"
657 <<"tracklet.=" << fTracklet
662 if(!track) return NULL;
663 // special care for EVE usage
664 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
665 return H->Projection(kYrez);
669 //________________________________________________________
670 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
672 // Store resolution/pulls of Kalman after last update with the TRD information
673 // at the radial position of the first tracklet. The following points are used
675 // - the (y,z,snp) of the first TRD tracklet
676 // - the (y, z, snp, tgl, pt) of the MC track reference
678 // Additionally the momentum resolution/pulls are calculated for usage in the
681 if(track) fkTrack = track;
685 //________________________________________________________
686 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
689 // Plot MC distributions
692 if(!HasMCdata()) return NULL;
693 if(track) fkTrack = track;
695 AliDebug(4, "No Track defined.");
698 Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
702 AliWarning("No output container defined.");
705 // retriev track characteristics
706 Int_t pdg = fkMC->GetPDG(),
707 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
710 label(fkMC->GetLabel());
712 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
713 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
714 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
717 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
718 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
720 Double_t x, y, z, pt, dydx, dzdx, dzdl;
721 Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
722 Double_t covR[7]/*, cov[3]*/;
724 /* if(DebugLevel()>=3){
725 // get first detector
727 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
728 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
729 det = fTracklet->GetDetector();
733 TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
735 m = fkTrack->GetMass();
736 if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
737 (*DebugStream()) << "MCkalman"
754 AliTRDcluster *c(NULL);
755 Double_t val[kNdim+1];
756 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
757 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
758 !fTracklet->IsOK())*/ continue;
760 x= x0 = fTracklet->GetX();
761 Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
762 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, phi, s)) continue;
764 // MC track position at reference radial position
766 Float_t ymc = y0 - dx*dydx0;
767 Float_t zmc = z0 - dx*dzdx0;
768 //phi -= TMath::Pi();
773 val[kSpeciesChgRC]= rc?0.:sign*sIdx;
774 val[kPt] = pt0<0.8?0:(pt0<1.5?1:2);//GetPtBin(pt0);
775 Double_t tilt(fTracklet->GetTilt());
777 // ,corr(1./(1. + t2))
778 // ,cost(TMath::Sqrt(corr));
780 AliExternalTrackParam *tin(fkTrack->GetTrackIn());
781 if(ily==0 && tin){ // trackIn residuals
782 // check radial position
783 if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
785 val[kBC] = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
786 val[kYrez] = tin->GetY()-ymc;
787 val[kZrez] = tin->GetZ()-zmc;
788 val[kPrez] = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
789 if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
792 //if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
795 dydx = fTracklet->GetYref(1);
796 dzdx = fTracklet->GetZref(1);
797 dzdl = fTracklet->GetTgl();
798 y = fTracklet->GetYref(0);
800 z = fTracklet->GetZref(0);
802 pt = TMath::Abs(fTracklet->GetPt());
803 fTracklet->GetCovRef(covR);
806 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
808 val[kNdim] = 1.e2*(pt/pt0-1.);
809 if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
810 /* // theta resolution/ tgl pulls
811 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
812 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
813 ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
814 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
815 // pt resolution \\ 1/pt pulls \\ p resolution for PID
816 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
817 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
818 ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
819 ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
820 ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
822 // Fill Debug stream for MC track
824 Int_t det(fTracklet->GetDetector());
825 (*DebugStream()) << "MC"
837 // Fill Debug stream for Kalman track
838 (*DebugStream()) << "MCtrack"
850 // tracklet residuals
851 dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
852 dzdx = fTracklet->GetZfit(1);
853 y = fTracklet->GetYfit(0);
855 z = fTracklet->GetZfit(0);
857 val[kYrez] = dy - dz*tilt;
858 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
859 val[kZrez] = dz + dy*tilt;
860 // val[kNdim] = pt/pt0-1.;
861 if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
864 // Fill Debug stream for tracklet
866 Float_t s2y = fTracklet->GetS2Y();
867 Float_t s2z = fTracklet->GetS2Z();
868 (*DebugStream()) << "MCtracklet"
879 AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
880 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
881 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
883 H = (THnSparse*)fContainer->At(kMCcluster);
884 val[kPt] = TMath::ATan(dydx0)*TMath::RadToDeg();
885 //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
887 Float_t padCorr(tilt*fTracklet->GetPadLength());
888 fTracklet->ResetClusterIter(kTRUE);
889 while((c = fTracklet->NextCluster())){
890 if(row0<0) row0 = c->GetPadRow();
891 x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
892 y = c->GetY() + padCorr*(c->GetPadRow() - row0);
899 val[kYrez] = dy - dz*tilt;
901 val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t q(TMath::Abs(c->GetQ()));
902 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
903 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
904 if(ic) val[kZrez] /= (ic*q);
908 // Fill calibration container
909 Float_t d = zr0 - zmc;
910 d -= ((Int_t)(2 * d)) / 2.0;
911 if (d > 0.25) d = 0.5 - d;
912 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
913 clInfo->SetCluster(c);
914 clInfo->SetMC(pdg, label);
915 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
916 clInfo->SetResolution(dy);
917 clInfo->SetAnisochronity(d);
918 clInfo->SetDriftLength(dx);
919 clInfo->SetTilt(tilt);
920 if(fMCcl) fMCcl->Add(clInfo);
921 else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
924 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
925 clInfoArr->SetOwner(kFALSE);
927 clInfoArr->Add(clInfo);
931 if(DebugLevel()>=5 && clInfoArr){
932 (*DebugStream()) << "MCcluster"
933 <<"clInfo.=" << clInfoArr
938 if(clInfoArr) delete clInfoArr;
943 //__________________________________________________________________________
944 Int_t AliTRDresolution::GetPtBin(Float_t pt)
946 // Find pt bin according to local pt segmentation
948 while(ipt<AliTRDresolution::kNpt){
949 if(pt<fgPtBin[ipt+1]) break;
955 //________________________________________________________
956 Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)
958 // return mean number of entries/bin of histogram "h"
959 // if option "opt" is given the following values are accepted:
960 // "<" : consider only entries less than "cut"
961 // ">" : consider only entries greater than "cut"
963 //Int_t dim(h->GetDimension());
964 Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());
965 Double_t sum(0.); Int_t n(0);
966 for(Int_t ix(1); ix<=nbx; ix++)
967 for(Int_t iy(1); iy<=nby; iy++)
968 for(Int_t iz(1); iz<=nbz; iz++){
969 if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}
971 if(strcmp(opt, "<")==0) {
972 if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
973 } else if(strcmp(opt, ">")==0){
974 if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
975 } else {sum += h->GetBinContent(ix, iy, iz); n++;}
981 //________________________________________________________
982 void AliTRDresolution::GetRangeZ(TH2 *h2, Float_t &min, Float_t &max)
984 // Get range on Z axis such to avoid outliers
986 Double_t cnt[20000], c, m, s;
987 Int_t nx(h2->GetXaxis()->GetNbins()), ny(h2->GetYaxis()->GetNbins()), nc(0);
988 for(Int_t ix(1); ix<=nx; ix++){
989 for(Int_t iy(1); iy<=ny; iy++){
990 if((c = h2->GetBinContent(ix, iy))<10) continue;
996 AliMathBase::EvaluateUni(nc, cnt, m, s, 0);
997 min = m-s; max = m+2.*s;
1000 //________________________________________________________
1001 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
1004 // Get the reference figures
1008 AliWarning("Please provide a canvas to draw results.");
1011 /* Int_t selection[100], n(0), selStart(0); //
1012 Int_t ly0(0), dly(5);
1013 TList *l(NULL); TVirtualPad *pad(NULL); */
1018 AliWarning(Form("Reference plot [%d] missing result", ifig));
1023 //________________________________________________________
1024 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
1026 // Build pt segments
1027 for(Int_t j(0); j<=kNpt; j++){
1028 pt0+=(TMath::Exp(j*j*dpt)-1.);
1033 //________________________________________________________
1034 void AliTRDresolution::MakeSummary()
1036 // Build summary plots
1039 AliError("Missing results");
1042 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
1043 TObjArray *arr(NULL); TH2 *h2(NULL);
1045 // cluster resolution
1047 gStyle->SetPalette(1);
1048 const Int_t nClViews(9);
1049 const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh"};
1050 const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1};
1051 const Int_t nTrkltViews(19);
1052 const Char_t *vTrkltName[nTrkltViews] = {
1053 "TrkltY", "TrkltYn", "TrkltYp", "TrkltY", "TrkltYn", "TrkltYp",
1054 "TrkltPh", "TrkltPhn", "TrkltPhp",
1055 "TrkltQ", "TrkltQn", "TrkltQp",
1056 "TrkltQS", "TrkltQSn", "TrkltQSp",
1057 "TrkltTBn", "TrkltTBp", "TrkltTBn", "TrkltTBp"
1058 // "TrkltYnl0", "TrkltYpl0", "TrkltPhnl0", "TrkltPhpl0", "TrkltQnl0", "TrkltQpl0", // electrons low pt
1059 /* "TrkltYnl1", "TrkltYpl1", "TrkltPhnl1", "TrkltPhpl1", "TrkltQnl1", "TrkltQpl1", // muons low pt
1060 "TrkltYnl2", "TrkltYpl2", "TrkltPhnl2", "TrkltPhpl2", "TrkltQnl2", "TrkltQpl2" // pions low pt*/
1062 const UChar_t vTrkltOpt[nTrkltViews] = {
1069 const Int_t nTrkInViews(5);
1070 const Char_t *vTrkInName[nTrkInViews][6] = {
1071 {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1072 {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1073 {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"},
1074 {"TrkInXnl", "TrkInXpl", "TrkInXl", "TrkInYnh", "TrkInYph", "TrkInYh"},
1075 {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"},
1076 //{"TrkInRCX", "TrkInRCY", "TrkInRCPh", "TrkInRCZl", "TrkInRCZi", "TrkInRCZh"}
1078 const UChar_t vTrkInOpt[nTrkInViews] = {0, 1, 0, 0, 0};
1079 const Float_t min[6] = {0.15, 0.15, 0.15, 0.15, 0.5, 0.5};
1080 const Float_t max[6] = {0.6, 0.6, 0.6, 0.6, 2.3, 2.3};
1081 const Char_t *ttt[6] = {"#sigma(#Deltay) [cm]", "#sigma(#Deltay) [cm]", "#sigma(#Deltay) [cm]", "#sigma(#Deltaz) [cm]", "#sigma(#Delta#phi) [deg]", "#sigma(#Delta#phi) [deg]"};
1083 const Int_t nTrkViews(27);
1084 const Char_t *vTrkName[nTrkViews] = {
1085 "TrkY", "TrkYn", "TrkYp",
1086 "TrkPh", "TrkPhn", "TrkPhp",
1087 "TrkDPt", "TrkDPtn", "TrkDPtp",
1088 "TrkYnl0", "TrkYpl0", "TrkPhnl0", "TrkPhpl0", "TrkDPtnl0", "TrkDPtpl0", // electrons low pt
1089 "TrkYnl1", "TrkYpl1", "TrkPhnl1", "TrkPhpl1", "TrkDPtnl1", "TrkDPtpl1", // muons low pt
1090 "TrkYnl2", "TrkYpl2", "TrkPhnl2", "TrkPhpl2", "TrkDPtnl2", "TrkDPtpl2" // pions low pt
1092 const Char_t *typName[] = {"", "MC"};
1093 const Int_t nx(2048), ny(1536);
1095 if((arr = (TObjArray*)fProj->At(kDetector))){
1096 cOut = new TCanvas(Form("%s_Det", GetName()), "Detector performance", 2*nx, 2*ny);
1097 cOut->Divide(AliTRDgeometry::kNlayer,AliTRDeventInfo::kCentralityClasses, 1.e-5, 1.e-5);
1098 for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1099 Float_t zmin(1.e10), zmax(0.);
1100 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1101 if(!(h2 = (TH2*)arr->FindObject(Form("HDetQ%d%d_yx", ily, icen)))) continue;
1102 Float_t m, M; GetRangeZ(h2, m, M);
1106 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1107 p=cOut->cd(icen*AliTRDgeometry::kNlayer+ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1108 if(!(h2 = (TH2*)arr->FindObject(Form("HDetQ%d%d_yx", ily, icen)))) continue;
1109 SetRangeZ(h2, zmin, zmax);
1111 MakeDetectorPlot(ily);
1114 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1116 for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
1117 if((arr = (TObjArray*)fProj->At(ityp?kMCcluster:kCluster))){
1118 for(Int_t iview(0); iview<nClViews; iview++){
1119 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vClName[iview], vClOpt[iview]), "Cluster Resolution", nx, ny);
1120 cOut->Divide(3,2, 1.e-5, 1.e-5);
1122 for(Int_t iplot(0); iplot<6; iplot++){
1123 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1124 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], iplot)))) continue;
1126 if(vClOpt[iview]==0) h2->Draw("colz");
1127 else if(vClOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [#mum]", 3.2e2, 5.e2, 1.e4);
1128 MakeDetectorPlot(iplot);
1130 if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1134 // tracklet systematic
1135 if((arr = (TObjArray*)fProj->At(ityp?kMCtracklet:kTracklet))){
1136 for(Int_t iview(0); iview<nTrkltViews; iview++){
1137 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkltName[iview], vTrkltOpt[iview]), "Tracklet Resolution", nx, ny);
1138 cOut->Divide(3,2, 1.e-5, 1.e-5);
1140 for(Int_t iplot(0); iplot<6; iplot++){
1141 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1142 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
1144 if(vTrkltOpt[iview]==0) h2->Draw("colz");
1145 else DrawSigma(h2, "#sigma(#Deltay) [cm]", .15, .6);
1146 MakeDetectorPlot(iplot);
1149 cOut->Modified();cOut->Update();
1150 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1155 // trackIn systematic
1156 if((arr = (TObjArray*)fProj->At(ityp?kMCtrackIn:kTrackIn))){
1157 for(Int_t iview(0); iview<nTrkInViews; iview++){
1158 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkInName[iview][0], vTrkInOpt[iview]), "Track IN Resolution", nx, ny);
1159 cOut->Divide(3,2, 1.e-5, 1.e-5);
1161 for(Int_t iplot(0); iplot<6; iplot++){
1162 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1163 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])))){
1164 AliInfo(Form("Missing H%s%s_2D", typName[ityp], vTrkInName[iview][iplot]));
1168 if(vTrkInOpt[iview]==0) h2->Draw("colz");
1169 else DrawSigma(h2, ttt[iplot], min[iplot], max[iplot]);
1170 MakeDetectorPlot(0);
1172 if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1176 for(Int_t iview(0); iview<2; iview++){
1177 const Char_t *plot = iview?"Ph":"Y";
1178 cOut = new TCanvas(Form("%s_%sSpec%s", GetName(), typName[ityp], plot), "Track IN Resolution", Int_t(1.5*nx), Int_t(1.5*ny));
1179 cOut->Divide(5,3, 1.e-5, 1.e-5);
1180 Int_t nplot(0); const Char_t *chName[] = {"p", "n", ""};
1181 for(Int_t iplot(0); iplot<3; iplot++){
1182 for(Int_t ispec(0); ispec<AliPID::kSPECIES; ispec++){
1183 p=cOut->cd(iplot*AliPID::kSPECIES+ispec+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1184 if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkIn%s%s%d_2D", typName[ityp], plot, chName[iplot], ispec)))) {
1185 AliInfo(Form("Missing H%sTrkIn%s%s%d_2D", typName[ityp], plot, chName[iplot], ispec));
1190 MakeDetectorPlot(0);
1193 if(nplot==15) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1198 // track MC systematic
1199 if((arr = (TObjArray*)fProj->At(kMCtrack))) {
1200 for(Int_t iview(0); iview<nTrkViews; iview++){
1201 cOut = new TCanvas(Form("%s_MC%s", GetName(), vTrkName[iview]), "Track Resolution", nx, ny);
1202 cOut->Divide(3,2, 1.e-5, 1.e-5);
1204 for(Int_t iplot(0); iplot<6; iplot++){
1205 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1206 if(!(h2 = (TH2*)arr->FindObject(Form("HMC%s%d_2D", vTrkName[iview], iplot)))) continue;
1207 h2->Draw("colz"); nplot++;
1209 if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1215 gStyle->SetPalette(1);
1218 //________________________________________________________
1219 void AliTRDresolution::DrawSigma(TH2 *h2, const Char_t *title, Float_t m, Float_t M, Float_t scale)
1221 // Draw error bars scaled with "scale" instead of content values
1222 //use range [m,M] if limits are specified
1225 TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1226 TH2F *h2e = new TH2F(Form("%s_E", h2->GetName()),
1227 Form("%s;%s;%s;%s", h2->GetTitle(), ax->GetTitle(), ay->GetTitle(), title),
1228 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
1229 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1231 TAxis *az(h2e->GetZaxis());
1232 if(M>m) az->SetRangeUser(m, M);
1234 az->SetTitleOffset(1.5);
1235 for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1236 for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1237 if(h2->GetBinContent(ix, iy)<-100.) continue;
1238 Float_t v(scale*h2->GetBinError(ix, iy));
1239 if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1240 h2e->SetBinContent(ix, iy, v);
1246 //________________________________________________________
1247 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1249 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1250 // The "range" vector should be initialized with 2 elements
1251 // Option "mod" can be any of
1252 // - 0 : gaussian like distribution
1253 // - 1 : tailed distribution
1255 Int_t nx(h2->GetNbinsX())
1256 , ny(h2->GetNbinsY())
1258 Double_t *data=new Double_t[n];
1259 for(Int_t ix(1), in(0); ix<=nx; ix++){
1260 for(Int_t iy(1); iy<=ny; iy++)
1261 data[in++] = h2->GetBinContent(ix, iy);
1263 Double_t mean, sigm;
1264 AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1266 range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1267 if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1268 AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1269 TH1S h1("h1SF0", "", 100, range[0], range[1]);
1274 case 0:// gaussian distribution
1276 TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1278 mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1279 range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1280 AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1]));
1283 case 1:// tailed distribution
1285 Int_t bmax(h1.GetMaximumBin());
1286 Int_t jBinMin(1), jBinMax(100);
1287 for(Int_t ibin(bmax); ibin--;){
1288 if(h1.GetBinContent(ibin)<1.){
1289 jBinMin=ibin; break;
1292 for(Int_t ibin(bmax); ibin++;){
1293 if(h1.GetBinContent(ibin)<1.){
1294 jBinMax=ibin; break;
1297 range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1298 AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1]));
1306 //________________________________________________________
1307 Bool_t AliTRDresolution::MakeProjectionDetector()
1310 const Int_t kNcontours(9);
1311 const Int_t kNstat(100);
1312 if(fProj && fProj->At(kDetector)) return kTRUE;
1314 AliError("Missing data container.");
1318 if(!(H = (THnSparse*)fContainer->FindObject("hDet2Cluster"))){
1319 AliInfo(Form("Missing/Wrong data @ hDet2Cluster."));
1322 Int_t ndim(H->GetNdimensions());
1323 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1324 TAxis *aa[kNdim], *ac(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1325 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1329 nCen = AliTRDeventInfo::kCentralityClasses;
1331 // build list of projections
1332 const Int_t nsel(AliTRDgeometry::kNlayer*AliTRDeventInfo::kCentralityClasses);
1333 // define rebinning strategy
1334 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 2, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1335 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1336 AliTRDresolutionProjection hp[kDetNproj]; TObjArray php(kDetNproj);
1337 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1338 for(Int_t icen(0); icen<nCen; icen++){
1339 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1340 isel++; // new selection
1341 hp[ih].Build(Form("HDetQ%d%d", ily, icen),
1342 Form("Detectors :: Q Ly[%d] Cen[%s]", ily, cenName[icen]),
1343 kEta, kPhi, kYrez, aa);
1344 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1345 hp[ih].SetShowRange(10., 55.);
1346 php.AddLast(&hp[ih++]); np[isel]++;
1349 AliInfo(Form("Build %3d 3D projections.", ih));
1351 Int_t ly(0), cen(0);
1352 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1353 v = H->GetBinContent(ib, coord); if(v<1.) continue;
1356 // centrality selection
1358 if(ac) cen = coord[4]-1;
1360 isel = cen*AliTRDgeometry::kNlayer+ly; Int_t ioff=isel;
1361 // AliDebug(4, Form("SELECTION[%d] :: ch[%c] pt[%c] sp[%d] ly[%d]\n", np[isel], ch==2?'Z':chName[ch], ptName[pt], sp, ly));
1362 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1364 TObjArray *arr(NULL);
1365 fProj->AddAt(arr = new TObjArray(kDetNproj), kDetector);
1367 TH2 *h2(NULL); Int_t jh(0);
1369 if(!hp[ih].fH) continue;
1370 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1371 arr->AddAt(h2, jh++);
1372 if(!(h2 = (TH2*)gDirectory->Get(Form("%s_yx", hp[ih].fH->GetName())))) continue;
1373 arr->AddAt(h2, jh++);
1375 AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1376 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1377 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HDetQ%d%d", ily, 0)))){
1378 for(Int_t icen(1); icen<nCen; icen++){
1379 if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HDetQ%d%d", ily, icen)))){
1383 pr0->fH->SetNameTitle(Form("HDetQ%d", ily), Form("Detectors :: Q Ly[%d]", ily));
1384 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1387 AliInfo(Form("Done %3d 2D projections.", jh));
1391 //________________________________________________________
1392 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1395 const Int_t kNcontours(9);
1396 const Int_t kNstat(100);
1397 Int_t cidx=mc?kMCcluster:kCluster;
1398 if(fProj && fProj->At(cidx)) return kTRUE;
1400 AliError("Missing data container.");
1403 const Char_t *projName[] = {"hCluster2Track", "hCluster2MC"};
1405 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1406 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1409 Int_t ndim(H->GetNdimensions());
1410 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1411 TAxis *aa[kNdim], *as(NULL), *apt(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1412 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1413 if(ndim > kPt) apt = H->GetAxis(kPt);
1414 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1415 // build list of projections
1416 const Int_t nsel(12);
1417 // define rebinning strategy
1418 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1419 AliTRDresolutionProjection hp[kClNproj]; TObjArray php(kClNproj);
1420 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1421 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1422 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1423 for(Int_t ich(0); ich<kNcharge; ich++){
1424 isel++; // new selection
1425 hp[ih].Build(Form("H%sClY%c%d", mc?"MC":"", chName[ich], ily),
1426 Form("Clusters[%c] :: #Deltay Ly[%d]", chSgn[ich], ily),
1427 kEta, kPhi, kYrez, aa);
1428 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1429 php.AddLast(&hp[ih++]); np[isel]++;
1430 hp[ih].Build(Form("H%sClQ%c%d", mc?"MC":"", chName[ich], ily),
1431 Form("Clusters[%c] :: Q Ly[%d]", chSgn[ich], ily),
1432 kEta, kPhi, kSpeciesChgRC, aa);
1433 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1434 hp[ih].SetShowRange(20., 40.);
1435 php.AddLast(&hp[ih++]); np[isel]++;
1436 hp[ih].Build(Form("H%sClYXTC%c%d", mc?"MC":"", chName[ich], ily),
1437 Form("Clusters[%c] :: #Deltay(x,TC) Ly[%d]", chSgn[ich], ily),
1438 kPrez, kZrez, kYrez, aa);
1439 php.AddLast(&hp[ih++]); np[isel]++;
1440 hp[ih].Build(Form("H%sClYXPh%c%d", mc?"MC":"", chName[ich], ily),
1441 Form("Clusters[%c] :: #Deltay(x,#Phi) Ly[%d]", chSgn[ich], ily),
1442 kPrez, kPt, kYrez, aa);
1443 php.AddLast(&hp[ih++]); np[isel]++;
1446 AliInfo(Form("Build %3d 3D projections.", ih));
1448 Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1);
1449 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1450 v = H->GetBinContent(ib, coord); if(v<1.) continue;
1453 if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1456 ch = 0; // [-] track
1457 if(chBin>0 && coord[kPt] > chBin) ch = 1; // [+] track
1459 isel = ly*2+ch; Int_t ioff=isel*4;
1460 // AliDebug(4, Form("SELECTION[%d] :: ch[%c] pt[%c] sp[%d] ly[%d]\n", np[isel], ch==2?'Z':chName[ch], ptName[pt], sp, ly));
1461 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1463 TObjArray *arr(NULL);
1464 fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1466 TH2 *h2(NULL); Int_t jh(0);
1468 if(!hp[ih].fH) continue;
1469 //if(hp[ih].fH->GetEntries()<100) continue;
1470 Int_t mid(1), nstat(kNstat);
1471 if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=300;}
1472 if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1473 arr->AddAt(h2, jh++);
1475 AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1476 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1477 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d", mc?"MC":"", chName[0], ily)))){
1478 if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d", mc?"MC":"", chName[1], ily)))){
1480 pr0->fH->SetNameTitle(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1481 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1484 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d", mc?"MC":"", chName[0], ily)))){
1485 if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d", mc?"MC":"", chName[1], ily)))){
1487 pr0->fH->SetNameTitle(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: #Deltay(x,#Phi) Ly[%d]", ily));
1488 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1492 AliInfo(Form("Done %3d 2D projections.", jh));
1496 //________________________________________________________
1497 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1500 const Int_t kNcontours(9);
1501 const Int_t kNstat(30);
1502 const Int_t kNstatQ(30);
1503 Int_t cidx=mc?kMCtracklet:kTracklet;
1504 if(fProj && fProj->At(cidx)) return kTRUE;
1506 AliError("Missing data container.");
1509 const Char_t *projName[] = {"hTracklet2Track", "hTracklet2MC"};
1511 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1512 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1515 const Int_t mdim(kNdim+8);
1516 Int_t ndim(H->GetNdimensions());
1517 Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
1518 TAxis *aa[mdim], *as(NULL), *ap(NULL), *ac(NULL); memset(aa, 0, sizeof(TAxis*) * mdim);
1519 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1520 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC); // init species/charge selection
1521 if(ndim > kPt) ap = H->GetAxis(kPt); // init pt selection
1522 if(ndim > kNdim+1) ac = H->GetAxis(kNdim+1); // init centrality selection
1523 // calculate size depending on debug level
1524 const Int_t nCen(ndim>kNdimTrklt?AliTRDeventInfo::kCentralityClasses:1);
1525 const Int_t nPt(ndim>kNdimTrklt?kNpt:1);
1526 const Int_t nSpc(1);//ndim>kNdimTrklt?fgkNbins[kSpeciesChgRC]:1);
1527 const Int_t nCh(ndim>kNdimTrklt?kNcharge:1);
1529 // build list of projections
1530 const Int_t nsel(AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer*kNpt*(AliPID::kSPECIES*kNcharge + 1));
1531 // define rebinning strategy
1532 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1533 AliTRDresolutionProjection hp[kTrkltNproj]; TObjArray php(kTrkltNproj);
1534 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1535 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1536 const Char_t ptName[kNpt] = {'l', 'i', 'h'};
1537 const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"};
1538 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1539 for(Int_t icen(0); icen<nCen; icen++){
1540 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1541 for(Int_t ipt(0); ipt<nPt; ipt++){
1542 for(Int_t isp(0); isp<nSpc; isp++){
1543 for(Int_t ich(0); ich<nCh; ich++){
1544 isel++; // new selection
1545 hp[ih].Build(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1546 Form("Tracklets[%s%c]:: #Deltay{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1547 kEta, kPhi, kYrez, aa);
1548 //hp[ih].SetShowRange(-0.1,0.1);
1549 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1550 php.AddLast(&hp[ih++]); np[isel]++;
1551 hp[ih].Build(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1552 Form("Tracklets[%s%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1553 kEta, kPhi, kPrez, aa);
1554 //hp[ih].SetShowRange(-0.5,0.5);
1555 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1556 php.AddLast(&hp[ih++]); np[isel]++;
1557 hp[ih].Build(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1558 Form("Tracklets[%s%c]:: dQdl{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1559 kEta, kPhi, kNdim, aa);
1560 hp[ih].SetShowRange(1.,2.3);
1561 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1562 php.AddLast(&hp[ih++]); np[isel]++;
1563 hp[ih].Build(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1564 Form("Tracklets[%s%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1565 kEta, kPhi, kNdim+2, aa);
1566 hp[ih].SetShowRange(30., 70.);
1567 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1568 php.AddLast(&hp[ih++]); np[isel]++;
1571 if(ndim==kNdimTrklt) continue;
1573 isel++; // new selection
1574 hp[ih].Build(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1575 Form("Tracklets[RC]:: #Deltaz{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1576 kEta, kPhi, kZrez, aa);
1577 // hp[ih].SetShowRange(-0.1,0.1);
1578 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1579 php.AddLast(&hp[ih++]); np[isel]++;
1580 hp[ih].Build(Form("H%sTrkltRCY%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1581 Form("Tracklets[RC]:: #Deltay{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1582 kEta, kPhi, kYrez, aa);
1583 //hp[ih].SetShowRange(-0.1,0.1);
1584 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1585 php.AddLast(&hp[ih++]); np[isel]++;
1586 hp[ih].Build(Form("H%sTrkltRCPh%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1587 Form("Tracklets[RC]:: #Delta#phi{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1588 kEta, kPhi, kPrez, aa);
1589 //hp[ih].SetShowRange(-0.1,0.1);
1590 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1591 php.AddLast(&hp[ih++]); np[isel]++;
1592 hp[ih].Build(Form("H%sTrkltRCQ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1593 Form("Tracklets[RC]:: dQdl{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1594 kEta, kPhi, kNdim, aa);
1595 //hp[ih].SetShowRange(-0.1,0.1);
1596 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1597 php.AddLast(&hp[ih++]); np[isel]++;
1601 AliInfo(Form("Build %3d 3D projections.", ih));
1603 Int_t ly(0), ch(0), sp(2), rcBin(as?as->FindBin(0.):-1), pt(0), cen(0), ioff(0), ispc(0);
1604 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1605 v = H->GetBinContent(ib, coord);
1607 ly = coord[kBC]-1; // layer selection
1609 ch = 0; sp=0;// [e-] track [dafault]
1610 if(rcBin>0){ // debug mode in which species are also saved
1611 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1612 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
1613 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
1617 if(ap) pt = TMath::Min(coord[kPt]-1, Int_t(kNpt)-1);
1618 // centrality selection
1620 if(ac) cen = coord[kNdim+1]-1;
1622 if(ndim==kNdimTrklt){
1627 isel = cen*AliTRDgeometry::kNlayer*kNpt*ispc+ly*kNpt*ispc+pt*ispc; isel+=sp<0?(nSpc*nCh):(sp*nCh+ch);
1630 AliDebug(4, Form("SELECTION[%d] :: ch[%c] pt[%c] sp[%d] ly[%d] cen[%d]", np[isel], ch==2?'Z':chName[ch], ptName[pt], sp, ly, cen));
1631 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1633 TObjArray *arr(NULL);
1634 fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1636 TH2 *h2(NULL); Int_t jh(0);
1638 if(!hp[ih].fH) continue;
1639 Int_t mid(0), nstat(kNstat);
1640 if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=kNstatQ;}
1641 if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1642 arr->AddAt(h2, jh++);
1644 // build combined performance plots
1645 AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1646 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1647 for(Int_t ich(0); ich<nCh; ich++){
1648 for(Int_t icen(0); icen<nCen; icen++){
1649 for(Int_t ipt(0); ipt<nPt; ipt++){
1651 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1652 for(Int_t isp(1); isp<nSpc; isp++){
1653 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1656 pr0->fH->SetNameTitle(Form("H%sTrkltY%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1657 Form("Tracklets[%c]:: #Deltay{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1658 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1659 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1662 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1663 for(Int_t isp(1); isp<nSpc; isp++){
1664 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1667 pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1668 Form("Tracklets[%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1669 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1670 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1673 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1674 for(Int_t isp(1); isp<nSpc; isp++){
1675 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1678 pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1679 Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1680 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1681 pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1682 Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1683 pr0->SetShowRange(2.4, 5.1);
1684 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1685 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1688 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1689 for(Int_t isp(1); isp<nSpc; isp++){
1690 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1693 pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1694 Form("Tracklets[%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1695 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1696 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1698 } // end pt integration
1700 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1701 pr0->fH->SetNameTitle(Form("H%sTrkltY%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1702 Form("Tracklets[%c]:: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1703 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1704 if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1707 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1708 pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1709 Form("Tracklets[%c]:: #Delta#phi Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1710 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1711 if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1714 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1715 pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1716 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1717 pr0->SetShowRange(1.,2.3);
1718 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1719 pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1720 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1721 pr0->SetShowRange(2.4,5.1);
1722 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1723 if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1726 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1727 pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1728 Form("Tracklets[%c]:: OccupancyTB Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1729 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1730 if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1732 } // end centrality integration
1734 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1735 pr0->fH->SetNameTitle(Form("H%sTrkltY%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Deltay Ly[%d]", chSgn[ich], ily));
1736 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1737 if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1740 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1741 pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Delta#phi Ly[%d]", chSgn[ich], ily));
1742 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1743 if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1746 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1747 pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1748 pr0->SetShowRange(1.,2.3);
1749 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1750 pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1751 pr0->SetShowRange(2.4,5.1);
1752 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1753 if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1756 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1757 pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: OccupancyTB Ly[%d]", chSgn[ich], ily));
1758 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1759 if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1761 } // end charge integration
1763 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1764 pr0->fH->SetNameTitle(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets :: #Deltay Ly[%d]", ily));
1765 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1768 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1769 pr0->fH->SetNameTitle(Form("H%sTrkltPh%d", mc?"MC":"", ily), Form("Tracklets :: #Delta#phi Ly[%d]", ily));
1770 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1773 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1774 pr0->fH->SetNameTitle(Form("H%sTrkltQ%d", mc?"MC":"", ily), Form("Tracklets :: dQdl Ly[%d]", ily));
1775 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1776 pr0->fH->SetNameTitle(Form("H%sTrkltQS%d", mc?"MC":"", ily), Form("Tracklets :: dQdl Ly[%d]", ily));
1777 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1780 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1781 pr0->fH->SetNameTitle(Form("H%sTrkltTB%d", mc?"MC":"", ily), Form("Tracklets :: OccupancyTB Ly[%d]", ily));
1782 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1785 /*! Row Cross processing*/
1786 for(Int_t icen(0); icen<nCen; icen++){
1788 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, icen)))){
1789 for(Int_t ipt(0); ipt<kNpt; ipt++){
1790 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen)))) continue;
1793 pr0->fH->SetNameTitle(Form("H%sTrkltRCZ%d%d", mc?"MC":"", ily, icen), Form("Tracklets[RC]:: #Deltaz Ly[%d] Cen[%s]", ily, cenName[icen]));
1794 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1795 if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))) (*pr1)+=(*pr0);
1797 } // end centrality integration for row cross
1799 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))){
1800 pr0->fH->SetNameTitle(Form("H%sTrkltRCZ%d", mc?"MC":"", ily), Form("Tracklets[RC] :: #Deltaz Ly[%d]", ily));
1801 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1804 AliInfo(Form("Done %3d 2D projections.", jh));
1808 //________________________________________________________
1809 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
1812 const Int_t kNcontours(9);
1813 const Int_t kNstat(30);
1814 Int_t cidx=mc?kMCtrackIn:kTrackIn;
1815 if(fProj && fProj->At(cidx)) return kTRUE;
1817 AliError("Missing data container.");
1820 const Char_t *projName[] = {"hTracklet2TRDin", "hTRDin2MC"};
1822 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1823 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1827 const Int_t mdim(kNdim+3);
1828 Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
1829 Int_t ndim(H->GetNdimensions());
1830 TAxis *aa[mdim], *as(NULL), *ap(NULL), *abf(NULL); memset(aa, 0, sizeof(TAxis*) * (mdim));
1831 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1832 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1833 if(ndim > kPt) ap = H->GetAxis(kPt);
1834 if(ndim > (kNdim+2)) abf = H->GetAxis(kNdim+2);
1835 AliInfo(Form("Using : Species[%c] Pt[%c] BunchFill[%c]", as?'y':'n', ap?'y':'n', abf?'y':'n'));
1836 const Int_t nPt(ndim>kNdimTrkIn?kNpt:1);
1838 // build list of projections
1839 const Int_t nsel(kNpt*(AliPID::kSPECIES*kNcharge + 1));
1840 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1841 const Char_t ptName[kNpt] = {'l', 'i', 'h'};
1842 const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"};
1843 // define rebinning strategy
1844 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1845 AliTRDresolutionProjection hp[kMCTrkInNproj]; TObjArray php(kMCTrkInNproj+2);
1846 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1847 // define list of projections
1848 for(Int_t ipt(0); ipt<nPt; ipt++){
1849 for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
1850 for(Int_t ich(0); ich<kNcharge; ich++){
1851 isel++; // new selection
1852 hp[ih].Build(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
1853 Form("TrackIn[%s%c]:: #Deltay{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
1854 kEta, kPhi, kYrez, aa);
1855 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1856 php.AddLast(&hp[ih++]); np[isel]++;
1857 hp[ih].Build(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
1858 Form("TrackIn[%s%c]:: #Delta#phi{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
1859 kEta, kPhi, kPrez, aa);
1860 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1861 php.AddLast(&hp[ih++]); np[isel]++;
1862 hp[ih].Build(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
1863 Form("TrackIn[%s%c]:: #Deltax{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
1864 kEta, kPhi, kNdim+1, aa);
1865 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1866 php.AddLast(&hp[ih++]); np[isel]++;
1869 isel++; // RC tracks
1870 hp[ih].Build(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt]),
1871 Form("TrackIn[RC]:: #Deltaz{%s}", ptCut[ipt]),
1872 kEta, kPhi, kZrez, aa);
1873 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1874 php.AddLast(&hp[ih++]); np[isel]++;
1875 hp[ih].Build(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt]),
1876 Form("TrackIn[RC]:: #Deltay{%s}", ptCut[ipt]),
1877 kEta, kPhi, kYrez, aa);
1878 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1879 php.AddLast(&hp[ih++]); np[isel]++;
1880 hp[ih].Build(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt]),
1881 Form("TrackIn[RC]:: #Delta#phi{%s}", ptCut[ipt]),
1882 kEta, kPhi, kPrez, aa);
1883 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1884 php.AddLast(&hp[ih++]); np[isel]++;
1885 hp[ih].Build(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt]),
1886 Form("TrackIn[RC]:: #Deltax{%s}", ptCut[ipt]),
1887 kEta, kPhi, kNdim+1, aa);
1888 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1889 php.AddLast(&hp[ih++]); np[isel]++;
1891 AliInfo(Form("Build %3d 3D projections.", ih));
1894 Int_t ch(0), pt(0), sp(1), rcBin(as?as->FindBin(0.):-1), ioff(0);
1895 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1896 v = H->GetBinContent(ib, coord);
1898 if(fBCbinTOF>0 && coord[kBC]!=fBCbinTOF) continue; // TOF bunch cross cut
1899 if(fBCbinFill>0 && abf && coord[kNdim+2]!=fBCbinTOF) continue; // Fill bunch cut
1901 ch = 0; sp=1;// [-] track
1902 if(rcBin>0){ // debug mode in which species are also saved
1903 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1904 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
1905 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
1909 if(ap) pt = TMath::Min(coord[kPt]-1, Int_t(kNpt)-1);
1911 isel = pt*11; isel+=sp<0?10:(sp*kNcharge+ch);
1912 ioff = pt*34; ioff+=3*(sp<0?10:(sp*kNcharge+ch));
1913 AliDebug(4, Form("SELECTION[%d] :: ch[%c] pt[%c] sp[%d]\n", np[isel], ch==2?'Z':chName[ch], ptName[pt], sp));
1914 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1916 TObjArray *arr(NULL);
1917 fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx);
1919 TH2 *h2(NULL); Int_t jh(0);
1921 if(!hp[ih].fH) continue;
1922 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
1923 arr->AddAt(h2, jh++);
1925 // build combined performance plots
1926 // combine up the tree of projections
1927 AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1928 AliTRDresolutionProjection xlow[2], specY[kNcharge*AliPID::kSPECIES], specPh[kNcharge*AliPID::kSPECIES];
1929 for(Int_t ich(0); ich<kNcharge; ich++){
1930 // PID dependency - summation over pt
1931 for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
1933 Int_t idx(ich*AliPID::kSPECIES+isp);
1934 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
1935 specY[idx] = (*pr0);
1936 specY[idx].SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
1937 specY[idx].fH->SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp),
1938 Form("TrackIn[%s%c]:: #Deltay", AliPID::ParticleLatexName(isp), chSgn[ich]));
1939 for(Int_t ipt(1); ipt<nPt; ipt++){
1940 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
1943 php.AddLast(&specY[idx]);
1944 if((h2 = specY[idx].Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
1945 if((h2 = (TH2*)gDirectory->Get(Form("%s_yx", specY[idx].fH->GetName())))) arr->AddAt(h2, jh++);
1946 if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specY[idx];
1949 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
1950 specPh[idx] = (*pr0);
1951 specPh[idx].SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
1952 specPh[idx].fH->SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp),
1953 Form("TrackIn[%s%c]:: #Delta#phi", AliPID::ParticleLatexName(isp), chSgn[ich]));
1954 specPh[idx].SetShowRange(-1.5, 1.5);
1955 for(Int_t ipt(1); ipt<nPt; ipt++){
1956 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
1957 specPh[idx]+=(*pr1);
1959 php.AddLast(&specPh[idx]);
1960 if((h2 = specPh[idx].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1961 if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specPh[idx];
1964 // pt dependency - summation over PID
1965 for(Int_t ipt(0); ipt<nPt; ipt++){
1967 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
1968 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
1969 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
1972 pr0->fH->SetNameTitle(Form("H%sTrkInY%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
1973 Form("TrackIn[%c]:: #Deltay{%s}", chSgn[ich], ptCut[ipt]));
1974 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1975 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
1978 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
1979 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
1980 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
1983 pr0->fH->SetNameTitle(Form("H%sTrkInPh%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
1984 Form("TrackIn[%c]:: #Delta#phi{%s}", chSgn[ich], ptCut[ipt]));
1985 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1986 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
1989 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
1990 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
1991 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
1994 pr0->fH->SetNameTitle(Form("H%sTrkInX%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
1995 Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[ipt]));
1996 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1999 xlow[ich].SetNameTitle(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5),
2000 Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[0]));
2001 php.AddLast(&xlow[ich]);
2003 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2007 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2008 pr0->fH->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", chName[ich]),
2009 Form("TrackIn[%c]:: #Deltay", chSgn[ich]));
2010 pr0->SetShowRange(-0.3, 0.3);
2011 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2012 if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2015 if(ich && (pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[2], 0)))){
2016 if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[2], 0)))){
2018 pr0->fH->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", ptName[2]), Form("TrackIn :: #Deltay{%s}", ptCut[2]));
2019 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2023 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2024 pr0->fH->SetNameTitle(Form("H%sTrkInPh%c", mc?"MC":"", chName[ich]),
2025 Form("TrackIn[%c]:: #Delta#phi", chSgn[ich]));
2026 pr0->SetShowRange(-1., 1.);
2027 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2028 if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2031 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2032 pr0->fH->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", chName[ich]),
2033 Form("TrackIn[%c]:: #Deltax", chSgn[ich]));
2034 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2035 if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2038 if(ich && (pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 5)))){
2039 if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5)))){
2041 pr0->fH->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", ptName[0]), Form("TrackIn :: #Deltax{%s}", ptCut[0]));
2042 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2046 for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2048 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))){
2049 pr0->fH->SetNameTitle(Form("H%sTrkInY%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Deltay", AliPID::ParticleLatexName(isp)));
2050 pr0->SetShowRange(-0.3, 0.3);
2051 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2052 if((h2 = (TH2*)gDirectory->Get(Form("%s_yx", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
2055 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))){
2056 pr0->fH->SetNameTitle(Form("H%sTrkInPh%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Delta#phi", AliPID::ParticleLatexName(isp)));
2057 pr0->SetShowRange(-1., 1.);
2058 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2062 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2063 pr0->fH->SetNameTitle(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: #Deltay");
2064 pr0->SetShowRange(-0.3, 0.3);
2065 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2068 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2069 pr0->fH->SetNameTitle(Form("H%sTrkInPh", mc?"MC":""), "TrackIn :: #Delta#phi");
2070 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2073 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2074 pr0->fH->SetNameTitle(Form("H%sTrkInX", mc?"MC":""), "TrackIn :: #Deltax");
2075 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2078 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[0])))){
2079 for(Int_t ipt(0); ipt<kNpt; ipt++){
2080 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt])))) continue;
2083 pr0->fH->SetNameTitle(Form("H%sTrkInRCZ", mc?"MC":""), "TrackIn[RC]:: #Deltaz");
2084 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2087 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[0])))){
2088 for(Int_t ipt(0); ipt<kNpt; ipt++){
2089 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt])))) continue;
2092 pr0->fH->SetNameTitle(Form("H%sTrkInRCY", mc?"MC":""), "TrackIn[RC]:: #Deltay");
2093 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2096 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[0])))){
2097 for(Int_t ipt(0); ipt<kNpt; ipt++){
2098 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt])))) continue;
2101 pr0->fH->SetNameTitle(Form("H%sTrkInRCPh", mc?"MC":""), "TrackIn[RC]:: #Delta#phi");
2102 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2105 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[0])))){
2106 for(Int_t ipt(0); ipt<kNpt; ipt++){
2107 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt])))) continue;
2110 pr0->fH->SetNameTitle(Form("H%sTrkInRCX", mc?"MC":""), "TrackIn[RC]:: #Deltax");
2111 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2113 AliInfo(Form("Done %3d 2D projections.", jh));
2118 //________________________________________________________
2119 Bool_t AliTRDresolution::MakeProjectionTrack()
2122 const Int_t kNcontours(9);
2123 const Int_t kNstat(30);
2124 Int_t cidx(kMCtrack);
2125 if(fProj && fProj->At(cidx)) return kTRUE;
2127 AliError("Missing data container.");
2131 if(!(H = (THnSparse*)fContainer->FindObject("hTRD2MC"))){
2132 AliError("Missing/Wrong data @ hTRD2MC.");
2135 Int_t ndim(H->GetNdimensions());
2136 Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
2137 TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
2138 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2139 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
2140 if(ndim > kPt) ap = H->GetAxis(kPt);
2142 // build list of projections
2143 const Int_t nsel(AliTRDgeometry::kNlayer*kNpt*AliPID::kSPECIES*7);//, npsel(3);
2144 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2145 const Char_t ptName[kNpt] = {'l', 'i', 'h'};
2146 const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"};
2147 // define rebinning strategy
2148 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
2149 AliTRDresolutionProjection hp[kTrkNproj]; TObjArray php(kTrkNproj);
2150 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2151 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2152 for(Int_t ipt(0); ipt<kNpt; ipt++){
2153 for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2154 for(Int_t ich(0); ich<kNcharge; ich++){
2155 isel++; // new selection
2156 hp[ih].Build(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2157 Form("Tracks[%s%c]:: #Deltay{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2158 kEta, kPhi, kYrez, aa);
2159 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2160 php.AddLast(&hp[ih++]); np[isel]++;
2161 hp[ih].Build(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2162 Form("Tracks[%s%c]:: #Delta#phi{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2163 kEta, kPhi, kPrez, aa);
2164 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2165 php.AddLast(&hp[ih++]); np[isel]++;
2166 hp[ih].Build(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2167 Form("Tracks[%s%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2168 kEta, kPhi, kNdim, aa);
2169 hp[ih].SetShowRange(0.,10.);
2170 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2171 php.AddLast(&hp[ih++]); np[isel]++;
2174 isel++; // new selection
2175 hp[ih].Build(Form("HMCTrkZ%c%d", ptName[ipt], ily),
2176 Form("Tracks[RC]:: #Deltaz{%s} Ly[%d]", ptCut[ipt], ily),
2177 kEta, kPhi, kZrez, aa);
2178 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2179 php.AddLast(&hp[ih++]); np[isel]++;
2183 Int_t ly(0), ch(0), pt(0), sp(2), rcBin(as?as->FindBin(0.):-1);
2184 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2185 v = H->GetBinContent(ib, coord);
2187 ly = coord[kBC]-1; // layer selection
2189 ch=0; sp=2;// [pi-] track [dafault]
2190 if(rcBin>0){ // debug mode in which species are also saved
2191 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2192 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
2193 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
2196 pt = 0; // low pt [default]
2197 if(ap) pt = coord[kPt]-1;
2199 Int_t ioff = ly*kNpt*31+pt*31; ioff+=3*(sp<0?10:(sp*kNcharge+ch));
2200 isel = ly*kNpt*11+pt*11; isel+=sp<0?10:(sp*kNcharge+ch);
2201 AliDebug(4, Form("SELECTION[%d] :: ch[%c] pt[%c] sp[%d] ly[%d]\n", np[isel], ch==2?'Z':chName[ch], ptName[pt], sp, ly));
2202 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
2204 TObjArray *arr(NULL);
2205 fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
2207 TH2 *h2(NULL); Int_t jh(0);
2209 if(!hp[ih].fH) continue;
2210 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2211 arr->AddAt(h2, jh++);
2214 // combine up the tree of projections
2215 AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
2216 //Int_t iproj(0), jproj(0);
2217 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2218 for(Int_t ich(0); ich<kNcharge; ich++){
2219 for(Int_t ipt(0); ipt<kNpt; ipt++){
2221 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2222 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2223 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2226 AliDebug(2, Form("Rename %s to HMCTrkY%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2227 pr0->fH->SetNameTitle(Form("HMCTrkY%c%c%d", chName[ich], ptName[ipt], ily),
2228 Form("Tracks[%c]:: #Deltay{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2229 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2230 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2233 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2234 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2235 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2238 AliDebug(2, Form("Rename %s to HMCTrkPh%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2239 pr0->fH->SetNameTitle(Form("HMCTrkPh%c%c%d", chName[ich], ptName[ipt], ily),
2240 Form("Tracks[%c]:: #Delta#phi{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2241 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2242 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2246 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2247 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2248 if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2251 AliDebug(2, Form("Rename %s to HMCTrkDPt%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2252 pr0->fH->SetNameTitle(Form("HMCTrkDPt%c%c%d", chName[ich], ptName[ipt], ily),
2253 Form("Tracks[%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2254 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2255 if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2259 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2260 pr0->fH->SetNameTitle(Form("HMCTrkY%c%d", chName[ich], ily),
2261 Form("Tracks[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
2262 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2263 if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2266 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2267 pr0->fH->SetNameTitle(Form("HMCTrkPh%c%d", chName[ich], ily),
2268 Form("Tracks[%c]:: #Delta#phi Ly[%d]", chSgn[ich], ily));
2269 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2270 if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2273 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2274 pr0->fH->SetNameTitle(Form("HMCTrkDPt%c%d", chName[ich], ily),
2275 Form("Tracks[%c]:: #Deltap_{t}/p_{t} Ly[%d]", chSgn[ich], ily));
2276 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2277 if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2281 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2282 pr0->fH->SetNameTitle(Form("HMCTrkY%d", ily), Form("Tracks :: #Deltay Ly[%d]", ily));
2283 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2286 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2287 pr0->fH->SetNameTitle(Form("HMCTrkPh%d", ily), Form("Tracks :: #Delta#phi Ly[%d]", ily));
2288 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2291 if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2292 pr0->fH->SetNameTitle(Form("HMCTrkDPt%d", ily), Form("Tracks :: #Deltap_{t}/p_{t} Ly[%d]", ily));
2293 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2296 AliInfo(Form("Done %3d 2D projections.", jh));
2300 //________________________________________________________
2301 Bool_t AliTRDresolution::PostProcess()
2303 // Fit, Project, Combine, Extract values from the containers filled during execution
2306 AliError("ERROR: list not available");
2310 AliInfo("Building array of projections ...");
2311 fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
2314 //PROCESS EXPERIMENTAL DISTRIBUTIONS
2315 // Clusters detector
2316 if(!MakeProjectionDetector()) return kFALSE;
2317 // Clusters residuals
2318 if(!MakeProjectionCluster()) return kFALSE;
2320 // Tracklet residual/pulls
2321 if(!MakeProjectionTracklet()) return kFALSE;
2323 // TRDin residual/pulls
2324 if(!MakeProjectionTrackIn()) return kFALSE;
2327 if(!HasMCdata()) return kTRUE;
2328 //PROCESS MC RESIDUAL DISTRIBUTIONS
2330 // CLUSTER Y RESOLUTION/PULLS
2331 if(!MakeProjectionCluster(kTRUE)) return kFALSE;
2334 // TRACKLET RESOLUTION/PULLS
2335 if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
2338 // TRACK RESOLUTION/PULLS
2339 if(!MakeProjectionTrack()) return kFALSE;
2342 // TRACK TRDin RESOLUTION/PULLS
2343 if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
2350 //________________________________________________________
2351 void AliTRDresolution::Terminate(Option_t *opt)
2353 AliTRDrecoTask::Terminate(opt);
2354 if(HasPostProcess()) PostProcess();
2357 //________________________________________________________
2358 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2360 // Helper function to avoid duplication of code
2361 // Make first guesses on the fit parameters
2363 // find the intial parameters of the fit !! (thanks George)
2364 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2366 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2367 f->SetParLimits(0, 0., 3.*sum);
2368 f->SetParameter(0, .9*sum);
2369 Double_t rms = h->GetRMS();
2370 f->SetParLimits(1, -rms, rms);
2371 f->SetParameter(1, h->GetMean());
2373 f->SetParLimits(2, 0., 2.*rms);
2374 f->SetParameter(2, rms);
2375 if(f->GetNpar() <= 4) return;
2377 f->SetParLimits(3, 0., sum);
2378 f->SetParameter(3, .1*sum);
2380 f->SetParLimits(4, -.3, .3);
2381 f->SetParameter(4, 0.);
2383 f->SetParLimits(5, 0., 1.e2);
2384 f->SetParameter(5, 2.e-1);
2387 //________________________________________________________
2388 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
2390 // Build performance histograms for AliTRDcluster.vs TRD track or MC
2391 // - y reziduals/pulls
2393 TObjArray *arr = new TObjArray(2);
2394 arr->SetName(name); arr->SetOwner();
2395 TH1 *h(NULL); char hname[100], htitle[300];
2397 // tracklet resolution/pull in y direction
2398 snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
2399 snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
2400 Float_t rr = range<0.?fDyRange:range;
2401 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2403 if(expand) nybins*=2;
2404 h = new TH3S(hname, htitle,
2405 48, -.48, .48, // phi
2407 nybins, -0.5, nybins-0.5);// segment
2410 snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
2411 snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, "Detector");
2412 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2413 h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
2420 //________________________________________________________
2421 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
2423 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
2424 // - y reziduals/pulls
2425 // - z reziduals/pulls
2427 TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
2429 TH1 *h(NULL); char hname[100], htitle[300];
2431 // tracklet resolution/pull in z direction
2432 snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
2433 snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
2434 if(!(h = (TH2S*)gROOT->FindObject(hname))){
2435 h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
2438 snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
2439 snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
2440 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2441 h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
2442 h->GetZaxis()->SetBinLabel(1, "no RC");
2443 h->GetZaxis()->SetBinLabel(2, "RC");
2447 // tracklet to track phi resolution
2448 snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
2449 snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
2451 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2452 h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
2459 //________________________________________________________
2460 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
2462 // Build performance histograms for AliExternalTrackParam.vs MC
2463 // - y resolution/pulls
2464 // - z resolution/pulls
2465 // - phi resolution, snp pulls
2466 // - theta resolution, tgl pulls
2467 // - pt resolution, 1/pt pulls, p resolution
2469 TObjArray *arr = BuildMonitorContainerTracklet(name);
2471 TH1 *h(NULL); char hname[100], htitle[300];
2475 snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
2476 snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
2477 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2478 h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
2483 snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
2484 snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
2485 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2486 h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
2490 snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
2491 snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
2492 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2493 h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
2497 const Int_t kNdpt(150);
2498 const Int_t kNspc = 2*AliPID::kSPECIES+1;
2499 Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
2500 Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
2501 for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
2502 for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
2503 for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
2506 snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
2507 snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
2508 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2509 h = new TH3S(hname, htitle,
2510 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2511 //ax = h->GetZaxis();
2512 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2516 snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
2517 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);
2518 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2519 h = new TH3S(hname, htitle,
2520 kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
2521 //ax = h->GetZaxis();
2522 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2526 snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
2527 snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
2528 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2529 h = new TH3S(hname, htitle,
2530 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2531 //ax = h->GetZaxis();
2532 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2540 //________________________________________________________
2541 TObjArray* AliTRDresolution::Histos()
2544 // Define histograms
2547 if(fContainer) return fContainer;
2549 fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
2551 const Int_t nhn(100); Char_t hn[nhn]; TString st;
2553 //++++++++++++++++++++++
2554 // cluster to detector
2555 snprintf(hn, nhn, "h%s", fgPerformanceName[kDetector]);
2556 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2557 const Int_t mdim(5);
2558 const Char_t *cldTitle[mdim] = {"layer", fgkTitle[kPhi], "pad row", "centrality", "q [a.u.]"/*, "occupancy", "tb [10 Hz]"*/};
2559 const Int_t cldNbins[mdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], 76, AliTRDeventInfo::kCentralityClasses, 100};
2560 const Double_t cldMin[mdim] = {-0.5, fgkMin[kPhi], -0.5, -0.5, 0.},
2561 cldMax[mdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], 75.5, AliTRDeventInfo::kCentralityClasses - 0.5, 1200.};
2562 st = "cluster proprieties;";
2563 // define minimum info to be saved in non debug mode
2564 Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimDet);
2565 for(Int_t idim(0); idim<ndim; idim++){ st += cldTitle[idim]; st+=";";}
2566 H = new THnSparseI(hn, st.Data(), ndim, cldNbins, cldMin, cldMax);
2568 fContainer->AddAt(H, kDetector);
2569 //++++++++++++++++++++++
2570 // cluster to tracklet residuals/pulls
2571 snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
2572 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2573 const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]"};
2574 const Int_t clNbins[kNdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 30, 30, 15};
2575 const Double_t clMin[kNdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45},
2576 clMax[kNdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45};
2577 st = "cluster spatial&charge resolution;";
2578 // define minimum info to be saved in non debug mode
2579 Int_t ndim=DebugLevel()>=1?Int_t(kNdim):Int_t(kNdimCl);
2580 for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2581 H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2583 fContainer->AddAt(H, kCluster);
2584 //++++++++++++++++++++++
2585 // tracklet to TRD track
2586 snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
2587 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2588 const Int_t mdim(kNdim+8);
2589 Char_t *trTitle[mdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2590 Int_t trNbins[mdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2591 Double_t trMin[mdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2592 Double_t trMax[mdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2593 // set specific fields
2594 trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
2595 trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
2596 trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
2597 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2598 Int_t jdim(kNdim); trTitle[jdim]=StrDup("dq/dl [a.u.]"); trNbins[jdim] = 100; trMin[jdim] = 0.; trMax[jdim] = 20;
2599 jdim++; trTitle[jdim]=StrDup("centrality"); trNbins[jdim] = AliTRDeventInfo::kCentralityClasses; trMin[jdim] = -.5; trMax[jdim] = AliTRDeventInfo::kCentralityClasses - 0.5;
2600 jdim++; trTitle[jdim]=StrDup("occupancy [%]"); trNbins[jdim] = 12; trMin[jdim] = 25.; trMax[jdim] = 85.;
2601 jdim++; trTitle[jdim]=StrDup("gap [cm]"); trNbins[jdim] = 80; trMin[jdim] = 0.1; trMax[jdim] = 2.1;
2602 jdim++; trTitle[jdim]=StrDup("size_{0} [cm]"); trNbins[jdim] = 16; trMin[jdim] = 0.15; trMax[jdim] = 1.75;
2603 jdim++; trTitle[jdim]=StrDup("pos_{0} [cm]"); trNbins[jdim] = 10; trMin[jdim] = 0.; trMax[jdim] = 3.5;
2604 jdim++; trTitle[jdim]=StrDup("size_{1} [cm]"); trNbins[jdim] = 16; trMin[jdim] = 0.15; trMax[jdim] = 1.75;
2605 jdim++; trTitle[jdim]=StrDup("pos_{1} [cm]"); trNbins[jdim] = 10; trMin[jdim] = 0.; trMax[jdim] = 3.5;
2607 st = "tracklet spatial&charge resolution;";
2608 // define minimum info to be saved in non debug mode
2609 Int_t ndim=DebugLevel()>=1?mdim:kNdimTrklt;
2610 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2611 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2613 fContainer->AddAt(H, kTracklet);
2614 //++++++++++++++++++++++
2615 // tracklet to TRDin
2616 snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
2617 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2618 // set specific fields
2619 const Int_t mdim(kNdim+3);
2620 Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*));
2621 Int_t trinNbins[mdim]; memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t));
2622 Double_t trinMin[mdim]; memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t));
2623 Double_t trinMax[mdim]; memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t));
2624 trinNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1; trinMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trinMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2625 trinTitle[kNdim]=StrDup("detector"); trinNbins[kNdim] = 540; trinMin[kNdim] = -0.5; trinMax[kNdim] = 539.5;
2626 trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4;
2627 trinTitle[kNdim+2]=StrDup("Fill Bunch"); trinNbins[kNdim+2]=3500; trinMin[kNdim+2]=-0.5; trinMax[kNdim+2]=3499.5;
2628 st = "r-#phi/z/angular residuals @ TRD entry;";
2629 // define minimum info to be saved in non debug mode
2630 Int_t ndim=DebugLevel()>=1?mdim:kNdimTrkIn;
2631 for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";}
2632 H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax);
2634 fContainer->AddAt(H, kTrackIn);
2635 // tracklet to TRDout
2636 // fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2639 // Resolution histos
2640 if(!HasMCdata()) return fContainer;
2642 //++++++++++++++++++++++
2643 // cluster to TrackRef residuals/pulls
2644 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
2645 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2646 const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", fgkTitle[kSpeciesChgRC], "#Phi [deg]"};
2647 const Int_t clNbins[kNdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 20, 10, Int_t(kNcharge)*AliPID::kSPECIES+1, 15};
2648 const Double_t clMin[kNdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, -AliPID::kSPECIES-0.5, -45},
2649 clMax[kNdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, AliPID::kSPECIES+0.5, 45};
2650 st = "MC cluster spatial resolution;";
2651 // define minimum info to be saved in non debug mode
2652 Int_t ndim=DebugLevel()>=1?kNdim:4;
2653 for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2654 H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2656 fContainer->AddAt(H, kMCcluster);
2657 //++++++++++++++++++++++
2658 // tracklet to TrackRef
2659 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
2660 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2661 Char_t *trTitle[kNdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2662 Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2663 Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2664 Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2665 // set specific fields
2666 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2667 trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
2668 trMin[kPrez] = -4.5; trMax[kPrez] = -trMin[kPrez];
2669 trMin[kZrez] = -1.5; trMax[kZrez] = -trMin[kZrez];
2670 trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2672 st = "MC tracklet spatial resolution;";
2673 // define minimum info to be saved in non debug mode
2674 Int_t ndim=DebugLevel()>=1?kNdim:4;
2675 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2676 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2678 fContainer->AddAt(H, kMCtracklet);
2679 //++++++++++++++++++++++
2680 // TRDin to TrackRef
2681 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
2682 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2683 st = "MC r-#phi/z/angular residuals @ TRD entry;";
2684 // set specific fields
2685 Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2686 Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2687 Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2688 trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
2689 trMin[kPrez] = -2.4; trMax[kPrez] = -trMin[kPrez];
2690 trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2691 trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2692 // define minimum info to be saved in non debug mode
2693 Int_t ndim=DebugLevel()>=1?kNdim:7;
2694 for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
2695 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2697 fContainer->AddAt(H, kMCtrackIn);
2698 //++++++++++++++++++++++
2699 // track to TrackRef
2700 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
2701 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2702 Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2703 Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2704 Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2705 Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2706 // set specific fields
2707 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2708 trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
2709 trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
2710 trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2711 trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2712 trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5;
2714 st = "MC track spatial&p_{t} resolution;";
2715 // define minimum info to be saved in non debug mode
2716 Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
2717 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2718 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2720 fContainer->AddAt(H, kMCtrack);
2725 //________________________________________________________
2726 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
2728 // Robust function to process sigma/mean for 2D plot dy(x)
2729 // For each x bin a gauss fit is performed on the y projection and the range
2730 // with the minimum chi2/ndf is choosen
2733 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
2736 if(!Int_t(h2->GetEntries())){
2737 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
2740 if(!g || !g[0]|| !g[1]) {
2741 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
2745 TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
2746 Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
2747 TF1 f("f", "gaus", ymin, ymax);
2749 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2750 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2752 if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
2753 Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
2757 Float_t chi2OverNdf(0.);
2758 for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
2759 x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
2760 ymin = ay->GetXmin(); ymax = ay->GetXmax();
2762 h = h2->ProjectionY("py", ix, ix);
2763 if((n=(Int_t)h->GetEntries())<stat){
2764 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
2767 // looking for a first order mean value
2768 f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
2770 chi2OverNdf = f.GetChisquare()/f.GetNDF();
2771 printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
2772 y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
2773 ey = f.GetParError(1);
2774 sy = f.GetParameter(2); esy = f.GetParError(2);
2775 // // looking for the best chi2/ndf value
2776 // while(ymin<y0 && ymax>y1){
2777 // f.SetParameter(1, y);
2778 // f.SetParameter(2, sy);
2779 // h->Fit(&f, "QNW", "", y0, y1);
2780 // printf(" range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
2781 // if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
2782 // chi2OverNdf = f.GetChisquare()/f.GetNDF();
2783 // y = f.GetParameter(1); ey = f.GetParError(1);
2784 // sy = f.GetParameter(2); esy = f.GetParError(2);
2785 // printf(" set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
2789 g[0]->SetPoint(np, x, y);
2790 g[0]->SetPointError(np, ex, ey);
2791 g[1]->SetPoint(np, x, sy);
2792 g[1]->SetPointError(np, ex, esy);
2799 //________________________________________________________
2800 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
2803 // Do the processing
2806 Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
2808 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2809 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2810 if(Int_t(h2->GetEntries())){
2811 AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
2813 AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
2818 const Int_t kINTEGRAL=1;
2819 for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
2820 Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
2821 Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
2822 TH1D *h = h2->ProjectionY(pn, abin, bbin);
2823 if((n=(Int_t)h->GetEntries())<150){
2824 AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
2828 Int_t ip = g[0]->GetN();
2829 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)));
2830 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
2831 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
2832 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
2833 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
2835 g[0]->SetPoint(ip, x, k*h->GetMean());
2836 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
2837 g[1]->SetPoint(ip, x, k*h->GetRMS());
2838 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2845 //____________________________________________________________________
2846 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2849 // Fit track with a staight line using the "np" clusters stored in the array "points".
2850 // The following particularities are stored in the clusters from points:
2851 // 1. pad tilt as cluster charge
2852 // 2. pad row cross or vertex constrain as fake cluster with cluster type 1
2853 // The parameters of the straight line fit are stored in the array "param" in the following order :
2854 // param[0] - x0 reference radial position
2855 // param[1] - y0 reference r-phi position @ x0
2856 // param[2] - z0 reference z position @ x0
2857 // param[3] - slope dy/dx
2858 // param[4] - slope dz/dx
2861 // Function should be used to refit tracks for B=0T
2865 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2868 TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2871 for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
2874 Double_t x, y, z, dx, tilt(0.);
2875 for(Int_t ip(0); ip<np; ip++){
2876 x = points[ip].GetX(); z = points[ip].GetZ();
2878 zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
2880 if(zfitter.Eval() != 0) return kFALSE;
2882 Double_t z0 = zfitter.GetParameter(0);
2883 Double_t dzdx = zfitter.GetParameter(1);
2884 for(Int_t ip(0); ip<np; ip++){
2885 if(points[ip].GetClusterType()) continue;
2886 x = points[ip].GetX();
2888 y = points[ip].GetY();
2889 z = points[ip].GetZ();
2890 tilt = points[ip].GetCharge();
2891 y -= tilt*(-dzdx*dx + z - z0);
2892 Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
2893 yfitter.AddPoint(&dx, y, 1.);
2895 if(yfitter.Eval() != 0) return kFALSE;
2896 Double_t y0 = yfitter.GetParameter(0);
2897 Double_t dydx = yfitter.GetParameter(1);
2899 param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
2900 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);
2904 //____________________________________________________________________
2905 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
2908 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
2909 // See function FitTrack for the data stored in the "clusters" array
2911 // The parameters of the straight line fit are stored in the array "param" in the following order :
2912 // par[0] - x0 reference radial position
2913 // par[1] - y0 reference r-phi position @ x0
2914 // par[2] - slope dy/dx
2917 // Function should be used to refit tracks for B=0T
2920 TLinearFitter yfitter(2, "pol1");
2922 // grep data for tracklet
2923 Double_t x0(0.), x[60], y[60], dy[60];
2925 for(Int_t ip(0); ip<np; ip++){
2926 if(points[ip].GetClusterType()) continue;
2927 if(points[ip].GetVolumeID() != ly) continue;
2928 Float_t xt(points[ip].GetX())
2929 ,yt(param[1] + param[3] * (xt - param[0]));
2931 y[nly] = points[ip].GetY();
2937 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
2940 // set radial reference for fit
2943 // find tracklet core
2944 Double_t mean(0.), sig(1.e3);
2945 AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
2947 // simple cluster error parameterization
2948 Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
2950 // fit tracklet core
2951 for(Int_t jly(0); jly<nly; jly++){
2952 if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
2953 Double_t dx(x[jly]-x0);
2954 yfitter.AddPoint(&dx, y[jly], 1.);
2956 if(yfitter.Eval() != 0) return kFALSE;
2958 par[1] = yfitter.GetParameter(0);
2959 par[2] = yfitter.GetParameter(1);
2963 //____________________________________________________________________
2964 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
2967 // Global selection mechanism of tracksbased on cluster to fit residuals
2968 // The parameters are the same as used ni function FitTrack().
2970 const Float_t kS(0.6), kM(0.2);
2971 TH1S h("h1", "", 100, -5.*kS, 5.*kS);
2972 Float_t dy, dz, s, m;
2973 for(Int_t ip(0); ip<np; ip++){
2974 if(points[ip].GetClusterType()) continue;
2975 Float_t x0(points[ip].GetX())
2976 ,y0(param[1] + param[3] * (x0 - param[0]))
2977 ,z0(param[2] + param[4] * (x0 - param[0]));
2978 dy=points[ip].GetY() - y0; h.Fill(dy);
2979 dz=points[ip].GetZ() - z0;
2981 TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
2982 fg.SetParameter(1, 0.);
2983 fg.SetParameter(2, 2.e-2);
2985 m=fg.GetParameter(1); s=fg.GetParameter(2);
2986 if(s>kS || TMath::Abs(m)>kM) return kFALSE;
2990 //________________________________________________________
2991 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2994 // Get the most probable value and the full width half mean
2995 // of a Landau distribution
2998 const Float_t dx = 1.;
2999 mpv = f->GetParameter(1);
3000 Float_t fx, max = f->Eval(mpv);
3003 while((fx = f->Eval(xm))>.5*max){
3012 while((fx = f->Eval(xM))>.5*max) xM += dx;
3016 // #include "TFile.h"
3017 // //________________________________________________________
3018 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
3021 // AliWarning("Use cluster position as in reconstruction.");
3022 // SetLoadCorrection();
3025 // TDirectory *cwd(gDirectory);
3026 // TString fileList;
3027 // FILE *filePtr = fopen(file, "rt");
3029 // AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
3030 // SetLoadCorrection();
3033 // TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
3034 // while(fileList.Gets(filePtr)){
3035 // if(!TFile::Open(fileList.Data())) {
3036 // AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
3038 // } else AliInfo(Form("\"%s\"", fileList.Data()));
3040 // TTree *tSys = (TTree*)gFile->Get("tSys");
3041 // h2->SetDirectory(gDirectory); h2->Reset("ICE");
3042 // tSys->Draw("det:t>>h2", "dx", "goff");
3043 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3044 // for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
3046 // h2->SetDirectory(cwd);
3051 // if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
3052 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3053 // printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
3054 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3055 // FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
3056 // fprintf(fout, " static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
3057 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3058 // printf("%03d|", idet);
3059 // fprintf(fout, " {");
3060 // for(Int_t it(0); it<30; it++){
3061 // printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
3062 // fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
3065 // fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
3067 // fprintf(fout, " };\n");
3069 // SetLoadCorrection();
3073 //________________________________________________________
3074 AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection()
3082 memset(fAx, 0, 3*sizeof(Int_t));
3083 memset(fRange, 0, 4*sizeof(Float_t));
3086 //________________________________________________________
3087 AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection()
3093 //________________________________________________________
3094 void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])
3096 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"
3097 if(!aa[ix] || !aa[iy] || !aa[iz]) return;
3098 TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);
3099 // check ax definiton to protect against older versions of the data
3100 if(ax->GetNbins()<=0 || (ax->GetXmax()-ax->GetXmin())<=0.){
3101 AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax()));
3104 if(ay->GetNbins()<=0 || (ay->GetXmax()-ay->GetXmin())<=0.){
3105 AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax()));
3108 if(az->GetNbins()<=0 || (az->GetXmax()-az->GetXmin())<=0.){
3109 AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));
3113 fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),
3114 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
3115 ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
3116 az->GetNbins(), az->GetXmin(), az->GetXmax());
3117 fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;
3118 fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;
3121 //________________________________________________________
3122 AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator=(const AliTRDresolutionProjection& rhs)
3125 if(this == &rhs) return *this;
3127 TNamed::operator=(rhs);
3128 if(fNrebin){fNrebin=0; delete [] fRebinX; delete [] fRebinY;}
3129 if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebinX, rhs.fRebinY);
3130 memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));
3131 memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));
3133 if(rhs.fH) fH=(TH3I*)rhs.fH->Clone(Form("%s_CLONE", rhs.fH->GetName()));
3137 //________________________________________________________
3138 AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator+=(const AliTRDresolutionProjection& other)
3140 // increment projections
3141 if(!fH || !other.fH) return *this;
3142 AliDebug(2, Form("%s+=%s [%s+=%s]", GetName(), other.GetName(), fH->GetName(), (other.fH)->GetName()));
3147 //________________________________________________________
3148 void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v)
3150 // increment bin with value "v" pointed by general coord in "bin"
3152 AliDebug(4, Form(" %s", fH->GetName()));
3153 fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));
3156 //________________________________________________________
3157 TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid, Bool_t del)
3159 // build the 2D projection and adjust binning
3161 const Char_t *title[] = {"Mean", "#mu", "MPV"};
3162 if(!fH) return NULL;
3163 TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());
3165 if(!(h2s = (TH2*)fH->Project3D("yx"))) return NULL;
3166 Int_t irebin(0), dxBin(1), dyBin(1);
3167 while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){
3168 h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);
3169 dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];
3172 Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
3173 if(h2s && del) delete h2s;
3176 TH1 *h(NULL); Int_t n(0);
3177 Float_t dz=(fRange[1]-fRange[1])/ncol;
3178 TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));
3179 TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),
3180 Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], (*tokenTitle)[0]->GetName(), tokenTitle->GetEntriesFast()>1?(*tokenTitle)[1]->GetName():""),
3181 nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());
3182 h2->SetContour(ncol);
3183 h2->GetZaxis()->CenterTitle();
3184 h2->GetZaxis()->SetTitleOffset(1.4);
3185 h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);
3186 //printf("%s[%s] nx[%d] ny[%d]\n", h2->GetName(), h2->GetTitle(), nx, ny);
3187 for(Int_t iy(0); iy<ny; iy++){
3188 for(Int_t ix(0); ix<nx; ix++){
3189 h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);
3190 Int_t ne((Int_t)h->Integral());
3192 h2->SetBinContent(ix+1, iy+1, -999);
3193 h2->SetBinError(ix+1, iy+1, 1.);
3196 Float_t v(h->GetMean()), ve(h->GetRMS());
3198 TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());
3199 fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);
3201 v = fg.GetParameter(1); ve = fg.GetParameter(2);
3202 } else if (mid==2) {
3203 TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());
3204 fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);
3206 v = fl.GetMaximumX(); ve = fl.GetParameter(2);
3207 /* TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());
3208 fgle.SetParameter(0, fl.GetParameter(0));
3209 fgle.SetParameter(1, fl.GetParameter(1));
3210 fgle.SetParameter(2, fl.GetParameter(2));
3211 fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);
3212 h->Fit(&fgle, "WQ");
3213 v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/
3215 if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);
3216 else h2->SetBinContent(ix+1, iy+1, v);
3217 h2->SetBinError(ix+1, iy+1, ve);
3222 if(n==nx*ny){delete h2; h2=NULL;} // clean empty projections
3226 //________________________________________________________
3227 void AliTRDresolution::SetRangeZ(TH2 *h2, Float_t min, Float_t max)
3229 // Set range on Z axis such to avoid outliers
3231 Float_t c(0.), dz(1.e-3*(max-min));
3232 for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){
3233 for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){
3234 if((c = h2->GetBinContent(ix, iy))<10) continue;
3235 if(c<=min) h2->SetBinContent(ix, iy, min+dz);
3238 h2->GetZaxis()->SetRangeUser(min, max);
3241 //________________________________________________________
3242 void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])
3244 // define rebinning strategy for this projection
3246 fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));
3247 fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));