]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDcheckTRK.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- bunch crossing selection in PWG1/TRD
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckTRK.cxx
CommitLineData
3ceb45ae 1/**************************************************************************
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercialf purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15
16
17////////////////////////////////////////////////////////////////////////////
18// //
19// TRD tracker systematic //
20//
21//
22// Authors: //
23// Alexandru Bercuci <A.Bercuci@gsi.de> //
24// //
25////////////////////////////////////////////////////////////////////////////
26
27#include "TROOT.h"
28#include "TAxis.h"
29#include "TH1.h"
30#include "TH2.h"
31#include "TH3.h"
32#include "TObjArray.h"
33#include "THnSparse.h"
34#include <TVectorT.h>
35
36#include "AliLog.h"
37
38#include "AliTRDcluster.h"
39#include "AliTRDseedV1.h"
40#include "AliTRDtrackV1.h"
41#include "AliTRDtrackerV1.h"
42#include "AliTRDtransform.h"
43
44#include "AliTRDcheckTRK.h"
45
46ClassImp(AliTRDcheckTRK)
47
48Bool_t AliTRDcheckTRK::fgKalmanUpdate = kTRUE;
49Bool_t AliTRDcheckTRK::fgClRecalibrate = kFALSE;
50Float_t AliTRDcheckTRK::fgKalmanStep = 2.;
51//__________________________________________________________________________
52AliTRDcheckTRK::AliTRDcheckTRK()
53 : AliTRDrecoTask()
54{
55// Default constructor
56 SetNameTitle("TRDtrackerSys", "TRD Tracker Systematic");
57 Float_t pt0(0.2);
58 for(Int_t j(0); j<=kNpt; j++){
59 pt0+=(TMath::Exp(j*j*.002)-1.);
60 fPtBins[j]=pt0;
61 }
62 memset(fProj, 0, 10*sizeof(TH1*));
63}
64
65//__________________________________________________________________________
66AliTRDcheckTRK::AliTRDcheckTRK(char* name)
67 : AliTRDrecoTask(name, "TRD Tracker Systematic")
68{
69// User constructor
70 InitFunctorList();
71 Float_t pt0(0.2);
72 for(Int_t j(0); j<=kNpt; j++){
73 pt0+=(TMath::Exp(j*j*.002)-1.);
74 fPtBins[j]=pt0;
75 }
76 memset(fProj, 0, 10*sizeof(TH1*));
77}
78
79//__________________________________________________________________________
80AliTRDcheckTRK::~AliTRDcheckTRK()
81{
82// Destructor
83}
84
85//__________________________________________________________________________
86Int_t AliTRDcheckTRK::GetPtBin(Float_t pt)
87{
88// Find pt bin according to local pt segmentation
89 Int_t ipt(0);
90 while(ipt<kNpt){
91 if(pt<fPtBins[ipt]) break;
92 ipt++;
93 }
94 return TMath::Max(0,ipt);
95}
96
97//__________________________________________________________________________
98Int_t AliTRDcheckTRK::GetSpeciesByMass(Float_t m)
99{
100// Find particle index by mass
101// 0 electron
102// 1 muon
103// 2 pion
104// 3 kaon
105// 4 proton
106
107 for(Int_t is(0); is<AliPID::kSPECIES; is++) if(TMath::Abs(m-AliPID::ParticleMass(is))<1.e-4) return is;
108 return -1;
109}
110
111
112//__________________________________________________________________________
113Bool_t AliTRDcheckTRK::GetRefFigure(Int_t ifig)
114{
115 switch(ifig){
116 case 0:
117 if(!MakeProjectionEtaPhi()) return kFALSE;
118 break;
119 }
120 return kTRUE;
121}
122
123//__________________________________________________________________________
124TObjArray* AliTRDcheckTRK::Histos()
125{
126 //
127 // Create QA histogram tree
128 //
129
130 if(fContainer) return fContainer;
131
132 THnSparseI *h(NULL);
133 fContainer = new TObjArray(kNclasses);
134 fContainer->SetOwner(kTRUE);
135
136 const Char_t *ccn[kNclasses] = {"Entry", "Propag"};
137 const Char_t *ctt[kNclasses] = {"r-#phi/z/angular residuals @ entry", "r-#phi/z/angular residuals each ly"};
138
139 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140 Int_t bins[kNdim+1] = {Int_t(kNbunchCross)/*bc*/, 180/*phi*/, 50/*eta*/, Int_t(kNcharge)*AliPID::kSPECIES+1/*chg*species*/, kNpt/*pt*/, 50/*dy*/, 50/*dz*/, 40/*dphi*/, AliTRDgeometry::kNlayer};
141 Double_t min[kNdim+1] = {-0.5, -TMath::Pi(), -1., -AliPID::kSPECIES-0.5, -0.5, -1.5, -2.5, -10., -0.5},
142 max[kNdim+1] = {Int_t(kNbunchCross)-0.5, TMath::Pi(), 1., AliPID::kSPECIES+0.5, kNpt-0.5, 1.5, 2.5, 10., AliTRDgeometry::kNlayer-0.5};
143 Char_t hn[100], ht[700];
144 snprintf(hn, 100, "h%s", ccn[kEntry]);
145 if(!(h = (THnSparseI*)gROOT->FindObject(hn))){
146 snprintf(ht, 700, "%s;bunch cross;#phi [rad];#eta;chg*spec*rc;bin_p_{t};#Deltay [cm];#Deltaz [cm];#Delta#phi [deg];", ctt[kEntry]);
147 h = new THnSparseI(hn, ht, kNdim, bins, min, max);
148 } else h->Reset();
149 fContainer->AddAt(h, kEntry);
150
151 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
152 min[5] = -0.15; max[5] = 0.15;
153 snprintf(hn, 100, "h%s", ccn[kPropagation]);
154 if(!(h = (THnSparseI*)gROOT->FindObject(hn))){
155 snprintf(ht, 700, "%s;bunch cross;#phi [rad];#eta;chg*spec*rc;bin_p_{t};#Deltay [cm];#Deltaz [cm];#Delta#phi [deg];layer;", ctt[kPropagation]);
156 h = new THnSparseI(hn, ht, kNdim+1, bins, min, max);
157 } else h->Reset();
158 fContainer->AddAt(h, kPropagation);
159
160 return fContainer;
161}
162
163//__________________________________________________________________________
164TH1* AliTRDcheckTRK::PlotEntry(const AliTRDtrackV1 *track)
165{
166// comment needed
167 if(track) fkTrack = track;
168 if(!fkTrack){
169 AliDebug(4, "No Track defined.");
170 return NULL;
171 }
172 // check container
173 THnSparseI *h=(THnSparseI*)fContainer->At(kEntry);
174 if(!h){
175 AliError(Form("Missing container @ %d", Int_t(kEntry)));
176 return NULL;
177 }
178 // check input track status
179 AliExternalTrackParam *tin(NULL);
180 if(!(tin = fkTrack->GetTrackIn())){
181 AliError("Track did not entered TRD fiducial volume.");
182 return NULL;
183 }
184 // check first tracklet
185 AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
186 if(!fTracklet){
187 AliDebug(3, "No Tracklet in ly[0]. Skip track.");
188 return NULL;
189 }
190 // check radial position
191 Double_t x = tin->GetX();
192 if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
193 AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
194 return NULL;
195 }
196
197 Double_t xyz[3];
198 if(!tin->GetXYZ(xyz)){
199 AliDebug(1, "Failed getting global track postion");
200 xyz[0]=x; xyz[1]=0.; xyz[2]=0.;
201 }
202 Float_t mass(fkTrack->GetMass()),
203 eta(tin->Eta()),
204 phi(TMath::ATan2(xyz[1], xyz[0]));
205 Int_t charge(fkTrack->Charge()),
206 species(GetSpeciesByMass(mass)),
207 bc(fkESD->GetTOFbc());
208 const Double_t *parR(tin->GetParameter());
209 Double_t dyt(parR[0] - fTracklet->GetYfit(0)), dzt(parR[1] - fTracklet->GetZfit(0)),
210 phit(fTracklet->GetYfit(1)),
211 tilt(fTracklet->GetTilt());
212
213 // correct for tilt rotation
214 Double_t dy = dyt - dzt*tilt,
215 dz = dzt + dyt*tilt;
216 phit += tilt*parR[3];
217 Double_t dphi = TMath::ASin(parR[2])-TMath::ATan(phit);
218
219 Double_t val[kNdim];
220 val[kBC] = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
221 val[kPhi] = phi;
222 val[kEta] = eta;
223 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:charge*(species+1);
224 val[kPt] = GetPtBin(tin->Pt());
225 val[kYrez] = dy;
226 val[kZrez] = dz;
227 val[kPrez] = dphi*TMath::RadToDeg();
228 h->Fill(val);
229 return NULL;
230}
231
232//__________________________________________________________________________
233TH1* AliTRDcheckTRK::PlotPropagation(const AliTRDtrackV1 *track)
234{
235// comment needed
236
237 if(track) fkTrack = track;
238 if(!fkTrack){
239 AliDebug(4, "No Track defined.");
240 return NULL;
241 }
242 // check container
243 THnSparseI *h=(THnSparseI*)fContainer->At(kPropagation);
244 if(!h){
245 AliError(Form("Missing container @ %d", Int_t(kPropagation)));
246 return NULL;
247 }
248 // check input track status
249 AliExternalTrackParam *tin(NULL);
250 if(!(tin = fkTrack->GetTrackIn())){
251 AliError("Track did not entered TRD fiducial volume.");
252 return NULL;
253 }
254
255 Float_t mass(fkTrack->GetMass());
256 Int_t charge(fkTrack->Charge()),
257 species(GetSpeciesByMass(mass)),
258 bc(fkESD->GetTOFbc()+1);
259 TVectorD dX(AliTRDgeometry::kNlayer),
260 dY(AliTRDgeometry::kNlayer),
261 dZ(AliTRDgeometry::kNlayer),
262 dPhi(AliTRDgeometry::kNlayer),
263 vPt(AliTRDgeometry::kNlayer),
264 vPhi(AliTRDgeometry::kNlayer),
265 vEta(AliTRDgeometry::kNlayer),
266 budget(AliTRDgeometry::kNlayer),
267 cCOV(AliTRDgeometry::kNlayer*15);
268 const Int_t nopt(1000); Char_t opt[nopt];
269
270 // propagation
271 //snprintf(opt, nopt, "");
272
273 AliTRDseedV1 *tr(NULL);
274 Double_t val[kNdim+1];
275 if(PropagateKalman(fkTrack, tin, &dX, &dY, &dZ, &dPhi, &vPt, &vPhi, &vEta, &budget, &cCOV, opt)){
276 val[kBC]=(bc>=kNbunchCross)?(kNbunchCross-1):bc;
277 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
278 if(dX[ily]<0.) continue;
279 if(!(tr = fkTrack->GetTracklet(ily))) continue;
280 val[kPhi] = vPhi[ily];
281 val[kEta] = vEta[ily];
282 val[kSpeciesChgRC]= tr->IsRowCross()?0:charge*(species+1);
283 val[kPt] = GetPtBin(vPt[ily]);
284 val[kYrez] = dY[ily];
285 val[kZrez] = dZ[ily];
286 val[kPrez] = dPhi[ily]*TMath::RadToDeg();
287 val[kNdim] = ily;
288 h->Fill(val);
289 }
290 }
291
292 return NULL;
293}
294
295
296//___________________________________________________
297Bool_t AliTRDcheckTRK::PropagateKalman(const AliTRDtrackV1 *t, AliExternalTrackParam *ref, TVectorD *dx, TVectorD *dy, TVectorD *dz, TVectorD *dphi, TVectorD *pt, TVectorD *eta, TVectorD *phi, TVectorD *budget, TVectorD *cov, Option_t */*opt*/)
298{
299// Propagate Kalman from the first TRD tracklet to
300// last one and save residuals in the y, z and pt.
301// The track is intialized from the reference "ref".
302//
303// This is to calibrate the dEdx and MS corrections
304
305 Int_t ntracklets(t->GetNumberOfTracklets());
306 if(!ntracklets) return kFALSE;
307 if(ref->Pt()<1.e-3) return kFALSE;
308
309 AliTRDseedV1 *tr(NULL);
310 for(Int_t itr(AliTRDgeometry::kNlayer); itr--;){
311 (*dx)[itr] = -1.; (*dy)[itr] = 100.; (*dz)[itr] = 100.; (*dphi)[itr] = 100.;
312 }
313
314 // Initialize TRD track to the reference
315 AliTRDtrackV1 tt;
316 tt.Set(ref->GetX(), ref->GetAlpha(), ref->GetParameter(), ref->GetCovariance());
317 tt.SetMass(t->GetMass());
318
319 Double_t x0(ref->GetX()), xyz[3] = {0., 0., 0.};
320 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
321 if(!(tr = t->GetTracklet(ily))) continue;
322 if(fgClRecalibrate){
323 AliTRDtransform trans(tr->GetDetector());
324 AliTRDcluster *c(NULL);
325 tr->ResetClusterIter(kFALSE);
326 while((c = tr->PrevCluster())) trans.Recalibrate(c, kFALSE);
327 tr->FitRobust();
328 }
329 if(!AliTRDtrackerV1::PropagateToX(tt, tr->GetX0(), fgKalmanStep)) continue;
330 if(fgKalmanUpdate){
331 Double_t x(tr->GetX0()),
332 p[2] = { tr->GetYfit(0), tr->GetZfit(0)},
333 covTrklt[3];
334 tr->GetCovAt(x, covTrklt);
335 if(!((AliExternalTrackParam&)tt).Update(p, covTrklt)) continue;
336 }
337 if(!tt.GetXYZ(xyz)) continue;
338 (*phi)[ily] = TMath::ATan2(xyz[1], xyz[0]);
339 (*eta)[ily] = tt.Eta();
340 (*dx)[ily] = tt.GetX() - x0;
341 (*pt)[ily] = tt.Pt();
342 if(budget) (*budget)[ily] = tt.GetBudget(0);
343 if(cov){
344 const Double_t *cc(tt.GetCovariance());
345 for(Int_t ic(0), jp(ily*15); ic<15; ic++, jp++) (*cov)[jp]=cc[ic];
346 }
347 const Double_t *parR(tt.GetParameter());
348 Double_t dyt(parR[0] - tr->GetYfit(0)), dzt(parR[1] - tr->GetZfit(0)),
349 dydx(tr->GetYfit(1)),
350 tilt(tr->GetTilt());
351 // correct for tilt rotation
352 (*dy)[ily] = dyt - dzt*tilt;
353 (*dz)[ily] = dzt + dyt*tilt;
354 dydx += tilt*parR[3];
355 (*dphi)[ily]= TMath::ASin(parR[2])-TMath::ATan(dydx);
356 }
357 return kTRUE;
358}
359
360
361//__________________________________________________________________________
362Bool_t AliTRDcheckTRK::MakeProjectionEtaPhi()
363{
364// Get dy residual projection as a function of eta and phi
365
366 if(fProj[0] && fProj[1]) return kTRUE;
367 if(!fContainer || fContainer->GetEntries() != kNclasses){
368 AliError("Missing/Wrong data container.");
369 return kFALSE;
370 }
371 THnSparse *H(NULL);
372 if(!(H = (THnSparse*)fContainer->At(kEntry))){
373 AliError(Form("Missing/Wrong data @ %d.", Int_t(kEntry)));
374 return kFALSE;
375 }
376
377 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
378 TAxis *abc(H->GetAxis(kBC)),
379 *aphi(H->GetAxis(kPhi)),
380 *aeta(H->GetAxis(kEta)),
381 *as(H->GetAxis(kSpeciesChgRC)),
382 *apt(H->GetAxis(kPt)),
383 *ay(H->GetAxis(kYrez)),
384 *az(H->GetAxis(kZrez)),
385 *ap(H->GetAxis(kPrez));
386 Int_t neta(aeta->GetNbins()), nphi(aphi->GetNbins());
387 TH3I *h3[3];
388 h3[0] = new TH3I("h30", Form("r-#phi residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
389 neta, aeta->GetXmin(), aeta->GetXmax(),
390 nphi, aphi->GetXmin(), aphi->GetXmax(),
391 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
392 h3[1] = new TH3I("h31", Form("z residuals for row cross;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), az->GetTitle()),
393 neta, aeta->GetXmin(), aeta->GetXmax(),
394 nphi, aphi->GetXmin(), aphi->GetXmax(),
395 az->GetNbins(), az->GetXmin(), az->GetXmax());
396 h3[2] = new TH3I("h32", Form("r-#phi residuals for pos tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
397 neta, aeta->GetXmin(), aeta->GetXmax(),
398 nphi, aphi->GetXmin(), aphi->GetXmax(),
399 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
400 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
401 v = H->GetBinContent(ib, coord);
402 if(coord[kBC]>1) continue; // bunch cross cut
403 // species selection
404 if(coord[kSpeciesChgRC]<6){
405 h3[0]->AddBinContent(
406 h3[0]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
407 } else if(coord[kSpeciesChgRC]==6) {
408 h3[1]->AddBinContent(
409 h3[1]->GetBin(coord[kEta], coord[kPhi], coord[kZrez]), v);
410 } else if(coord[kSpeciesChgRC]>6) {
411 h3[2]->AddBinContent(
412 h3[2]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
413 }
414 }
415 printf("start Z projection\n");
416 TH2F *h2[3];
417 for(Int_t iproj(0); iproj<3; iproj++){
418 h2[iproj] = new TH2F(Form("h2%d", iproj),
419 Form("%s;%s;%s;%s", h3[iproj]->GetTitle(), aeta->GetTitle(), aphi->GetTitle(), h3[iproj]->GetZaxis()->GetTitle()),
420 neta, aeta->GetXmin(), aeta->GetXmax(),
421 nphi, aphi->GetXmin(), aphi->GetXmax());
422 for(Int_t iphi(1); iphi<=nphi; iphi++){
423 for(Int_t ieta(1); ieta<=neta; ieta++){
424 TH1 *h = h3[iproj]->ProjectionZ(Form("hy%d", iproj), ieta, ieta, iphi, iphi);
425 h2[iproj]->SetBinContent(ieta, iphi, h->GetEntries()>20?h->GetMean():-999.);
426 }
427 }
428 }
429
430 return kTRUE;
431}