]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/macros/AliTrackerTest.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
This test macro to test AliTracker functionality.
[u/mrichter/AliRoot.git] / PWG1 / macros / AliTrackerTest.C
CommitLineData
8296a418 1
2/*
3 This test macro to test AliTracker functionality.
4 Additional focus on the tracking of cosmic tracks .
5
6 Test:
7 1. Simulate random tracks - AliTracker used to propagate track
8 2. Missalign and smear space points
9 3. Fit the tracks
10
11 4. Visualize results/residual/pulls using tree draw
12
13 The input MC paramters are stored together with reconstructed output in the tree.
14
15
16 Usage:
17 .L AliTrackerTest.C+g
18 AliTrackerTest(10000)
19 TFile f("fit.root");
20
21 Make a simple plots:
22 //Pulls in P3 parameter - tangent lambda
23 fit->Draw("(seedF0.fP[3]-mcU.fP[3])/sqrt(seedF0.fC[9])>>hisSP3(100,-5,5)","")
24 //Pulls of seeding for seed
25 fit->Draw("(seedF0.fP[4]-mcU.fP[4])/sqrt(seedF0.fC[14])>>hisSP4(100,-5,5)","")
26 //
27 //Pulls in P3 parameter - tangent lambda
28 fit->Draw("(trackF0.fP[3]-mcU.fP[3])/sqrt(trackF02.fC[9])>>hisP3(100,-5,5)","")
29 //Pulls in P4 parameter
30 fit->Draw("(trackF0.fP[4]-mcU.fP[4])/sqrt(trackF02.fC[14])>>hisP4(100,-5,5)","")
31
32 //
33 //
34 // seeding procedure test (errors scaled 0.01 cm)
35 fit->Draw("(seedF0.Pt()-mcD.Pt())/mcD.Pt()^2:mcD.Pt()>>hisSPt(10,0,400,100,-0.005,0.005)","abs(seedF0.fP[4])<3&&abs(pointsU.fNPoints+pointsD.fNPoints)>300","colz");
36 //"standart track
37 fit->Draw("(trackF0.Pt()-mcU.Pt())/mcU.Pt()^2:mcU.Pt()>>his0Pt(10,0,400,100,-0.005,0.005)","abs(trackF0.fP[4])<3&&abs(pointsU.fNPoints+pointsD.fNPoints)>300","colz");
38 // misaligned track - up down
39 fit->Draw("(trackF02.Pt()-mcU.Pt())/mcU.Pt()^2:mcU.Pt()>>his2Pt(10,0,400,100,-0.005,0.005)","abs(trackF02.fP[4])<3&&abs(pointsU.fNPoints+pointsD.fNPoints)>300","colz");
40 hisSPt->FitSlicesY();
41 his0Pt->FitSlicesY();
42 his2Pt->FitSlicesY();
43
44*/
45
46//ROOT includes
47#include <iostream>
48#include "TRandom.h"
49#include "TRandom3.h"
50#include "TNtuple.h"
51#include "TStopwatch.h"
52#include "TDatabasePDG.h"
53#include "TMath.h"
54#include "TGeoManager.h"
55#include "TClonesArray.h"
56#include "TTree.h"
57#include "TFile.h"
58//AliRoot includes
59#include "AliGeomManager.h"
60#include "AliMagF.h"
61#include "AliESDVertex.h"
62#include "AliExternalTrackParam.h"
63#include "TTreeStream.h"
64#include "AliTrackPointArray.h"
65#include "AliTrackerBase.h"
66#include "AliTracker.h"
67
68//
69//
70#include "refitTrack.C"
71
72void TestRotateMI(AliExternalTrackParam *paramIn, TTreeSRedirector *pcstream);
73
74
75Int_t simulateCosmicTrack(AliExternalTrackParam &paramU, AliExternalTrackParam &paramD,
76 AliTrackPointArray &pointsU, AliTrackPointArray &pointsD,
77 AliTrackPointArray &pointsF){
78 //
79 // Toy - Simulate cosmic muon track
80 // space points genearted with the step 1 cm - at given radius
81 //
82 // Return value - number of points
83 // paramU - parameters at the beginning of track
84 // paramD - parameters at the end of track
85 // pointsU - points in upper half of
86 // pointsD - points in lower part
87 // pointsF - all points along trajectory
88 //
89 Double_t rTPC1=250;
90 Double_t rTPC0=80;
91 Double_t kMuon = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
92 Double_t kMaxSnp = 0.85;
93 Double_t xyz[3]={0,0,0};
94 Double_t xyzdw[3]={0,0,0};
95 Double_t pxyzup[3]={0,0,0};
96 Double_t pxyzdw[3]={0,0,0};
97 Double_t pxyz[3]={0,0,0};
98 Double_t cv[21];
99 for (Int_t i=0; i<21; i++) cv[i]=0;
100 Float_t covPoint[6]={0,0,0,0,0,0};
101 UShort_t volid=0;
102 //
103 static TClonesArray arrup("AliTrackPoint",500);
104 static TClonesArray arrdw("AliTrackPoint",500);
105 static TClonesArray arr("AliTrackPoint",500);
106 if (arrup[0]==0){
107 // init point arrays
108 for (Int_t i=0;i<500;i++){
109 new (arrup[i]) AliTrackPoint;
110 new (arrdw[i]) AliTrackPoint;
111 new (arr[i]) AliTrackPoint;
112 }
113 }
114
115 //
116 xyz[0]=(gRandom->Rndm()-0.5)*250;
117 xyz[1]=250;
118 xyz[2]=(gRandom->Rndm()-0.5)*100;
119 //
120 Double_t pt = gRandom->Exp(8.);
121 if (gRandom->Rndm()>0.3) pt= gRandom->Exp(50.);
122 if (gRandom->Rndm()>0.6) pt= gRandom->Rndm()*400;
123 // Double_t pt = gRandom->Exp(60.);
124 Double_t pz = gRandom->Gaus(0,0.3)*pt;
125 Double_t phi = gRandom->Gaus(0,0.3);
126 pxyz[0] = pt*TMath::Sin(phi);
127 pxyz[1] = pt*TMath::Cos(phi);
128 pxyz[2] = pz;
129 Short_t sign= (gRandom->Rndm()>0.5)? -1:1;
130 //
131 pxyzup[0]=pxyz[0];
132 pxyzup[1]=pxyz[1];
133 pxyzup[2]=pxyz[2];
134 //
135 AliExternalTrackParam lparamU(xyz, pxyzup,cv,sign);
136 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
137 lparamU.Rotate(alpha);
138 paramU=lparamU;
139 //
140 //
141 //
142 Double_t txyzup[3]={0,0,0};
143 Double_t txyzdw[3]={0,0,0};
144 //
145 Int_t npointsUp=0;
146 AliTrackerBase::PropagateTrackTo(&lparamU,rTPC1,kMuon,3,kTRUE,kMaxSnp);
147 paramU=lparamU;
148 for (Double_t r=rTPC1; r>0; r-=1){
149 Bool_t status = AliTrackerBase::PropagateTrackTo(&lparamU,r,kMuon,3,kTRUE,kMaxSnp);
150 if (!status) break;
151 if (TMath::Abs(lparamU.GetSnp())>kMaxSnp) break;
152 if (r<rTPC0) continue;
153 lparamU.GetXYZ(txyzup);
154 new (arrup[npointsUp]) AliTrackPoint(txyzup[0],txyzup[1],txyzup[2],covPoint,volid,0.,0.,0.,0);
155 npointsUp++;
156 }
157 //
158 lparamU.GetXYZ(txyzup);
159 Double_t bz=AliTrackerBase::GetBz(txyzup);
160 Double_t maxd=100000;
161 Double_t pvtx[3]={0,0,0};
162 Double_t sigmavtx[3]={0.01,0.01,3000};
163 AliESDVertex *vtx= new AliESDVertex(pvtx,sigmavtx,"vertex");
164 Double_t dz[2]={0,0};
165 Double_t cvtx[3]={0.0,0.0,0};
166 lparamU.PropagateToDCA(vtx,bz,maxd,dz,cvtx);
167 //
168 // make step to other side
169 lparamU.PropagateTo(-30,bz);
170 lparamU.GetXYZ(xyzdw);
171 lparamU.GetPxPyPz(pxyzdw);
172 // invert the sign of the momentum
173 pxyzdw[0]=-pxyzdw[0];
174 pxyzdw[1]=-pxyzdw[1];
175 pxyzdw[2]=-pxyzdw[2];
176 Short_t sign2=-sign;
177 AliExternalTrackParam lparamD(xyzdw,pxyzdw,cv,sign2);
178 lparamU.GetXYZ(xyzdw);
179 lparamU.GetPxPyPz(pxyzdw);
180 Double_t alphadw = TMath::ATan2(xyzdw[1],xyzdw[0]);
181 lparamD.Rotate(alphadw);//I have to rotate gobal to local coordenate externalparam
182 Double_t radius0=TMath::Sqrt(xyzdw[1]*xyzdw[1]+xyzdw[0]*xyzdw[0]);
183 Int_t npointsDown=0;
184 for (Double_t r=radius0; r<rTPC1; r+=1){
185 Bool_t status = AliTrackerBase::PropagateTrackTo(&lparamD,r,kMuon,3,kTRUE,0.99);
186 if (!status) continue;
187 if (TMath::Abs(lparamD.GetSnp())>kMaxSnp) continue;
188 if(r>rTPC0){
189 lparamD.GetXYZ(txyzdw);
190 new (arrdw[npointsDown]) AliTrackPoint(txyzdw[0],txyzdw[1],txyzdw[2],covPoint,volid,0.,0.,0.,0);
191 npointsDown++;
192 }
193 }
194 //
195 // Fill MC point arrays
196 //
197 Int_t npoints=npointsUp+npointsDown;
198 AliTrackPointArray lpointsF(npoints);
199 AliTrackPointArray lpointsU(npointsUp);
200 AliTrackPointArray lpointsD(npointsDown);
201 for (Int_t i=0; i<npointsUp;i++){
202 AliTrackPoint *point = (AliTrackPoint*)arrup[i];
203 lpointsF.AddPoint(i,point);
204 lpointsU.AddPoint(i,point);
205 }
206 //
207 for (Int_t i=0; i<npointsDown;i++){
208 AliTrackPoint *point = (AliTrackPoint*)arrdw[i];
209 lpointsF.AddPoint(i+npointsUp,point);
210 lpointsD.AddPoint(i,point);
211 }
212 // export points
213 AliTrackerBase::PropagateTrackTo(&lparamD,rTPC1,kMuon,3,kTRUE,kMaxSnp);
214 paramD=lparamD;
215 //
216 pointsU=lpointsU;
217 pointsD=lpointsD;
218 pointsF=lpointsF;
219 return npoints;
220}
221
222
223AliTrackPointArray * SmearPoints(AliTrackPointArray &pointArray, Double_t sigmaY, Double_t sigmaZ, Double_t shiftY,Double_t shiftZ){
224 //
225 // Smear ideal points form the simulation
226 // 1. Smear the input points (in the "local frame")
227 // 2. Assing corresponding covariance matrix
228 // 3. Add systematic shift - shift y and shift z
229
230 Int_t npoints=pointArray.GetNPoints();
231 AliTrackPointArray *outputArray = new AliTrackPointArray(npoints);
232 Double_t xyz[3]={0,0,0};
233 Float_t covPoint[6]={0,0,0, sigmaY*sigmaY,0,sigmaZ*sigmaZ}; //covaraince at the local frame
234 //
235 //
236 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
237 AliTrackPoint pointIn;
238 pointArray.GetPoint(pointIn,ipoint);
239 Double_t alpha = TMath::ATan2(pointIn.GetY(),pointIn.GetX());
240 AliTrackPoint pr = pointIn.Rotate(alpha);
241 xyz[0]=pr.GetX(); //local x
242 xyz[1]=pr.GetY()+gRandom->Gaus(0,sigmaY); //local y
243 xyz[2]=pr.GetZ()+gRandom->Gaus(0,sigmaZ); //local z
244 if (pointIn.GetY()>0) xyz[1]+=shiftY;
245 if (pointIn.GetY()>0) xyz[2]+=shiftZ;
246 //
247 pr.SetXYZ(xyz[0],xyz[1],xyz[2],covPoint); // set covariance matrix
248 AliTrackPoint pg= pr.Rotate(-alpha);
249 AliTrackPoint prCheck= pg.Rotate(alpha);
250 outputArray->AddPoint(ipoint,&pg);
251 }
252 return outputArray;
253}
254
255
256
257AliExternalTrackParam * MakeSeed(AliTrackPointArray &pointArray, Int_t seedDelta){
258 //
259 // Example: creation of seed
260 // Make seed for array of track points
261 // seedDelta - gap between seeding points
262 //
263 Int_t npoints=pointArray.GetNPoints();
264 if(npoints<=3) return 0; //not enough points to make a trac
265 if (npoints-2*seedDelta-1<0) return 0;
266 AliTrackPoint point1;
267 AliTrackPoint point2;
268 AliTrackPoint point3;
269 pointArray.GetPoint(point1,npoints-1);
270 pointArray.GetPoint(point2,npoints-seedDelta-1);
271 pointArray.GetPoint(point3,npoints-2*seedDelta-1);
272 //
273 AliExternalTrackParam * trackParam = AliTrackerBase::MakeSeed(point1, point2, point3);
274 return trackParam;
275}
276
277
278void AliTrackerTest(Int_t ntracks) {
279 //
280 //
281 //
282 TGeoManager::Import("./geometry.root");
283 AliGeomManager::LoadGeometry("./geometry.root");
284 AliMagF::BMap_t smag = AliMagF::k5kG;
285 Double_t kMuon = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
286 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., smag));
287 //
288
289 TTreeSRedirector *pcstream = new TTreeSRedirector("fit.root");
290 AliExternalTrackParam paramU;
291 AliExternalTrackParam paramD;
292 AliTrackPointArray pointsU;
293 AliTrackPointArray pointsD;
294 AliTrackPointArray pointsF;
295
296 for (Int_t i=0;i<ntracks;i++) {
297 if (i%10==0) printf("Track\t%d\n",i);
298 Int_t npoints = simulateCosmicTrack(paramU, paramD, pointsU, pointsD, pointsF);
299 if (npoints<10) continue;
300 //
301 // point for test of seeding - with scaled errors
302 AliTrackPointArray* pointsFS=SmearPoints(pointsF, 0.01,0.01,0,0);
303 AliTrackPointArray* pointsDS=SmearPoints(pointsD, 0.01,0.01,0,0);
304 AliTrackPointArray* pointsUS=SmearPoints(pointsU, 0.01,0.01,0,0);
305 //test of seeding routines
306 AliExternalTrackParam *seedF0 = MakeSeed(*pointsFS,50);
307 AliExternalTrackParam *seedD0 = MakeSeed(*pointsDS,50);
308 AliExternalTrackParam *seedU0 = MakeSeed(*pointsUS,50);
309 // points smeeared according TPC resolution
310 AliTrackPointArray* pointsF0=SmearPoints(pointsF, 0.1,0.1,0,0);
311 AliTrackPointArray* pointsD0=SmearPoints(pointsD, 0.1,0.1,0,0);
312 AliTrackPointArray* pointsU0=SmearPoints(pointsU, 0.1,0.1,0,0);
313 // points smeared according TPCresolution - + systematic shift
314 AliTrackPointArray* pointsF02=SmearPoints(pointsF, 0.1,0.1,0.2,0.2);
315 AliTrackPointArray* pointsD02=SmearPoints(pointsD, 0.1,0.1,0.2,0.2);
316 AliTrackPointArray* pointsU02=SmearPoints(pointsU, 0.1,0.1,0.2,0.2);
317 //test of tracking
318 AliExternalTrackParam *trackF0 = MakeSeed(*pointsF0,10);
319 AliExternalTrackParam *trackD0 = MakeSeed(*pointsD0,10);
320 AliExternalTrackParam *trackU0 = MakeSeed(*pointsU0,10);
321 AliExternalTrackParam *trackF02 = MakeSeed(*pointsF02,10);
322 AliExternalTrackParam *trackD02 = MakeSeed(*pointsD02,10);
323 AliExternalTrackParam *trackU02 = MakeSeed(*pointsU02,10);
324 //
325 //
326 if (!trackF0) continue;
327 if (!trackD0) continue;
328 if (!trackU0) continue;
329 if (!seedF0) continue;
330 if (!seedD0) continue;
331 if (!seedU0) continue;
332 AliTrackerBase::FitTrack(trackF0, pointsF0,kMuon,3);
333 AliTrackerBase::FitTrack(trackU0, pointsU0,kMuon,3);
334 AliTrackerBase::FitTrack(trackD0, pointsD0,kMuon,3);
335 AliTrackerBase::FitTrack(trackF02, pointsF02,kMuon,3);
336 AliTrackerBase::FitTrack(trackD02, pointsD02,kMuon,3);
337 AliTrackerBase::FitTrack(trackU02, pointsU02,kMuon,3);
338 //
339 TestRotateMI(trackF0, pcstream);
340 if (trackF0&&trackD0&&trackU0){
341 (*pcstream)<<"fit"<<
342 "pointsD.="<<&pointsD<<
343 "pointsD0.="<<pointsD0<<
344 "pointsD02.="<<pointsD02<<
345 //
346 "pointsU.="<<&pointsU<<
347 "pointsU0.="<<pointsU0<<
348 "pointsU02.="<<pointsU02<<
349 //
350 "pointsF.="<<&pointsF<<
351 "pointsF0.="<<pointsF0<<
352 "pointsF02.="<<pointsF02<<
353 //
354 "mcU.="<<&paramU<< //MC track down
355 "mcD.="<<&paramD<< //MC track top
356 // track fit - 0 misalignemnt
357 "seedF0.="<<seedF0<< //full track
358 "seedD0.="<<seedD0<< //down track
359 "seedU0.="<<seedU0<< //up track
360 //
361 "trackF0.="<<trackF0<< //full track
362 "trackD0.="<<trackD0<< //down track
363 "trackU0.="<<trackU0<< //up track
364 // track fit - 0.2 cm misalignemnt
365 "trackF02.="<<trackF02<< //full track
366 "trackD02.="<<trackD02<< //down track
367 "trackU02.="<<trackU02<< //up track
368 "\n";
369 }
370 }
371 delete pcstream;
372
373}
374
375void TestRotateMI(AliExternalTrackParam *paramIn, TTreeSRedirector *pcstream){
376 //
377 // test of rotation function
378 // rotate by 360 degrees
379 // dump the state vector to the tree after each rotation
380 AliExternalTrackParam param(*paramIn);
381 AliExternalTrackParam paramMI(*paramIn);
382 for (Int_t idiv=0; idiv<=18; idiv++){
383 Double_t alphaRot = paramIn->GetAlpha()+2*TMath::Pi()*idiv/18.;
384 param.Rotate(alphaRot);
385 paramMI.RotateMI(alphaRot);
386 (*pcstream)<<"rotateTest"<<
387 "pIn.="<<paramIn<<
388 "pRot.="<<&param<<
389 "pRotMI.="<<&paramMI<<
390 "idiv="<<idiv<<
391 "\n";
392 }
393 /*
394 to check:
395
396 */
397}