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 ////////////////////////////////////////////////////////////////////////////
55 #include <TObjArray.h>
59 #include <THnSparse.h>
65 #include <TGraphErrors.h>
66 #include <TGraphAsymmErrors.h>
67 #include <TLinearFitter.h>
71 #include <TTreeStream.h>
72 #include <TGeoManager.h>
73 #include <TDatabasePDG.h>
77 #include "AliESDtrack.h"
78 #include "AliMathBase.h"
79 #include "AliTrackPointArray.h"
81 #include "AliTRDresolution.h"
82 #include "AliTRDgeometry.h"
83 #include "AliTRDtransform.h"
84 #include "AliTRDpadPlane.h"
85 #include "AliTRDcluster.h"
86 #include "AliTRDseedV1.h"
87 #include "AliTRDtrackV1.h"
88 #include "AliTRDReconstructor.h"
89 #include "AliTRDrecoParam.h"
90 #include "AliTRDpidUtil.h"
91 #include "AliTRDinfoGen.h"
92 #include "AliAnalysisManager.h"
93 #include "info/AliTRDclusterInfo.h"
94 #include "info/AliTRDeventInfo.h"
96 ClassImp(AliTRDresolution)
97 ClassImp(AliTRDresolution::AliTRDrecoProjection)
99 Int_t const AliTRDresolution::fgkNbins[kNdim] = {
100 Int_t(kNbunchCross)/*bc*/,
108 }; //! no of bins/projection
109 Double_t const AliTRDresolution::fgkMin[kNdim] = {
118 }; //! low limits for projections
119 Double_t const AliTRDresolution::fgkMax[kNdim] = {
128 }; //! high limits for projections
129 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
138 }; //! title of projection
140 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
149 // ,"Tracklet2TRDout"
153 //________________________________________________________
154 AliTRDresolution::AliTRDresolution()
167 // Default constructor
169 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
170 SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
173 //________________________________________________________
174 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
175 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
187 // Default constructor
191 SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
193 SetUseExchangeContainers();
194 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
195 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
199 //________________________________________________________
200 AliTRDresolution::~AliTRDresolution()
205 if (AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
206 if(fProj){fProj->Delete(); delete fProj;}
207 if(fCl){fCl->Delete(); delete fCl;}
208 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
212 //________________________________________________________
213 void AliTRDresolution::UserCreateOutputObjects()
215 // spatial resolution
217 AliTRDrecoTask::UserCreateOutputObjects();
218 if(UseExchangeContainers()) InitExchangeContainers();
221 //________________________________________________________
222 void AliTRDresolution::InitExchangeContainers()
224 // Init containers for subsequent tasks (AliTRDclusterResolution)
226 fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
227 fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
228 PostData(kClToTrk, fCl);
229 PostData(kClToMC, fMCcl);
232 //________________________________________________________
233 void AliTRDresolution::UserExec(Option_t *opt)
239 if(fCl) fCl->Delete();
240 if(fMCcl) fMCcl->Delete();
241 AliTRDrecoTask::UserExec(opt);
244 //________________________________________________________
245 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
247 // Helper function to calculate pulls in the yz plane
248 // using proper tilt rotation
249 // Uses functionality defined by AliTRDseedV1.
253 Double_t t2(tilt*tilt);
254 // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
258 cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
259 cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
260 cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
262 Double_t sqr[3]={0., 0., 0.};
263 if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
264 Double_t invsqr[3]={0., 0., 0.};
265 if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
266 Double_t tmp(dyz[0]);
267 dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
268 dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
273 //________________________________________________________
274 TH1* AliTRDresolution::DetCluster(const TObjArray *cls)
277 // Plot the cluster distributions
280 if(cls) fkClusters = cls;
282 AliDebug(4, "No Clusters defined.");
286 if(!(ncl = fkClusters->GetEntriesFast())){
287 AliDebug(1, "No RecPoints defined.");
291 AliTRDcluster *cl(NULL);
292 for(Int_t icl(0); icl<ncl; icl++){
293 if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])) continue;
294 det = cl->GetDetector(); break;
297 AliDebug(1, "No useful clusters defined.");
301 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kDetector)))){
302 AliWarning("No output container defined.");
305 Int_t ly(AliTRDgeometry::GetLayer(det)),
306 stk(AliTRDgeometry::GetStack(det));
308 alpha((0.5+AliTRDgeometry::GetSector(det))*AliTRDgeometry::GetAlpha()),
309 cs(TMath::Cos(alpha)),
310 sn(TMath::Sin(alpha));
311 for(Int_t icl(0); icl<ncl; icl++){
312 if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])){
313 AliDebug(2, Form("Missing clusters @ %d", icl));
316 if(cl->IsMasked()) printf("Mask %03d[%02d_%d_%d] %d[%d]", det, AliTRDgeometry::GetSector(det), stk, ly, cl->GetPadMaskedPosition(), cl->GetPadMaskedStatus());
318 val[kPhi] = TMath::ATan2(cl->GetX()*sn + cl->GetY()*cs, cl->GetX()*cs - cl->GetY()*sn);
319 val[kEta] = (5-stk)*16-cl->GetPadRow()-1-(stk<3?4:0);
320 val[kYrez]= fEvent->GetMultiplicity();
321 val[4] = TMath::Abs(cl->GetQ());
322 val[5] = fTriggerList?fTriggerSlot:(cl->IsFivePad()?1:cl->GetNPads());
328 //________________________________________________________
329 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
332 // Plot the cluster distributions
335 if(track) fkTrack = track;
337 AliDebug(4, "No Track defined.");
340 if(fkESD && TMath::Abs(fkESD->GetTOFbc()) > 1){
341 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
344 if(fkESD && fPt<fPtThreshold){
345 AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
349 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
350 AliWarning("No output container defined.");
354 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
355 Double_t val[kNdim+2],
356 alpha(0.), cs(-2.), sn(0.); //Float_t exb, vd, t0, s2, dl, dt;
357 TObjArray *clInfoArr(NULL);
358 AliTRDseedV1 *fTracklet(NULL);
359 AliTRDcluster *c(NULL), *cc(NULL);
360 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
361 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
362 if(!fTracklet->IsOK()) continue;
363 //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
366 alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
367 cs = TMath::Cos(alpha);
368 sn = TMath::Sin(alpha);
370 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
371 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
372 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
373 val[kPt] = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
374 Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
376 Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
377 fTracklet->ResetClusterIter(kTRUE);
378 while((c = fTracklet->NextCluster())){
379 Float_t xc(c->GetX()),
380 q(TMath::Abs(c->GetQ()));
381 if(row0<0) row0 = c->GetPadRow();
383 val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
384 val[kPrez] = fTracklet->GetX0()-xc;
385 val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
386 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
387 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
388 if(ic) val[kZrez] /= (ic*q);
389 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
390 val[kNdim] = fEvent?fEvent->GetMultiplicity():0.;
391 val[kNdim+1] = fTriggerList?fTriggerSlot:(c->IsFivePad()?1:c->GetNPads());
393 /* // tilt rotation of covariance for clusters
394 Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
395 cov[0] = (sy2+t2*sz2)*corr;
396 cov[1] = tilt*(sz2 - sy2)*corr;
397 cov[2] = (t2*sy2+sz2)*corr;
398 // sum with track covariance
399 cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
400 Double_t dyz[2]= {dy[1], dz[1]};
401 Pulls(dyz, cov, tilt);*/
403 // Get z-position with respect to anode wire
404 Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
405 zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
406 Int_t istk = geo->GetStack(c->GetDetector());
407 AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
408 Float_t rowZ = pp->GetRow0();
409 Float_t d = rowZ - zt + pp->GetAnodeWireOffset();
410 d -= ((Int_t)(2 * d)) / 2.0;
411 if (d > 0.25) d = 0.5 - d;
413 AliTRDclusterInfo *clInfo(NULL);
414 clInfo = new AliTRDclusterInfo;
415 clInfo->SetCluster(c);
416 //Float_t covF[] = {cov[0], cov[1], cov[2]};
417 clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
418 clInfo->SetResolution(val[kYrez]);
419 clInfo->SetAnisochronity(d);
420 clInfo->SetDriftLength(val[kZrez]);
421 clInfo->SetTilt(fTracklet->GetTilt());
422 if(fCl) fCl->Add(clInfo);
423 //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
427 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
428 clInfoArr->SetOwner(kFALSE);
430 clInfoArr->Add(clInfo);
433 if(DebugLevel()>=2 && clInfoArr){
434 ULong_t status = fkESD->GetStatus();
435 (*DebugStream()) << "cluster"
436 <<"status=" << status
437 <<"clInfo.=" << clInfoArr
442 if(clInfoArr) delete clInfoArr;
444 if(!track) return NULL;
445 // special care for EVE usage
447 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
448 return H->Projection(kYrez);
452 //________________________________________________________
453 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
455 // Plot normalized residuals for tracklets to track.
457 // We start from the result that if X=N(|m|, |Cov|)
459 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
461 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
462 // reference position.
463 if(track) fkTrack = track;
465 AliDebug(4, "No Track defined.");
468 if(fkESD && TMath::Abs(fkESD->GetTOFbc())>1){
469 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
473 if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
474 AliWarning("No output container defined.");
478 // down scale PID resolution
479 Int_t spc(fSpecies); if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
480 const Int_t ndim(kNdim+8);
482 alpha(0.), cs(-2.), sn(0.);
483 // Float_t sz[AliTRDseedV1::kNtb], pos[AliTRDseedV1::kNtb]; Int_t ntbGap[AliTRDseedV1::kNtb];
484 AliTRDseedV1 *fTracklet(NULL);
485 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
486 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
487 if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()) continue;
490 alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
491 cs = TMath::Cos(alpha);
492 sn = TMath::Sin(alpha);
494 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
495 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());//fTracklet->GetTgl();
496 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
498 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:spc;// fSpecies;
499 val[kPt] = GetPtBin(fPt); //fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fTracklet->GetMomentum());
500 Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)),
501 dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)),
502 dydx(fTracklet->GetYfit(1)),
503 tilt(fTracklet->GetTilt());
504 // correct for tilt rotation
505 val[kYrez] = dyt - dzt*tilt;
506 val[kZrez] = fTracklet->IsRowCross()?(dzt + dyt*tilt):(fTracklet->GetdQdl()*3.e-4-1.5);
507 dydx+= tilt*fTracklet->GetZref(1);
508 val[kPrez] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
509 val[kNdim] = fEvent?fEvent->GetMultiplicity():0;
510 val[kNdim+1] = fTriggerList?fTriggerSlot:(1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb);
511 /* Int_t n = fTracklet->GetChargeGaps(sz, pos, ntbGap);
512 val[kNdim+2] = 0.; for(Int_t igap(0); igap<n; igap++) val[kNdim+2] += sz[igap];
513 for(Int_t ifill(0); ifill<3; ifill++){ val[kNdim+3+ifill]=0.;val[kNdim+4+ifill]=0.;}
514 for(Int_t igap(0), ifill(0); igap<n&&ifill<2; igap++){
515 if(ntbGap[igap]<2) continue;
516 val[kNdim+3+ifill] = sz[igap];
517 val[kNdim+4+ifill] = pos[igap];
522 // // compute covariance matrix
523 // fTracklet->GetCovAt(x, cov);
524 // fTracklet->GetCovRef(covR);
525 // cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
526 // Double_t dyz[2]= {dy[1], dz[1]};
527 // Pulls(dyz, cov, tilt);
528 // ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
529 // ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
532 Bool_t rc(fTracklet->IsRowCross());
533 UChar_t err(fTracklet->GetErrorMsg());
534 Double_t x(fTracklet->GetX()),
535 pt(fTracklet->GetPt()),
536 yt(fTracklet->GetYref(0)),
537 zt(fTracklet->GetZref(0)),
538 phi(fTracklet->GetYref(1)),
539 tht(fTracklet->GetZref(1));
540 Int_t ncl(fTracklet->GetN()), det(fTracklet->GetDetector());
541 (*DebugStream()) << "tracklet"
552 <<"dy=" << val[kYrez]
553 <<"dz=" << val[kZrez]
554 <<"dphi="<< val[kPrez]
555 <<"dQ ="<< val[kNdim]
561 if(!track) return NULL;
562 // special care for EVE usage
564 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
565 return H->Projection(kYrez);
569 //________________________________________________________
570 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
572 // Store resolution/pulls of Kalman before updating with the TRD information
573 // at the radial position of the first tracklet. The following points are used
575 // - the (y,z,snp) of the first TRD tracklet
576 // - the (y, z, snp, tgl, pt) of the MC track reference
578 // Additionally the momentum resolution/pulls are calculated for usage in the
581 if(track) fkTrack = track;
583 AliDebug(4, "No Track defined.");
587 TH1 *h(NULL); // EVE projection
589 THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
591 AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
594 // check input track status
595 AliExternalTrackParam *tin(NULL);
596 if(!(tin = fkTrack->GetTrackIn())){
597 AliError("Track did not entered TRD fiducial volume.");
600 // check first tracklet
601 AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
603 AliDebug(3, "No Tracklet in ly[0]. Skip track.");
604 if(!track) return NULL;
605 // special care for EVE usage
606 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
607 return H->Projection(kYrez);
609 if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()){
610 AliDebug(3, "Tracklet or Chamber not OK. Skip track.");
611 if(!track) return NULL;
612 // special care for EVE usage
613 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
614 return H->Projection(kYrez);
616 // check radial position
617 Double_t x = tin->GetX();
618 if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
619 AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
620 if(!track) return NULL;
621 // special care for EVE usage
622 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
623 return H->Projection(kYrez);
626 // down scale PID resolution
627 Int_t spc(fSpecies); if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
630 Float_t tpc(fkESD->GetTPCdedx());
631 Float_t tof(fkESD->GetTOFbeta());
632 (*DebugStream()) << "trackIn"
636 <<"tracklet.=" << fTracklet
641 //Int_t bc(fkESD?fkESD->GetTOFbc()/2:0);
642 const Double_t *parR(tin->GetParameter());
643 Double_t dyt(fTracklet->GetYfit(0)-parR[0]), dzt(fTracklet->GetZfit(0)-parR[1]),
644 phit(fTracklet->GetYfit(1)),
645 tilt(fTracklet->GetTilt()),
646 norm(1./TMath::Sqrt((1.-parR[2])*(1.+parR[2])));
648 // correct for tilt rotation
649 Double_t dy = dyt - dzt*tilt,
651 dx = dy/(parR[2]*norm-parR[3]*norm*tilt);
652 phit += tilt*parR[3];
653 Double_t dphi = TMath::ATan(phit) - TMath::ASin(parR[2]);
655 Double_t val[kNdim+2];
656 val[kBC] = fTriggerSlot;//bc==0?0:(bc<0?-1.:1.);
657 Double_t alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha(),
658 cs = TMath::Cos(alpha),
659 sn = TMath::Sin(alpha);
660 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
661 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
662 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
664 val[kZrez] = fTracklet->IsRowCross()?dz:(fTracklet->GetdQdl()*5.e-4 - 2.5);
665 val[kPrez] = dphi*TMath::RadToDeg();
666 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:spc;
667 val[kPt] = GetPtBin(fPt);
668 val[kNdim] = GetPtBin(fTracklet->GetMomentum());
670 //val[kNdim+2] = fEvent?fEvent->GetBunchFill():0;
673 if(!track) return NULL;
674 // special care for EVE usage
675 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
676 return H->Projection(kYrez);
680 //________________________________________________________
681 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
683 // Store resolution/pulls of Kalman after last update with the TRD information
684 // at the radial position of the first tracklet. The following points are used
686 // - the (y,z,snp) of the first TRD tracklet
687 // - the (y, z, snp, tgl, pt) of the MC track reference
689 // Additionally the momentum resolution/pulls are calculated for usage in the
692 if(track) fkTrack = track;
696 //________________________________________________________
697 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
700 // Plot MC distributions
703 if(!HasMCdata()) return NULL;
704 if(track) fkTrack = track;
706 AliDebug(4, "No Track defined.");
709 Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
713 AliWarning("No output container defined.");
716 // retriev track characteristics
717 Int_t pdg = fkMC->GetPDG(),
718 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
721 label(fkMC->GetLabel());
723 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
724 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
725 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
726 if(sIdx==3) sIdx=2; if(sIdx>3) sIdx=3; // downscale PID
729 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
730 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
732 Double_t x, y, z, pt, dydx, dzdx, dzdl;
733 Float_t pt0, p0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
734 Double_t covR[7]/*, cov[3]*/;
736 /* if(DebugLevel()>=3){
737 // get first detector
739 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
740 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
741 det = fTracklet->GetDetector();
745 TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
747 m = fkTrack->GetMass();
748 if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
749 (*DebugStream()) << "MCkalman"
766 AliTRDcluster *c(NULL);
767 Double_t val[kNdim+2];
768 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
769 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
770 !fTracklet->IsOK())*/ continue;
772 x= x0 = fTracklet->GetX();
773 Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
774 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, p0, eta, phi, s)) continue;
776 // MC track position at reference radial position
778 Float_t ymc = y0 - dx*dydx0;
779 Float_t zmc = z0 - dx*dzdx0;
780 //phi -= TMath::Pi();
785 val[kSpeciesChgRC]= rc?0.:sign*sIdx;
786 val[kPt] = GetPtBin(pt0);
787 val[kNdim]= GetPtBin(p0);
788 Double_t tilt(fTracklet->GetTilt());
790 // ,corr(1./(1. + t2))
791 // ,cost(TMath::Sqrt(corr));
793 AliExternalTrackParam *tin(fkTrack->GetTrackIn());
794 if(ily==0 && tin){ // trackIn residuals
795 // check radial position
796 if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
798 val[kBC] = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
799 val[kYrez] = tin->GetY()-ymc;
800 val[kZrez] = rc?(tin->GetZ()-zmc):(fTracklet->GetdQdl()*5e-4 - 2.5);
801 val[kPrez] = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
802 val[kNdim+1] = 0.;//dx;
803 if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
806 //if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
809 dydx = fTracklet->GetYref(1);
810 dzdx = fTracklet->GetZref(1);
811 dzdl = fTracklet->GetTgl();
812 y = fTracklet->GetYref(0);
814 z = fTracklet->GetZref(0);
816 pt = TMath::Abs(fTracklet->GetPt());
817 fTracklet->GetCovRef(covR);
820 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
822 val[kPt] = GetPtBin(pt0);
823 val[kNdim] = 1.e2*(pt/pt0-1.);
824 if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
825 /* // theta resolution/ tgl pulls
826 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
827 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
828 ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
829 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
830 // pt resolution \\ 1/pt pulls \\ p resolution for PID
831 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
832 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
833 ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
834 ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
835 ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
837 // Fill Debug stream for MC track
839 Int_t det(fTracklet->GetDetector());
840 (*DebugStream()) << "MC"
852 // Fill Debug stream for Kalman track
853 (*DebugStream()) << "MCtrack"
865 // tracklet residuals
866 dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
867 dzdx = fTracklet->GetZfit(1);
868 y = fTracklet->GetYfit(0);
870 z = fTracklet->GetZfit(0);
872 val[kYrez] = dy - dz*tilt;
873 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
874 val[kZrez] = rc?(dz + dy*tilt):(fTracklet->GetdQdl()*3.e-4 - 1.5);
875 val[kNdim] = fEvent?fEvent->GetMultiplicity():0;
876 val[kNdim+1] = 1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb;
877 // val[kNdim] = pt/pt0-1.;
878 if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
881 // Fill Debug stream for tracklet
883 Float_t s2y = fTracklet->GetS2Y();
884 Float_t s2z = fTracklet->GetS2Z();
885 (*DebugStream()) << "MCtracklet"
896 AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
897 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
898 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
900 H = (THnSparse*)fContainer->At(kMCcluster);
901 val[kPt] = TMath::ATan(dydx0)*TMath::RadToDeg();
902 //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
904 Float_t padCorr(tilt*fTracklet->GetPadLength()),
905 corr(1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0));
906 fTracklet->ResetClusterIter(kTRUE);
907 while((c = fTracklet->NextCluster())){
908 if(row0<0) row0 = c->GetPadRow();
909 x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
910 y = c->GetY() + padCorr*(c->GetPadRow() - row0);
917 val[kYrez] = dy - dz*tilt;
919 val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t q(TMath::Abs(c->GetQ()));
920 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
921 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
922 if(ic) val[kZrez] /= (ic*q);
923 val[kSpeciesChgRC]= rc?0.:(TMath::Max(q*corr, Float_t(3.)));
924 val[kNdim+1] = c->IsFivePad()?1:c->GetNPads();
928 // Fill calibration container
929 Float_t d = zr0 - zmc;
930 d -= ((Int_t)(2 * d)) / 2.0;
931 if (d > 0.25) d = 0.5 - d;
932 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
933 clInfo->SetCluster(c);
934 clInfo->SetMC(pdg, label);
935 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
936 clInfo->SetResolution(dy);
937 clInfo->SetAnisochronity(d);
938 clInfo->SetDriftLength(dx);
939 clInfo->SetTilt(tilt);
940 if(fMCcl) fMCcl->Add(clInfo);
941 else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
944 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
945 clInfoArr->SetOwner(kFALSE);
947 clInfoArr->Add(clInfo);
951 if(DebugLevel()>=5 && clInfoArr){
952 (*DebugStream()) << "MCcluster"
953 <<"clInfo.=" << clInfoArr
958 if(clInfoArr) delete clInfoArr;
959 if(!track) return NULL;
960 // special care for EVE usage
961 if(H && (h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
962 return H?H->Projection(kYrez):NULL;
965 //________________________________________________________
966 void AliTRDresolution::GetRangeZ(TH2 *h2, Float_t &min, Float_t &max)
968 // Get range on Z axis such to avoid outliers
970 Double_t cnt[20000], c, m, s;
971 Int_t nx(h2->GetXaxis()->GetNbins()), ny(h2->GetYaxis()->GetNbins()), nc(0);
972 for(Int_t ix(1); ix<=nx; ix++){
973 for(Int_t iy(1); iy<=ny; iy++){
974 if((c = h2->GetBinContent(ix, iy))<10) continue;
980 AliMathBase::EvaluateUni(nc, cnt, m, s, 0);
981 min = m-s; max = m+2.*s;
984 //________________________________________________________
985 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
988 // Get the reference figures
992 AliWarning("Please provide a canvas to draw results.");
995 /* Int_t selection[100], n(0), selStart(0); //
996 Int_t ly0(0), dly(5);
997 TList *l(NULL); TVirtualPad *pad(NULL); */
1002 AliWarning(Form("Reference plot [%d] missing result", ifig));
1006 //________________________________________________________
1007 void AliTRDresolution::MakeSummary()
1009 // Build summary plots
1012 AliError("Missing results");
1015 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
1016 TObjArray *arr(NULL); TH2 *h2(NULL);
1018 // cluster resolution
1020 gStyle->SetPalette(1);
1021 const Int_t nClViews(9);
1022 const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh"};
1023 const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1};
1024 const Float_t vClSignMin[2] = {2.6e2, 4.4e2},
1025 vClSignMax[2] = {4.4e2, 6.2e2},
1026 vClMin[nClViews] = {3.2e2, vClSignMin[Int_t(fBsign)], vClSignMin[Int_t(!fBsign)], 0., 0., 0., 0., 0., 3.2e2},
1027 vClMax[nClViews] = {5.e2, vClSignMax[Int_t(fBsign)], vClSignMax[Int_t(!fBsign)], 0., 0., 0., 0., 0., 5.e2};
1028 const Int_t nTrkltViews(20);
1029 const Char_t *vTrkltName[nTrkltViews] = {
1030 "TrkltY", "TrkltYn", "TrkltYp", "TrkltY", "TrkltYn", "TrkltYp",
1032 "TrkltPh", "TrkltPhn", "TrkltPhp",
1033 "TrkltQ", "TrkltQn", "TrkltQp",
1034 "TrkltQS", "TrkltQSn", "TrkltQSp",
1035 "TrkltTBn", "TrkltTBp", "TrkltTBn", "TrkltTBp"
1036 // "TrkltYnl0", "TrkltYpl0", "TrkltPhnl0", "TrkltPhpl0", "TrkltQnl0", "TrkltQpl0", // electrons low pt
1037 /* "TrkltYnl1", "TrkltYpl1", "TrkltPhnl1", "TrkltPhpl1", "TrkltQnl1", "TrkltQpl1", // muons low pt
1038 "TrkltYnl2", "TrkltYpl2", "TrkltPhnl2", "TrkltPhpl2", "TrkltQnl2", "TrkltQpl2" // pions low pt*/
1040 const UChar_t vTrkltOpt[nTrkltViews] = {
1048 const Int_t nTrkInViews(5);
1049 const Char_t *vTrkInName[nTrkInViews][6] = {
1050 {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1051 {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1052 {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"},
1053 {"TrkInXnl", "TrkInXpl", "TrkInXl", "TrkInYnh", "TrkInYph", "TrkInYh"},
1054 {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"},
1055 //{"TrkInRCX", "TrkInRCY", "TrkInRCPh", "TrkInRCZl", "TrkInRCZi", "TrkInRCZh"}
1057 const UChar_t vTrkInOpt[nTrkInViews] = {0, 1, 0, 0, 0};
1058 const Float_t min[6] = {0.15, 0.15, 0.15, 0.15, 0.5, 0.5};
1059 const Float_t max[6] = {0.6, 0.6, 0.6, 0.6, 2.3, 2.3};
1060 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]"};
1062 const Int_t nTrkViews(27);
1063 const Char_t *vTrkName[nTrkViews] = {
1064 "TrkY", "TrkYn", "TrkYp",
1065 "TrkPh", "TrkPhn", "TrkPhp",
1066 "TrkDPt", "TrkDPtn", "TrkDPtp",
1067 "TrkYnl0", "TrkYpl0", "TrkPhnl0", "TrkPhpl0", "TrkDPtnl0", "TrkDPtpl0", // electrons low pt
1068 "TrkYnl1", "TrkYpl1", "TrkPhnl1", "TrkPhpl1", "TrkDPtnl1", "TrkDPtpl1", // muons low pt
1069 "TrkYnl2", "TrkYpl2", "TrkPhnl2", "TrkPhpl2", "TrkDPtnl2", "TrkDPtpl2" // pions low pt
1071 const Char_t *typName[] = {"", "MC"};
1072 const Int_t nx(2048), ny(1536);
1074 if((arr = (TObjArray*)fProj->FindObject("hDet2Cluster"))){
1075 cOut = new TCanvas(Form("%s_DetOccupancy", GetName()), "Detector performance", 2*nx, 2*ny);
1076 cOut->Divide(AliTRDgeometry::kNlayer,AliTRDeventInfo::kCentralityClasses, 1.e-5, 1.e-5);
1078 for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1079 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1080 p=cOut->cd(icen*AliTRDgeometry::kNlayer+ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1081 if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d%dEn", ily, icen)))) continue;
1082 if(SetNormZ(h2, 1, -1, 1, -1, 10.) < 1.e3) continue; // cut all bins with < 10 entries
1083 SetRangeZ(h2, -90., 90, -200.);
1084 //h2->GetXaxis()->SetNdivisions(1010);
1085 h2->GetZaxis()->SetTitle("Rel. Det. Occup. [%]");
1086 h2->GetZaxis()->CenterTitle();
1087 h2->SetContour(9); h2->Draw("colz"); n++;
1088 MakeDetectorPlot(ily, "pad");
1091 if(n>=AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1093 cOut = new TCanvas(Form("%s_DetCharge", GetName()), "Detector performance", nx, ny);
1094 cOut->Divide(3,2, 1.e-5, 1.e-5);
1095 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1096 p=cOut->cd(ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1097 if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d_2D", ily)))) continue;
1098 SetNormZ(h2, 1, -1, 1, -1, 10.); // cut all <q> < 10
1099 SetRangeZ(h2, -30., 30., -200.);
1100 //h2->GetXaxis()->SetNdivisions(1010);
1101 h2->GetZaxis()->SetTitle("Rel. Mean(q) [%]");
1102 h2->GetZaxis()->CenterTitle();
1104 MakeDetectorPlot(ily, "pad");
1106 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1108 for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
1109 if((arr = (TObjArray*)fProj->FindObject(ityp?"hCluster2MC":"hCluster2Track"))){
1110 for(Int_t iview(0); iview<nClViews; iview++){
1111 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vClName[iview], vClOpt[iview]), "Cluster Resolution", nx, ny);
1112 cOut->Divide(3,2, 1.e-5, 1.e-5);
1114 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1115 p=cOut->cd(ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1116 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], ily)))) continue;
1118 if(vClOpt[iview]==0) h2->Draw("colz");
1119 else if(vClOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [#mum]", vClMin[iview], vClMax[iview], 1.e4);
1120 if(iview<5) MakeDetectorPlot(ily);
1122 if(nplot==AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1126 // tracklet systematic
1127 if((arr = (TObjArray*)fProj->FindObject(ityp?"hTracklet2MC":"hTracklet2Track"))){
1128 for(Int_t iview(0); iview<nTrkltViews; iview++){
1129 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkltName[iview], vTrkltOpt[iview]), "Tracklet Resolution", nx, ny);
1130 cOut->Divide(3,2, 1.e-5, 1.e-5);
1132 for(Int_t iplot(0); iplot<6; iplot++){
1133 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1134 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
1136 if(vTrkltOpt[iview]==0) h2->Draw("colz");
1137 else if (vTrkltOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [cm]", .15, .4);
1138 else if (vTrkltOpt[iview]==2) DrawSigma(h2, "#sigma(occupancy) [%]", 10.5, 15.);
1139 MakeDetectorPlot(iplot);
1142 cOut->Modified();cOut->Update();
1143 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1148 // trackIn systematic
1149 if((arr = (TObjArray*)fProj->FindObject(ityp?"hTRDin2MC":"hTracklet2TRDin"))){
1150 for(Int_t iview(0); iview<nTrkInViews; iview++){
1151 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkInName[iview][0], vTrkInOpt[iview]), "Track IN Resolution", nx, ny);
1152 cOut->Divide(3,2, 1.e-5, 1.e-5);
1154 for(Int_t iplot(0); iplot<6; iplot++){
1155 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1156 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])))){
1157 AliInfo(Form("Missing H%s%s_2D", typName[ityp], vTrkInName[iview][iplot]));
1161 if(vTrkInOpt[iview]==0) h2->Draw("colz");
1162 else DrawSigma(h2, ttt[iplot], min[iplot], max[iplot]);
1163 MakeDetectorPlot(0);
1165 if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1169 const Float_t zmin[] = {1., 61., 15.},
1170 zmax[] = {10., 79., 33.};
1171 cOut = new TCanvas(Form("%s_%sTrkInSpc", GetName(), typName[ityp]), "Track IN PID", Int_t(1.5*ny), Int_t(1.5*ny));
1172 cOut->Divide(3,3, 1.e-5, 1.e-5);
1173 Int_t nplot(0); const Char_t *chName[] = {"p", "n", ""};
1174 for(Int_t ich(0), ipad(1); ich<3; ich++){
1176 if(!(h2s = (TH2*)arr->FindObject(Form("H%sTrkInY%sEn", typName[ityp], chName[ich])))) {
1177 AliInfo(Form("Missing H%sTrkIn%sEn", typName[ityp], chName[ich]));
1180 Int_t irebin(0), dxBin(1), dyBin(1);
1181 const Int_t nrebin(5); Int_t rebinX[nrebin] = {1, 2, 1, 2, 1}, rebinY[nrebin] = {2, 1, 2, 1, 2};
1182 if(h2s->GetNbinsY()==180){ // old binning
1183 rebinX[0] = 1; rebinY[0] = 2;
1184 rebinX[1] = 2; rebinY[1] = 1;
1185 rebinX[2] = 2; rebinY[2] = 1;
1186 rebinX[3] = 1; rebinY[3] = 5;
1187 rebinX[4] = 1; rebinY[4] = 1; // dummy
1189 while(irebin<nrebin && (AliTRDrecoTask::GetMeanStat(h2s, .5, ">")<200)){
1190 h2s->Rebin2D(rebinX[irebin], rebinY[irebin]);
1191 dxBin*=rebinX[irebin];dyBin*=rebinY[irebin];irebin++;
1193 AliDebug(2, Form("Rebin level[%d] @ chg[%d]. Binning[%2dx%2d]", irebin, ich, h2s->GetNbinsX(), h2s->GetNbinsY()));
1194 for(Int_t ispec(0); ispec<kNspc; ispec++){
1195 p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1196 if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkInY%s%dEn", typName[ityp], chName[ich], ispec)))) {
1197 AliInfo(Form("Missing H%sTrkIn%s%dEn", typName[ityp], chName[ich], ispec));
1201 h2->Rebin2D(dxBin,dyBin);
1202 h2->Divide(h2, h2s, 1.e2);
1204 h2->GetZaxis()->SetRangeUser(zmin[ispec], zmax[ispec]);
1205 h2->GetZaxis()->SetTitle("Rel. Abundancy [%]");h2->GetZaxis()->CenterTitle();
1206 TString tit(h2->GetTitle()); TObjArray *atit = tit.Tokenize("::");
1207 h2->SetTitle(Form("%s :: Relative Abundancy", ((*atit)[0])->GetName()));
1208 atit->Delete(); delete atit;
1210 MakeDetectorPlot(0);
1213 if(nplot==9) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1216 const char *chQ[] = {"Q", "QS"};
1217 for(Int_t iq(0); iq<2; iq++){
1218 cOut = new TCanvas(Form("%s_%sTrkIn%s", GetName(), typName[ityp], chQ[iq]), "Track IN PID", Int_t(1.5*ny), Int_t(1.5*ny));
1219 cOut->Divide(3,3, 1.e-5, 1.e-5);
1221 for(Int_t ich(0), ipad(1); ich<3; ich++){
1222 for(Int_t ispec(0); ispec<kNspc; ispec++){
1223 p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1224 if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkIn%s%s%d_2D", typName[ityp], chQ[iq], chName[ich], ispec)))) {
1225 AliInfo(Form("Missing H%sTrkIn%s%s%d_2D", typName[ityp], chQ[iq], chName[ich], ispec));
1230 MakeDetectorPlot(0);
1233 if(nplot==9) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1238 // track MC systematic
1239 if((arr = (TObjArray*)fProj->FindObject("hTRD2MC"))) {
1240 for(Int_t iview(0); iview<nTrkViews; iview++){
1241 cOut = new TCanvas(Form("%s_MC%s", GetName(), vTrkName[iview]), "Track Resolution", nx, ny);
1242 cOut->Divide(3,2, 1.e-5, 1.e-5);
1244 for(Int_t iplot(0); iplot<6; iplot++){
1245 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1246 if(!(h2 = (TH2*)arr->FindObject(Form("HMC%s%d_2D", vTrkName[iview], iplot)))) continue;
1247 h2->Draw("colz"); nplot++;
1249 if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1255 gStyle->SetPalette(1);
1258 //________________________________________________________
1259 void AliTRDresolution::DrawSigma(TH2 *h2, const Char_t *title, Float_t m, Float_t M, Float_t scale)
1261 // Draw error bars scaled with "scale" instead of content values
1262 //use range [m,M] if limits are specified
1265 TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1266 TH2F *h2e = new TH2F(Form("%s_E", h2->GetName()),
1267 Form("%s;%s;%s;%s", h2->GetTitle(), ax->GetTitle(), ay->GetTitle(), title),
1268 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
1269 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1271 TAxis *az(h2e->GetZaxis());
1272 if(M>m) az->SetRangeUser(m, M);
1274 az->SetTitleOffset(1.5);
1275 for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1276 for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1277 if(h2->GetBinContent(ix, iy)<-100.) continue;
1278 Float_t v(scale*h2->GetBinError(ix, iy));
1279 if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1280 h2e->SetBinContent(ix, iy, v);
1286 //________________________________________________________
1287 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1289 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1290 // The "range" vector should be initialized with 2 elements
1291 // Option "mod" can be any of
1292 // - 0 : gaussian like distribution
1293 // - 1 : tailed distribution
1295 Int_t nx(h2->GetNbinsX())
1296 , ny(h2->GetNbinsY())
1298 Double_t *data=new Double_t[n];
1299 for(Int_t ix(1), in(0); ix<=nx; ix++){
1300 for(Int_t iy(1); iy<=ny; iy++)
1301 data[in++] = h2->GetBinContent(ix, iy);
1303 Double_t mean, sigm;
1304 AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1306 range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1307 if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1308 AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1309 TH1S h1("h1SF0", "", 100, range[0], range[1]);
1314 case 0:// gaussian distribution
1316 TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1318 mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1319 range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1320 AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1]));
1323 case 1:// tailed distribution
1325 Int_t bmax(h1.GetMaximumBin());
1326 Int_t jBinMin(1), jBinMax(100);
1327 for(Int_t ibin(bmax); ibin--;){
1328 if(h1.GetBinContent(ibin)<1.){
1329 jBinMin=ibin; break;
1332 for(Int_t ibin(bmax); ibin++;){
1333 if(h1.GetBinContent(ibin)<1.){
1334 jBinMax=ibin; break;
1337 range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1338 AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1]));
1346 //________________________________________________________
1347 Bool_t AliTRDresolution::MakeProjectionDetector()
1350 const Int_t kNcontours(9);
1351 const Int_t kNstat(100);
1352 if(fProj && fProj->At(kDetector)) return kTRUE;
1354 AliError("Missing data container.");
1358 if(!(H = (THnSparse*)fContainer->FindObject("hDet2Cluster"))){
1359 AliInfo(Form("Missing/Wrong data @ hDet2Cluster."));
1362 Int_t ndim(H->GetNdimensions());
1363 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1364 TAxis *aa[kNdim], *an(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1365 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1366 if(ndim < 5) aa[4] = new TAxis(1, -0.5, 0.5);
1372 // build list of projections
1373 const Int_t nsel(8*AliTRDgeometry::kNlayer*AliTRDeventInfo::kCentralityClasses);
1374 // define rebinning strategy
1375 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 1, 2, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1376 if(aa[1]->GetNbins()==180){ // old binning
1377 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1378 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1379 rebinEtaPhiX[2] = 2; rebinEtaPhiY[2] = 1;
1380 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1381 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1383 //const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1384 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1385 AliTRDrecoProjection hp[kDetNproj]; TObjArray php(kDetNproj);
1386 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1387 for(Int_t ipad(0); ipad<nPad; ipad++){
1388 for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1389 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1390 isel++; // new selection
1391 hp[ih].Build(Form("HDet%d%d%d", ily, icen, ipad),
1392 Form("Detectors :: Ly[%d] Cen[%s] Pads[%s]", ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1394 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1395 hp[ih].SetShowRange(10., 55.);
1396 php.AddLast(&hp[ih++]); np[isel]++;
1400 AliInfo(Form("Build %3d 3D projections.", ih));
1402 Int_t ly(0), cen(0), npad(0);
1403 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1404 v = H->GetBinContent(ib, coord); if(v<1.) continue;
1405 ly = coord[kBC]-1; // layer selection
1406 cen = coord[kYrez]-1; // centrality selection
1407 npad = 0; // no. of pads selection
1408 if(an) npad = TMath::Min(coord[5]-1, Int_t(kNpads));
1409 isel = npad*AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer+cen*AliTRDgeometry::kNlayer+ly;
1410 ((AliTRDrecoProjection*)php.At(isel))->Increment(coord, v);
1411 //Int_t ioff=isel;for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1413 TObjArray *arr(NULL);
1414 fProj->AddAt(arr = new TObjArray(kDetNproj), kDetector);
1415 arr->SetName("hDet2Cluster");
1417 TH2 *h2(NULL); Int_t jh(0);
1419 if(!hp[ih].H()) continue;
1420 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1421 arr->AddAt(h2, jh++);
1422 if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].H()->GetName())))) continue;
1423 arr->AddAt(h2, jh++);
1425 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1426 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1427 for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1428 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, 0)))){
1429 for(Int_t ipad(1); ipad<nPad; ipad++){
1430 if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, ipad)))){
1434 pr0->H()->SetNameTitle(Form("HDet%d%d", ily, icen), Form("Detectors :: Ly[%d] Cen[%s]", ily, cenName[icen]));
1435 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0, kFALSE))) arr->AddAt(h2, jh++);
1436 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
1437 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))) (*pr1)+=(*pr0);
1440 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))){
1441 pr0->H()->SetNameTitle(Form("HDet%d", ily), Form("Detectors :: Ly[%d]", ily));
1442 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1445 AliInfo(Form("Done %3d 2D projections.", jh));
1449 //________________________________________________________
1450 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1453 const Int_t kNcontours(9);
1454 const Int_t kNstat(100);
1455 Int_t cidx=mc?kMCcluster:kCluster;
1456 if(fProj && fProj->At(cidx)) return kTRUE;
1458 AliError("Missing data container.");
1461 const Char_t *projName[] = {"hCluster2Track", "hCluster2MC"};
1463 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1464 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1467 Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimCl));
1468 Int_t coord[10]; memset(coord, 0, sizeof(Int_t) * 10); Double_t v = 0.;
1469 TAxis *aa[kNdim], *as(NULL), *apt(NULL), *acen(NULL), *anp(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1470 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1471 if(ndim > Int_t(kPt)) apt = H->GetAxis(kPt);
1472 if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC);
1473 if(ndim > Int_t(kNdim)) acen = H->GetAxis(kNdim);
1474 if(ndim > Int_t(kNdim)+1) anp = H->GetAxis(kNdim+1);
1475 // calculate size depending on debug level
1476 const Int_t nCh(apt?2:1);
1477 const Int_t nCen(acen?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1478 const Int_t nNpad(anp?(Int_t(kNpads)+1):1);
1480 // build list of projections
1481 const Int_t nsel(AliTRDgeometry::kNlayer*kNcharge*(kNpads+1)*AliTRDeventInfo::kCentralityClasses);
1482 // define rebinning strategy
1483 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1484 if(aa[1]->GetNbins()==180){ // old binning
1485 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1486 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1487 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
1488 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1489 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1491 AliTRDrecoProjection hp[kClNproj]; TObjArray php(kClNproj);
1492 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1493 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1494 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1495 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1496 for(Int_t ich(0); ich<nCh; ich++){
1497 for(Int_t icen(0); icen<nCen; icen++){
1498 for(Int_t ipad(0); ipad<nNpad; ipad++){
1499 isel++; // new selection
1500 hp[ih].Build(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1501 Form("Clusters[%c] :: #Deltay Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1502 kEta, kPhi, kYrez, aa);
1503 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1504 php.AddLast(&hp[ih++]); np[isel]++;
1506 hp[ih].Build(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1507 Form("Clusters[%c] :: Q Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1508 kEta, kPhi, kSpeciesChgRC, aa);
1509 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1510 hp[ih].SetShowRange(24., 33.);
1511 php.AddLast(&hp[ih++]); np[isel]++;
1512 hp[ih].Build(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1513 Form("Clusters[%c] :: #Deltay(x,TC) Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1514 kPrez, kZrez, kYrez, aa);
1515 php.AddLast(&hp[ih++]); np[isel]++;
1516 hp[ih].Build(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1517 Form("Clusters[%c] :: #Deltay(x,#Phi) Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1518 kPrez, kPt, kYrez, aa);
1519 php.AddLast(&hp[ih++]); np[isel]++;
1524 AliInfo(Form("Build %3d 3D %s projections.", ih, mc?"MC":""));
1526 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1527 Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1), ioff(0), cen(0), npad(0);
1528 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1529 v = H->GetBinContent(ib, coord); if(v<1.) continue;
1532 if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1535 ch = 0; // [-] track
1536 if(chBin>0 && coord[kPt] > chBin) ch = 1; // [+] track
1537 cen = 0; // centrality selection
1538 if(acen) cen = coord[kNdim]-1;
1539 npad = 0; // no. of pads selection
1540 if(anp) npad = TMath::Min(coord[kNdim+1]-1, Int_t(kNpads));
1543 isel = ly*nCh*nCen*nNpad
1549 isel = ly; ioff = isel;
1552 AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1555 if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
1556 AliError(Form("Missing projection %d", ioff));
1559 if(strcmp(pr0->H()->GetName(), Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ch], ly, cen, npad))!=0){
1560 AliError(Form("Projection mismatch :: request[H%sClY%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ly, cen, npad, pr0->H()->GetName()));
1563 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1565 if(HasDump3DFor(kCluster)){
1566 TDirectory *cwd = gDirectory;
1567 TFile::Open(Form("ResDump_%s.root", H->GetName()), "RECREATE");
1568 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
1569 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
1570 if(!pr0->H()) continue;
1571 TH3 *h3=(TH3*)pr0->H()->Clone();
1578 TObjArray *arr(NULL);
1579 fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1580 arr->SetName(projName[Int_t(mc)]);
1582 TH2 *h2(NULL); Int_t jh(0);
1584 if(!hp[ih].H()) continue;
1585 if(strchr(hp[ih].H()->GetName(), 'Q')){
1586 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1587 arr->AddAt(h2, jh++);
1588 if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].H()->GetName())))) continue;
1589 arr->AddAt(h2, jh++);
1591 if((h2 = hp[ih].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1594 Double_t m(0.), s(0.), trend(0.);
1595 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1596 for(Int_t ich(0); ich<nCh; ich++){
1597 for(Int_t icen(0); icen<nCen; icen++){
1599 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1600 for(Int_t ipad(1); ipad<nNpad; ipad++){
1601 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1604 pr0->H()->SetNameTitle(Form("H%sClY%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1605 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1606 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1609 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1610 for(Int_t ipad(1); ipad<nNpad; ipad++){
1611 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1614 pr0->H()->SetNameTitle(Form("H%sClQ%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: Q Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1615 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2, kFALSE))) arr->AddAt(h2, jh++);
1616 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
1617 pr0->H()->SetName(Form("H%sClQS%c%d%d", mc?"MC":"", chName[ich], ily, icen));
1618 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1619 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1622 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1623 for(Int_t ipad(1); ipad<nNpad; ipad++){
1624 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1627 pr0->H()->SetNameTitle(Form("H%sClYXTC%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay(x,TC) Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1628 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1629 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1632 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1633 for(Int_t ipad(1); ipad<nNpad; ipad++){
1634 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1637 pr0->H()->SetNameTitle(Form("H%sClYXPh%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay(x,#Phi) Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1638 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1639 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1641 } // end centrality integration
1643 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1644 pr0->H()->SetNameTitle(Form("H%sClY%c%d", mc?"MC":"", chName[ich], ily),
1645 Form("Clusters[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
1646 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1647 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1650 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1651 pr0->H()->SetNameTitle(Form("H%sClQ%c%d", mc?"MC":"", chName[ich], ily),
1652 Form("Clusters[%c]:: Q Ly[%d]", chSgn[ich], ily));
1653 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
1654 pr0->H()->SetName(Form("H%sClQS%c%d", mc?"MC":"", chName[ich], ily));
1655 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1656 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1659 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1660 pr0->H()->SetNameTitle(Form("H%sClYXTC%c%d", mc?"MC":"", chName[ich], ily),
1661 Form("Clusters[%c]:: #Deltay(x,TC) Ly[%d]", chSgn[ich], ily));
1662 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1663 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1666 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1667 pr0->H()->SetNameTitle(Form("H%sClYXPh%c%d", mc?"MC":"", chName[ich], ily),
1668 Form("Clusters[%c]:: #Deltay(x,#Phi) Ly[%d]", chSgn[ich], ily));
1669 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1670 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1672 } // end charge integration
1674 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1675 pr0->H()->SetNameTitle(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1676 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1677 if((trend=pr0->GetTrendValue(1,&m,&s))>-100.) PutTrendValue(Form("%sClS%d", mc?"MC":"", ily), s);
1680 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1681 pr0->H()->SetNameTitle(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1682 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1686 AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
1690 //________________________________________________________
1691 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1694 const Int_t kNcontours(9);
1695 const Int_t kNstat(30);
1696 const Int_t kNstatQ(30);
1697 Int_t cidx=mc?kMCtracklet:kTracklet;
1698 if(fProj && fProj->At(cidx)) return kTRUE;
1700 AliError("Missing data container.");
1703 const Char_t *projName[] = {"hTracklet2Track", "hTracklet2MC"};
1705 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1706 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1709 const Int_t mdim(kNdim+8);
1710 Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimTrklt));
1711 Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
1712 TAxis *aa[mdim], *as(NULL), *ap(NULL), *ac(NULL); memset(aa, 0, sizeof(TAxis*) * mdim);
1713 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1714 if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC); // init species/charge selection
1715 if(ndim > Int_t(kPt)) ap = H->GetAxis(kPt); // init pt selection
1716 if(ndim > Int_t(kNdim)) ac = H->GetAxis(kNdim); // init centrality selection
1717 // calculate size depending on debug level
1718 const Int_t nCen(debug?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1719 const Int_t nPt(debug?Int_t(kNpt+2):1);
1720 const Int_t nSpc(1);//ndim>kNdimTrklt?fgkNbins[kSpeciesChgRC]:1);
1721 const Int_t nCh(debug?Int_t(kNcharge):1);
1723 // build list of projections
1724 const Int_t nsel(AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer*(kNpt+2)*(AliPID::kSPECIES*kNcharge + 1));
1725 // define rebinning strategy
1726 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1727 if(aa[1]->GetNbins()==180){ // old binning
1728 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1729 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1730 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
1731 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1732 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1734 AliTRDrecoProjection hp[kTrkltNproj]; TObjArray php(kTrkltNproj);
1735 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1736 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1737 const Char_t ptName[kNpt+2] = {'L', 'l', 'i', 'h', 'H'};
1738 const Char_t *ptCut[kNpt+2] = {"p_{t}[GeV/c]<0.5", "0.5<=p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "1.5<=p_{t}[GeV/c]<5.0", "p_{t}[GeV/c]>=5.0"};
1739 // const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1740 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1741 for(Int_t icen(0); icen<nCen; icen++){
1742 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1743 for(Int_t ipt(0); ipt<nPt; ipt++){
1744 for(Int_t isp(0); isp<nSpc; isp++){
1745 for(Int_t ich(0); ich<nCh; ich++){
1746 isel++; // new selection
1747 hp[ih].Build(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1748 Form("Tracklets[%s%c]:: #Deltay{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1749 kEta, kPhi, kYrez, aa);
1750 //hp[ih].SetShowRange(-0.1,0.1);
1751 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1752 php.AddLast(&hp[ih++]); np[isel]++;
1753 hp[ih].Build(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1754 Form("Tracklets[%s%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1755 kEta, kPhi, kPrez, aa);
1756 //hp[ih].SetShowRange(-0.5,0.5);
1757 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1758 php.AddLast(&hp[ih++]); np[isel]++;
1759 hp[ih].Build(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1760 Form("Tracklets[%s%c]:: dQdl{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1761 kEta, kPhi, kZrez, aa);
1762 hp[ih].SetShowRange(1.,2.3);
1763 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1764 php.AddLast(&hp[ih++]); np[isel]++;
1765 hp[ih].Build(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1766 Form("Tracklets[%s%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1767 kEta, kPhi, kNdim+1, aa);
1768 hp[ih].SetShowRange(30., 70.);
1769 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1770 php.AddLast(&hp[ih++]); np[isel]++;
1773 if(ndim==kNdimTrklt) continue;
1775 isel++; // new selection
1776 hp[ih].Build(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1777 Form("Tracklets[RC]:: #Deltaz{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1778 kEta, kPhi, kZrez, aa);
1779 // hp[ih].SetShowRange(-0.1,0.1);
1780 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1781 php.AddLast(&hp[ih++]); np[isel]++;
1782 /* hp[ih].Build(Form("H%sTrkltRCY%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1783 Form("Tracklets[RC]:: #Deltay{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1784 kEta, kPhi, kYrez, aa);
1785 //hp[ih].SetShowRange(-0.1,0.1);
1786 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1787 php.AddLast(&hp[ih++]); np[isel]++;
1788 hp[ih].Build(Form("H%sTrkltRCPh%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1789 Form("Tracklets[RC]:: #Delta#phi{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1790 kEta, kPhi, kPrez, aa);
1791 //hp[ih].SetShowRange(-0.1,0.1);
1792 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1793 php.AddLast(&hp[ih++]); np[isel]++;
1794 hp[ih].Build(Form("H%sTrkltRCQ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1795 Form("Tracklets[RC]:: dQdl{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1796 kEta, kPhi, kNdim, aa);
1797 //hp[ih].SetShowRange(-0.1,0.1);
1798 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1799 php.AddLast(&hp[ih++]); np[isel]++;*/
1803 AliInfo(Form("Build %3d 3D %s projections.", ih, mc?"MC":""));
1805 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1806 Int_t ly(0), ch(0), sp(2), rcBin(as?as->FindBin(0.):-1), pt(0), cen(0), ioff(0), jspc(nSpc*nCh+1), kspc(nSpc*nCh*4+1);
1807 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1808 v = H->GetBinContent(ib, coord);
1810 ly = coord[kBC]-1; // layer selection
1812 ch = 0; sp=0;// [e-] track [dafault]
1813 if(rcBin>0){ // debug mode in which species/charge are also saved
1814 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1815 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
1816 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
1820 if(ap) pt = TMath::Min(coord[kPt], Int_t(kNpt)+1);
1821 // centrality selection
1823 if(ac) cen = coord[kNdim]-1;
1825 if(ndim==kNdimTrklt){
1829 isel = cen*AliTRDgeometry::kNlayer*nPt*jspc+ly*nPt*jspc+pt*jspc; isel+=sp<0?(nSpc*nCh):ch;
1830 ioff = cen*AliTRDgeometry::kNlayer*nPt*kspc+ly*nPt*kspc+pt*kspc; ioff+=sp<0?((nSpc*nCh)*4):(ch*4);
1833 AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1836 if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
1837 AliError(Form("Missing projection %d", ioff));
1841 if(strcmp(pr0->H()->GetName(), Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ch], ptName[pt], sp, ly, cen))!=0){
1842 AliError(Form("Projection mismatch :: request[H%sTrkltY%c%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, ly, cen, pr0->H()->GetName()));
1846 if(strcmp(pr0->H()->GetName(), Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[pt], ly, cen))!=0){
1847 AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c%d%d] found[%s]", mc?"MC":"", ptName[pt], ly, cen, pr0->H()->GetName()));
1851 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1853 if(HasDump3DFor(kTracklet)){
1854 TDirectory *cwd = gDirectory;
1855 TFile::Open(Form("ResDump_%s.root", H->GetName()), "RECREATE");
1856 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
1857 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
1858 if(!pr0->H()) continue;
1859 TH3 *h3=(TH3*)pr0->H()->Clone();
1866 TObjArray *arr(NULL);
1867 fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1868 arr->SetName(projName[Int_t(mc)]);
1870 TH2 *h2(NULL); Int_t jh(0);
1872 if(!hp[ih].H()) continue;
1873 Int_t mid(0), nstat(kNstat);
1874 if(strchr(hp[ih].H()->GetName(), 'Q')){ mid=2; nstat=kNstatQ;}
1875 if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1876 arr->AddAt(h2, jh++);
1878 // build combined performance plots
1879 Double_t m(0.), s(0.), trend(0.);
1880 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1881 for(Int_t ich(0); ich<nCh; ich++){
1882 for(Int_t icen(0); icen<nCen; icen++){
1883 for(Int_t ipt(0); ipt<nPt; ipt++){
1885 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1886 for(Int_t isp(1); isp<nSpc; isp++){
1887 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1890 pr0->H()->SetNameTitle(Form("H%sTrkltY%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1891 Form("Tracklets[%c]:: #Deltay{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1892 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1893 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1896 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1897 for(Int_t isp(1); isp<nSpc; isp++){
1898 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1901 pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1902 Form("Tracklets[%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1903 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1904 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1907 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1908 for(Int_t isp(1); isp<nSpc; isp++){
1909 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1912 pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1913 Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1914 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1915 pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1916 Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1917 pr0->SetShowRange(2.4, 5.1);
1918 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1919 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1922 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1923 for(Int_t isp(1); isp<nSpc; isp++){
1924 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1927 pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1928 Form("Tracklets[%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1929 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1930 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1932 } // end pt integration
1934 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1935 pr0->H()->SetNameTitle(Form("H%sTrkltY%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1936 Form("Tracklets[%c]:: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1937 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1938 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1941 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1942 pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1943 Form("Tracklets[%c]:: #Delta#phi Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1944 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1945 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1948 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1949 pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1950 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1951 pr0->SetShowRange(1.,2.3);
1952 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1953 pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1954 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1955 pr0->SetShowRange(2.4,5.1);
1956 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1957 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1960 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1961 pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1962 Form("Tracklets[%c]:: OccupancyTB Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1963 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1964 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1966 } // end centrality integration
1968 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1969 pr0->H()->SetNameTitle(Form("H%sTrkltY%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Deltay Ly[%d]", chSgn[ich], ily));
1970 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1971 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1974 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1975 pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Delta#phi Ly[%d]", chSgn[ich], ily));
1976 pr0->SetShowRange(-.9,.9);
1977 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1978 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1981 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1982 pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1983 pr0->SetShowRange(1.,2.3);
1984 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1985 pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1986 pr0->SetShowRange(2.4,5.1);
1987 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1988 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1991 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1992 pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: OccupancyTB Ly[%d]", chSgn[ich], ily));
1993 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1994 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1996 } // end charge integration
1998 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1999 pr0->H()->SetNameTitle(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets :: #Deltay Ly[%d]", ily));
2000 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2001 if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2002 PutTrendValue(Form("%sTrkltY%d", mc?"MC":"", ily), trend);
2003 PutTrendValue(Form("%sTrkltYS%d", mc?"MC":"", ily), s);
2007 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
2008 pr0->H()->SetNameTitle(Form("H%sTrkltPh%d", mc?"MC":"", ily), Form("Tracklets :: #Delta#phi Ly[%d]", ily));
2009 pr0->SetShowRange(-.45,.45);
2010 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2013 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
2014 pr0->H()->SetNameTitle(Form("H%sTrkltQ%d", mc?"MC":"", ily), Form("Tracklets :: dQdl Ly[%d]", ily));
2015 pr0->SetShowRange(1.,2.3);
2016 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2017 pr0->H()->SetName(Form("H%sTrkltQS%d", mc?"MC":"", ily));
2018 pr0->SetShowRange(2.4,5.1);
2019 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2022 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
2023 pr0->H()->SetNameTitle(Form("H%sTrkltTB%d", mc?"MC":"", ily), Form("Tracklets :: OccupancyTB Ly[%d]", ily));
2024 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
2027 /*! Row Cross processing*/
2028 for(Int_t icen(0); icen<nCen; icen++){
2030 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, icen)))){
2031 for(Int_t ipt(0); ipt<kNpt; ipt++){
2032 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen)))) continue;
2035 pr0->H()->SetNameTitle(Form("H%sTrkltRCZ%d%d", mc?"MC":"", ily, icen), Form("Tracklets[RC]:: #Deltaz Ly[%d] Cen[%s]", ily, cenName[icen]));
2036 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2037 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))) (*pr1)+=(*pr0);
2039 } // end centrality integration for row cross
2041 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))){
2042 pr0->H()->SetNameTitle(Form("H%sTrkltRCZ%d", mc?"MC":"", ily), Form("Tracklets[RC] :: #Deltaz Ly[%d]", ily));
2043 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2046 AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
2050 //________________________________________________________
2051 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
2054 const Int_t kNcontours(9);
2055 const Int_t kNstat(30);
2056 Int_t cidx=mc?kMCtrackIn:kTrackIn;
2057 if(fProj && fProj->At(cidx)) return kTRUE;
2059 AliError("Missing data container.");
2062 const Char_t *projName[] = {"hTracklet2TRDin", "hTRDin2MC"};
2064 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
2065 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
2069 const Int_t mdim(kNdim+3);
2070 Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
2071 Int_t ndim(H->GetNdimensions());
2072 TAxis *aa[mdim], *as(NULL), *ap(NULL), *apt(NULL), *ax(NULL), *abf(NULL); memset(aa, 0, sizeof(TAxis*) * (mdim));
2073 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2074 if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC);
2075 if(ndim > Int_t(kPt)) apt = H->GetAxis(kPt);
2076 if(ndim > Int_t(kNdim)) ap = H->GetAxis(kNdim);
2077 if(ndim > Int_t(kNdim)+1) ax = H->GetAxis(kNdim+1);
2078 if(ndim > Int_t(kNdim)+2) abf = H->GetAxis(kNdim+2);
2079 //AliInfo(Form("Using : Species[%c] Pt[%c] BunchFill[%c]", as?'y':'n', ap?'y':'n', abf?'y':'n'));
2080 const Int_t nPt(ndim>Int_t(kNdimTrkIn)?Int_t(kNpt+2):1);
2082 // build list of projections
2083 const Int_t nsel((kNpt+2)*(kNspc*kNcharge + 1));
2084 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2085 const Char_t *spcName[kNspc] = {"e", "#mu#pi", "Kp"};
2086 const Char_t ptName[kNpt+2] = {'L', 'l', 'i', 'h', 'H'};
2087 const Char_t *ptCut[kNpt+2] = {"p_{t}[GeV/c]<0.5", "0.5<=p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "1.5<=p_{t}[GeV/c]<5.0", "p_{t}[GeV/c]>=5.0"};
2088 const Char_t *pCut[kNpt+2] = {"p[GeV/c]<0.5", "0.5<=p[GeV/c]<0.8", "0.8<=p[GeV/c]<1.5", "1.5<=p[GeV/c]<5.0", "p[GeV/c]>=5.0"};
2089 // define rebinning strategy
2090 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
2091 if(aa[1]->GetNbins()==180){ // old binning
2092 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
2093 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
2094 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
2095 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
2096 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
2098 AliTRDrecoProjection hp[kMCTrkInNproj]; TObjArray php(kMCTrkInNproj+2);
2099 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2100 // define list of projections
2101 for(Int_t ipt(0); ipt<nPt; ipt++){
2102 for(Int_t isp(0); isp<kNspc; isp++){
2103 for(Int_t ich(0); ich<kNcharge; ich++){
2104 isel++; // new selection
2105 hp[ih].Build(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2106 Form("TrackIn[%s%c]:: #Deltay{%s}", spcName[isp], chSgn[ich], ptCut[ipt]),
2107 kEta, kPhi, kYrez, aa);
2108 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2109 php.AddLast(&hp[ih++]); np[isel]++;
2110 hp[ih].Build(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2111 Form("TrackIn[%s%c]:: #Delta#phi{%s}", spcName[isp], chSgn[ich], ptCut[ipt]),
2112 kEta, kPhi, kPrez, aa);
2113 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2114 php.AddLast(&hp[ih++]); np[isel]++;
2115 hp[ih].Build(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2116 Form("TrackIn[%s%c]:: dQdl {%s}", spcName[isp], chSgn[ich], pCut[ipt]),
2117 kEta, kPhi, kZrez, aa);
2118 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2119 php.AddLast(&hp[ih++]); np[isel]++;
2121 hp[ih].Build(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2122 Form("TrackIn[%s%c]:: #Deltax{%s}", spcName[isp], chSgn[ich], ptCut[ipt]),
2123 kEta, kPhi, kNdim+1, aa);
2124 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2125 php.AddLast(&hp[ih++]); np[isel]++;
2128 isel++; // RC tracks
2129 hp[ih].Build(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt]),
2130 Form("TrackIn[RC]:: #Deltaz{%s}", ptCut[ipt]),
2131 kEta, kPhi, kZrez, aa);
2132 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2133 php.AddLast(&hp[ih++]); np[isel]++;
2134 hp[ih].Build(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt]),
2135 Form("TrackIn[RC]:: #Deltay{%s}", ptCut[ipt]),
2136 kEta, kPhi, kYrez, aa);
2137 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2138 php.AddLast(&hp[ih++]); np[isel]++;
2139 hp[ih].Build(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt]),
2140 Form("TrackIn[RC]:: #Delta#phi{%s}", ptCut[ipt]),
2141 kEta, kPhi, kPrez, aa);
2142 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2143 php.AddLast(&hp[ih++]); np[isel]++;
2145 hp[ih].Build(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt]),
2146 Form("TrackIn[RC]:: #Deltax{%s}", ptCut[ipt]),
2147 kEta, kPhi, kNdim+1, aa);
2148 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2149 php.AddLast(&hp[ih++]); np[isel]++;
2151 AliInfo(Form("Build %3d 3D %s projections.", ih, mc?"MC":""));
2154 Int_t ch(0), pt(0), p(0), sp(1), rcBin(as?as->FindBin(0.):-1), ioff(0), joff(0), jsel(0);
2155 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
2156 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2157 v = H->GetBinContent(ib, coord);
2159 if(fBCbinTOF>0 && coord[kBC]!=fBCbinTOF) continue; // TOF bunch cross cut
2160 if(fBCbinFill>0 && abf && coord[kNdim+2]!=fBCbinTOF) continue; // Fill bunch cut
2161 if(coord[kBC]==3) continue;
2163 ch = 0; sp=1;// [pi-] track
2164 if(rcBin>0){ // debug mode in which species are also saved
2165 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2166 // take care of old data format (2*AliPID::kSPECIES+1)
2167 if(as->GetNbins() == kNcharge*AliPID::kSPECIES+1){
2169 else if(sp==2) sp=1;
2171 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
2172 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
2175 pt = 0; p = 0; // low pt
2176 if(apt) pt = TMath::Min(coord[kPt], Int_t(kNpt)+1);
2177 if(ap ) p = TMath::Min(coord[kNdim], Int_t(kNpt)+1);
2179 isel = pt*(kNspc*kNcharge+1); isel+=(ch==2?(kNspc*kNcharge):(sp*kNcharge+ch)); ioff = isel*(ax?4:3);
2180 jsel = p *(kNspc*kNcharge+1); jsel+=(ch==2?(kNspc*kNcharge):(sp*kNcharge+ch)); joff = jsel*(ax?4:3);
2182 AliError(Form("Wrong selection %d [%3d]", ioff, ih));
2185 if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
2186 AliError(Form("Missing projection %d", ioff));
2190 if(strcmp(pr0->H()->GetName(), Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ch], ptName[pt], sp))!=0){
2191 AliError(Form("Projection mismatch :: request[H%sTrkInY%c%c%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, pr0->H()->GetName()));
2195 if(strcmp(pr0->H()->GetName(), Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[pt]))!=0){
2196 AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c] found[%s]", mc?"MC":"", ptName[pt], pr0->H()->GetName()));
2200 AliDebug(2, Form("Found %s for selection sp[%d] ch[%d] pt[%d]", pr0->H()->GetName(), sp, ch, pt));
2201 for(Int_t jh(0); jh<np[isel]; jh++){
2202 if(ch<2 && jh==2) ((AliTRDrecoProjection*)php.At(joff+jh))->Increment(coord, v); // special care for dQdl=f(p)
2203 else ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v); // all = f(p_t)
2206 if(HasDump3DFor(kTrackIn)){
2207 TDirectory *cwd = gDirectory;
2208 TFile::Open(Form("ResDump_%s.root", H->GetName()), "RECREATE");
2209 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
2210 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
2211 if(!pr0->H()) continue;
2212 TH3 *h3=(TH3*)pr0->H()->Clone();
2219 TObjArray *arr(NULL);
2220 fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx);
2221 arr->SetName(projName[Int_t(mc)]);
2223 TH2 *h2(NULL); Int_t jh(0);
2225 if(!hp[ih].H()) continue;
2226 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2227 arr->AddAt(h2, jh++);
2229 // build combined performance plots
2230 // combine up the tree of projections
2231 Double_t m(0.), s(0.), trend(0.);
2232 AliTRDrecoProjection xlow[2], specY[kNcharge*kNspc], specPh[kNcharge*kNspc], specQ[kNcharge*kNspc];
2233 for(Int_t ich(0); ich<kNcharge; ich++){
2234 // PID dependency - summation over pt
2235 for(Int_t isp(0); isp<kNspc; isp++){
2237 Int_t idx(ich*kNspc+isp);
2238 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2239 specY[idx] = (*pr0);
2240 specY[idx].SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2241 specY[idx].H()->SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp),
2242 Form("TrackIn[%s%c]:: #Deltay", spcName[isp], chSgn[ich]));
2243 if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2244 PutTrendValue(Form("%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), trend);
2245 PutTrendValue(Form("%sTrkInYS%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), s);
2247 for(Int_t ipt(1); ipt<nPt; ipt++){
2248 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2249 if((trend=pr1->GetTrendValue(1,&m,&s))>-100.){
2250 PutTrendValue(Form("%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), trend);
2251 PutTrendValue(Form("%sTrkInYS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), s);
2255 php.AddLast(&specY[idx]);
2256 if((h2 = specY[idx].Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2257 if((trend=specY[idx].GetTrendValue(1,&m,&s))>-100.){
2258 PutTrendValue(Form("%sTrkInY%c%d", mc?"MC":"", chName[ich], isp), trend);
2259 PutTrendValue(Form("%sTrkInYS%c%d", mc?"MC":"", chName[ich], isp), s);
2261 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", specY[idx].H()->GetName())))) arr->AddAt(h2, jh++);
2262 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specY[idx];
2265 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2266 specPh[idx] = (*pr0);
2267 specPh[idx].SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2268 specPh[idx].H()->SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp),
2269 Form("TrackIn[%s%c]:: #Delta#phi", spcName[isp], chSgn[ich]));
2270 specPh[idx].SetShowRange(-1.5, 1.5);
2271 if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2272 PutTrendValue(Form("%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), trend);
2273 PutTrendValue(Form("%sTrkInPhS%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), s);
2275 for(Int_t ipt(1); ipt<nPt; ipt++){
2276 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2277 if((trend=pr1->GetTrendValue(1,&m,&s))>-100.){
2278 PutTrendValue(Form("%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), trend);
2279 PutTrendValue(Form("%sTrkInPhS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), s);
2281 specPh[idx]+=(*pr1);
2283 php.AddLast(&specPh[idx]);
2284 if((h2 = specPh[idx].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2285 if((trend=specPh[idx].GetTrendValue(1,&m,&s))>-100.){
2286 PutTrendValue(Form("%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp), trend);
2287 PutTrendValue(Form("%sTrkInPhS%c%d", mc?"MC":"", chName[ich], isp), s);
2289 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specPh[idx];
2292 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2293 specQ[idx] = (*pr0);
2294 specQ[idx].SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp), "Sum over p");
2295 specQ[idx].H()->SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp),
2296 Form("TrackIn[%s%c]:: dQdl", spcName[isp], chSgn[ich]));
2297 specQ[idx].SetShowRange(-2.2, -1.75);
2298 specQ[idx].H()->GetZaxis()->SetTitle("dQdl [a.u.]");
2299 if((trend = pr0->GetTrendValue(2, &m))>-100.){
2300 PutTrendValue(Form("%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), trend);
2301 PutTrendValue(Form("%sTrkInQS%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp), m);
2303 for(Int_t ipt(1); ipt<nPt; ipt++){
2304 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2305 if((trend=pr1->GetTrendValue(2, &m))>-100.){
2306 PutTrendValue(Form("%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), trend);
2307 PutTrendValue(Form("%sTrkInQS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp), m);
2311 php.AddLast(&specQ[idx]);
2312 if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2313 specQ[idx].H()->SetName(Form("H%sTrkInQS%c%d", mc?"MC":"", chName[ich], isp));
2314 specQ[idx].SetShowRange(-1.85, -1.4);
2315 if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2316 if((trend=specQ[idx].GetTrendValue(2, &m))>-100.){
2317 PutTrendValue(Form("%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp), trend);
2318 PutTrendValue(Form("%sTrkInQS%c%d", mc?"MC":"", chName[ich], isp), m);
2320 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specQ[idx];
2322 } // end PID loop for pt integration
2324 // pt dependency - summation over PID
2325 for(Int_t ipt(0); ipt<nPt; ipt++){
2327 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2328 for(Int_t isp(1); isp<kNspc; isp++){
2329 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2332 pr0->H()->SetNameTitle(Form("H%sTrkInY%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2333 Form("TrackIn[%c]:: #Deltay{%s}", chSgn[ich], ptCut[ipt]));
2334 pr0->SetShowRange(-0.3, 0.3);
2335 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2336 if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2337 PutTrendValue(Form("%sTrkInY%c%c", mc?"MC":"", chName[ich], ptName[ipt]), trend);
2338 PutTrendValue(Form("%sTrkInYS%c%c", mc?"MC":"", chName[ich], ptName[ipt]), s);
2340 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2343 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2344 for(Int_t isp(1); isp<kNspc; isp++){
2345 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2348 pr0->H()->SetNameTitle(Form("H%sTrkInPh%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2349 Form("TrackIn[%c]:: #Delta#phi{%s}", chSgn[ich], ptCut[ipt]));
2350 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2351 if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2352 PutTrendValue(Form("%sTrkInPh%c%c", mc?"MC":"", chName[ich], ptName[ipt]), trend);
2353 PutTrendValue(Form("%sTrkInPhS%c%c", mc?"MC":"", chName[ich], ptName[ipt]), s);
2355 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2358 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2359 for(Int_t isp(1); isp<kNspc; isp++){
2360 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2363 pr0->H()->SetNameTitle(Form("H%sTrkInX%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2364 Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[ipt]));
2365 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2366 PutTrendValue(Form("%sTrkInX%c%c", mc?"MC":"", chName[ich], ptName[ipt]), pr0->GetTrendValue(1));
2369 xlow[ich].SetNameTitle(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5),
2370 Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[0]));
2371 php.AddLast(&xlow[ich]);
2373 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2375 } // end pt loop for PID integration
2378 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2379 pr0->H()->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", chName[ich]),
2380 Form("TrackIn[%c]:: #Deltay", chSgn[ich]));
2381 pr0->SetShowRange(-0.3, 0.3);
2382 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2383 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2384 if((trend=pr0->GetTrendValue(1, &m, &s))>-100.){
2385 PutTrendValue(Form("%sTrkInY%c", mc?"MC":"", chName[ich]), trend);
2386 PutTrendValue(Form("%sTrkInYS%c", mc?"MC":"", chName[ich]), s);
2388 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2391 if(ich && (pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[3], 0)))){
2392 if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[3], 0)))){
2394 pr0->H()->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", ptName[3]), Form("TrackIn :: #Deltay{%s}", ptCut[3]));
2395 pr0->SetShowRange(-0.3, 0.3);
2396 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2400 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2401 pr0->H()->SetNameTitle(Form("H%sTrkInPh%c", mc?"MC":"", chName[ich]),
2402 Form("TrackIn[%c]:: #Delta#phi", chSgn[ich]));
2403 pr0->SetShowRange(-1., 1.);
2404 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2405 if((trend=pr0->GetTrendValue(1, &m, &s))>-100.){
2406 PutTrendValue(Form("%sTrkInPh%c", mc?"MC":"", chName[ich]), trend);
2407 PutTrendValue(Form("%sTrkInPhS%c", mc?"MC":"", chName[ich]), s);
2409 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2412 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2413 pr0->H()->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", chName[ich]),
2414 Form("TrackIn[%c]:: #Deltax", chSgn[ich]));
2415 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2416 PutTrendValue(Form("%sTrkInX%c", mc?"MC":"", chName[ich]), pr0->GetTrendValue(1));
2417 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2420 if(ich && (pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[1], 5)))){
2421 if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[1], 5)))){
2423 pr0->H()->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", ptName[1]), Form("TrackIn :: #Deltax{%s}", ptCut[1]));
2424 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2427 } // end charge loop
2429 for(Int_t isp(0); isp<kNspc; isp++){
2431 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))){
2432 pr0->H()->SetNameTitle(Form("H%sTrkInY%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Deltay", spcName[isp]));
2433 pr0->SetShowRange(-0.3, 0.3);
2434 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2435 PutTrendValue(Form("%sTrkInY%d", mc?"MC":"", isp), pr0->GetTrendValue(1));
2436 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2439 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))){
2440 pr0->H()->SetNameTitle(Form("H%sTrkInPh%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Delta#phi", spcName[isp]));
2441 pr0->SetShowRange(-1., 1.);
2442 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2443 PutTrendValue(Form("%sTrkInPh%d", mc?"MC":"", isp), pr0->GetTrendValue(1));
2446 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))){
2447 pr0->H()->SetNameTitle(Form("H%sTrkInQ%d", mc?"MC":"", isp), Form("TrackIn[%s] :: dQdl", spcName[isp]));
2448 pr0->SetShowRange(-2.2, -1.75);
2449 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2450 pr0->H()->SetName(Form("H%sTrkInQS%d", mc?"MC":"", isp));
2451 pr0->SetShowRange(-1.85, -1.4);
2452 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2453 if((trend=pr0->GetTrendValue(2, &m))>-100.){
2454 PutTrendValue(Form("%sTrkInQ%d", mc?"MC":"", isp), trend);
2455 PutTrendValue(Form("%sTrkInQS%d", mc?"MC":"", isp), m);
2458 } // end PID processing
2461 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2462 pr0->H()->SetNameTitle(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: #Deltay");
2463 pr0->SetShowRange(-0.3, 0.3);
2464 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2465 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2468 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2469 pr0->H()->SetNameTitle(Form("H%sTrkInPh", mc?"MC":""), "TrackIn :: #Delta#phi");
2470 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2473 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2474 pr0->H()->SetNameTitle(Form("H%sTrkInX", mc?"MC":""), "TrackIn :: #Deltax");
2475 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2478 // Row Cross processing
2480 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[0])))){
2481 for(Int_t ipt(0); ipt<kNpt; ipt++){
2482 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt])))) continue;
2485 pr0->H()->SetNameTitle(Form("H%sTrkInRCZ", mc?"MC":""), "TrackIn[RC]:: #Deltaz");
2486 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2487 if((trend=pr0->GetTrendValue(1,&m,&s))>-100.){
2488 PutTrendValue(Form("%sTrkInRCZ", mc?"MC":""), trend);
2489 PutTrendValue(Form("%sTrkInRCZS", mc?"MC":""), s);
2493 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[0])))){
2494 for(Int_t ipt(0); ipt<kNpt; ipt++){
2495 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt])))) continue;
2498 pr0->H()->SetNameTitle(Form("H%sTrkInRCY", mc?"MC":""), "TrackIn[RC]:: #Deltay");
2499 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2502 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[0])))){
2503 for(Int_t ipt(0); ipt<kNpt; ipt++){
2504 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt])))) continue;
2507 pr0->H()->SetNameTitle(Form("H%sTrkInRCPh", mc?"MC":""), "TrackIn[RC]:: #Delta#phi");
2508 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2511 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[0])))){
2512 for(Int_t ipt(0); ipt<kNpt; ipt++){
2513 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt])))) continue;
2516 pr0->H()->SetNameTitle(Form("H%sTrkInRCX", mc?"MC":""), "TrackIn[RC]:: #Deltax");
2517 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2519 AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
2524 //________________________________________________________
2525 Bool_t AliTRDresolution::MakeProjectionTrack()
2528 const Int_t kNcontours(9);
2529 const Int_t kNstat(30);
2530 Int_t cidx(kMCtrack);
2531 if(fProj && fProj->At(cidx)) return kTRUE;
2533 AliError("Missing data container.");
2537 if(!(H = (THnSparse*)fContainer->FindObject("hTRD2MC"))){
2538 AliError("Missing/Wrong data @ hTRD2MC.");
2541 Int_t ndim(H->GetNdimensions());
2542 Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
2543 TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
2544 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2545 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
2546 if(ndim > kPt) ap = H->GetAxis(kPt);
2548 // build list of projections
2549 const Int_t nsel(AliTRDgeometry::kNlayer*kNpt*AliPID::kSPECIES*7);//, npsel(3);
2550 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2551 const Char_t ptName[kNpt] = {'l', 'i', 'h'};
2552 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"};
2553 // define rebinning strategy
2554 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
2555 if(aa[1]->GetNbins()==180){ // old binning
2556 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
2557 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
2558 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
2559 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
2560 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
2562 AliTRDrecoProjection hp[kTrkNproj]; TObjArray php(kTrkNproj);
2563 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2564 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2565 for(Int_t ipt(0); ipt<kNpt; ipt++){
2566 for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2567 for(Int_t ich(0); ich<kNcharge; ich++){
2568 isel++; // new selection
2569 hp[ih].Build(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2570 Form("Tracks[%s%c]:: #Deltay{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2571 kEta, kPhi, kYrez, aa);
2572 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2573 php.AddLast(&hp[ih++]); np[isel]++;
2574 hp[ih].Build(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2575 Form("Tracks[%s%c]:: #Delta#phi{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2576 kEta, kPhi, kPrez, aa);
2577 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2578 php.AddLast(&hp[ih++]); np[isel]++;
2579 hp[ih].Build(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2580 Form("Tracks[%s%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2581 kEta, kPhi, kNdim, aa);
2582 hp[ih].SetShowRange(0.,10.);
2583 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2584 php.AddLast(&hp[ih++]); np[isel]++;
2587 isel++; // new selection
2588 hp[ih].Build(Form("HMCTrkZ%c%d", ptName[ipt], ily),
2589 Form("Tracks[RC]:: #Deltaz{%s} Ly[%d]", ptCut[ipt], ily),
2590 kEta, kPhi, kZrez, aa);
2591 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2592 php.AddLast(&hp[ih++]); np[isel]++;
2595 AliInfo(Form("Build %3d 3D projections.", ih));
2597 Int_t ly(0), ch(0), pt(0), sp(2), rcBin(as?as->FindBin(0.):-1);
2598 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2599 v = H->GetBinContent(ib, coord);
2601 ly = coord[kBC]-1; // layer selection
2603 ch=0; sp=2;// [pi-] track [dafault]
2604 if(rcBin>0){ // debug mode in which species are also saved
2605 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2606 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
2607 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
2610 pt = 0; // low pt [default]
2611 if(ap) pt = coord[kPt]-1;
2613 Int_t ioff = ly*kNpt*31+pt*31; ioff+=3*(sp<0?10:(sp*kNcharge+ch));
2614 isel = ly*kNpt*11+pt*11; isel+=sp<0?10:(sp*kNcharge+ch);
2615 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));
2616 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
2618 TObjArray *arr(NULL);
2619 fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
2620 arr->SetName("hTRD2MC");
2622 TH2 *h2(NULL); Int_t jh(0);
2624 if(!hp[ih].H()) continue;
2625 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2626 arr->AddAt(h2, jh++);
2629 // combine up the tree of projections
2630 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
2631 //Int_t iproj(0), jproj(0);
2632 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2633 for(Int_t ich(0); ich<kNcharge; ich++){
2634 for(Int_t ipt(0); ipt<kNpt; ipt++){
2636 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2637 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2638 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2641 AliDebug(2, Form("Rename %s to HMCTrkY%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], ily));
2642 pr0->H()->SetNameTitle(Form("HMCTrkY%c%c%d", chName[ich], ptName[ipt], ily),
2643 Form("Tracks[%c]:: #Deltay{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2644 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2645 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2648 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2649 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2650 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2653 AliDebug(2, Form("Rename %s to HMCTrkPh%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], ily));
2654 pr0->H()->SetNameTitle(Form("HMCTrkPh%c%c%d", chName[ich], ptName[ipt], ily),
2655 Form("Tracks[%c]:: #Delta#phi{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2656 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2657 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2661 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2662 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2663 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2666 AliDebug(2, Form("Rename %s to HMCTrkDPt%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], ily));
2667 pr0->H()->SetNameTitle(Form("HMCTrkDPt%c%c%d", chName[ich], ptName[ipt], ily),
2668 Form("Tracks[%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2669 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2670 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2674 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2675 pr0->H()->SetNameTitle(Form("HMCTrkY%c%d", chName[ich], ily),
2676 Form("Tracks[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
2677 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2678 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2681 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2682 pr0->H()->SetNameTitle(Form("HMCTrkPh%c%d", chName[ich], ily),
2683 Form("Tracks[%c]:: #Delta#phi Ly[%d]", chSgn[ich], ily));
2684 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2685 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2688 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2689 pr0->H()->SetNameTitle(Form("HMCTrkDPt%c%d", chName[ich], ily),
2690 Form("Tracks[%c]:: #Deltap_{t}/p_{t} Ly[%d]", chSgn[ich], ily));
2691 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2692 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2696 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2697 pr0->H()->SetNameTitle(Form("HMCTrkY%d", ily), Form("Tracks :: #Deltay Ly[%d]", ily));
2698 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2701 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2702 pr0->H()->SetNameTitle(Form("HMCTrkPh%d", ily), Form("Tracks :: #Delta#phi Ly[%d]", ily));
2703 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2706 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2707 pr0->H()->SetNameTitle(Form("HMCTrkDPt%d", ily), Form("Tracks :: #Deltap_{t}/p_{t} Ly[%d]", ily));
2708 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2711 AliInfo(Form("Done %3d 2D projections.", jh));
2715 //________________________________________________________
2716 Bool_t AliTRDresolution::PostProcess()
2718 // Fit, Project, Combine, Extract values from the containers filled during execution
2721 AliError("ERROR: list not available");
2725 AliInfo("Building array of projections ...");
2726 fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
2728 // set pt/p segmentation. guess from data
2730 if(!(H = (THnSparse*)fContainer->FindObject("hTracklet2TRDin"))){
2731 AliError("Missing/Wrong data @ hTracklet2TRDin.");
2734 fNpt=H->GetAxis(kPt)->GetNbins()+1;
2735 if(!MakeMomSegmentation()) return kFALSE;
2737 //PROCESS EXPERIMENTAL DISTRIBUTIONS
2738 // Clusters detector
2739 if(HasProcess(kDetector)) if(!MakeProjectionDetector()) return kFALSE;
2740 // Clusters residuals
2741 if(HasProcess(kCluster)) if(!MakeProjectionCluster()) return kFALSE;
2743 // Tracklet residual/pulls
2744 if(HasProcess(kTracklet)) if(!MakeProjectionTracklet()) return kFALSE;
2746 // TRDin residual/pulls
2747 if(HasProcess(kTrackIn)) if(!MakeProjectionTrackIn()) return kFALSE;
2750 if(!HasMCdata()) return kTRUE;
2751 //PROCESS MC RESIDUAL DISTRIBUTIONS
2753 // CLUSTER Y RESOLUTION/PULLS
2754 if(HasProcess(kCluster)) if(!MakeProjectionCluster(kTRUE)) return kFALSE;
2757 // TRACKLET RESOLUTION/PULLS
2758 if(HasProcess(kTracklet)) if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
2761 // TRACK RESOLUTION/PULLS
2762 if(HasProcess(kTracklet)) if(!MakeProjectionTrack()) return kFALSE;
2765 // TRACK TRDin RESOLUTION/PULLS
2766 if(HasProcess(kTrackIn)) if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
2773 //________________________________________________________
2774 void AliTRDresolution::Terminate(Option_t *opt)
2776 AliTRDrecoTask::Terminate(opt);
2777 if(HasPostProcess()) PostProcess();
2780 //________________________________________________________
2781 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2783 // Helper function to avoid duplication of code
2784 // Make first guesses on the fit parameters
2786 // find the intial parameters of the fit !! (thanks George)
2787 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2789 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2790 f->SetParLimits(0, 0., 3.*sum);
2791 f->SetParameter(0, .9*sum);
2792 Double_t rms = h->GetRMS();
2793 f->SetParLimits(1, -rms, rms);
2794 f->SetParameter(1, h->GetMean());
2796 f->SetParLimits(2, 0., 2.*rms);
2797 f->SetParameter(2, rms);
2798 if(f->GetNpar() <= 4) return;
2800 f->SetParLimits(3, 0., sum);
2801 f->SetParameter(3, .1*sum);
2803 f->SetParLimits(4, -.3, .3);
2804 f->SetParameter(4, 0.);
2806 f->SetParLimits(5, 0., 1.e2);
2807 f->SetParameter(5, 2.e-1);
2810 //________________________________________________________
2811 TObjArray* AliTRDresolution::Histos()
2814 // Define histograms
2817 if(fContainer) return fContainer;
2819 fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
2821 const Int_t nhn(100); Char_t hn[nhn]; TString st;
2822 Int_t triggerBins(fTriggerList?(TMath::Power(2.,fTriggerList->GetEntriesFast())-1):0);
2824 //++++++++++++++++++++++
2825 // cluster to detector
2826 snprintf(hn, nhn, "h%s", fgPerformanceName[kDetector]);
2827 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2828 const Int_t mdim(6);
2829 const Char_t *cldTitle[mdim] = {"layer", fgkTitle[kPhi], "pad row", "centrality", "q [a.u.]", triggerBins?"trigger":"no. of pads"/*, "tb [10 Hz]"*/};
2830 const Int_t cldNbins[mdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], 76, AliTRDeventInfo::kCentralityClasses, 50, triggerBins?triggerBins:kNpads};
2831 const Double_t cldMin[mdim] = {-0.5, fgkMin[kPhi], -0.5, -0.5, 0., 0.5},
2832 cldMax[mdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], 75.5, AliTRDeventInfo::kCentralityClasses - 0.5, 1200., triggerBins?(triggerBins+.5):(kNpads+.5)};
2833 st = "cluster proprieties;";
2834 // define minimum info to be saved in non debug mode
2835 Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimDet);
2836 for(Int_t idim(0); idim<ndim; idim++){ st += cldTitle[idim]; st+=";";}
2837 H = new THnSparseI(hn, st.Data(), ndim, cldNbins, cldMin, cldMax);
2839 fContainer->AddAt(H, kDetector);
2840 //++++++++++++++++++++++
2841 // cluster to tracklet residuals/pulls
2842 snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
2843 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2844 const Int_t mdim(10);
2845 const Char_t *clTitle[mdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]", "centrality", triggerBins?"trigger":"no. of pads"};
2846 const Int_t clNbins[mdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 30, 30, 15, AliTRDeventInfo::kCentralityClasses, triggerBins?triggerBins:kNpads};
2847 const Double_t clMin[mdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45, -0.5, 0.5},
2848 clMax[mdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45, AliTRDeventInfo::kCentralityClasses - 0.5, triggerBins?(triggerBins+.5):(kNpads+.5)};
2849 st = "cluster spatial&charge resolution;";
2850 // define minimum info to be saved in non debug mode
2851 Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimCl);
2852 for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2853 H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2855 fContainer->AddAt(H, kCluster);
2856 //++++++++++++++++++++++
2857 // tracklet to TRD track
2858 snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
2859 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2860 const Int_t mdim(kNdim+8);
2861 Char_t *trTitle[mdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2862 Int_t trNbins[mdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2863 Double_t trMin[mdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2864 Double_t trMax[mdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2865 // set specific fields
2866 trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
2867 trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
2868 trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
2869 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2871 trTitle[jdim]=StrDup("centrality"); trNbins[jdim] = AliTRDeventInfo::kCentralityClasses; trMin[jdim] = -.5; trMax[jdim] = AliTRDeventInfo::kCentralityClasses - 0.5;
2874 trTitle[jdim]=StrDup("trigger"); trNbins[jdim] = triggerBins;
2875 trMin[jdim] = 0.5; trMax[jdim] = triggerBins+.5;
2877 trTitle[jdim]=StrDup("occupancy [%]"); trNbins[jdim] = 12;
2878 trMin[jdim] = 25.; trMax[jdim] = 85.;
2881 st = "tracklet spatial&charge resolution;";
2882 // define minimum info to be saved in non debug mode
2883 Int_t ndim=DebugLevel()>=1?(jdim+1):kNdimTrklt;
2884 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2885 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2887 fContainer->AddAt(H, kTracklet);
2888 //++++++++++++++++++++++
2889 // tracklet to TRDin
2890 snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
2891 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2892 // set specific fields
2893 const Int_t mdim(kNdim+2);
2894 Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*));
2895 Int_t trinNbins[mdim]; memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t));
2896 Double_t trinMin[mdim]; memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t));
2897 Double_t trinMax[mdim]; memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t));
2899 trinTitle[kBC]=StrDup("trigger"); trinNbins[kBC] = triggerBins;
2900 trinMin[kBC] = 0.5; trinMax[kBC] = triggerBins+.5;
2902 // trinNbins[kSpeciesChgRC] = Int_t(kNcharge)*(kNspc-1)+1; trinMin[kSpeciesChgRC] = -kNspc+0.5; trinMax[kSpeciesChgRC] = kNspc-0.5;
2903 trinNbins[kPt]=fNpt-1; trinMax[kPt] = trinNbins[kPt]-.5;
2904 trinTitle[kNdim]=StrDup("bin_p"); trinNbins[kNdim] = fNpt-1; trinMin[kNdim] = -0.5; trinMax[kNdim] = trinNbins[kNdim]-.5;
2905 trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4;
2906 //trinTitle[kNdim+2]=StrDup("Fill Bunch"); trinNbins[kNdim+2]=3500; trinMin[kNdim+2]=-0.5; trinMax[kNdim+2]=3499.5;
2907 st = "r-#phi/z/angular residuals @ TRD entry;";
2908 // define minimum info to be saved in non debug mode
2909 Int_t ndim=(DebugLevel()>=1?mdim:(kNdim+1));//kNdimTrkIn;
2910 for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";}
2911 H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax);
2913 fContainer->AddAt(H, kTrackIn);
2914 // tracklet to TRDout
2915 // fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2918 // Resolution histos
2919 if(!HasMCdata()) return fContainer;
2921 //++++++++++++++++++++++
2922 // cluster to TrackRef residuals/pulls
2923 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
2924 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2925 H = (THnSparseI*)((THnSparseI*)fContainer->At(kCluster))->Clone(hn);
2926 H->SetTitle("MC cluster spatial resolution");
2928 fContainer->AddAt(H, kMCcluster);
2929 //++++++++++++++++++++++
2930 // tracklet to TrackRef
2931 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
2932 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2933 H = (THnSparseI*)((THnSparseI*)fContainer->At(kTracklet))->Clone(hn);
2934 H->SetTitle("MC tracklet spatial resolution");
2936 fContainer->AddAt(H, kMCtracklet);
2937 //++++++++++++++++++++++
2938 // TRDin to TrackRef
2939 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
2940 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2941 H = (THnSparseI*)((THnSparseI*)fContainer->At(kTrackIn))->Clone(hn);
2942 H->SetTitle("MC r-#phi/z/angular residuals @ TRD entry");
2944 fContainer->AddAt(H, kMCtrackIn);
2945 //++++++++++++++++++++++
2946 // track to TrackRef
2947 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
2948 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2949 Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2950 Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2951 Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2952 Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2953 // set specific fields
2954 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2955 trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
2956 trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
2957 trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2958 trNbins[kPt] =fNpt-1;trMax[kPt] = trNbins[kPt]-.5;
2959 //trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2960 trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5;
2962 st = "MC track spatial&p_{t} resolution;";
2963 // define minimum info to be saved in non debug mode
2964 Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
2965 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2966 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2968 fContainer->AddAt(H, kMCtrack);
2973 //____________________________________________________________________
2974 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2977 // Fit track with a staight line using the "np" clusters stored in the array "points".
2978 // The following particularities are stored in the clusters from points:
2979 // 1. pad tilt as cluster charge
2980 // 2. pad row cross or vertex constrain as fake cluster with cluster type 1
2981 // The parameters of the straight line fit are stored in the array "param" in the following order :
2982 // param[0] - x0 reference radial position
2983 // param[1] - y0 reference r-phi position @ x0
2984 // param[2] - z0 reference z position @ x0
2985 // param[3] - slope dy/dx
2986 // param[4] - slope dz/dx
2989 // Function should be used to refit tracks for B=0T
2993 if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2996 TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2999 for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
3002 Double_t x, y, z, dx, tilt(0.);
3003 for(Int_t ip(0); ip<np; ip++){
3004 x = points[ip].GetX(); z = points[ip].GetZ();
3006 zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
3008 if(zfitter.Eval() != 0) return kFALSE;
3010 Double_t z0 = zfitter.GetParameter(0);
3011 Double_t dzdx = zfitter.GetParameter(1);
3012 for(Int_t ip(0); ip<np; ip++){
3013 if(points[ip].GetClusterType()) continue;
3014 x = points[ip].GetX();
3016 y = points[ip].GetY();
3017 z = points[ip].GetZ();
3018 tilt = points[ip].GetCharge();
3019 y -= tilt*(-dzdx*dx + z - z0);
3020 Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
3021 yfitter.AddPoint(&dx, y, 1.);
3023 if(yfitter.Eval() != 0) return kFALSE;
3024 Double_t y0 = yfitter.GetParameter(0);
3025 Double_t dydx = yfitter.GetParameter(1);
3027 param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
3028 if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
3032 //____________________________________________________________________
3033 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
3036 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
3037 // See function FitTrack for the data stored in the "clusters" array
3039 // The parameters of the straight line fit are stored in the array "param" in the following order :
3040 // par[0] - x0 reference radial position
3041 // par[1] - y0 reference r-phi position @ x0
3042 // par[2] - slope dy/dx
3045 // Function should be used to refit tracks for B=0T
3048 TLinearFitter yfitter(2, "pol1");
3050 // grep data for tracklet
3051 Double_t x0(0.), x[60], y[60], dy[60];
3053 for(Int_t ip(0); ip<np; ip++){
3054 if(points[ip].GetClusterType()) continue;
3055 if(points[ip].GetVolumeID() != ly) continue;
3056 Float_t xt(points[ip].GetX())
3057 ,yt(param[1] + param[3] * (xt - param[0]));
3059 y[nly] = points[ip].GetY();
3065 if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
3068 // set radial reference for fit
3071 // find tracklet core
3072 Double_t mean(0.), sig(1.e3);
3073 AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
3075 // simple cluster error parameterization
3076 Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
3078 // fit tracklet core
3079 for(Int_t jly(0); jly<nly; jly++){
3080 if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
3081 Double_t dx(x[jly]-x0);
3082 yfitter.AddPoint(&dx, y[jly], 1.);
3084 if(yfitter.Eval() != 0) return kFALSE;
3086 par[1] = yfitter.GetParameter(0);
3087 par[2] = yfitter.GetParameter(1);
3091 //____________________________________________________________________
3092 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
3095 // Global selection mechanism of tracksbased on cluster to fit residuals
3096 // The parameters are the same as used ni function FitTrack().
3098 const Float_t kS(0.6), kM(0.2);
3099 TH1S h("h1", "", 100, -5.*kS, 5.*kS);
3100 Float_t dy, dz, s, m;
3101 for(Int_t ip(0); ip<np; ip++){
3102 if(points[ip].GetClusterType()) continue;
3103 Float_t x0(points[ip].GetX())
3104 ,y0(param[1] + param[3] * (x0 - param[0]))
3105 ,z0(param[2] + param[4] * (x0 - param[0]));
3106 dy=points[ip].GetY() - y0; h.Fill(dy);
3107 dz=points[ip].GetZ() - z0;
3109 TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
3110 fg.SetParameter(1, 0.);
3111 fg.SetParameter(2, 2.e-2);
3113 m=fg.GetParameter(1); s=fg.GetParameter(2);
3114 if(s>kS || TMath::Abs(m)>kM) return kFALSE;
3118 //________________________________________________________
3119 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
3122 // Get the most probable value and the full width half mean
3123 // of a Landau distribution
3126 const Float_t dx = 1.;
3127 mpv = f->GetParameter(1);
3128 Float_t fx, max = f->Eval(mpv);
3131 while((fx = f->Eval(xm))>.5*max){
3140 while((fx = f->Eval(xM))>.5*max) xM += dx;
3144 // #include "TFile.h"
3145 // //________________________________________________________
3146 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
3149 // AliWarning("Use cluster position as in reconstruction.");
3150 // SetLoadCorrection();
3153 // TDirectory *cwd(gDirectory);
3154 // TString fileList;
3155 // FILE *filePtr = fopen(file, "rt");
3157 // AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
3158 // SetLoadCorrection();
3161 // TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
3162 // while(fileList.Gets(filePtr)){
3163 // if(!TFile::Open(fileList.Data())) {
3164 // AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
3166 // } else AliInfo(Form("\"%s\"", fileList.Data()));
3168 // TTree *tSys = (TTree*)gFile->Get("tSys");
3169 // h2->SetDirectory(gDirectory); h2->Reset("ICE");
3170 // tSys->Draw("det:t>>h2", "dx", "goff");
3171 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3172 // for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
3174 // h2->SetDirectory(cwd);
3179 // if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>=2){
3180 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3181 // printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
3182 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3183 // FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
3184 // fprintf(fout, " static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
3185 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3186 // printf("%03d|", idet);
3187 // fprintf(fout, " {");
3188 // for(Int_t it(0); it<30; it++){
3189 // printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
3190 // fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
3193 // fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
3195 // fprintf(fout, " };\n");
3197 // SetLoadCorrection();
3201 //________________________________________________________
3202 void AliTRDresolution::SetProcesses(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3205 if(det) SETBIT(fSteer, kDetector); else CLRBIT(fSteer, kDetector);
3206 if(cl) SETBIT(fSteer, kCluster); else CLRBIT(fSteer, kCluster);
3207 if(trklt) SETBIT(fSteer, kTracklet); else CLRBIT(fSteer, kTracklet);
3208 if(trkin) SETBIT(fSteer, kTrackIn); else CLRBIT(fSteer, kTrackIn);
3212 //________________________________________________________
3213 void AliTRDresolution::SetDump3D(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3216 if(det) SETBIT(fSteer, 4+kDetector); else CLRBIT(fSteer, 4+kDetector);
3217 if(cl) SETBIT(fSteer, 4+kCluster); else CLRBIT(fSteer, 4+kCluster);
3218 if(trklt) SETBIT(fSteer, 4+kTracklet); else CLRBIT(fSteer, 4+kTracklet);
3219 if(trkin) SETBIT(fSteer, 4+kTrackIn); else CLRBIT(fSteer, 4+kTrackIn);