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");
35 // gSystem->Load("libTRDqaRec");
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 "AliESDHeader.h"
79 #include "AliMathBase.h"
80 #include "AliTrackPointArray.h"
82 #include "AliTRDresolution.h"
83 #include "AliTRDgeometry.h"
84 #include "AliTRDtransform.h"
85 #include "AliTRDpadPlane.h"
86 #include "AliTRDcluster.h"
87 #include "AliTRDseedV1.h"
88 #include "AliTRDtrackV1.h"
89 #include "AliTRDReconstructor.h"
90 #include "AliTRDrecoParam.h"
91 #include "AliTRDpidUtil.h"
92 #include "AliTRDinfoGen.h"
93 #include "AliAnalysisManager.h"
94 #include "info/AliTRDclusterInfo.h"
95 #include "info/AliTRDeventInfo.h"
97 ClassImp(AliTRDresolution)
98 ClassImp(AliTRDresolution::AliTRDrecoProjection)
100 Int_t const AliTRDresolution::fgkNbins[kNdim] = {
101 Int_t(kNbunchCross)/*bc*/,
109 }; //! no of bins/projection
110 Double_t const AliTRDresolution::fgkMin[kNdim] = {
117 -2.5, /*chg*species*/
119 }; //! low limits for projections
120 Double_t const AliTRDresolution::fgkMax[kNdim] = {
129 }; //! high limits for projections
130 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
139 }; //! title of projection
141 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
150 // ,"Tracklet2TRDout"
154 //________________________________________________________
155 AliTRDresolution::AliTRDresolution()
169 // Default constructor
171 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
172 SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
175 //________________________________________________________
176 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
177 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
190 // Default constructor
194 SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
196 SetUseExchangeContainers();
197 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
198 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
202 //________________________________________________________
203 AliTRDresolution::~AliTRDresolution()
208 if (AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
213 if(fCl){fCl->Delete(); delete fCl;}
214 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
218 //________________________________________________________
219 void AliTRDresolution::UserCreateOutputObjects()
221 // spatial resolution
223 AliTRDrecoTask::UserCreateOutputObjects();
224 if(UseExchangeContainers()) InitExchangeContainers();
227 //________________________________________________________
228 void AliTRDresolution::InitExchangeContainers()
230 // Init containers for subsequent tasks (AliTRDclusterResolution)
232 fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
233 fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
234 PostData(kClToTrk, fCl);
235 PostData(kClToMC, fMCcl);
238 //________________________________________________________
239 void AliTRDresolution::UserExec(Option_t *opt)
245 if(fCl) fCl->Delete();
246 if(fMCcl) fMCcl->Delete();
247 AliTRDrecoTask::UserExec(opt);
250 //________________________________________________________
251 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
253 // Helper function to calculate pulls in the yz plane
254 // using proper tilt rotation
255 // Uses functionality defined by AliTRDseedV1.
259 Double_t t2(tilt*tilt);
260 // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
264 cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
265 cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
266 cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
268 Double_t sqr[3]={0., 0., 0.};
269 if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
270 Double_t invsqr[3]={0., 0., 0.};
271 if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
272 Double_t tmp(dyz[0]);
273 dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
274 dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
279 //________________________________________________________
280 TH1* AliTRDresolution::DetCluster(const TObjArray *cls)
283 // Plot the cluster distributions
286 if(cls) fkClusters = cls;
288 AliDebug(4, "No Clusters defined.");
292 if(!(ncl = fkClusters->GetEntriesFast())){
293 AliDebug(1, "No RecPoints defined.");
297 AliTRDcluster *cl(NULL);
298 for(Int_t icl(0); icl<ncl; icl++){
299 if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])) continue;
300 det = cl->GetDetector(); break;
303 AliDebug(1, "No useful clusters defined.");
307 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kDetector)))){
308 AliWarning("No output container defined.");
311 Int_t ly(AliTRDgeometry::GetLayer(det)),
312 stk(AliTRDgeometry::GetStack(det));
314 alpha((0.5+AliTRDgeometry::GetSector(det))*AliTRDgeometry::GetAlpha()),
315 cs(TMath::Cos(alpha)),
316 sn(TMath::Sin(alpha));
317 for(Int_t icl(0); icl<ncl; icl++){
318 if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])){
319 AliDebug(2, Form("Missing clusters @ %d", icl));
322 if(cl->IsMasked()) printf("Mask %03d[%02d_%d_%d] %d[%d]", det, AliTRDgeometry::GetSector(det), stk, ly, cl->GetPadMaskedPosition(), cl->GetPadMaskedStatus());
324 val[kPhi] = TMath::ATan2(cl->GetX()*sn + cl->GetY()*cs, cl->GetX()*cs - cl->GetY()*sn);
325 val[kEta] = (5-stk)*16-cl->GetPadRow()-1-(stk<3?4:0);
326 val[kYrez]= fEvent->GetMultiplicity();
327 val[4] = TMath::Abs(cl->GetQ());
328 val[5] = fTriggerList?fTriggerSlot:(cl->IsFivePad()?1:cl->GetNPads());
334 //________________________________________________________
335 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
338 // Plot the cluster distributions
341 if(track) fkTrack = track;
343 AliDebug(4, "No Track defined.");
346 if(fkESD && TMath::Abs(fkESD->GetTOFbc()) > 1){
347 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
350 if(fkESD && fPt<fPtThreshold){
351 AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
355 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
356 AliWarning("No output container defined.");
360 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
361 Double_t val[kNdim+2],
362 alpha(0.), cs(-2.), sn(0.); //Float_t exb, vd, t0, s2, dl, dt;
363 TObjArray *clInfoArr(NULL);
364 AliTRDseedV1 *fTracklet(NULL);
365 AliTRDcluster *c(NULL), *cc(NULL);
366 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
367 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
368 if(!fTracklet->IsOK()) continue;
369 //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
372 alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
373 cs = TMath::Cos(alpha);
374 sn = TMath::Sin(alpha);
376 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
377 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
378 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
379 val[kPt] = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
380 Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
382 Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
383 fTracklet->ResetClusterIter(kTRUE);
384 while((c = fTracklet->NextCluster())){
385 Float_t xc(c->GetX()),
386 q(TMath::Abs(c->GetQ()));
387 if(row0<0) row0 = c->GetPadRow();
389 val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
390 val[kPrez] = fTracklet->GetX0()-xc;
391 val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
392 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
393 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
394 if(ic) val[kZrez] /= (ic*q);
395 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
396 val[kNdim] = fEvent?fEvent->GetMultiplicity():0.;
397 val[kNdim+1] = fTriggerList?fTriggerSlot:(c->IsFivePad()?1:c->GetNPads());
399 /* // tilt rotation of covariance for clusters
400 Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
401 cov[0] = (sy2+t2*sz2)*corr;
402 cov[1] = tilt*(sz2 - sy2)*corr;
403 cov[2] = (t2*sy2+sz2)*corr;
404 // sum with track covariance
405 cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
406 Double_t dyz[2]= {dy[1], dz[1]};
407 Pulls(dyz, cov, tilt);*/
409 // Get z-position with respect to anode wire
410 Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
411 zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
412 Int_t istk = geo->GetStack(c->GetDetector());
413 AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
414 Float_t rowZ = pp->GetRow0();
415 Float_t d = rowZ - zt + pp->GetAnodeWireOffset();
416 d -= ((Int_t)(2 * d)) / 2.0;
417 if (d > 0.25) d = 0.5 - d;
419 AliTRDclusterInfo *clInfo(NULL);
420 clInfo = new AliTRDclusterInfo;
421 clInfo->SetCluster(c);
422 //Float_t covF[] = {cov[0], cov[1], cov[2]};
423 clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
424 clInfo->SetResolution(val[kYrez]);
425 clInfo->SetAnisochronity(d);
426 clInfo->SetDriftLength(val[kZrez]);
427 clInfo->SetTilt(fTracklet->GetTilt());
428 if(fCl) fCl->Add(clInfo);
429 //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
433 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
434 clInfoArr->SetOwner(kFALSE);
436 clInfoArr->Add(clInfo);
439 if(DebugLevel()>=2 && clInfoArr){
440 ULong_t status = fkESD->GetStatus();
441 (*DebugStream()) << "cluster"
442 <<"status=" << status
443 <<"clInfo.=" << clInfoArr
448 if(clInfoArr) delete clInfoArr;
450 if(!track) return NULL;
451 // special care for EVE usage
453 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
454 return H->Projection(kYrez);
457 //________________________________________________________
458 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
460 // Plot normalized residuals for tracklets to track.
462 // We start from the result that if X=N(|m|, |Cov|)
464 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
466 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
467 // reference position.
468 if(track) fkTrack = track;
470 AliDebug(4, "No Track defined.");
473 if(fkESD && TMath::Abs(fkESD->GetTOFbc())>1){
474 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
478 if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
479 AliWarning("No output container defined.");
483 // down scale PID resolution
484 Int_t spc(fSpecies); if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
485 const Int_t ndim(kNdim+8);
487 alpha(0.), cs(-2.), sn(0.);
488 // Float_t sz[AliTRDseedV1::kNtb], pos[AliTRDseedV1::kNtb]; Int_t ntbGap[AliTRDseedV1::kNtb];
489 AliTRDseedV1 *fTracklet(NULL);
490 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
491 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
492 if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()) continue;
495 alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
496 cs = TMath::Cos(alpha);
497 sn = TMath::Sin(alpha);
499 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
500 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());//fTracklet->GetTgl();
501 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
503 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:spc;// fSpecies;
504 val[kPt] = GetPtBin(fPt); //fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fTracklet->GetMomentum());
505 Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)),
506 dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)),
507 dydx(fTracklet->GetYfit(1)),
508 tilt(fTracklet->GetTilt());
509 // correct for tilt rotation
510 val[kYrez] = dyt - dzt*tilt;
511 val[kZrez] = fTracklet->IsRowCross()?(dzt + dyt*tilt):(fTracklet->GetdQdl()*3.e-4-1.5);
512 dydx+= tilt*fTracklet->GetZref(1);
513 val[kPrez] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
514 val[kNdim] = fEvent?fEvent->GetMultiplicity():0;
515 val[kNdim+1] = fTriggerList?fTriggerSlot:(1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb);
516 /* Int_t n = fTracklet->GetChargeGaps(sz, pos, ntbGap);
517 val[kNdim+2] = 0.; for(Int_t igap(0); igap<n; igap++) val[kNdim+2] += sz[igap];
518 for(Int_t ifill(0); ifill<3; ifill++){ val[kNdim+3+ifill]=0.;val[kNdim+4+ifill]=0.;}
519 for(Int_t igap(0), ifill(0); igap<n&&ifill<2; igap++){
520 if(ntbGap[igap]<2) continue;
521 val[kNdim+3+ifill] = sz[igap];
522 val[kNdim+4+ifill] = pos[igap];
527 // // compute covariance matrix
528 // fTracklet->GetCovAt(x, cov);
529 // fTracklet->GetCovRef(covR);
530 // cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
531 // Double_t dyz[2]= {dy[1], dz[1]};
532 // Pulls(dyz, cov, tilt);
533 // ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
534 // ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
537 Bool_t rc(fTracklet->IsRowCross());
538 UChar_t err(fTracklet->GetErrorMsg());
539 Double_t x(fTracklet->GetX()),
540 pt(fTracklet->GetPt()),
541 yt(fTracklet->GetYref(0)),
542 zt(fTracklet->GetZref(0)),
543 phi(fTracklet->GetYref(1)),
544 tht(fTracklet->GetZref(1));
545 Int_t ncl(fTracklet->GetN()), det(fTracklet->GetDetector());
546 (*DebugStream()) << "tracklet"
557 <<"dy=" << val[kYrez]
558 <<"dz=" << val[kZrez]
559 <<"dphi="<< val[kPrez]
560 <<"dQ ="<< val[kNdim]
566 if(!track) return NULL;
567 // special care for EVE usage
569 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
570 return H->Projection(kYrez);
574 //________________________________________________________
575 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
577 // Store resolution/pulls of Kalman before updating with the TRD information
578 // at the radial position of the first tracklet. The following points are used
580 // - the (y,z,snp) of the first TRD tracklet
581 // - the (y, z, snp, tgl, pt) of the MC track reference
583 // Additionally the momentum resolution/pulls are calculated for usage in the
586 if(track) fkTrack = track;
588 AliDebug(4, "No Track defined.");
592 TH1 *h(NULL); // EVE projection
594 THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
596 AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
599 // check input track status
600 AliExternalTrackParam *tin(NULL);
601 if(!(tin = fkTrack->GetTrackIn())){
602 AliError("Track did not entered TRD fiducial volume.");
605 // check first tracklet
606 AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
608 AliDebug(3, "No Tracklet in ly[0]. Skip track.");
609 if(!track) return NULL;
610 // special care for EVE usage
611 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
612 return H->Projection(kYrez);
614 if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()){
615 AliDebug(3, "Tracklet or Chamber not OK. Skip track.");
616 if(!track) return NULL;
617 // special care for EVE usage
618 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
619 return H->Projection(kYrez);
621 // check radial position
622 Double_t x = tin->GetX();
623 if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
624 AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
625 if(!track) return NULL;
626 // special care for EVE usage
627 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
628 return H->Projection(kYrez);
631 // down scale PID resolution
632 Int_t spc(fSpecies); if(spc==3) spc=2; if(spc==-3) spc=-2; if(spc>3) spc=3; if(spc<-3) spc=-3;
633 // read V0 PID if available
635 if(fkESD->IsElectron()) v0pid = -1;
636 else if(fkESD->IsPion()) v0pid = 0;
637 else if(fkESD->IsProton()) v0pid = 1;
638 if(DebugLevel()>=1/* && v0pid>-2*/){
639 Float_t tpc(fkESD->GetTPCdedx());
640 Float_t tof(fkESD->GetTOFbeta());
641 Int_t ev(fEvent->GetEventHeader()->GetEventNumberInFile());
642 AliTRDtrackV1 t(*fkTrack); t.SetOwner();
643 (*DebugStream()) << "trackIn"
653 //Int_t bc(fkESD?fkESD->GetTOFbc()/2:0);
654 const Double_t *parR(tin->GetParameter());
655 Double_t dyt(fTracklet->GetYfit(0)-parR[0]), dzt(fTracklet->GetZfit(0)-parR[1]),
656 phit(fTracklet->GetYfit(1)),
657 tilt(fTracklet->GetTilt()),
658 norm(1./TMath::Sqrt((1.-parR[2])*(1.+parR[2])));
660 // correct for tilt rotation
661 Double_t dy = dyt - dzt*tilt,
663 dx = dy/(parR[2]*norm-parR[3]*norm*tilt);
664 phit += tilt*parR[3];
665 Double_t dphi = TMath::ATan(phit) - TMath::ASin(parR[2]);
667 Double_t val[kNdim+2];
668 val[kBC] = fTriggerSlot?fTriggerSlot:v0pid;//bc==0?0:(bc<0?-1.:1.);
669 Double_t alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha(),
670 cs = TMath::Cos(alpha),
671 sn = TMath::Sin(alpha);
672 val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
673 Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
674 val[kEta] = -TMath::Log(TMath::Tan(0.5 * (0.5*TMath::Pi() - TMath::ATan(tgl))));
676 val[kZrez] = fTracklet->IsRowCross()?dz:(fTracklet->GetdQdl()*5.e-4 - 2.5);
677 val[kPrez] = dphi*TMath::RadToDeg();
678 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:spc;
679 Float_t pz(0.); if(fTracklet->GetMomentum()-fPt>1.e-5) pz = TMath::Sqrt((fTracklet->GetMomentum()-fPt)*(fTracklet->GetMomentum()+fPt));
680 val[kPt] = fTracklet->IsRowCross()?GetPtBin(pz):GetPtBin(fPt);
681 val[kNdim] = GetPtBin(fTracklet->GetMomentum());
683 //val[kNdim+2] = fEvent?fEvent->GetBunchFill():0;
686 if(!track) return NULL;
687 // special care for EVE usage
688 if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
689 return H->Projection(kYrez);
693 //________________________________________________________
694 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
696 // Store resolution/pulls of Kalman after last update with the TRD information
697 // at the radial position of the first tracklet. The following points are used
699 // - the (y,z,snp) of the first TRD tracklet
700 // - the (y, z, snp, tgl, pt) of the MC track reference
702 // Additionally the momentum resolution/pulls are calculated for usage in the
705 if(track) fkTrack = track;
709 //________________________________________________________
710 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
713 // Plot MC distributions
716 if(!HasMCdata()) return NULL;
717 if(track) fkTrack = track;
719 AliDebug(4, "No Track defined.");
722 // Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
726 AliWarning("No output container defined.");
729 // retriev track characteristics
730 Int_t pdg = fkMC->GetPDG(),
731 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
734 label(fkMC->GetLabel());
736 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
737 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
738 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
739 if(sIdx==3) sIdx=2; if(sIdx>3) sIdx=3; // downscale PID
742 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
743 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
745 Double_t x, y, z, pt, dydx, dzdx/*, dzdl*/;
746 Float_t pt0, p0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
747 Double_t covR[7]/*, cov[3]*/;
749 /* if(DebugLevel()>=3){
750 // get first detector
752 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
753 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
754 det = fTracklet->GetDetector();
758 TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
760 m = fkTrack->GetMass();
761 if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
762 (*DebugStream()) << "MCkalman"
779 AliTRDcluster *c(NULL);
780 Double_t val[kNdim+2];
781 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
782 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
783 !fTracklet->IsOK())*/ continue;
785 x= x0 = fTracklet->GetX();
786 Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
787 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, p0, eta, phi, s)) continue;
789 // MC track position at reference radial position
791 Float_t ymc = y0 - dx*dydx0;
792 Float_t zmc = z0 - dx*dzdx0;
793 //phi -= TMath::Pi();
798 val[kSpeciesChgRC]= rc?0.:sign*sIdx;
799 val[kPt] = GetPtBin(pt0);
800 val[kNdim]= GetPtBin(p0);
801 Double_t tilt(fTracklet->GetTilt());
803 // ,corr(1./(1. + t2))
804 // ,cost(TMath::Sqrt(corr));
806 AliExternalTrackParam *tin(fkTrack->GetTrackIn());
807 if(ily==0 && tin){ // trackIn residuals
808 // check radial position
809 if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
812 switch(TMath::Abs(pdg)){
813 case 11: exactPID = -1;break;
814 case 211: exactPID = 0;break;
815 case 2212: exactPID = 1;break;
818 val[kYrez] = tin->GetY()-ymc;
819 val[kZrez] = rc?(tin->GetZ()-zmc):(fTracklet->GetdQdl()*5e-4 - 2.5);
820 val[kPrez] = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
821 val[kNdim+1] = 0.;//dx;
822 if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
823 val[kBC] = ily; // reset for subsequent components
825 if(DebugLevel()>=1 && exactPID>-2){
826 Float_t tpc(fkESD->GetTPCdedx());
827 Float_t tof(fkESD->GetTOFbeta());
828 AliTRDtrackV1 t(*fkTrack); t.SetOwner();
829 (*DebugStream()) << "MCtrackIn"
839 //if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
842 dydx = fTracklet->GetYref(1);
843 dzdx = fTracklet->GetZref(1);
844 //dzdl = fTracklet->GetTgl();
845 y = fTracklet->GetYref(0);
847 z = fTracklet->GetZref(0);
849 pt = TMath::Abs(fTracklet->GetPt());
850 fTracklet->GetCovRef(covR);
853 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
855 val[kNdim] = 1.e2*(pt/pt0-1.);
856 if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
857 /* // theta resolution/ tgl pulls
858 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
859 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
860 ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
861 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
862 // pt resolution \\ 1/pt pulls \\ p resolution for PID
863 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
864 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
865 ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
866 ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
867 ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
869 // Fill Debug stream for MC track
871 Int_t det(fTracklet->GetDetector());
872 (*DebugStream()) << "MC"
884 // Fill Debug stream for Kalman track
885 (*DebugStream()) << "MCtrack"
897 // tracklet residuals
898 dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
899 dzdx = fTracklet->GetZfit(1);
900 y = fTracklet->GetYfit(0);
902 z = fTracklet->GetZfit(0);
904 val[kYrez] = dy - dz*tilt;
905 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
906 val[kZrez] = rc?(dz + dy*tilt):(fTracklet->GetdQdl()*3.e-4 - 1.5);
907 val[kNdim] = fEvent?fEvent->GetMultiplicity():0;
908 val[kNdim+1] = 1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb;
909 // val[kNdim] = pt/pt0-1.;
910 if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
913 // Fill Debug stream for tracklet
915 Float_t s2y = fTracklet->GetS2Y();
916 Float_t s2z = fTracklet->GetS2Z();
917 (*DebugStream()) << "MCtracklet"
928 AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
929 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
930 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
932 H = (THnSparse*)fContainer->At(kMCcluster);
933 val[kPt] = TMath::ATan(dydx0)*TMath::RadToDeg();
934 //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
936 Float_t padCorr(tilt*fTracklet->GetPadLength()),
937 corr(1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0));
938 fTracklet->ResetClusterIter(kTRUE);
939 while((c = fTracklet->NextCluster())){
940 if(row0<0) row0 = c->GetPadRow();
941 x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
942 y = c->GetY() + padCorr*(c->GetPadRow() - row0);
949 val[kYrez] = dy - dz*tilt;
951 val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t q(TMath::Abs(c->GetQ()));
952 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
953 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
954 if(ic) val[kZrez] /= (ic*q);
955 val[kSpeciesChgRC]= rc?0.:(TMath::Max(q*corr, Float_t(3.)));
956 val[kNdim+1] = c->IsFivePad()?1:c->GetNPads();
960 // Fill calibration container
961 Float_t d = zr0 - zmc;
962 d -= ((Int_t)(2 * d)) / 2.0;
963 if (d > 0.25) d = 0.5 - d;
964 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
965 clInfo->SetCluster(c);
966 clInfo->SetMC(pdg, label);
967 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
968 clInfo->SetResolution(dy);
969 clInfo->SetAnisochronity(d);
970 clInfo->SetDriftLength(dx);
971 clInfo->SetTilt(tilt);
972 if(fMCcl) fMCcl->Add(clInfo);
973 else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
976 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
977 clInfoArr->SetOwner(kFALSE);
979 clInfoArr->Add(clInfo);
983 if(DebugLevel()>=5 && clInfoArr){
984 (*DebugStream()) << "MCcluster"
985 <<"clInfo.=" << clInfoArr
990 if(clInfoArr) delete clInfoArr;
991 if(!track) return NULL;
992 // special care for EVE usage
993 if(H && (h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
994 return H?H->Projection(kYrez):NULL;
997 //________________________________________________________
998 void AliTRDresolution::GetRangeZ(TH2 *h2, Float_t &min, Float_t &max)
1000 // Get range on Z axis such to avoid outliers
1002 Double_t cnt[20000], c, m, s;
1003 Int_t nx(h2->GetXaxis()->GetNbins()), ny(h2->GetYaxis()->GetNbins()), nc(0);
1004 for(Int_t ix(1); ix<=nx; ix++){
1005 for(Int_t iy(1); iy<=ny; iy++){
1006 if((c = h2->GetBinContent(ix, iy))<10) continue;
1008 if(nc==20000) break;
1010 if(nc==20000) break;
1012 AliMathBase::EvaluateUni(nc, cnt, m, s, 0);
1013 min = m-s; max = m+2.*s;
1016 //________________________________________________________
1017 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
1020 // Get the reference figures
1024 AliWarning("Please provide a canvas to draw results.");
1027 /* Int_t selection[100], n(0), selStart(0); //
1028 Int_t ly0(0), dly(5);
1029 TList *l(NULL); TVirtualPad *pad(NULL); */
1034 AliWarning(Form("Reference plot [%d] missing result", ifig));
1038 //________________________________________________________
1039 void AliTRDresolution::MakeSummary()
1041 // Build summary plots
1044 AliError("Missing results");
1047 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
1048 TObjArray *arr(NULL); TH2 *h2(NULL);
1049 TH2 *h2e[100] = {NULL}; Int_t ih2e(0); // save sigma histos for later deletion
1051 // cluster resolution
1053 gStyle->SetPalette(1);
1054 const Int_t nClViews(9);
1055 const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh"};
1056 const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1};
1057 const Float_t vClSignMin[2] = {2.6e2, 4.4e2},
1058 vClSignMax[2] = {4.4e2, 6.2e2},
1059 vClMin[nClViews] = {3.2e2, vClSignMin[Int_t(fBsign)], vClSignMin[Int_t(!fBsign)], 0., 0., 0., 0., 0., 3.2e2},
1060 vClMax[nClViews] = {5.e2, vClSignMax[Int_t(fBsign)], vClSignMax[Int_t(!fBsign)], 0., 0., 0., 0., 0., 5.e2};
1061 const Int_t nTrkltViews(20);
1062 const Char_t *vTrkltName[nTrkltViews] = {
1063 "TrkltY", "TrkltYn", "TrkltYp", "TrkltY", "TrkltYn", "TrkltYp",
1065 "TrkltPh", "TrkltPhn", "TrkltPhp",
1066 "TrkltQ", "TrkltQn", "TrkltQp",
1067 "TrkltQS", "TrkltQSn", "TrkltQSp",
1068 "TrkltTBn", "TrkltTBp", "TrkltTBn", "TrkltTBp"
1069 // "TrkltYnl0", "TrkltYpl0", "TrkltPhnl0", "TrkltPhpl0", "TrkltQnl0", "TrkltQpl0", // electrons low pt
1070 /* "TrkltYnl1", "TrkltYpl1", "TrkltPhnl1", "TrkltPhpl1", "TrkltQnl1", "TrkltQpl1", // muons low pt
1071 "TrkltYnl2", "TrkltYpl2", "TrkltPhnl2", "TrkltPhpl2", "TrkltQnl2", "TrkltQpl2" // pions low pt*/
1073 const UChar_t vTrkltOpt[nTrkltViews] = {
1081 const Int_t nTrkInViews(5);
1082 const Char_t *vTrkInName[nTrkInViews][6] = {
1083 {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1084 {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1085 {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"},
1086 {"TrkInXnl", "TrkInXpl", "TrkInXl", "TrkInYnh", "TrkInYph", "TrkInYh"},
1087 {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"},
1088 //{"TrkInRCX", "TrkInRCY", "TrkInRCPh", "TrkInRCZl", "TrkInRCZi", "TrkInRCZh"}
1090 const UChar_t vTrkInOpt[nTrkInViews] = {0, 1, 0, 0, 0};
1091 const Float_t min[6] = {0.15, 0.15, 0.15, 0.15, 0.5, 0.5};
1092 const Float_t max[6] = {0.6, 0.6, 0.6, 0.6, 2.3, 2.3};
1093 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]"};
1095 const Int_t nTrkViews(27);
1096 const Char_t *vTrkName[nTrkViews] = {
1097 "TrkY", "TrkYn", "TrkYp",
1098 "TrkPh", "TrkPhn", "TrkPhp",
1099 "TrkDPt", "TrkDPtn", "TrkDPtp",
1100 "TrkYnl0", "TrkYpl0", "TrkPhnl0", "TrkPhpl0", "TrkDPtnl0", "TrkDPtpl0", // electrons low pt
1101 "TrkYnl1", "TrkYpl1", "TrkPhnl1", "TrkPhpl1", "TrkDPtnl1", "TrkDPtpl1", // muons low pt
1102 "TrkYnl2", "TrkYpl2", "TrkPhnl2", "TrkPhpl2", "TrkDPtnl2", "TrkDPtpl2" // pions low pt
1104 const Char_t *typName[] = {"", "MC"};
1105 const Int_t nx(2048), ny(1536);
1107 if((arr = (TObjArray*)fProj->FindObject("hDet2Cluster"))){
1108 cOut = new TCanvas(Form("%s_DetOccupancy", GetName()), "Detector performance", 2*nx, 2*ny);
1109 cOut->Divide(AliTRDgeometry::kNlayer,AliTRDeventInfo::kCentralityClasses, 1.e-5, 1.e-5);
1111 for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1112 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1113 p=cOut->cd(icen*AliTRDgeometry::kNlayer+ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1114 if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d%dEn", ily, icen)))) continue;
1115 if(SetNormZ(h2, 1, -1, 1, -1, 10.) < 1.e3) continue; // cut all bins with < 10 entries
1116 SetRangeZ(h2, -90., 90, -200.);
1117 //h2->GetXaxis()->SetNdivisions(1010);
1118 h2->GetZaxis()->SetTitle("Rel. Det. Occup. [%]");
1119 h2->GetZaxis()->CenterTitle();
1120 h2->SetContour(9); h2->Draw("colz"); n++;
1121 MakeDetectorPlot(ily, "p");
1124 if(n>=AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1126 cOut = new TCanvas(Form("%s_DetCharge", GetName()), "Detector performance", nx, ny);
1127 cOut->Divide(3,2, 1.e-5, 1.e-5);
1128 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1129 p=cOut->cd(ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1130 if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d_2D", UseLYselectTrklt()?fLYselect:ily)))) continue;
1131 SetNormZ(h2, 1, -1, 1, -1, 10.); // cut all <q> < 10
1132 SetRangeZ(h2, -30., 30., -200.);
1133 //h2->GetXaxis()->SetNdivisions(1010);
1134 h2->GetZaxis()->SetTitle("Rel. Mean(q) [%]");
1135 h2->GetZaxis()->CenterTitle();
1137 MakeDetectorPlot(ily, "p");
1139 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1141 for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
1142 if((arr = (TObjArray*)fProj->FindObject(ityp?"hCluster2MC":"hCluster2Track"))){
1143 for(Int_t iview(0); iview<nClViews; iview++){
1144 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vClName[iview], vClOpt[iview]), "Cluster Resolution", nx, ny);
1145 cOut->Divide(3,2, 1.e-5, 1.e-5);
1147 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1148 p=cOut->cd(ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1149 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], UseLYselectTrklt()?fLYselect:ily)))) continue;
1151 if(vClOpt[iview]==0) h2->Draw("colz");
1152 else if(vClOpt[iview]==1) h2e[ih2e++] = DrawSigma(h2, "#sigma(#Deltay) [#mum]", vClMin[iview], vClMax[iview], 1.e4);
1153 if(iview<5) MakeDetectorPlot(ily);
1155 if(nplot==AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1159 // tracklet systematic
1160 if((arr = (TObjArray*)fProj->FindObject(ityp?"hTracklet2MC":"hTracklet2Track"))){
1161 for(Int_t iview(0); iview<nTrkltViews; iview++){
1162 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkltName[iview], vTrkltOpt[iview]), "Tracklet Resolution", nx, ny);
1163 cOut->Divide(3,2, 1.e-5, 1.e-5);
1165 for(Int_t iplot(0); iplot<6; iplot++){
1166 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1167 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
1169 if(vTrkltOpt[iview]==0) h2->Draw("colz");
1170 else if (vTrkltOpt[iview]==1) h2e[ih2e++] = DrawSigma(h2, "#sigma(#Deltay) [cm]", .15, .4);
1171 else if (vTrkltOpt[iview]==2) h2e[ih2e++] = DrawSigma(h2, "#sigma(occupancy) [%]", 10.5, 15.);
1172 MakeDetectorPlot(iplot);
1175 cOut->Modified();cOut->Update();
1176 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1181 // trackIn systematic
1182 for(Int_t iv0(0); iv0<2; iv0++){
1183 const Char_t *prefix = ityp?"MC":(iv0?"V0":"");
1184 if(ityp && !iv0) arr = (TObjArray*)fProj->FindObject("hTRDin2MC");
1185 else if(!ityp && !iv0) arr = (TObjArray*)fProj->FindObject("hTracklet2TRDin");
1186 else if(!ityp && iv0) arr = (TObjArray*)fProj->FindObject("hTracklet2TRDinV0");
1189 for(Int_t iview(0); iview<nTrkInViews; iview++){
1190 cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), prefix, vTrkInName[iview][0], vTrkInOpt[iview]), "Track IN Resolution", nx, ny);
1191 cOut->Divide(3,2, 1.e-5, 1.e-5);
1193 for(Int_t iplot(0); iplot<6; iplot++){
1194 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1195 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", prefix, vTrkInName[iview][iplot])))){
1196 AliInfo(Form("Missing H%s%s_2D", prefix, vTrkInName[iview][iplot]));
1200 if(vTrkInOpt[iview]==0) h2->Draw("colz");
1201 else h2e[ih2e++] = DrawSigma(h2, ttt[iplot], min[iplot], max[iplot]);
1202 MakeDetectorPlot(0);
1204 if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1208 const Float_t zmin[] = {1., 61., 15.},
1209 zmax[] = {10., 79., 33.};
1210 cOut = new TCanvas(Form("%s_%sTrkInSpc", GetName(), prefix), "Track IN PID", Int_t(1.5*ny), Int_t(1.5*ny));
1211 cOut->Divide(3,3, 1.e-5, 1.e-5);
1212 Int_t nplot(0); const Char_t *chName[] = {"p", "n", ""};
1213 for(Int_t ich(0), ipad(1); ich<3; ich++){
1215 if(!(h2s = (TH2*)arr->FindObject(Form("H%sTrkInY%sEn", prefix, chName[ich])))) {
1216 AliInfo(Form("Missing H%sTrkIn%sEn", prefix, chName[ich]));
1219 Int_t irebin(0), dxBin(1), dyBin(1);
1220 const Int_t nrebin(5); Int_t rebinX[nrebin] = {1, 2, 1, 2, 1}, rebinY[nrebin] = {2, 1, 2, 1, 2};
1221 if(h2s->GetNbinsY()==180){ // old binning
1222 rebinX[0] = 1; rebinY[0] = 2;
1223 rebinX[1] = 2; rebinY[1] = 1;
1224 rebinX[2] = 2; rebinY[2] = 1;
1225 rebinX[3] = 1; rebinY[3] = 5;
1226 rebinX[4] = 1; rebinY[4] = 1; // dummy
1228 while(irebin<nrebin && (AliTRDrecoTask::GetMeanStat(h2s, .5, 1)<200)){
1229 h2s->Rebin2D(rebinX[irebin], rebinY[irebin]);
1230 dxBin*=rebinX[irebin];dyBin*=rebinY[irebin];irebin++;
1232 AliDebug(2, Form("Rebin level[%d] @ chg[%d]. Binning[%2dx%2d]", irebin, ich, h2s->GetNbinsX(), h2s->GetNbinsY()));
1233 for(Int_t ispec(0); ispec<kNspc; ispec++){
1234 p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1235 if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkInY%s%dEn", prefix, chName[ich], ispec)))) {
1236 AliInfo(Form("Missing H%sTrkIn%s%dEn", prefix, chName[ich], ispec));
1240 h2->Rebin2D(dxBin,dyBin);
1241 h2->Divide(h2, h2s, 1.e2);
1243 h2->GetZaxis()->SetRangeUser(zmin[ispec], zmax[ispec]);
1244 h2->GetZaxis()->SetTitle("Rel. Abundancy [%]");h2->GetZaxis()->CenterTitle();
1245 TString tit(h2->GetTitle()); TObjArray *atit = tit.Tokenize("::");
1246 h2->SetTitle(Form("%s :: Relative Abundancy", ((*atit)[0])->GetName()));
1247 atit->Delete(); delete atit;
1249 MakeDetectorPlot(0);
1252 if(nplot==9) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1254 // pt resolution plot
1255 cOut = new TCanvas(Form("%s_%sTrkInPt", GetName(), prefix), "TrackIn Pt", Int_t(1.5*ny), Int_t(1.5*ny));
1256 cOut->Divide(3,2, 1.e-5, 1.e-5);
1257 for(Int_t ich(0), ipad(1); ich<2; ich++){
1258 for(Int_t ispc(0); ispc<kNspc; ispc++){
1259 p=cOut->cd(ipad++); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1260 if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkInPt%s%d_2D", prefix, chName[ich], ispc)))) continue;
1261 h2->GetZaxis()->SetRangeUser(0.5, 1.2); h2->SetContour(7); h2->Draw("colz");
1262 MakeDetectorPlot(0);
1265 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1266 // MPV(Q) & <Q plots>
1267 const char *chQ[] = {"Q", "QS"};
1268 for(Int_t iq(0); iq<2; iq++){
1269 cOut = new TCanvas(Form("%s_%sTrkIn%s", GetName(), prefix, chQ[iq]), "Track IN PID", Int_t(1.5*ny), Int_t(1.5*ny));
1270 cOut->Divide(3,3, 1.e-5, 1.e-5);
1272 for(Int_t ich(0), ipad(1); ich<3; ich++){
1273 for(Int_t ispec(0); ispec<kNspc; ispec++){
1274 p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1275 if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkIn%s%s%d_2D", prefix, chQ[iq], chName[ich], ispec)))) {
1276 AliInfo(Form("Missing H%sTrkIn%s%s%d_2D", prefix, chQ[iq], chName[ich], ispec));
1281 MakeDetectorPlot(0);
1284 if(nplot==9) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1290 // track MC systematic
1291 if((arr = (TObjArray*)fProj->FindObject("hTRD2MC"))) {
1292 for(Int_t iview(0); iview<nTrkViews; iview++){
1293 cOut = new TCanvas(Form("%s_MC%s", GetName(), vTrkName[iview]), "Track Resolution", nx, ny);
1294 cOut->Divide(3,2, 1.e-5, 1.e-5);
1296 for(Int_t iplot(0); iplot<6; iplot++){
1297 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1298 if(!(h2 = (TH2*)arr->FindObject(Form("HMC%s%d_2D", vTrkName[iview], iplot)))) continue;
1299 h2->Draw("colz"); nplot++;
1301 if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1305 gStyle->SetPalette(1);
1308 for(Int_t ih(ih2e); ih--;) delete h2e[ih];
1311 //________________________________________________________
1312 TH2* AliTRDresolution::DrawSigma(TH2 *h2, const Char_t *title, Float_t m, Float_t M, Float_t scale)
1314 // Draw error bars scaled with "scale" instead of content values
1315 //use range [m,M] if limits are specified
1317 if(!h2) return NULL;
1318 TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1319 TH2F *h2e = new TH2F(Form("%s_E", h2->GetName()),
1320 Form("%s;%s;%s;%s", h2->GetTitle(), ax->GetTitle(), ay->GetTitle(), title),
1321 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
1322 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1324 TAxis *az(h2e->GetZaxis());
1325 if(M>m) az->SetRangeUser(m, M);
1327 az->SetTitleOffset(1.5);
1328 for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1329 for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1330 if(h2->GetBinContent(ix, iy)<-100.) continue;
1331 Float_t v(scale*h2->GetBinError(ix, iy));
1332 if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1333 h2e->SetBinContent(ix, iy, v);
1340 //________________________________________________________
1341 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1343 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1344 // The "range" vector should be initialized with 2 elements
1345 // Option "mod" can be any of
1346 // - 0 : gaussian like distribution
1347 // - 1 : tailed distribution
1349 Int_t nx(h2->GetNbinsX())
1350 , ny(h2->GetNbinsY())
1352 Double_t *data=new Double_t[n];
1353 for(Int_t ix(1), in(0); ix<=nx; ix++){
1354 for(Int_t iy(1); iy<=ny; iy++)
1355 data[in++] = h2->GetBinContent(ix, iy);
1357 Double_t mean, sigm;
1358 AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1360 range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1361 if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1362 AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1363 TH1S h1("h1SF0", "", 100, range[0], range[1]);
1368 case 0:// gaussian distribution
1370 TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1372 mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1373 range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1374 AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1]));
1377 case 1:// tailed distribution
1379 Int_t bmax(h1.GetMaximumBin());
1380 Int_t jBinMin(1), jBinMax(100);
1381 for(Int_t ibin(bmax); ibin--;){
1382 if(h1.GetBinContent(ibin)<1.){
1383 jBinMin=ibin; break;
1386 for(Int_t ibin(bmax); ibin++;){
1387 if(h1.GetBinContent(ibin)<1.){
1388 jBinMax=ibin; break;
1391 range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1392 AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1]));
1400 //________________________________________________________
1401 Bool_t AliTRDresolution::MakeProjectionDetector()
1404 const Int_t kNcontours(9);
1405 const Int_t kNstat(100);
1406 if(fProj && fProj->At(kDetector)) return kTRUE;
1408 AliError("Missing data container.");
1412 if(!(H = (THnSparse*)fContainer->FindObject("hDet2Cluster"))){
1413 AliInfo("Missing/Wrong data @ hDet2Cluster.");
1416 Int_t ndim(H->GetNdimensions());
1417 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1418 TAxis *aa[kNdim], *an(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1419 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1420 if(ndim < 5) aa[4] = new TAxis(1, -0.5, 0.5);
1426 // build list of projections
1427 const Int_t nsel(8*AliTRDgeometry::kNlayer*AliTRDeventInfo::kCentralityClasses);
1428 // define rebinning strategy
1429 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 1, 2, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1430 if(aa[1]->GetNbins()==180){ // old binning
1431 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1432 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1433 rebinEtaPhiX[2] = 2; rebinEtaPhiY[2] = 1;
1434 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1435 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1437 //const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1438 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1439 AliTRDrecoProjection hp[kDetNproj]; TObjArray php(kDetNproj);
1440 Int_t ih(0), isel(-1), np[nsel]={0};
1441 for(Int_t ipad(0); ipad<nPad; ipad++){
1442 for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1443 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1444 isel++; // new selection
1445 hp[ih].Build(Form("HDet%d%d%d", ily, icen, ipad),
1446 Form("Detectors :: Ly[%d] Cen[%s] Pads[%s]", ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1448 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1449 hp[ih].SetShowRange(10., 55.);
1450 php.AddLast(&hp[ih++]); np[isel]++;
1454 AliInfo(Form("Build %3d 3D projections.", ih));
1456 Int_t ly(0), cen(0), npad(0);
1457 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1458 v = H->GetBinContent(ib, coord); if(v<1.) continue;
1459 ly = coord[kBC]-1; // layer selection
1460 cen = coord[kYrez]-1; // centrality selection
1461 npad = 0; // no. of pads selection
1462 if(an) npad = TMath::Min(coord[5]-1, Int_t(kNpads));
1463 isel = npad*AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer+cen*AliTRDgeometry::kNlayer+ly;
1464 ((AliTRDrecoProjection*)php.At(isel))->Increment(coord, v);
1465 //Int_t ioff=isel;for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1467 TObjArray *arr(NULL);
1468 fProj->AddAt(arr = new TObjArray(kDetNproj), kDetector);
1469 arr->SetName("hDet2Cluster"); arr->SetOwner();
1471 TH2 *h2(NULL); Int_t jh(0);
1473 if(!hp[ih].H()) continue;
1474 if((h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) arr->AddAt(h2, jh++);
1475 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].H()->GetName())))) arr->AddAt(h2, jh++);
1477 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1478 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1479 for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1480 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, 0)))){
1481 for(Int_t ipad(1); ipad<nPad; ipad++){
1482 if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, ipad)))){
1486 pr0->H()->SetNameTitle(Form("HDet%d%d", ily, icen), Form("Detectors :: Ly[%d] Cen[%s]", ily, cenName[icen]));
1487 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0, kFALSE))) arr->AddAt(h2, jh++);
1488 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
1489 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))) (*pr1)+=(*pr0);
1492 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))){
1493 pr0->H()->SetNameTitle(Form("HDet%d", UseLYselectTrklt()?fLYselect:ily), Form("Detectors :: Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
1494 if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1497 AliInfo(Form("Done %3d 2D projections.", jh));
1501 //________________________________________________________
1502 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1505 const Int_t kNcontours(9);
1506 const Int_t kNstat(100);
1507 Int_t cidx=mc?kMCcluster:kCluster;
1508 if(fProj && fProj->At(cidx)) return kTRUE;
1510 AliError("Missing data container.");
1513 const Char_t *projName[] = {"hCluster2Track", "hCluster2MC"};
1515 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1516 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1519 Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimCl));
1520 Int_t coord[10]; memset(coord, 0, sizeof(Int_t) * 10); Double_t v = 0.;
1521 TAxis *aa[kNdim], *as(NULL), *apt(NULL), *acen(NULL), *anp(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1522 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1523 if(ndim > Int_t(kPt)) apt = H->GetAxis(kPt);
1524 if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC);
1525 if(ndim > Int_t(kNdim)) acen = H->GetAxis(kNdim);
1526 if(ndim > Int_t(kNdim)+1) anp = H->GetAxis(kNdim+1);
1527 // calculate size depending on debug level
1528 const Int_t nCh(apt?2:1);
1529 const Int_t nCen(acen?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1530 const Int_t nNpad(anp?(Int_t(kNpads)+1):1);
1532 // build list of projections
1533 const Int_t nsel(AliTRDgeometry::kNlayer*kNcharge*(kNpads+1)*AliTRDeventInfo::kCentralityClasses);
1534 // define rebinning strategy
1535 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1536 if(aa[1]->GetNbins()==180){ // old binning
1537 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1538 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1539 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
1540 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1541 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1543 AliTRDrecoProjection hp[kClNproj]; TObjArray php(kClNproj);
1544 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1545 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1546 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1547 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1548 for(Int_t ich(0); ich<nCh; ich++){
1549 for(Int_t icen(0); icen<nCen; icen++){
1550 for(Int_t ipad(0); ipad<nNpad; ipad++){
1551 isel++; // new selection
1552 hp[ih].Build(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1553 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"),
1554 kEta, kPhi, kYrez, aa);
1555 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1556 php.AddLast(&hp[ih++]); np[isel]++;
1558 hp[ih].Build(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1559 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"),
1560 kEta, kPhi, kSpeciesChgRC, aa);
1561 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1562 hp[ih].SetShowRange(24., 33.);
1563 php.AddLast(&hp[ih++]); np[isel]++;
1564 hp[ih].Build(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1565 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"),
1566 kPrez, kZrez, kYrez, aa);
1567 php.AddLast(&hp[ih++]); np[isel]++;
1568 hp[ih].Build(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1569 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"),
1570 kPrez, kPt, kYrez, aa);
1571 php.AddLast(&hp[ih++]); np[isel]++;
1576 AliInfo(Form("Build %3d 3D %s projections.", ih, mc?"MC":""));
1578 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1579 Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1), ioff(0), cen(0), npad(0);
1580 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1581 v = H->GetBinContent(ib, coord); if(v<1.) continue;
1584 if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1587 ch = 0; // [-] track
1588 if(chBin>0 && coord[kPt] > chBin) ch = 1; // [+] track
1589 cen = 0; // centrality selection
1590 if(acen) cen = coord[kNdim]-1;
1591 npad = 0; // no. of pads selection
1592 if(anp) npad = TMath::Min(coord[kNdim+1]-1, Int_t(kNpads));
1595 isel = ly*nCh*nCen*nNpad
1601 if(mc && ly<0) ly=0; // fix for bug in PlotMC
1602 isel = ly; ioff = isel;
1605 AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1608 if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
1609 AliError(Form("Missing projection %d", ioff));
1612 if(strcmp(pr0->H()->GetName(), Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ch], ly, cen, npad))!=0){
1613 AliError(Form("Projection mismatch :: request[H%sClY%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ly, cen, npad, pr0->H()->GetName()));
1616 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1618 if(HasDump3DFor(kCluster)){
1619 TDirectory *cwd = gDirectory;
1620 TFile::Open(Form("DumpRes_%s.root", H->GetName()), "RECREATE");
1621 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
1622 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
1623 if(!pr0->H()) continue;
1624 TH3 *h3=(TH3*)pr0->H()->Clone();
1631 TObjArray *arr(NULL);
1632 fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1633 arr->SetName(projName[Int_t(mc)]); arr->SetOwner();
1635 TH2 *h2(NULL); Int_t jh(0);
1637 if(!hp[ih].H()) continue;
1638 if(strchr(hp[ih].H()->GetName(), 'Q')){
1639 if((h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) arr->AddAt(h2, jh++);
1640 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].H()->GetName())))) arr->AddAt(h2, jh++);
1642 if((h2 = hp[ih].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1645 Double_t m(0.), s(0.), se(0.), trend(0.);
1646 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1647 for(Int_t ich(0); ich<nCh; ich++){
1648 for(Int_t icen(0); icen<nCen; icen++){
1650 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1651 for(Int_t ipad(1); ipad<nNpad; ipad++){
1652 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1655 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]));
1656 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1657 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1660 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1661 for(Int_t ipad(1); ipad<nNpad; ipad++){
1662 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1665 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]));
1666 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2, kFALSE))) arr->AddAt(h2, jh++);
1667 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
1668 pr0->H()->SetName(Form("H%sClQS%c%d%d", mc?"MC":"", chName[ich], ily, icen));
1669 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1670 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1673 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1674 for(Int_t ipad(1); ipad<nNpad; ipad++){
1675 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1678 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]));
1679 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1680 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1683 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1684 for(Int_t ipad(1); ipad<nNpad; ipad++){
1685 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad)))) continue;
1688 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]));
1689 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1690 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1692 } // end centrality integration
1694 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1695 pr0->H()->SetNameTitle(Form("H%sClY%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily),
1696 Form("Clusters[%c]:: #Deltay Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
1697 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1698 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1701 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1702 pr0->H()->SetNameTitle(Form("H%sClQ%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily),
1703 Form("Clusters[%c]:: Q Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
1704 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
1705 pr0->H()->SetName(Form("H%sClQS%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily));
1706 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1707 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1710 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1711 pr0->H()->SetNameTitle(Form("H%sClYXTC%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily),
1712 Form("Clusters[%c]:: #Deltay(x,TC) Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
1713 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1714 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1717 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1718 pr0->H()->SetNameTitle(Form("H%sClYXPh%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily),
1719 Form("Clusters[%c]:: #Deltay(x,#Phi) Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
1720 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1721 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1723 } // end charge integration
1725 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1726 pr0->H()->SetNameTitle(Form("H%sClY%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), Form("Clusters :: #Deltay Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
1727 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1728 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.) PutTrendValue(Form("%sClS%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), s, se);
1731 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1732 pr0->H()->SetNameTitle(Form("H%sClYXPh%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), Form("Clusters :: #Deltay Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
1733 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1737 AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
1741 //________________________________________________________
1742 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1745 const Int_t kNcontours(9);
1746 const Int_t kNstat(30);
1747 const Int_t kNstatQ(30);
1748 Int_t cidx=mc?kMCtracklet:kTracklet;
1750 AliError("Missing results container.");
1754 AliError("Missing data container.");
1757 const Char_t *projName[] = {"hTracklet2Track", "hTracklet2MC"};
1759 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1760 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1763 const Int_t mdim(kNdim+8);
1764 Int_t ndim(H->GetNdimensions()); //Bool_t debug(ndim>Int_t(kNdimTrklt));
1765 Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
1766 TAxis *aa[mdim], *as(NULL), *ap(NULL), *ac(NULL); memset(aa, 0, sizeof(TAxis*) * mdim);
1767 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1768 if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC); // init species/charge selection
1769 if(ndim > Int_t(kPt)) ap = H->GetAxis(kPt); // init pt selection
1770 // if(ndim > Int_t(kNdim)) ac = H->GetAxis(kNdim); // init centrality selection
1771 // calculate size depending on debug level
1772 const Int_t nCen(ac?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1773 const Int_t nPt(ap?Int_t(fNpt+2):1);
1774 const Int_t nSpc(as?Int_t(kNspc):1);
1775 const Int_t nCh(as?Int_t(kNcharge):1);
1776 const Int_t nLy(UseLYselectTrklt()?1:AliTRDgeometry::kNlayer);
1778 // build list of projections
1779 const Int_t nsel(AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer*(fgNPt+2)*(AliPID::kSPECIES*kNcharge + 1));
1780 // define rebinning strategy
1781 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
1782 if(aa[1]->GetNbins()==180){ // old binning
1783 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
1784 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
1785 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
1786 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
1787 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
1789 AliTRDrecoProjection hp[kTrkltNproj]; TObjArray php(kTrkltNproj);
1790 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1791 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1792 const Char_t ptShortName[5] = {'L', 'l', 'i', 'h', 'H'};
1793 Char_t ptName[fgNPt+2] = {0};
1794 Char_t *ptCut[fgNPt+2] = {NULL};
1796 // const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1797 const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1798 for(Int_t icen(0); icen<nCen; icen++){
1799 for(Int_t ily(0); ily<nLy; ily++){
1800 for(Int_t ipt(0); ipt<nPt; ipt++){
1802 ptName[ipt]= nPt>5?(64+ipt):ptShortName[ipt];
1803 ptCut[ipt] = StrDup(Form("#it{%4.2f<=p_{t}^{%s}[GeV/c]<%4.2f}",ipt?fgPt[ipt-1]:0., mc?"MC":"", ipt>fgNPt?99.99:fgPt[ipt]));
1805 for(Int_t isp(0); isp<nSpc; isp++){
1806 for(Int_t ich(0); ich<nCh; ich++){
1807 isel++; // new selection
1808 AliDebug(3, Form("Building sel[%3d|%4d]", isel, ih));
1809 hp[ih].Build(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily, icen),
1810 Form("Tracklets[%s%c]:: #Deltay{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]),
1811 kEta, kPhi, kYrez, aa);
1812 //hp[ih].SetShowRange(-0.1,0.1);
1813 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1814 php.AddLast(&hp[ih++]); np[isel]++;
1815 hp[ih].Build(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily, icen),
1816 Form("Tracklets[%s%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]),
1817 kEta, kPhi, kPrez, aa);
1818 //hp[ih].SetShowRange(-0.5,0.5);
1819 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1820 php.AddLast(&hp[ih++]); np[isel]++;
1821 hp[ih].Build(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily, icen),
1822 Form("Tracklets[%s%c]:: dQdl{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]),
1823 kEta, kPhi, kZrez, aa);
1824 hp[ih].SetShowRange(1.,2.3);
1825 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1826 php.AddLast(&hp[ih++]); np[isel]++;
1827 // hp[ih].Build(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1828 // Form("Tracklets[%s%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1829 // kEta, kPhi, kNdim+1, aa);
1830 // hp[ih].SetShowRange(30., 70.);
1831 // hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1832 // php.AddLast(&hp[ih++]); np[isel]++;
1835 if(ndim==kNdimTrklt) continue;
1837 isel++; // new selection
1838 AliDebug(3, Form("Building selRC[%3d|%4d]", isel, ih));
1839 hp[ih].Build(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], UseLYselectTrklt()?fLYselect:ily, icen),
1840 Form("Tracklets[RC]:: #Deltaz{%s} Ly[%d] Cen[%s]", ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]),
1841 kEta, kPhi, kZrez, aa);
1842 // hp[ih].SetShowRange(-0.1,0.1);
1843 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1844 php.AddLast(&hp[ih++]); np[isel]++;
1845 /* hp[ih].Build(Form("H%sTrkltRCY%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1846 Form("Tracklets[RC]:: #Deltay{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1847 kEta, kPhi, kYrez, aa);
1848 //hp[ih].SetShowRange(-0.1,0.1);
1849 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1850 php.AddLast(&hp[ih++]); np[isel]++;
1851 hp[ih].Build(Form("H%sTrkltRCPh%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1852 Form("Tracklets[RC]:: #Delta#phi{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1853 kEta, kPhi, kPrez, aa);
1854 //hp[ih].SetShowRange(-0.1,0.1);
1855 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1856 php.AddLast(&hp[ih++]); np[isel]++;
1857 hp[ih].Build(Form("H%sTrkltRCQ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1858 Form("Tracklets[RC]:: dQdl{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1859 kEta, kPhi, kNdim, aa);
1860 //hp[ih].SetShowRange(-0.1,0.1);
1861 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1862 php.AddLast(&hp[ih++]); np[isel]++;*/
1866 AliInfo(Form("Build %3d 3D %s projections%s", ih, mc?"MC":"", UseLYselectTrklt()?Form(" for Ly[%d].", fLYselect):"."));
1868 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
1869 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*3/*4*/+1);
1870 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1871 v = H->GetBinContent(ib, coord);
1873 ly = coord[kBC]-1; // layer selection
1874 if(UseLYselectTrklt()&& (ly!=fLYselect)) continue;
1876 ch = 0; sp=0;// [e-] track [dafault]
1877 if(rcBin>0){ // debug mode in which species/charge are also saved
1878 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1879 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
1880 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
1884 if(ap) pt = coord[kPt];//TMath::Min(coord[kPt], Int_t(kNpt)+1);
1885 // centrality selection
1887 if(ac) cen = coord[kNdim]-1;
1889 if(ndim==kNdimTrklt){
1890 if(mc && ly<0) ly=0; // fix for bug in PlotMC
1894 isel = cen*nLy*nPt*jspc+(!UseLYselectTrklt())*ly*nPt*jspc+pt*jspc; isel+=sp<0?(nSpc*nCh):(sp*nCh+ch);
1895 ioff = cen*nLy*nPt*kspc+(!UseLYselectTrklt())*ly*nPt*kspc+pt*kspc; ioff+=sp<0?((nSpc*nCh)*3/*4*/):(sp*nCh*3+ch*3/*4*/);
1898 AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1901 if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
1902 AliError(Form("Missing projection %d", ioff));
1906 if(strcmp(pr0->H()->GetName(), Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ch], ptName[pt], sp, UseLYselectTrklt()?fLYselect:ly, cen))!=0){
1907 AliError(Form("Projection mismatch :: request[H%sTrkltY%c%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, UseLYselectTrklt()?fLYselect:ly, cen, pr0->H()->GetName()));
1911 if(strcmp(pr0->H()->GetName(), Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[pt], UseLYselectTrklt()?fLYselect:ly, cen))!=0){
1912 AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c%d%d] found[%s]", mc?"MC":"", ptName[pt], UseLYselectTrklt()?fLYselect:ly, cen, pr0->H()->GetName()));
1916 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
1918 if(HasDump3DFor(kTracklet)){
1919 TDirectory *cwd = gDirectory;
1920 TFile::Open(Form("DumpRes_%s.root", H->GetName()), "RECREATE");
1921 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
1922 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
1923 if(!pr0->H()) continue;
1924 TH3 *h3=(TH3*)pr0->H()->Clone();
1931 TH2 *h2(NULL); Int_t jh(0);
1932 TObjArray *arr(NULL);
1933 if(!(arr = (TObjArray*)fProj->At(cidx))){
1934 fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1935 arr->SetName(projName[Int_t(mc)]); arr->SetOwner();
1936 } else jh = arr->GetEntriesFast();
1939 if(!hp[ih].H()) continue;
1940 Int_t mid(0), nstat(kNstat);
1941 if(strchr(hp[ih].H()->GetName(), 'Q')){ mid=2; nstat=kNstatQ;}
1942 if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1943 arr->AddAt(h2, jh++);
1945 // build combined performance plots
1946 Double_t m(0.), s(0.), se(0.), trend(0.);
1947 for(Int_t ily(0); ily<nLy; ily++){
1948 for(Int_t ich(0); ich<nCh; ich++){
1949 for(Int_t ipt(0); ipt<nPt; ipt++){
1950 for(Int_t icen(0); icen<nCen; icen++){
1952 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, icen)))){
1953 for(Int_t isp(1); isp<nSpc; isp++){
1954 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily, icen)))) continue;
1957 pr0->H()->SetNameTitle(Form("H%sTrkltY%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily, icen),
1958 Form("Tracklets[%c]:: #Deltay{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]));
1959 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1960 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
1963 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, icen)))){
1964 for(Int_t isp(1); isp<nSpc; isp++){
1965 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily, icen)))) continue;
1968 pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily, icen),
1969 Form("Tracklets[%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]));
1970 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1971 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
1974 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, icen)))){
1975 for(Int_t isp(1); isp<nSpc; isp++){
1976 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily, icen)))) continue;
1979 pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily, icen),
1980 Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]));
1981 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1982 pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily, icen),
1983 Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]));
1984 pr0->SetShowRange(2.4, 5.1);
1985 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1986 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
1989 // if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, icen)))){
1990 // for(Int_t isp(1); isp<nSpc; isp++){
1991 // if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily, icen)))) continue;
1994 // pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily, icen),
1995 // Form("Tracklets[%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily, cenName[icen]));
1996 // if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1997 // if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
1999 } // end centrality integration
2001 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2002 pr0->H()->SetNameTitle(Form("H%sTrkltY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2003 Form("Tracklets[%c]:: #Deltay Ly[%d] %s", chSgn[ich], UseLYselectTrklt()?fLYselect:ily, ptCut[ipt]));
2004 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2005 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2006 PutTrendValue(Form("%sTrkltY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily), trend, m);
2007 PutTrendValue(Form("%sTrkltYS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily), s, se);
2009 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2012 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2013 pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2014 Form("Tracklets[%c]:: #Delta#phi Ly[%d] %s", chSgn[ich], UseLYselectTrklt()?fLYselect:ily, ptCut[ipt]));
2015 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2016 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2017 PutTrendValue(Form("%sTrkltPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily), trend, m);
2018 PutTrendValue(Form("%sTrkltPhS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily), s, se);
2020 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2023 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2024 pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2025 Form("Tracklets[%c]:: dQdl Ly[%d] %s", chSgn[ich], UseLYselectTrklt()?fLYselect:ily, ptCut[ipt]));
2026 pr0->SetShowRange(1.,2.3);
2027 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
2028 pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2029 Form("Tracklets[%c]:: dQdl Ly[%d] %s", chSgn[ich], UseLYselectTrklt()?fLYselect:ily, ptCut[ipt]));
2030 pr0->SetShowRange(2.4,5.1);
2031 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2032 if((trend=pr0->GetTrendValue(2,&m,&s,&se))>-100.){
2033 PutTrendValue(Form("%sTrkltQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily), trend, m);
2034 PutTrendValue(Form("%sTrkltQS%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily), s, se);
2036 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2039 // if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2040 // pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2041 // Form("Tracklets[%c]:: OccupancyTB Ly[%d] %s", chSgn[ich], UseLYselectTrklt()?fLYselect:ily, ptCut[ipt]));
2042 // if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
2043 // if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2045 } // end pt integration
2047 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2048 pr0->H()->SetNameTitle(Form("H%sTrkltY%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), Form("Tracklets[%c] :: #Deltay Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2049 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2050 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2051 PutTrendValue(Form("%sTrkltY%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), trend, m);
2052 PutTrendValue(Form("%sTrkltYS%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), s, se);
2054 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2057 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2058 pr0->H()->SetNameTitle(Form("H%sTrkltPh%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), Form("Tracklets[%c] :: #Delta#phi Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2059 pr0->SetShowRange(-.9,.9);
2060 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2061 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2062 PutTrendValue(Form("%sTrkltPh%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), trend, m);
2063 PutTrendValue(Form("%sTrkltPhS%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), s, se);
2065 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2068 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2069 pr0->H()->SetNameTitle(Form("H%sTrkltQ%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2070 pr0->SetShowRange(1.,2.3);
2071 if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
2072 pr0->H()->SetNameTitle(Form("H%sTrkltQS%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2073 pr0->SetShowRange(2.4,5.1);
2074 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2075 if((trend=pr0->GetTrendValue(2,&m,&s,&se))>-100.){
2076 PutTrendValue(Form("%sTrkltQ%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), trend, m);
2077 PutTrendValue(Form("%sTrkltQS%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), s, se);
2079 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2082 // if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2083 // pr0->H()->SetNameTitle(Form("H%sTrkltTB%c%d", mc?"MC":"", chName[ich], UseLYselectTrklt()?fLYselect:ily), Form("Tracklets[%c] :: OccupancyTB Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2084 // if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
2085 // if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2087 } // end charge integration
2089 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2090 pr0->H()->SetNameTitle(Form("H%sTrkltY%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), Form("Tracklets :: #Deltay Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2091 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2092 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2093 PutTrendValue(Form("%sTrkltY%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), trend, m);
2094 PutTrendValue(Form("%sTrkltYS%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), s, se);
2098 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2099 pr0->H()->SetNameTitle(Form("H%sTrkltPh%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), Form("Tracklets :: #Delta#phi Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2100 pr0->SetShowRange(-.45,.45);
2101 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2104 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2105 pr0->H()->SetNameTitle(Form("H%sTrkltQ%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), Form("Tracklets :: dQdl Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2106 pr0->SetShowRange(1.,2.3);
2107 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2108 pr0->H()->SetName(Form("H%sTrkltQS%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily));
2109 pr0->SetShowRange(2.4,5.1);
2110 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2113 // if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily, 0)))){
2114 // pr0->H()->SetNameTitle(Form("H%sTrkltTB%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), Form("Tracklets :: OccupancyTB Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2115 // if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
2118 /*! Row Cross processing*/
2119 for(Int_t icen(0); icen<nCen; icen++){
2121 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], UseLYselectTrklt()?fLYselect:ily, icen)))){
2122 for(Int_t ipt(0); ipt<nPt; ipt++){
2123 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], UseLYselectTrklt()?fLYselect:ily, icen)))) continue;
2126 pr0->H()->SetNameTitle(Form("H%sTrkltRCZ%d%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily, icen), Form("Tracklets[RC]:: #Deltaz Ly[%d] Cen[%s]", UseLYselectTrklt()?fLYselect:ily, cenName[icen]));
2127 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2128 if(icen && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], UseLYselectTrklt()?fLYselect:ily, 0)))) (*pr1)+=(*pr0);
2130 } // end centrality integration for row cross
2132 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], UseLYselectTrklt()?fLYselect:ily, 0)))){
2133 pr0->H()->SetNameTitle(Form("H%sTrkltRCZ%d", mc?"MC":"", UseLYselectTrklt()?fLYselect:ily), Form("Tracklets[RC] :: #Deltaz Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2134 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2137 AliInfo(Form("Done %3d 2D %s projections.", jh, mc?"MC":""));
2139 // clean local memory allocation
2140 for(Int_t i(fgNPt+2);i--;) if(ptCut[i]) delete [] ptCut[i];
2144 //________________________________________________________
2145 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc, Bool_t v0)
2148 const Int_t kNcontours(9);
2149 const Int_t kNstat(30);
2150 Int_t cidx=mc?kMCtrackIn:(v0?(kV0TrackIn-1):kTrackIn);
2151 if(fProj && fProj->At(cidx)){
2152 AliInfo(Form("Nothing to do for container %s.", ((TObjArray*)fProj->At(cidx))->GetName() ));
2156 AliError("Missing data container.");
2159 const Char_t *projName[] = {"hTracklet2TRDin", "hTRDin2MC", "hTracklet2TRDinV0"};
2161 if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
2162 AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
2165 const Char_t *prefix = mc?"MC":(v0?"V0":"");
2167 const Int_t mdim(kNdim+3);
2168 Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
2169 Int_t ndim(H->GetNdimensions());
2170 TAxis *aa[mdim], *as(NULL), *ap(NULL), *apt(NULL), *ax(NULL), *abf(NULL); memset(aa, 0, sizeof(TAxis*) * (mdim));
2171 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2172 if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC);
2173 if(ndim > Int_t(kPt)) apt = H->GetAxis(kPt);
2174 if(ndim > Int_t(kNdim)) ap = H->GetAxis(kNdim);
2175 if(ndim > Int_t(kNdim)+1) ax = H->GetAxis(kNdim+1);
2176 if(ndim > Int_t(kNdim)+2) abf = H->GetAxis(kNdim+2);
2177 //AliInfo(Form("Using : Species[%c] Pt[%c] BunchFill[%c]", as?'y':'n', ap?'y':'n', abf?'y':'n'));
2178 const Int_t nPt(apt?(apt->GetNbins()+2):1);
2180 // build list of projections
2181 const Int_t nsel((fgNPt+2)*(kNspc*kNcharge + 1) + kNspc*kNcharge);
2182 const Int_t nprj((fgNPt+2)*(kNspc*kNcharge + 1)*4 + kNspc*kNcharge);
2183 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2184 const Char_t *spcName[2][kNspc] = {{"e", "#mu#pi", "Kp"},
2186 const Char_t ptShortName[5] = {'L', 'l', 'i', 'h', 'H'};
2187 Char_t ptName[fgNPt+2] = {0};
2188 Char_t *ptCut[fgNPt+2] = {NULL};
2189 Char_t *pCut[fgNPt+2] = {NULL};
2190 // define rebinning strategy
2191 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
2192 if(aa[1]->GetNbins()==180){ // old binning
2193 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
2194 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
2195 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
2196 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
2197 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
2199 AliTRDrecoProjection hp[nprj]; TObjArray php(nprj+2);
2200 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2201 // define list of projections
2202 for(Int_t ipt(0); ipt<nPt; ipt++){
2203 ptName[ipt]= nPt>5?(64+ipt):ptShortName[ipt];
2204 ptCut[ipt] = StrDup(Form("#it{%4.2f<=p_{t}^{%s}[GeV/c]<%4.2f}",ipt?fgPt[ipt-1]:0., mc?"MC":"", ipt>fgNPt?99.99:fgPt[ipt]));
2205 pCut[ipt] = StrDup(Form("#it{%4.2f<=p^{%s}[GeV/c]<%4.2f}",ipt?fgPt[ipt-1]:0., mc?"MC":"", ipt>fgNPt?99.99:fgPt[ipt]));
2206 for(Int_t isp(0); isp<kNspc; isp++){
2207 for(Int_t ich(0); ich<kNcharge; ich++){
2208 isel++; // new selection
2209 AliDebug(3, Form("Building sel[%3d|%4d] spc[%s%c] pt[%s]", isel, ih, spcName[v0][isp], chSgn[ich], ptCut[ipt]));
2210 hp[ih].Build(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[ipt], isp),
2211 Form("TrackIn[%s%c]:: #Deltay{%s}", spcName[v0][isp], chSgn[ich], ptCut[ipt]),
2212 kEta, kPhi, kYrez, aa);
2213 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2214 php.AddLast(&hp[ih++]); np[isel]++;
2215 hp[ih].Build(Form("H%sTrkInPh%c%c%d", prefix, chName[ich], ptName[ipt], isp),
2216 Form("TrackIn[%s%c]:: #Delta#phi{%s}", spcName[v0][isp], chSgn[ich], ptCut[ipt]),
2217 kEta, kPhi, kPrez, aa);
2218 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2219 php.AddLast(&hp[ih++]); np[isel]++;
2220 hp[ih].Build(Form("H%sTrkInQ%c%c%d", prefix, chName[ich], ptName[ipt], isp),
2221 Form("TrackIn[%s%c]:: dQdl {%s}", spcName[v0][isp], chSgn[ich], pCut[ipt]),
2222 kEta, kPhi, kZrez, aa);
2223 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2224 php.AddLast(&hp[ih++]); np[isel]++;
2226 hp[ih].Build(Form("H%sTrkInX%c%c%d", prefix, chName[ich], ptName[ipt], isp),
2227 Form("TrackIn[%s%c]:: #Deltax{%s}", spcName[v0][isp], chSgn[ich], ptCut[ipt]),
2228 kEta, kPhi, kNdim+1, aa);
2229 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2230 php.AddLast(&hp[ih++]); np[isel]++;
2233 isel++; // RC projections
2234 AliDebug(3, Form("Building RCsel[%3d] pt[%s]", isel, ptCut[ipt]));
2235 hp[ih].Build(Form("H%sTrkInRCZ%c", prefix, ptName[ipt]),
2236 Form("TrackIn[RC]:: #Deltaz{%s}", ptCut[ipt]),
2237 kEta, kPhi, kZrez, aa);
2238 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2239 php.AddLast(&hp[ih++]); np[isel]++;
2240 hp[ih].Build(Form("H%sTrkInRCY%c", prefix, ptName[ipt]),
2241 Form("TrackIn[RC]:: #Deltay{%s}", ptCut[ipt]),
2242 kEta, kPhi, kYrez, aa);
2243 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2244 php.AddLast(&hp[ih++]); np[isel]++;
2245 hp[ih].Build(Form("H%sTrkInRCPh%c", prefix, ptName[ipt]),
2246 Form("TrackIn[RC]:: #Delta#phi{%s}", ptCut[ipt]),
2247 kEta, kPhi, kPrez, aa);
2248 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2249 php.AddLast(&hp[ih++]); np[isel]++;
2251 hp[ih].Build(Form("H%sTrkInRCX%c", prefix, ptName[ipt]),
2252 Form("TrackIn[RC]:: #Deltax{%s}", ptCut[ipt]),
2253 kEta, kPhi, kNdim+1, aa);
2254 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2255 php.AddLast(&hp[ih++]); np[isel]++;
2258 for(Int_t isp(0); isp<kNspc; isp++){
2259 for(Int_t ich(0); ich<kNcharge; ich++){
2260 isel++; // new selection
2261 AliDebug(3, Form("Building PTsel[%3d] spc[%s%c]", isel, spcName[v0][isp], chSgn[ich]));
2262 hp[ih].Build(Form("H%sTrkInPt%c%d", prefix, chName[ich], isp),
2263 Form("TrackIn[%s%c]:: P_{t}[GeV/c]", spcName[v0][isp], chSgn[ich]),
2264 kEta, kPhi, kPt, aa);
2265 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2266 php.AddLast(&hp[ih++]); //np[isel]++;
2269 AliInfo(Form("Build %3d 3D %s projections.", ih, prefix));
2272 Int_t ch(0), pt(0), p(0), sp(1), rcBin(as?as->FindBin(0.):-1), ioff(0), joff(0), jsel(0), ksel(0);
2273 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
2274 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2275 v = H->GetBinContent(ib, coord);
2277 if(fBCbinTOF>0 && coord[kBC]!=fBCbinTOF) continue; // TOF bunch cross cut
2278 if(fBCbinFill>0 && abf && coord[kNdim+2]!=fBCbinTOF) continue; // Fill bunch cut
2279 if(coord[kBC]==3) continue;
2281 ch = 0; sp=1;// [pi-] track
2283 if(v0){ // V0 PID for dQdl
2284 if(!coord[kBC]) continue;
2286 } else {// TPC PID for statistics
2287 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2290 // take care of old data format (2*AliPID::kSPECIES+1)
2291 if(as->GetNbins() == kNcharge*AliPID::kSPECIES+1){
2293 else if(sp==2) sp=1;
2295 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
2296 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
2299 pt = 0; p = 0; // low pt
2300 if(apt) pt = coord[kPt];//TMath::Min(coord[kPt], Int_t(kNpt)+1);
2301 if(ap ) p = coord[kNdim];//TMath::Min(coord[kNdim], Int_t(kNpt)+1);
2303 isel = pt*(kNspc*kNcharge+1); isel+=(ch==2?(kNspc*kNcharge):(sp*kNcharge+ch)); ioff = isel*(ax?4:3);
2304 jsel = p *(kNspc*kNcharge+1); jsel+=(ch==2?(kNspc*kNcharge):(sp*kNcharge+ch)); joff = jsel*(ax?4:3);
2305 ksel = ch==2?0:nPt*(kNspc*kNcharge+1); ksel *= (ax?4:3); ksel+=sp*kNcharge+ch;
2307 AliError(Form("Wrong selection %d [%3d]", ioff, ih));
2310 if(!(pr0=(AliTRDrecoProjection*)php.At(ioff))) {
2311 AliError(Form("Missing projection @ %d", ioff));
2315 if(strcmp(pr0->H()->GetName(), Form("H%sTrkInY%c%c%d", prefix, chName[ch], ptName[pt], sp))!=0){
2316 AliError(Form("Projection mismatch :: request[H%sTrkInY%c%c%d] found[%s]", prefix, chName[ch], ptName[pt], sp, pr0->H()->GetName()));
2320 if(strcmp(pr0->H()->GetName(), Form("H%sTrkInRCZ%c", prefix, ptName[pt]))!=0){
2321 AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c] found[%s]", prefix, ptName[pt], pr0->H()->GetName()));
2325 AliDebug(2, Form("Found %s for selection sp[%d] ch[%d] pt[%d]", pr0->H()->GetName(), sp, ch, pt));
2326 for(Int_t jh(0); jh<np[isel]; jh++){
2327 if(ch<2 && jh==2) ((AliTRDrecoProjection*)php.At(joff+jh))->Increment(coord, v); // special care for dQdl=f(p)
2328 else ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v); // all = f(p_t)
2331 if(!(pr0 = (AliTRDrecoProjection*)php.At(ksel))){
2332 AliError(Form("Missing P_t projection @ %d", ksel));
2335 AliDebug(2, Form("Found %s for selection[%d] sp[%d] ch[%d]", pr0->H()->GetName(), ksel, sp, ch));
2336 pr0->Increment(coord, v); // p_t spectra
2339 if(HasDump3DFor(kTrackIn)){
2340 TDirectory *cwd = gDirectory;
2341 TFile::Open(Form("Dump%s_%s.root", GetName(), H->GetName()), "UPDATE");
2342 for(Int_t ip(0); ip<php.GetEntriesFast(); ip++){
2343 if(!(pr0 = (AliTRDrecoProjection*)php.At(ip))) continue;
2344 if(!pr0->H()) continue;
2345 TH3 *h3=(TH3*)pr0->H()->Clone();
2352 TObjArray *arr(NULL);
2353 fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx);
2354 arr->SetName(mc?projName[1]:(v0?projName[2]:projName[0])); arr->SetOwner();
2356 TH2 *h2(NULL); Int_t jh(0);
2358 if(!hp[ih].H()) continue;
2359 if(strstr(hp[ih].H()->GetName(), "TrkInPt")){
2360 for(Int_t ipt(0); ipt<nPt; ipt++) arr->AddAt(hp[ih].Projection2Dbin(ipt), jh++);
2362 if((h2 = hp[ih].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
2365 // build combined performance plots
2366 // combine up the tree of projections
2367 Double_t m(0.), s(0.), se(0.), trend(0.);
2368 AliTRDrecoProjection xlow[2], specY[kNcharge*kNspc], specPh[kNcharge*kNspc], specQ[kNcharge*kNspc];
2369 for(Int_t ich(0); ich<kNcharge; ich++){
2370 // PID dependency - summation over pt
2371 for(Int_t isp(0); isp<kNspc; isp++){
2373 Int_t idx(ich*kNspc+isp);
2374 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[0], isp)))){
2375 specY[idx] = (*pr0);
2376 specY[idx].SetNameTitle(Form("H%sTrkInY%c%d", prefix, chName[ich], isp), "Sum over pt");
2377 specY[idx].H()->SetNameTitle(Form("H%sTrkInY%c%d", prefix, chName[ich], isp),
2378 Form("TrackIn[%s%c]:: #Deltay", spcName[v0][isp], chSgn[ich]));
2379 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2380 PutTrendValue(Form("%sTrkInY%c%c%d", prefix, chName[ich], ptName[0], isp), trend, m);
2381 PutTrendValue(Form("%sTrkInYS%c%c%d", prefix, chName[ich], ptName[0], isp), s, se);
2383 for(Int_t ipt(1); ipt<nPt; ipt++){
2384 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[ipt], isp)))) continue;
2385 if((trend=pr1->GetTrendValue(1,&m,&s, &se))>-100.){
2386 PutTrendValue(Form("%sTrkInY%c%c%d", prefix, chName[ich], ptName[ipt], isp), trend, m);
2387 PutTrendValue(Form("%sTrkInYS%c%c%d", prefix, chName[ich], ptName[ipt], isp), s, se);
2391 php.AddLast(&specY[idx]);
2392 if((h2 = specY[idx].Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2393 if((trend=specY[idx].GetTrendValue(1,&m,&s,&se))>-100.){
2394 PutTrendValue(Form("%sTrkInY%c%d", prefix, chName[ich], isp), trend, m);
2395 PutTrendValue(Form("%sTrkInYS%c%d", prefix, chName[ich], isp), s,se);
2397 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", specY[idx].H()->GetName())))) arr->AddAt(h2, jh++);
2398 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%d", prefix, chName[0], isp)))) (*pr1)+=specY[idx];
2401 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[ich], ptName[0], isp)))){
2402 specPh[idx] = (*pr0);
2403 specPh[idx].SetNameTitle(Form("H%sTrkInPh%c%d", prefix, chName[ich], isp), "Sum over pt");
2404 specPh[idx].H()->SetNameTitle(Form("H%sTrkInPh%c%d", prefix, chName[ich], isp),
2405 Form("TrackIn[%s%c]:: #Delta#phi", spcName[v0][isp], chSgn[ich]));
2406 specPh[idx].SetShowRange(-1.5, 1.5);
2407 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2408 PutTrendValue(Form("%sTrkInPh%c%c%d", prefix, chName[ich], ptName[0], isp), trend, m);
2409 PutTrendValue(Form("%sTrkInPhS%c%c%d", prefix, chName[ich], ptName[0], isp), s, se);
2411 for(Int_t ipt(1); ipt<nPt; ipt++){
2412 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[ich], ptName[ipt], isp)))) continue;
2413 if((trend=pr1->GetTrendValue(1,&m,&s,&se))>-100.){
2414 PutTrendValue(Form("%sTrkInPh%c%c%d", prefix, chName[ich], ptName[ipt], isp), trend, m);
2415 PutTrendValue(Form("%sTrkInPhS%c%c%d", prefix, chName[ich], ptName[ipt], isp), s, se);
2417 specPh[idx]+=(*pr1);
2419 php.AddLast(&specPh[idx]);
2420 if((h2 = specPh[idx].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2421 if((trend=specPh[idx].GetTrendValue(1,&m,&s,&se))>-100.){
2422 PutTrendValue(Form("%sTrkInPh%c%d", prefix, chName[ich], isp), trend, m);
2423 PutTrendValue(Form("%sTrkInPhS%c%d", prefix, chName[ich], isp), s,se);
2425 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%d", prefix, chName[0], isp)))) (*pr1)+=specPh[idx];
2428 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", prefix, chName[ich], ptName[0], isp)))){
2429 specQ[idx] = (*pr0);
2430 specQ[idx].SetNameTitle(Form("H%sTrkInQ%c%d", prefix, chName[ich], isp), "Sum over p");
2431 specQ[idx].H()->SetNameTitle(Form("H%sTrkInQ%c%d", prefix, chName[ich], isp),
2432 Form("TrackIn[%s%c]:: dQdl", spcName[v0][isp], chSgn[ich]));
2433 specQ[idx].SetShowRange(-2.2, -1.75);
2434 specQ[idx].H()->GetZaxis()->SetTitle("dQdl [a.u.]");
2435 if((trend = pr0->GetTrendValue(2, &m))>-100.) PutTrendValue(Form("%sTrkInQ%c%c%d", prefix, chName[ich], ptName[0], isp), trend, m);
2436 if((trend = pr0->GetTrendValue(0, &m))>-100.) PutTrendValue(Form("%sTrkInQS%c%c%d", prefix, chName[ich], ptName[0], isp), trend, m);
2437 for(Int_t ipt(1); ipt<nPt; ipt++){
2438 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", prefix, chName[ich], ptName[ipt], isp)))) continue;
2439 if((trend=pr1->GetTrendValue(2, &m))>-100.) PutTrendValue(Form("%sTrkInQ%c%c%d", prefix, chName[ich], ptName[ipt], isp), trend, m);
2440 if((trend=pr1->GetTrendValue(0, &m))>-100.) PutTrendValue(Form("%sTrkInQS%c%c%d", prefix, chName[ich], ptName[ipt], isp), trend, m);
2443 php.AddLast(&specQ[idx]);
2444 if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2445 specQ[idx].H()->SetName(Form("H%sTrkInQS%c%d", prefix, chName[ich], isp));
2446 specQ[idx].SetShowRange(-1.85, -1.4);
2447 if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2448 if((trend=specQ[idx].GetTrendValue(2, &m))>-100.) PutTrendValue(Form("%sTrkInQ%c%d", prefix, chName[ich], isp), trend, m);
2449 if((trend=specQ[idx].GetTrendValue(0, &m))>-100.) PutTrendValue(Form("%sTrkInQS%c%d", prefix, chName[ich], isp), trend, m);
2450 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%d", prefix, chName[0], isp)))) (*pr1)+=specQ[idx];
2452 } // end PID loop for pt integration
2454 // pt dependency - summation over PID
2455 for(Int_t ipt(0); ipt<nPt; ipt++){
2457 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[ipt], 0)))){
2458 for(Int_t isp(1); isp<kNspc; isp++){
2459 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[ipt], isp)))) continue;
2462 pr0->H()->SetNameTitle(Form("H%sTrkInY%c%c", prefix, chName[ich], ptName[ipt]),
2463 Form("TrackIn[%c]:: #Deltay{%s}", chSgn[ich], ptCut[ipt]));
2464 pr0->SetShowRange(-0.3, 0.3);
2465 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2466 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2467 PutTrendValue(Form("%sTrkInY%c%c", prefix, chName[ich], ptName[ipt]), trend, m);
2468 PutTrendValue(Form("%sTrkInYS%c%c", prefix, chName[ich], ptName[ipt]), s, se);
2470 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2473 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[ich], ptName[ipt], 0)))){
2474 for(Int_t isp(1); isp<kNspc; isp++){
2475 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[ich], ptName[ipt], isp)))) continue;
2478 pr0->H()->SetNameTitle(Form("H%sTrkInPh%c%c", prefix, chName[ich], ptName[ipt]),
2479 Form("TrackIn[%c]:: #Delta#phi{%s}", chSgn[ich], ptCut[ipt]));
2480 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2481 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2482 PutTrendValue(Form("%sTrkInPh%c%c", prefix, chName[ich], ptName[ipt]), trend, m);
2483 PutTrendValue(Form("%sTrkInPhS%c%c", prefix, chName[ich], ptName[ipt]), s, se);
2485 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2488 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[ich], ptName[ipt], 0)))){
2489 for(Int_t isp(1); isp<kNspc; isp++){
2490 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[ich], ptName[ipt], isp)))) continue;
2493 pr0->H()->SetNameTitle(Form("H%sTrkInX%c%c", prefix, chName[ich], ptName[ipt]),
2494 Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[ipt]));
2495 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2496 PutTrendValue(Form("%sTrkInX%c%c", prefix, chName[ich], ptName[ipt]), pr0->GetTrendValue(1));
2499 xlow[ich].SetNameTitle(Form("H%sTrkInX%c%c%d", prefix, chName[ich], ptName[0], 5),
2500 Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[0]));
2501 php.AddLast(&xlow[ich]);
2503 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2505 } // end pt loop for PID integration
2508 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[0], 0)))){
2509 pr0->H()->SetNameTitle(Form("H%sTrkInY%c", prefix, chName[ich]),
2510 Form("TrackIn[%c]:: #Deltay", chSgn[ich]));
2511 pr0->SetShowRange(-0.3, 0.3);
2512 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2513 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2514 if((trend=pr0->GetTrendValue(1, &m, &s, &se))>-100.){
2515 PutTrendValue(Form("%sTrkInY%c", prefix, chName[ich]), trend, m);
2516 PutTrendValue(Form("%sTrkInYS%c", prefix, chName[ich]), s, se);
2518 if(ich && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2521 if(ich && (pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[0], ptName[3], 0)))){
2522 if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[ich], ptName[3], 0)))){
2524 pr0->H()->SetNameTitle(Form("H%sTrkInY%c", prefix, ptName[3]), Form("TrackIn :: #Deltay{%s}", ptCut[3]));
2525 pr0->SetShowRange(-0.3, 0.3);
2526 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2530 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[ich], ptName[0], 0)))){
2531 pr0->H()->SetNameTitle(Form("H%sTrkInPh%c", prefix, chName[ich]),
2532 Form("TrackIn[%c]:: #Delta#phi", chSgn[ich]));
2533 pr0->SetShowRange(-1., 1.);
2534 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2535 if((trend=pr0->GetTrendValue(1, &m, &s, &se))>-100.){
2536 PutTrendValue(Form("%sTrkInPh%c", prefix, chName[ich]), trend, m);
2537 PutTrendValue(Form("%sTrkInPhS%c", prefix, chName[ich]), s, se);
2539 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2542 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[ich], ptName[0], 0)))){
2543 pr0->H()->SetNameTitle(Form("H%sTrkInX%c", prefix, chName[ich]),
2544 Form("TrackIn[%c]:: #Deltax", chSgn[ich]));
2545 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2546 PutTrendValue(Form("%sTrkInX%c", prefix, chName[ich]), pr0->GetTrendValue(1));
2547 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2550 if(ich && (pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[0], ptName[1], 5)))){
2551 if((pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[ich], ptName[1], 5)))){
2553 pr0->H()->SetNameTitle(Form("H%sTrkInX%c", prefix, ptName[1]), Form("TrackIn :: #Deltax{%s}", ptCut[1]));
2554 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2557 } // end charge loop
2559 for(Int_t isp(0); isp<kNspc; isp++){
2561 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%d", prefix, chName[0], isp)))){
2562 pr0->H()->SetNameTitle(Form("H%sTrkInY%d", prefix, isp), Form("TrackIn[%s] :: #Deltay", spcName[v0][isp]));
2563 pr0->SetShowRange(-0.3, 0.3);
2564 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2565 PutTrendValue(Form("%sTrkInY%d", prefix, isp), pr0->GetTrendValue(1));
2566 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2569 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%d", prefix, chName[0], isp)))){
2570 pr0->H()->SetNameTitle(Form("H%sTrkInPh%d", prefix, isp), Form("TrackIn[%s] :: #Delta#phi", spcName[v0][isp]));
2571 pr0->SetShowRange(-1., 1.);
2572 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2573 PutTrendValue(Form("%sTrkInPh%d", prefix, isp), pr0->GetTrendValue(1));
2576 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInQ%c%d", prefix, chName[0], isp)))){
2577 pr0->H()->SetNameTitle(Form("H%sTrkInQ%d", prefix, isp), Form("TrackIn[%s] :: dQdl", spcName[v0][isp]));
2578 pr0->SetShowRange(-2.2, -1.75);
2579 if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2580 pr0->H()->SetName(Form("H%sTrkInQS%d", prefix, isp));
2581 pr0->SetShowRange(-1.85, -1.4);
2582 if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2583 if((trend=pr0->GetTrendValue(2, &m))>-100.) PutTrendValue(Form("%sTrkInQ%d", prefix, isp), trend, m);
2584 if((trend=pr0->GetTrendValue(0, &m))>-100.) PutTrendValue(Form("%sTrkInQS%d", prefix, isp), trend, m);
2586 } // end PID processing
2589 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", prefix, chName[0], ptName[0], 0)))){
2590 pr0->H()->SetNameTitle(Form("H%sTrkInY", prefix), "TrackIn :: #Deltay");
2591 pr0->SetShowRange(-0.3, 0.3);
2592 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2593 if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->H()->GetName())))) arr->AddAt(h2, jh++);
2596 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", prefix, chName[0], ptName[0], 0)))){
2597 pr0->H()->SetNameTitle(Form("H%sTrkInPh", prefix), "TrackIn :: #Delta#phi");
2598 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2601 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", prefix, chName[0], ptName[0], 0)))){
2602 pr0->H()->SetNameTitle(Form("H%sTrkInX", prefix), "TrackIn :: #Deltax");
2603 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2606 // Row Cross processing
2608 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCZ%c", prefix, ptName[0])))){
2609 for(Int_t ipt(0); ipt<nPt; ipt++){
2610 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCZ%c", prefix, ptName[ipt])))) continue;
2613 pr0->H()->SetNameTitle(Form("H%sTrkInRCZ", prefix), "TrackIn[RC]:: #Deltaz");
2614 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2615 if((trend=pr0->GetTrendValue(1,&m,&s,&se))>-100.){
2616 PutTrendValue(Form("%sTrkInRCZ", prefix), trend, m);
2617 PutTrendValue(Form("%sTrkInRCZS", prefix), s, se);
2621 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCY%c", prefix, ptName[0])))){
2622 for(Int_t ipt(0); ipt<nPt; ipt++){
2623 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCY%c", prefix, ptName[ipt])))) continue;
2626 pr0->H()->SetNameTitle(Form("H%sTrkInRCY", prefix), "TrackIn[RC]:: #Deltay");
2627 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2630 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCPh%c", prefix, ptName[0])))){
2631 for(Int_t ipt(0); ipt<nPt; ipt++){
2632 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCPh%c", prefix, ptName[ipt])))) continue;
2635 pr0->H()->SetNameTitle(Form("H%sTrkInRCPh", prefix), "TrackIn[RC]:: #Delta#phi");
2636 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2639 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCX%c", prefix, ptName[0])))){
2640 for(Int_t ipt(0); ipt<nPt; ipt++){
2641 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("H%sTrkInRCX%c", prefix, ptName[ipt])))) continue;
2644 pr0->H()->SetNameTitle(Form("H%sTrkInRCX", prefix), "TrackIn[RC]:: #Deltax");
2645 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2649 for(Int_t isp(0); isp<kNspc; isp++){
2650 for(Int_t ich(0); ich<kNcharge; ich++){
2651 TH2 *h2pt[26]={NULL}, *h2s(NULL);
2652 for(Int_t ipt(0); ipt<nPt; ipt++){
2653 if(!(h2pt[ipt] = (TH2*)arr->FindObject(Form("H%sTrkInPt%c%d%d_2D", prefix, chName[ich], isp, ipt)))){
2654 AliWarning(Form("Missing \"H%sTrkInPt%c%d%d_2D\"", prefix, chName[ich], isp, ipt));
2657 if(!h2s) h2s = (TH2F*)h2pt[ipt]->Clone("h2s");
2658 else h2s->Add(h2pt[ipt]);
2661 Int_t irebin(0), dxBin(1), dyBin(1);
2662 while(irebin<nEtaPhi && (AliTRDrecoTask::GetMeanStat(h2s, .5, 1)<200)){
2663 h2s->Rebin2D(rebinEtaPhiX[irebin], rebinEtaPhiY[irebin]);
2664 dxBin*=rebinEtaPhiX[irebin];dyBin*=rebinEtaPhiY[irebin];irebin++;
2666 AliDebug(2, Form("Rebin level[%d] @ chg[%d] spc[%d]. Binning[%2dx%2d]", irebin, ich, isp, h2s->GetNbinsX(), h2s->GetNbinsY()));
2667 Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
2668 for(Int_t ipt(0); ipt<nPt; ipt++) h2pt[ipt]->Rebin2D(dxBin, dyBin);
2671 h2 = new TH2F(Form("H%sTrkInPt%c%d_2D", prefix, chName[ich], isp),
2672 Form("TrackIn[%s%c]:: <P_{t}>[GeV/c];%s;%s", spcName[v0][isp], chSgn[ich], aa[2]->GetTitle(), aa[1]->GetTitle()),
2673 nx, aa[2]->GetXmin(), aa[2]->GetXmax(),
2674 ny, aa[1]->GetXmin(), aa[1]->GetXmax());
2675 arr->AddAt(h2, jh++);
2676 for(Int_t ix(1); ix<=nx; ix++){ // eta
2677 for(Int_t iy(1); iy<=ny; iy++){ // phi
2678 Float_t w[26]={0.}, sw(0.);
2679 for(Int_t ipt(0); ipt<nPt; ipt++){
2680 if(!h2pt[ipt]) continue;
2681 w[ipt] = h2pt[ipt]->GetBinContent(ix, iy);
2684 if(sw<=0.) h2->SetBinContent(ix, iy, 0.);
2687 for(Int_t ipt(0); ipt<nPt; ipt++){
2689 ptm += w[ipt]*fgPt[ipt?ipt-1:0];
2691 h2->SetBinContent(ix, iy, ptm);
2695 PutTrendValue(Form("%sTrkInPt%c%d", prefix, chName[ich], isp), GetMeanStat(h2, 0.01, 1));
2698 AliInfo(Form("Done %3d 2D %s projections.", jh, prefix));
2700 // clean local memory allocation
2701 for(Int_t i(fgNPt+2);i--;){
2702 if(ptCut[i]) delete [] ptCut[i];
2703 if(pCut[i]) delete [] pCut[i];
2709 //________________________________________________________
2710 Bool_t AliTRDresolution::MakeProjectionTrack()
2713 const Int_t kNcontours(9);
2714 const Int_t kNstat(30);
2715 Int_t cidx(kMCtrack);
2716 if(fProj && fProj->At(cidx)) return kTRUE;
2718 AliError("Missing data container.");
2722 if(!(H = (THnSparse*)fContainer->FindObject("hTRD2MC"))){
2723 AliError("Missing/Wrong data @ hTRD2MC.");
2726 Int_t ndim(H->GetNdimensions());
2727 Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
2728 TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
2729 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2730 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
2731 if(ndim > kPt) ap = H->GetAxis(kPt);
2732 const Int_t nPt(ap?(ap->GetNbins()+2):1);
2734 // build list of projections
2735 const Int_t nsel(AliTRDgeometry::kNlayer*fgNPt*AliPID::kSPECIES*7);//, npsel(3);
2736 const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2737 const Char_t ptShortName[5] = {'L', 'l', 'i', 'h', 'H'};
2738 Char_t ptName[fgNPt+2] = {0};
2739 Char_t *ptCut[fgNPt+2] = {NULL};
2740 // define rebinning strategy
2741 const Int_t nEtaPhi(5); Int_t rebinEtaPhiX[nEtaPhi] = {1, 3, 1, 3, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 2, 1, 2};
2742 if(aa[1]->GetNbins()==180){ // old binning
2743 rebinEtaPhiX[0] = 1; rebinEtaPhiY[0] = 2;
2744 rebinEtaPhiX[1] = 2; rebinEtaPhiY[1] = 1;
2745 rebinEtaPhiX[2] = 5; rebinEtaPhiY[2] = 1;
2746 rebinEtaPhiX[3] = 1; rebinEtaPhiY[3] = 5;
2747 rebinEtaPhiX[4] = 1; rebinEtaPhiY[4] = 1; // dummy
2749 AliTRDrecoProjection hp[kTrkNproj]; TObjArray php(kTrkNproj);
2750 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2751 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2752 for(Int_t ipt(0); ipt<nPt; ipt++){
2754 ptName[ipt]= nPt>5?(64+ipt):ptShortName[ipt];
2755 ptCut[ipt] = StrDup(Form("#it{%4.2f<=p_{t}^{MC}[GeV/c]<%4.2f}",ipt?fgPt[ipt-1]:0., ipt>fgNPt?99.99:fgPt[ipt]));
2757 for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2758 for(Int_t ich(0); ich<kNcharge; ich++){
2759 isel++; // new selection
2760 hp[ih].Build(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily),
2761 Form("Tracks[%s%c]:: #Deltay{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily),
2762 kEta, kPhi, kYrez, aa);
2763 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2764 php.AddLast(&hp[ih++]); np[isel]++;
2765 hp[ih].Build(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily),
2766 Form("Tracks[%s%c]:: #Delta#phi{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily),
2767 kEta, kPhi, kPrez, aa);
2768 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2769 php.AddLast(&hp[ih++]); np[isel]++;
2770 hp[ih].Build(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily),
2771 Form("Tracks[%s%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily),
2772 kEta, kPhi, kNdim, aa);
2773 hp[ih].SetShowRange(0.,10.);
2774 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2775 php.AddLast(&hp[ih++]); np[isel]++;
2778 isel++; // new selection
2779 hp[ih].Build(Form("HMCTrkZ%c%d", ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2780 Form("Tracks[RC]:: #Deltaz{%s} Ly[%d]", ptCut[ipt], UseLYselectTrklt()?fLYselect:ily),
2781 kEta, kPhi, kZrez, aa);
2782 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2783 php.AddLast(&hp[ih++]); np[isel]++;
2786 AliInfo(Form("Build %3d 3D projections.", ih));
2788 Int_t ly(0), ch(0), pt(0), sp(2), rcBin(as?as->FindBin(0.):-1);
2789 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2790 v = H->GetBinContent(ib, coord);
2792 ly = coord[kBC]-1; // layer selection
2793 if(ly<0) ly=0; // fix for bug in PlotMC
2795 ch=0; sp=2;// [pi-] track [dafault]
2796 if(rcBin>0){ // debug mode in which species are also saved
2797 sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2798 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
2799 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
2802 pt = 0; // low pt [default]
2803 if(ap) pt = coord[kPt]-1;
2805 Int_t ioff = ly*nPt*31+pt*31; ioff+=3*(sp<0?10:(sp*kNcharge+ch));
2806 isel = ly*nPt*11+pt*11; isel+=sp<0?10:(sp*kNcharge+ch);
2807 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));
2808 for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDrecoProjection*)php.At(ioff+jh))->Increment(coord, v);
2810 TObjArray *arr(NULL);
2811 fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
2812 arr->SetName("hTRD2MC"); arr->SetOwner();
2814 TH2 *h2(NULL); Int_t jh(0);
2816 if(!hp[ih].H()) continue;
2817 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2818 arr->AddAt(h2, jh++);
2821 // combine up the tree of projections
2822 AliTRDrecoProjection *pr0(NULL), *pr1(NULL);
2823 //Int_t iproj(0), jproj(0);
2824 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2825 for(Int_t ich(0); ich<kNcharge; ich++){
2826 for(Int_t ipt(0); ipt<nPt; ipt++){
2828 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily)))){
2829 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2830 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily)))) continue;
2833 AliDebug(2, Form("Rename %s to HMCTrkY%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily));
2834 pr0->H()->SetNameTitle(Form("HMCTrkY%c%c%d", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2835 Form("Tracks[%c]:: #Deltay{%s} Ly[%d]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily));
2836 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2837 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))) (*pr1)+=(*pr0);
2840 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily)))){
2841 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2842 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily)))) continue;
2845 AliDebug(2, Form("Rename %s to HMCTrkPh%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily));
2846 pr0->H()->SetNameTitle(Form("HMCTrkPh%c%c%d", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2847 Form("Tracks[%c]:: #Delta#phi{%s} Ly[%d]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily));
2848 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2849 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))) (*pr1)+=(*pr0);
2853 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], 0, UseLYselectTrklt()?fLYselect:ily)))){
2854 for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2855 if(!(pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, UseLYselectTrklt()?fLYselect:ily)))) continue;
2858 AliDebug(2, Form("Rename %s to HMCTrkDPt%c%c%d", pr0->H()->GetName(), chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily));
2859 pr0->H()->SetNameTitle(Form("HMCTrkDPt%c%c%d", chName[ich], ptName[ipt], UseLYselectTrklt()?fLYselect:ily),
2860 Form("Tracks[%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", chSgn[ich], ptCut[ipt], UseLYselectTrklt()?fLYselect:ily));
2861 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2862 if(ipt && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))) (*pr1)+=(*pr0);
2866 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))){
2867 pr0->H()->SetNameTitle(Form("HMCTrkY%c%d", chName[ich], UseLYselectTrklt()?fLYselect:ily),
2868 Form("Tracks[%c]:: #Deltay Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2869 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2870 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))) (*pr1)+=(*pr0);
2873 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))){
2874 pr0->H()->SetNameTitle(Form("HMCTrkPh%c%d", chName[ich], UseLYselectTrklt()?fLYselect:ily),
2875 Form("Tracks[%c]:: #Delta#phi Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2876 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2877 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))) (*pr1)+=(*pr0);
2880 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))){
2881 pr0->H()->SetNameTitle(Form("HMCTrkDPt%c%d", chName[ich], UseLYselectTrklt()?fLYselect:ily),
2882 Form("Tracks[%c]:: #Deltap_{t}/p_{t} Ly[%d]", chSgn[ich], UseLYselectTrklt()?fLYselect:ily));
2883 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2884 if(ich==1 && (pr1 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))) (*pr1)+=(*pr0);
2888 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))){
2889 pr0->H()->SetNameTitle(Form("HMCTrkY%d", UseLYselectTrklt()?fLYselect:ily), Form("Tracks :: #Deltay Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2890 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2893 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))){
2894 pr0->H()->SetNameTitle(Form("HMCTrkPh%d", UseLYselectTrklt()?fLYselect:ily), Form("Tracks :: #Delta#phi Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2895 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2898 if((pr0 = (AliTRDrecoProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, UseLYselectTrklt()?fLYselect:ily)))){
2899 pr0->H()->SetNameTitle(Form("HMCTrkDPt%d", UseLYselectTrklt()?fLYselect:ily), Form("Tracks :: #Deltap_{t}/p_{t} Ly[%d]", UseLYselectTrklt()?fLYselect:ily));
2900 if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2903 AliInfo(Form("Done %3d 2D projections.", jh));
2907 //________________________________________________________
2908 Bool_t AliTRDresolution::PostProcess()
2910 // Fit, Project, Combine, Extract values from the containers filled during execution
2913 AliError("ERROR: list not available");
2917 AliInfo("Building array of projections ...");
2918 fProj = new TObjArray(kNclasses+1); fProj->SetOwner(kTRUE);
2920 // set pt/p segmentation. guess from data
2922 if(!(H = (THnSparse*)fContainer->FindObject("hTracklet2TRDin"))){
2923 AliWarning("Missing/Wrong data @ hTracklet2TRDin needed for infering pt/p segmentation.");
2926 // protect for backward compatibility
2927 fNpt=H->GetAxis(kPt)?(H->GetAxis(kPt)->GetNbins()+1):1;
2928 if(!MakeMomSegmentation()) return kFALSE;
2930 //PROCESS EXPERIMENTAL DISTRIBUTIONS
2931 // Clusters detector
2932 if(HasProcess(kDetector)) MakeProjectionDetector();
2933 // Clusters residuals
2934 if(HasProcess(kCluster)) MakeProjectionCluster();
2936 // Tracklet residual/pulls
2937 if(HasProcess(kTracklet)) MakeProjectionTracklet();
2939 // TRDin residual/pulls
2940 if(HasProcess(kTrackIn)){
2941 MakeProjectionTrackIn();
2942 MakeProjectionTrackIn(kFALSE, kTRUE);
2946 if(!HasMCdata()) return kTRUE;
2947 //PROCESS MC RESIDUAL DISTRIBUTIONS
2949 // CLUSTER Y RESOLUTION/PULLS
2950 if(HasProcess(kCluster)) MakeProjectionCluster(kTRUE);
2953 // TRACKLET RESOLUTION/PULLS
2954 if(HasProcess(kTracklet)) if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
2957 // TRACK RESOLUTION/PULLS
2958 // if(HasProcess(kTracklet)) if(!MakeProjectionTrack()) return kFALSE;
2959 // fNRefFigures+=16;
2961 // TRACK TRDin RESOLUTION/PULLS
2962 if(HasProcess(kTrackIn)) if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
2969 //________________________________________________________
2970 void AliTRDresolution::Terminate(Option_t *opt)
2972 AliTRDrecoTask::Terminate(opt);
2973 if(HasPostProcess()) PostProcess();
2976 //________________________________________________________
2977 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2979 // Helper function to avoid duplication of code
2980 // Make first guesses on the fit parameters
2982 // find the intial parameters of the fit !! (thanks George)
2983 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2985 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2986 f->SetParLimits(0, 0., 3.*sum);
2987 f->SetParameter(0, .9*sum);
2988 Double_t rms = h->GetRMS();
2989 f->SetParLimits(1, -rms, rms);
2990 f->SetParameter(1, h->GetMean());
2992 f->SetParLimits(2, 0., 2.*rms);
2993 f->SetParameter(2, rms);
2994 if(f->GetNpar() <= 4) return;
2996 f->SetParLimits(3, 0., sum);
2997 f->SetParameter(3, .1*sum);
2999 f->SetParLimits(4, -.3, .3);
3000 f->SetParameter(4, 0.);
3002 f->SetParLimits(5, 0., 1.e2);
3003 f->SetParameter(5, 2.e-1);
3006 //________________________________________________________
3007 TObjArray* AliTRDresolution::Histos()
3010 // Define histograms
3013 if(fContainer) return fContainer;
3015 fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
3017 const Int_t nhn(100); Char_t hn[nhn]; TString st;
3018 Int_t triggerBins(fTriggerList?(TMath::Power(2.,fTriggerList->GetEntriesFast())-1):0);
3020 //++++++++++++++++++++++
3021 // cluster to detector
3022 snprintf(hn, nhn, "h%s", fgPerformanceName[kDetector]);
3023 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3024 const Int_t mdim(6);
3025 const Char_t *cldTitle[mdim] = {"layer", fgkTitle[kPhi], "pad row", "centrality", "q [a.u.]", triggerBins?"trigger":"no. of pads"/*, "tb [10 Hz]"*/};
3026 const Int_t cldNbins[mdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], 76, AliTRDeventInfo::kCentralityClasses, 50, triggerBins?triggerBins:kNpads};
3027 const Double_t cldMin[mdim] = {-0.5, fgkMin[kPhi], -0.5, -0.5, 0., 0.5},
3028 cldMax[mdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], 75.5, AliTRDeventInfo::kCentralityClasses - 0.5, 1200., triggerBins?(triggerBins+.5):(kNpads+.5)};
3029 st = "cluster proprieties;";
3030 // define minimum info to be saved in non debug mode
3031 Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimDet);
3032 for(Int_t idim(0); idim<ndim; idim++){ st += cldTitle[idim]; st+=";";}
3033 H = new THnSparseI(hn, st.Data(), ndim, cldNbins, cldMin, cldMax);
3035 fContainer->AddAt(H, kDetector);
3036 //++++++++++++++++++++++
3037 // cluster to tracklet residuals/pulls
3038 snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
3039 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3040 const Int_t mdim(10);
3041 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"};
3042 const Int_t clNbins[mdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 30, 30, 15, AliTRDeventInfo::kCentralityClasses, triggerBins?triggerBins:kNpads};
3043 const Double_t clMin[mdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45, -0.5, 0.5},
3044 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)};
3045 st = "cluster spatial&charge resolution;";
3046 // define minimum info to be saved in non debug mode
3047 Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimCl);
3048 for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
3049 H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
3051 fContainer->AddAt(H, kCluster);
3052 //++++++++++++++++++++++
3053 // tracklet to TRD track
3054 snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
3055 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3056 const Int_t mdim(kNdim+8);
3057 Char_t *trTitle[mdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
3058 Int_t trNbins[mdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
3059 Double_t trMin[mdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
3060 Double_t trMax[mdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
3061 // set specific fields
3062 trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
3063 trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
3064 trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
3065 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
3066 trNbins[kPt]=fNpt-1; trMax[kPt] = trNbins[kPt]-.5;
3068 trTitle[jdim]=StrDup("centrality"); trNbins[jdim] = AliTRDeventInfo::kCentralityClasses; trMin[jdim] = -.5; trMax[jdim] = AliTRDeventInfo::kCentralityClasses - 0.5;
3071 trTitle[jdim]=StrDup("trigger"); trNbins[jdim] = triggerBins;
3072 trMin[jdim] = 0.5; trMax[jdim] = triggerBins+.5;
3074 trTitle[jdim]=StrDup("occupancy [%]"); trNbins[jdim] = 12;
3075 trMin[jdim] = 25.; trMax[jdim] = 85.;
3078 st = "tracklet spatial&charge resolution;";
3079 // define minimum info to be saved in non debug mode
3080 Int_t ndim=DebugLevel()>=1?(jdim+1):kNdimTrklt;
3081 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
3082 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
3084 fContainer->AddAt(H, kTracklet);
3085 //++++++++++++++++++++++
3086 // tracklet to TRDin
3087 snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
3088 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3089 // set specific fields
3090 const Int_t mdim(kNdim+2);
3091 Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*));
3092 Int_t trinNbins[mdim]; memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t));
3093 Double_t trinMin[mdim]; memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t));
3094 Double_t trinMax[mdim]; memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t));
3096 trinTitle[kBC]=StrDup("trigger"); trinNbins[kBC] = triggerBins;
3097 trinMin[kBC] = 0.5; trinMax[kBC] = triggerBins+.5;
3099 // trinNbins[kSpeciesChgRC] = Int_t(kNcharge)*(kNspc-1)+1; trinMin[kSpeciesChgRC] = -kNspc+0.5; trinMax[kSpeciesChgRC] = kNspc-0.5;
3100 trinNbins[kPt]=fNpt-1; trinMax[kPt] = trinNbins[kPt]-.5;
3101 trinTitle[kNdim]=StrDup("p [bin]"); trinNbins[kNdim] = fNpt-1; trinMin[kNdim] = -0.5; trinMax[kNdim] = trinNbins[kNdim]-.5;
3102 trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4;
3103 //trinTitle[kNdim+2]=StrDup("Fill Bunch"); trinNbins[kNdim+2]=3500; trinMin[kNdim+2]=-0.5; trinMax[kNdim+2]=3499.5;
3104 st = "r-#phi/z/angular residuals @ TRD entry;";
3105 // define minimum info to be saved in non debug mode
3106 Int_t ndim=(DebugLevel()>=1?mdim:(kNdim+1));//kNdimTrkIn;
3107 for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";}
3108 H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax);
3110 fContainer->AddAt(H, kTrackIn);
3111 // tracklet to TRDout
3112 // fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
3115 // Resolution histos
3116 if(!HasMCdata()) return fContainer;
3118 //++++++++++++++++++++++
3119 // cluster to TrackRef residuals/pulls
3120 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
3121 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3122 H = (THnSparseI*)((THnSparseI*)fContainer->At(kCluster))->Clone(hn);
3123 H->SetTitle("MC cluster spatial resolution");
3125 fContainer->AddAt(H, kMCcluster);
3126 //++++++++++++++++++++++
3127 // tracklet to TrackRef
3128 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
3129 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3130 H = (THnSparseI*)((THnSparseI*)fContainer->At(kTracklet))->Clone(hn);
3131 H->SetTitle("MC tracklet spatial resolution");
3133 fContainer->AddAt(H, kMCtracklet);
3134 //++++++++++++++++++++++
3135 // TRDin to TrackRef
3136 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
3137 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3138 H = (THnSparseI*)((THnSparseI*)fContainer->At(kTrackIn))->Clone(hn);
3139 H->SetTitle("MC r-#phi/z/angular residuals @ TRD entry");
3141 fContainer->AddAt(H, kMCtrackIn);
3142 //++++++++++++++++++++++
3143 // track to TrackRef
3144 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
3145 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
3146 Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
3147 Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
3148 Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
3149 Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
3150 // set specific fields
3151 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
3152 trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
3153 trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
3154 trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
3155 trNbins[kPt] =fNpt-1;trMax[kPt] = trNbins[kPt]-.5;
3156 //trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
3157 trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5;
3159 st = "MC track spatial&p_{t} resolution;";
3160 // define minimum info to be saved in non debug mode
3161 Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
3162 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
3163 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
3165 fContainer->AddAt(H, kMCtrack);
3170 //____________________________________________________________________
3171 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
3174 // Fit track with a staight line using the "np" clusters stored in the array "points".
3175 // The following particularities are stored in the clusters from points:
3176 // 1. pad tilt as cluster charge
3177 // 2. pad row cross or vertex constrain as fake cluster with cluster type 1
3178 // The parameters of the straight line fit are stored in the array "param" in the following order :
3179 // param[0] - x0 reference radial position
3180 // param[1] - y0 reference r-phi position @ x0
3181 // param[2] - z0 reference z position @ x0
3182 // param[3] - slope dy/dx
3183 // param[4] - slope dz/dx
3186 // Function should be used to refit tracks for B=0T
3190 if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
3193 TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
3196 for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
3199 Double_t x, y, z, dx, tilt(0.);
3200 for(Int_t ip(0); ip<np; ip++){
3201 x = points[ip].GetX(); z = points[ip].GetZ();
3203 zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
3205 if(zfitter.Eval() != 0) return kFALSE;
3207 Double_t z0 = zfitter.GetParameter(0);
3208 Double_t dzdx = zfitter.GetParameter(1);
3209 for(Int_t ip(0); ip<np; ip++){
3210 if(points[ip].GetClusterType()) continue;
3211 x = points[ip].GetX();
3213 y = points[ip].GetY();
3214 z = points[ip].GetZ();
3215 tilt = points[ip].GetCharge();
3216 y -= tilt*(-dzdx*dx + z - z0);
3217 Float_t xyz[3] = {(Float_t)x, (Float_t)y, (Float_t)z}; points[ip].SetXYZ(xyz);
3218 yfitter.AddPoint(&dx, y, 1.);
3220 if(yfitter.Eval() != 0) return kFALSE;
3221 Double_t y0 = yfitter.GetParameter(0);
3222 Double_t dydx = yfitter.GetParameter(1);
3224 param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
3225 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);
3229 //____________________________________________________________________
3230 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
3233 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
3234 // See function FitTrack for the data stored in the "clusters" array
3236 // The parameters of the straight line fit are stored in the array "param" in the following order :
3237 // par[0] - x0 reference radial position
3238 // par[1] - y0 reference r-phi position @ x0
3239 // par[2] - slope dy/dx
3242 // Function should be used to refit tracks for B=0T
3245 TLinearFitter yfitter(2, "pol1");
3247 // grep data for tracklet
3248 Double_t x0(0.), x[60], y[60], dy[60];
3250 for(Int_t ip(0); ip<np; ip++){
3251 if(points[ip].GetClusterType()) continue;
3252 if(points[ip].GetVolumeID() != ly) continue;
3253 Float_t xt(points[ip].GetX())
3254 ,yt(param[1] + param[3] * (xt - param[0]));
3256 y[nly] = points[ip].GetY();
3262 if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
3265 // set radial reference for fit
3268 // find tracklet core
3269 Double_t mean(0.), sig(1.e3);
3270 AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
3272 // simple cluster error parameterization
3273 Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
3275 // fit tracklet core
3276 for(Int_t jly(0); jly<nly; jly++){
3277 if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
3278 Double_t dx(x[jly]-x0);
3279 yfitter.AddPoint(&dx, y[jly], 1.);
3281 if(yfitter.Eval() != 0) return kFALSE;
3283 par[1] = yfitter.GetParameter(0);
3284 par[2] = yfitter.GetParameter(1);
3288 //____________________________________________________________________
3289 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
3292 // Global selection mechanism of tracksbased on cluster to fit residuals
3293 // The parameters are the same as used ni function FitTrack().
3295 const Float_t kS(0.6), kM(0.2);
3296 TH1S h("h1", "", 100, -5.*kS, 5.*kS);
3297 Float_t dy/*, dz*/, s, m;
3298 for(Int_t ip(0); ip<np; ip++){
3299 if(points[ip].GetClusterType()) continue;
3300 Float_t x0(points[ip].GetX())
3301 ,y0(param[1] + param[3] * (x0 - param[0]))
3302 /*,z0(param[2] + param[4] * (x0 - param[0]))*/;
3303 dy=points[ip].GetY() - y0; h.Fill(dy);
3304 //dz=points[ip].GetZ() - z0;
3306 TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
3307 fg.SetParameter(1, 0.);
3308 fg.SetParameter(2, 2.e-2);
3310 m=fg.GetParameter(1); s=fg.GetParameter(2);
3311 if(s>kS || TMath::Abs(m)>kM) return kFALSE;
3315 //________________________________________________________
3316 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
3319 // Get the most probable value and the full width half mean
3320 // of a Landau distribution
3323 const Float_t dx = 1.;
3324 mpv = f->GetParameter(1);
3325 Float_t fx, max = f->Eval(mpv);
3328 while((fx = f->Eval(xm))>.5*max){
3337 while((fx = f->Eval(xM))>.5*max) xM += dx;
3341 // #include "TFile.h"
3342 // //________________________________________________________
3343 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
3346 // AliWarning("Use cluster position as in reconstruction.");
3347 // SetLoadCorrection();
3350 // TDirectory *cwd(gDirectory);
3351 // TString fileList;
3352 // FILE *filePtr = fopen(file, "rt");
3354 // AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
3355 // SetLoadCorrection();
3358 // TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
3359 // while(fileList.Gets(filePtr)){
3360 // if(!TFile::Open(fileList.Data())) {
3361 // AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
3363 // } else AliInfo(Form("\"%s\"", fileList.Data()));
3365 // TTree *tSys = (TTree*)gFile->Get("tSys");
3366 // h2->SetDirectory(gDirectory); h2->Reset("ICE");
3367 // tSys->Draw("det:t>>h2", "dx", "goff");
3368 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3369 // for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
3371 // h2->SetDirectory(cwd);
3376 // if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>=2){
3377 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3378 // printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
3379 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3380 // FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
3381 // fprintf(fout, " static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
3382 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3383 // printf("%03d|", idet);
3384 // fprintf(fout, " {");
3385 // for(Int_t it(0); it<30; it++){
3386 // printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
3387 // fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
3390 // fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
3392 // fprintf(fout, " };\n");
3394 // SetLoadCorrection();
3398 //________________________________________________________
3399 void AliTRDresolution::SetProcesses(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3402 if(det) SETBIT(fSteer, kDetector); else CLRBIT(fSteer, kDetector);
3403 if(cl) SETBIT(fSteer, kCluster); else CLRBIT(fSteer, kCluster);
3404 if(trklt) SETBIT(fSteer, kTracklet); else CLRBIT(fSteer, kTracklet);
3405 if(trkin) SETBIT(fSteer, kTrackIn); else CLRBIT(fSteer, kTrackIn);
3409 //________________________________________________________
3410 void AliTRDresolution::SetDump3D(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3413 if(det) SETBIT(fSteer, 4+kDetector); else CLRBIT(fSteer, 4+kDetector);
3414 if(cl) SETBIT(fSteer, 4+kCluster); else CLRBIT(fSteer, 4+kCluster);
3415 if(trklt) SETBIT(fSteer, 4+kTracklet); else CLRBIT(fSteer, 4+kTracklet);
3416 if(trkin) SETBIT(fSteer, 4+kTrackIn); else CLRBIT(fSteer, 4+kTrackIn);