]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSTrackerV1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added function GetModuleType to return the detector type for a specific module.
[u/mrichter/AliRoot.git] / ITS / AliITSTrackerV1.cxx
CommitLineData
04b80329 1//The purpose of this class is to permorm the ITS tracking.
2//The constructor has the task to inizialize some private members.
3//The method DoTracking is written to be called by a macro. It gets the event number, the minimum and maximum
4//order number of TPC tracks that are to be tracked trough the ITS, and the file where the recpoints are
5//registered.
6//The method Recursivetracking is a recursive function that performs the tracking trough the ITS
7//The method Intersection found the layer, ladder and detector whre the intersection take place and caluclate
8//the cohordinates of this intersection. It returns an integer that is 0 if the intersection has been found
9//successfully.
10//The two mwthods Kalmanfilter and kalmanfiltervert operate the kalmanfilter without and with the vertex
11//imposition respectively.
12//The authors thank Mariana Bondila to have help them to resolve some problems. July-2000
13
14
15#include <fstream.h>
13888096 16#include <TMath.h>
13888096 17#include <TBranch.h>
18#include <TVector.h>
13888096 19#include <TFile.h>
20#include <TTree.h>
cbcfa38d 21#include <TStopwatch.h>
22
13888096 23#include "TParticle.h"
24#include "AliRun.h"
25#include "AliITS.h"
04b80329 26#include "AliITSsegmentationSSD.h"
13888096 27#include "AliITSgeomSPD.h"
28#include "AliITSgeomSDD.h"
29#include "AliITSgeomSSD.h"
30#include "AliITSgeom.h"
13888096 31#include "AliITSRecPoint.h"
04b80329 32#include "stdlib.h"
13888096 33#include "AliKalmanTrack.h"
34#include "AliMagF.h"
04b80329 35#include "AliITSTrackV1.h"
36#include "AliITSIOTrack.h"
37#include "AliITSRad.h"
13888096 38#include "../TPC/AliTPCtracker.h"
13888096 39#include "AliITSTrackerV1.h"
40
13888096 41ClassImp(AliITSTrackerV1)
42
43
44//________________________________________________________________
45
04b80329 46AliITSTrackerV1::AliITSTrackerV1(AliITS* IITTSS, Bool_t flag) {
47//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
48// Class constructor. It does some initializations.
13888096 49
04b80329 50 fITS = IITTSS;
51 fPtref=0.;
70a9314b 52 fChi2max=0.; //aggiunto il 31-7-2001
53 fepsphi=4.; //aggiunto il 1-8-2001
54 fepsz=4.; //aggiunto il 1-8-2001
04b80329 55 fflagvert=flag;
56 Int_t imax=200,jmax=450;
57 frl = new AliITSRad(imax,jmax);
13888096 58
a3894645 59/////////////////////////////////////// gets information on geometry ///////////////////////////////////
60
61 AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
62
63 Int_t ll=1, dd=1;
64 TVector det(9);
65
66 //cout<<" nlad ed ndet \n";
cbcfa38d 67 Int_t ia;
a3894645 68 for(ia=0; ia<6; ia++) {
69 fNlad[ia]=g1->GetNladders(ia+1);
70 fNdet[ia]=g1->GetNdetectors(ia+1);
71 //cout<<fNlad[i]<<" "<<fNdet[i]<<"\n";
72 }
73 //getchar();
74
75 //cout<<" raggio medio = ";
cbcfa38d 76 Int_t ib;
a3894645 77 for(ib=0; ib<6; ib++) {
78 g1->GetCenterThetaPhi(ib+1,ll,dd,det);
cbcfa38d 79 Double_t r1=TMath::Sqrt(det(0)*det(0)+det(1)*det(1));
80 g1->GetCenterThetaPhi(ib+1,ll,dd+1,det);
81 Double_t r2=TMath::Sqrt(det(0)*det(0)+det(1)*det(1));
82 fAvrad[ib]=(r1+r2)/2.;
a3894645 83 //cout<<fAvrad[ib]<<" ";
84 }
85 //cout<<"\n"; getchar();
86
87 fDetx[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDx();
88 fDetz[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDz();
89
90 fDetx[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDx();
91 fDetz[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDz();
92
93 fDetx[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDx();
94 fDetz[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDz();
95
96 fDetx[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDx();
97 fDetz[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDz();
98
99 fDetx[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDx();
100 fDetz[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDz();
101
102 fDetx[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDx();
103 fDetz[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDz();
104
105 //cout<<" Detx Detz\n";
106 //for(Int_t la=0; la<6; la++) cout<<" "<<fDetx[la]<<" "<<fDetz[la]<<"\n";
107 //getchar();
108//////////////////////////////////////////////////////////////////////////////////////////////////////////
109
cbcfa38d 110////////////////// allocate memory and define matrices fzmin, fzmax, fphimin and fphimax /////////////////////////////////
111
112 Double_t epsz=1.2;
113 Double_t epszdrift=0.05;
114
115 fzmin = new Double_t*[6]; fzmax = new Double_t*[6];
116 Int_t im1, im2, im2max;
117 for(im1=0; im1<6; im1++) {
118 im2max=fNdet[im1];
119 fzmin[im1] = new Double_t[im2max]; fzmax[im1] = new Double_t[im2max];
120 }
121
122 for(im1=0; im1<6; im1++) {
123 im2max=fNdet[im1];
124 for(im2=0; im2<im2max; im2++) {
125 g1->GetCenterThetaPhi(im1+1,1,im2+1,det);
126 if(im2!=0) fzmin[im1][im2]=det(2)-fDetz[im1];
127 else
128 fzmin[im1][im2]=det(2)-(fDetz[im1])*epsz;
129 if(im2!=(im2max-1)) fzmax[im1][im2]=det(2)+fDetz[im1];
130 else
131 fzmax[im1][im2]=det(2)+fDetz[im1]*epsz;
132 if(im1==2 || im1==3) {fzmin[im1][im2]-=epszdrift; fzmax[im1][im2]+=epszdrift;}
133 }
134 }
135
136 fphimin = new Double_t*[6]; fphimax = new Double_t*[6];
137 for(im1=0;im1<6;im1++) {
138 im2max=fNlad[im1];
139 fphimin[im1] = new Double_t[im2max]; fphimax[im1] = new Double_t[im2max];
140 }
141
142 fphidet = new Double_t*[6];
143 for(im1=0; im1<6; im1++) {
144 im2max=fNlad[im1];
145 fphidet[im1] = new Double_t[im2max];
146 }
147
148 Float_t global[3],local[3];
149 Double_t pigre=TMath::Pi();
150 Double_t xmin,ymin,xmax,ymax;
151
152 for(im1=0; im1<6; im1++) {
153 im2max=fNlad[im1];
154 for(im2=0; im2<im2max; im2++) {
155 Int_t idet=2;
156 g1->GetCenterThetaPhi(im1+1,im2+1,idet,det);
157 fphidet[im1][im2]= TMath::ATan2(Double_t(det(1)),Double_t(det(0)));
158 if(fphidet[im1][im2]<0.) fphidet[im1][im2]+=2.*pigre;
159 local[1]=local[2]=0.;
160 local[0]= -(fDetx[im1]);
161 if(im1==0) local[0]= (fDetx[im1]); //to take into account different reference system
162 g1->LtoG(im1+1,im2+1,idet,local,global);
163 xmax=global[0]; ymax=global[1];
164 local[0]= (fDetx[im1]);
165 if(im1==0) local[0]= -(fDetx[im1]); //take into account different reference system
166 g1->LtoG(im1+1,im2+1,idet,local,global);
167 xmin=global[0]; ymin=global[1];
168 fphimin[im1][im2]= TMath::ATan2(ymin,xmin); if(fphimin[im1][im2]<0.) fphimin[im1][im2]+=2.*pigre;
169 fphimax[im1][im2]= TMath::ATan2(ymax,xmax); if(fphimax[im1][im2]<0.) fphimax[im1][im2]+=2.*pigre;
70a9314b 170 //cout<<" layer detector phimin phimax= "<<im1<<" "<<im2<<" "<<fphimin[im1][im2]<<" "<<fphimax[im1][im2]<<"\n"; getchar();
cbcfa38d 171 }
172 }
173
174//////////////////////////////////////////////////////////////////////////////////////////////////////////
175
a3894645 176//////////////////////////////////////// gets magnetic field factor ////////////////////////////////
177
178 AliMagF * fieldPointer = gAlice->Field();
179 fFieldFactor = (Double_t)fieldPointer->Factor();
180 //cout<< " field factor = "<<fFieldFactor<<"\n"; getchar();
181
182/////////////////////////////////////////////////////////////////////////////////////////////////////////
183
04b80329 184}
13888096 185
04b80329 186AliITSTrackerV1::AliITSTrackerV1(const AliITSTrackerV1 &cobj) {
187//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
188// copy constructor
13888096 189
cbcfa38d 190 *fITS = *cobj.fITS;
191 *fresult = *cobj.fresult;
192 fPtref = cobj.fPtref;
70a9314b 193 fChi2max = cobj.fChi2max; //aggiunto il 31-7-2001
194 fepsphi = cobj.fepsphi; //aggiunto il 1-8-2001
195 fepsz = cobj.fepsz; //aggiunto il 1-8-2001
cbcfa38d 196 **fvettid = **cobj.fvettid;
197 fflagvert = cobj.fflagvert;
198 Int_t imax=200,jmax=450;
199 frl = new AliITSRad(imax,jmax);
200 *frl = *cobj.frl;
201 fFieldFactor = cobj.fFieldFactor;
202 Int_t i,im1,im2,im2max;
203 for(i=0; i<6; i++) {
204 fNlad[i] = cobj.fNlad[i];
205 fNdet[i] = cobj.fNdet[i];
206 fAvrad[i] = cobj.fAvrad[i];
207 fDetx[i] = cobj.fDetx[i];
208 fDetz[i] = cobj.fDetz[i];
209 }
210 fzmin = new Double_t*[6]; fzmax = new Double_t*[6];
211 for(im1=0; im1<6; im1++) {
212 im2max=fNdet[im1];
213 fzmin[im1] = new Double_t[im2max]; fzmax[im1] = new Double_t[im2max];
214 }
215 fphimin = new Double_t*[6]; fphimax = new Double_t*[6];
216 for(im1=0;im1<6;im1++) {
217 im2max=fNlad[im1];
218 fphimin[im1] = new Double_t[im2max]; fphimax[im1] = new Double_t[im2max];
219 }
220
221 fphidet = new Double_t*[6];
222 for(im1=0; im1<6; im1++) {
223 im2max=fNlad[im1];
224 fphidet[im1] = new Double_t[im2max];
225 }
226 for(im1=0; im1<6; im1++) {
227 im2max=fNdet[im1];
228 for(im2=0; im2<im2max; im2++) {
229 fzmin[im1][im2]=cobj.fzmin[im1][im2];
230 fzmax[im1][im2]=cobj.fzmax[im1][im2];
231 }
232 }
233 for(im1=0; im1<6; im1++) {
234 im2max=fNlad[im1];
235 for(im2=0; im2<im2max; im2++) {
236 fphimin[im1][im2]=cobj.fphimin[im1][im2];
237 fphimax[im1][im2]=cobj.fphimax[im1][im2];
238 fphidet[im1][im2]=cobj.fphidet[im1][im2];
239 }
240 }
241}
242
243AliITSTrackerV1::~AliITSTrackerV1(){
244//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
245// class destructor
246
247 //cout<<" CpuTimeKalman = "<<fTimerKalman->CpuTime()<<"\n";
248 //cout<<" CpuTimeIntersection = "<<fTimerIntersection->CpuTime()<<"\n";
249 //cout<<" CpuTimeIntersection = "<<TStopwatch::GetCPUTime()<<"\n";
250 //delete fTimerKalman;
251 //delete fTimerIntersection;
13888096 252
253}
254
cbcfa38d 255
04b80329 256AliITSTrackerV1 &AliITSTrackerV1::operator=(AliITSTrackerV1 obj) {
257//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
258// assignement operator
13888096 259
cbcfa38d 260 *fITS = *obj.fITS;
261 *fresult = *obj.fresult;
262 fPtref = obj.fPtref;
70a9314b 263 fChi2max = obj.fChi2max; //aggiunto il 31-7-2001
264 fepsphi = obj.fepsphi; //aggiunto il 1-8-2001
265 fepsz = obj.fepsz; //aggiunto io 1-8-2001
cbcfa38d 266 **fvettid = **obj.fvettid;
267 fflagvert = obj.fflagvert;
268 Int_t imax=200,jmax=450;
269 frl = new AliITSRad(imax,jmax);
270 *frl = *obj.frl;
271 fFieldFactor = obj.fFieldFactor;
272 Int_t i;
273 for(i=0; i<6; i++) {
274 fNlad[i] = obj.fNlad[i];
275 fNdet[i] = obj.fNdet[i];
276 fAvrad[i] = obj.fAvrad[i];
277 fDetx[i] = obj.fDetx[i];
278 fDetz[i] = obj.fDetz[i];
279 }
280 fzmin = new Double_t*[6]; fzmax = new Double_t*[6];
281 Int_t im1, im2, im2max;
282 for(im1=0; im1<6; im1++) {
283 im2max=fNdet[im1];
284 fzmin[im1] = new Double_t[im2max]; fzmax[im1] = new Double_t[im2max];
285 }
286 fphimin = new Double_t*[6]; fphimax = new Double_t*[6];
287 for(im1=0;im1<6;im1++) {
288 im2max=fNlad[im1];
289 fphimin[im1] = new Double_t[im2max]; fphimax[im1] = new Double_t[im2max];
290 }
291
292 fphidet = new Double_t*[6];
293 for(im1=0; im1<6; im1++) {
294 im2max=fNlad[im1];
295 fphidet[im1] = new Double_t[im2max];
296 }
297 for(im1=0; im1<6; im1++) {
298 im2max=fNdet[im1];
299 for(im2=0; im2<im2max; im2++) {
300 fzmin[im1][im2]=obj.fzmin[im1][im2];
301 fzmax[im1][im2]=obj.fzmax[im1][im2];
302 }
303 }
304 for(im1=0; im1<6; im1++) {
305 im2max=fNlad[im1];
306 for(im2=0; im2<im2max; im2++) {
307 fphimin[im1][im2]=obj.fphimin[im1][im2];
308 fphimax[im1][im2]=obj.fphimax[im1][im2];
309 fphidet[im1][im2]=obj.fphidet[im1][im2];
310 }
311 }
13888096 312
04b80329 313 return *this;
314
315}
13888096 316
317
04b80329 318//________________________________________________________________
13888096 319
13888096 320
04b80329 321void AliITSTrackerV1::DoTracking(Int_t evNumber, Int_t minTr, Int_t maxTr, TFile *file) {
322//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
323//The method needs the event number, the minimum and maximum order number of TPC tracks that
324//are to be tracked trough the ITS, and the file where the recpoints are registered.
325//The method can be called by a macro. It preforms the tracking for all good TPC tracks
583731c7 326
583731c7 327
13888096 328 printf("begin DoTracking - file %p\n",file);
a3894645 329
13888096 330 struct GoodTrack {
331 Int_t lab,code;
332 Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
333 Bool_t flag;
334 };
335
336
337 gAlice->GetEvent(0);
338
a3894645 339 AliKalmanTrack *kkprov;
340 kkprov->SetConvConst(100/0.299792458/0.2/fFieldFactor);
13888096 341
04b80329 342
13888096 343 TFile *cf=TFile::Open("AliTPCclusters.root");
344 AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
345 if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
346
347 AliTPCtracker *tracker = new AliTPCtracker(digp);
04b80329 348
13888096 349 // Load clusters
350 tracker->LoadInnerSectors();
351 tracker->LoadOuterSectors();
352
353
354 GoodTrack gt[15000];
355 Int_t ngood=0;
356 ifstream in("itsgood_tracks");
357
358 cerr<<"Reading itsgood tracks...\n";
359 while (in>>gt[ngood].lab>>gt[ngood].code
360 >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
361 >>gt[ngood].x >>gt[ngood].y >>gt[ngood].z
362 >>gt[ngood].pxg >>gt[ngood].pyg >>gt[ngood].pzg
363 >>gt[ngood].ptg >>gt[ngood].flag) {
364 ngood++;
365 cerr<<ngood<<'\r';
366 if (ngood==15000) {
367 cerr<<"Too many good tracks !\n";
368 break;
369 }
370 }
371 if (!in.eof()) cerr<<"Read error (itsgood_tracks) !\n";
372
373
374// Load tracks
375 TFile *tf=TFile::Open("AliTPCtracks.root");
376 if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return ;}
377 TObjArray tracks(200000);
378 TTree *tracktree=(TTree*)tf->Get("TPCf");
379 if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n";}
380 TBranch *tbranch=tracktree->GetBranch("tracks");
381 Int_t nentr=(Int_t)tracktree->GetEntries();
382 Int_t kk;
04b80329 383
384 AliTPCtrack *ioTrackTPC=0;
13888096 385 for (kk=0; kk<nentr; kk++) {
04b80329 386 ioTrackTPC=new AliTPCtrack;
387 tbranch->SetAddress(&ioTrackTPC);
13888096 388 tracktree->GetEvent(kk);
04b80329 389 tracker->CookLabel(ioTrackTPC,0.1);
390 tracks.AddLast(ioTrackTPC);
13888096 391 }
392 delete tracker;
393 tf->Close();
394
395
396 Int_t nt = tracks.GetEntriesFast();
397 cerr<<"Number of found tracks "<<nt<<endl;
398
04b80329 399 TVector dataOut(9);
13888096 400 Int_t kkk=0;
401
402 Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
403
404 ////////////////////////////// good tracks definition in TPC ////////////////////////////////
405
406 ofstream out1 ("AliITSTrag.out");
04b80329 407 Int_t i;
13888096 408 for (i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
409 out1.close();
410
411
412 TVector vec(5);
04b80329 413 TTree *tr=gAlice->TreeR();
414 Int_t nent=(Int_t)tr->GetEntries();
415 //TClonesArray *recPoints = RecPoints(); // nuova eliminata
416 //TClonesArray *recPoints = ITS->RecPoints(); // nuova
417 //TObjArray *
418 frecPoints = fITS->RecPoints(); // nuovissima tolta
13888096 419
420 Int_t numbpoints;
421 Int_t totalpoints=0;
422 Int_t *np = new Int_t[nent];
04b80329 423 fvettid = new Int_t* [nent];
13888096 424 Int_t mod;
425
426 for (mod=0; mod<nent; mod++) {
04b80329 427 fvettid[mod]=0;
428 fITS->ResetRecPoints(); // nuova
13888096 429 gAlice->TreeR()->GetEvent(mod); //first entry in TreeR is empty
04b80329 430 numbpoints = frecPoints->GetEntries();
13888096 431 totalpoints+=numbpoints;
432 np[mod] = numbpoints;
433 //cout<<" mod = "<<mod<<" numbpoints = "<<numbpoints<<"\n"; getchar();
04b80329 434 fvettid[mod] = new Int_t[numbpoints];
13888096 435 Int_t ii;
04b80329 436 for (ii=0;ii<numbpoints; ii++) *(fvettid[mod]+ii)=0;
13888096 437 }
438
cbcfa38d 439 AliTPCtrack *track=0;
13888096 440
441
04b80329 442 if(minTr < 0) {minTr = 0; maxTr = nt-1;}
13888096 443
444/*
445 ///////////////////////////////// Definition of vertex end its error ////////////////////////////
446 ////////////////////////// In the future it will be given by a method ///////////////////////////
447 Double_t Vx=0.;
448 Double_t Vy=0.;
449 Double_t Vz=0.;
450
451 Float_t sigmavx=0.0050; // 50 microns
452 Float_t sigmavy=0.0050; // 50 microns
453 Float_t sigmavz=0.010; // 100 microns
454
455 //Vx+=gRandom->Gaus(0,sigmavx); Vy+=gRandom->Gaus(0,sigmavy); Vz+=gRandom->Gaus(0,sigmavz);
456 TVector vertex(3), ervertex(3)
457 vertex(0)=Vx; vertex(1)=Vy; vertex(2)=Vz;
458 ervertex(0)=sigmavx; ervertex(1)=sigmavy; ervertex(2)=sigmavz;
459 /////////////////////////////////////////////////////////////////////////////////////////////////
460*/
461
462
463 TTree tracktree1("TreeT","Tree with ITS tracks");
04b80329 464 AliITSIOTrack *ioTrack=0;
465 tracktree1.Branch("ITStracks","AliITSIOTrack",&ioTrack,32000,0);
13888096 466
467 ofstream out ("AliITSTra.out");
04b80329 468
469
13888096 470 Int_t j;
04b80329 471 for (j=minTr; j<=maxTr; j++) {
13888096 472 track=(AliTPCtrack*)tracks.UncheckedAt(j);
473 Int_t flaglab=0;
474 if (!track) continue;
475 ////// elimination of not good tracks ////////////
476 Int_t ilab=TMath::Abs(track->GetLabel());
477 Int_t iii;
478 for (iii=0;iii<ngood;iii++) {
479 //cout<<" ilab, gt[iii].lab = "<<ilab<<" "<<gt[iii].lab<<"\n"; getchar();
480 if (ilab==gt[iii].lab) {
481 flaglab=1;
482 ptg=gt[iii].ptg;
483 pxg=gt[iii].pxg;
484 pyg=gt[iii].pyg;
485 pzg=gt[iii].pzg;
486 break;
487 }
488 }
489 //cout<<" j flaglab = " <<j<<" "<<flaglab<<"\n"; getchar();
490 if (!flaglab) continue;
491 //cout<<" j = " <<j<<"\n"; getchar();
04b80329 492
493 ////// propagation to the end of TPC //////////////
13888096 494 Double_t xk=77.415;
495 track->PropagateTo(xk, 28.94, 1.204e-3); //Ne
496 xk -=0.01;
497 track->PropagateTo(xk, 44.77, 1.71); //Tedlar
498 xk -=0.04;
04b80329 499 track->PropagateTo(xk, 44.86, 1.45); //kevlar
13888096 500 xk -=2.0;
501 track->PropagateTo(xk, 41.28, 0.029); //Nomex
502 xk-=16;
503 track->PropagateTo(xk,36.2,1.98e-3); //C02
504 xk -=0.01;
505 track->PropagateTo(xk, 24.01, 2.7); //Al
506 xk -=0.01;
507 track->PropagateTo(xk, 44.77, 1.71); //Tedlar
508 xk -=0.04;
04b80329 509 track->PropagateTo(xk, 44.86, 1.45); //kevlar
13888096 510 xk -=0.5;
511 track->PropagateTo(xk, 41.28, 0.029); //Nomex
512
04b80329 513 ///////////////////////////////////////////////////////////////
13888096 514
04b80329 515
516 AliITSTrackV1 trackITS(*track);
517
518 if(fresult) delete fresult;
519 fresult = new AliITSTrackV1(trackITS);
520
521 AliITSTrackV1 primaryTrack(trackITS);
522
523
524 TVector vgeant(3);
525 vgeant=(*fresult).GetVertex();
13888096 526
04b80329 527 // Definition of dv and zv for vertex constraint
13888096 528 Double_t sigmaDv=0.0050; Double_t sigmaZv=0.010;
529 //Double_t sigmaDv=0.0015; Double_t sigmaZv=0.0015;
530 Double_t uniform= gRandom->Uniform();
531 Double_t signdv;
532 if(uniform<=0.5) signdv=-1.;
533 else
534 signdv=1.;
535
04b80329 536 Double_t vr=TMath::Sqrt(vgeant(0)*vgeant(0)+ vgeant(1)*vgeant(1));
537 Double_t dv=gRandom->Gaus(signdv*vr,(Float_t)sigmaDv);
538 Double_t zv=gRandom->Gaus(vgeant(2),(Float_t)sigmaZv);
13888096 539
04b80329 540 //cout<<" Dv e Zv = "<<dv<<" "<<zv<<"\n";
541 trackITS.SetDv(dv); trackITS.SetZv(zv);
13888096 542 trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv);
04b80329 543 (*fresult).SetDv(dv); (*fresult).SetZv(zv);
544 (*fresult).SetsigmaDv(sigmaDv); (*fresult).SetsigmaZv(sigmaZv);
545 primaryTrack.SetDv(dv); primaryTrack.SetZv(zv);
546 primaryTrack.SetsigmaDv(sigmaDv); primaryTrack.SetsigmaZv(sigmaZv);
13888096 547
04b80329 548 primaryTrack.PrimaryTrack(frl);
549 TVector d2=primaryTrack.Getd2();
550 TVector tgl2=primaryTrack.Gettgl2();
551 TVector dtgl=primaryTrack.Getdtgl();
13888096 552 trackITS.Setd2(d2); trackITS.Settgl2(tgl2); trackITS.Setdtgl(dtgl);
04b80329 553 (*fresult).Setd2(d2); (*fresult).Settgl2(tgl2); (*fresult).Setdtgl(dtgl);
13888096 554 /*
555 trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
04b80329 556 (*result).SetVertex(vertex); (*result).SetErrorVertex(ervertex);
13888096 557 */
04b80329 558
559
560 TList *list= new TList();
561
562 list->AddLast(&trackITS);
563
564 fPtref=TMath::Abs( (trackITS).GetPt() );
70a9314b 565 if(fPtref>1.0) fChi2max=30.; //aggiunto il 31-7-2001
566 if(fPtref<=1.0) fChi2max=20.;
567 if(fPtref<0.4 ) fChi2max=20.;
568 if(fPtref<0.2 ) fChi2max=8.;
569 if(fPtref<0.1 ) fChi2max=5.;
570 if(fPtref<0.2){fepsphi=3.; fepsz=3.;}
571
572
573 cout << "\n Pt = " << fPtref <<"\n"; //stampa
04b80329 574
575 RecursiveTracking(list); // nuova ITS
576 list->Delete();
577 delete list;
578
579 Int_t itot=-1;
580 TVector vecTotLabRef(18);
581 Int_t lay, k;
582 for(lay=5; lay>=0; lay--) {
583 TVector vecLabRef(3);
584 vecLabRef=(*fresult).GetLabTrack(lay);
585 Float_t clustZ=(*fresult).GetZclusterTrack( lay);
586 for(k=0; k<3; k++){
587 Int_t lpp=(Int_t)vecLabRef(k);
588 if(lpp>=0) {
589 TParticle *p=(TParticle*) gAlice->Particle(lpp);
590 Int_t pcode=p->GetPdgCode();
591 if(pcode==11) vecLabRef(k)=p->GetFirstMother();
592 }
593 itot++; vecTotLabRef(itot)=vecLabRef(k);
594 if(vecLabRef(k)==0. && clustZ == 0.) vecTotLabRef(itot) =-3.; }
595 }
596 Long_t labref;
597 Int_t freq;
598 (*fresult).Search(vecTotLabRef, labref, freq);
599
600
601 //if(freq < 6) labref=-labref; // cinque - sei
602 if(freq < 5) labref=-labref; // cinque - sei
603 (*fresult).SetLabel(labref);
604
13888096 605 // cout<<" progressive track number = "<<j<<"\r";
606 // cout<<j<<"\r";
04b80329 607 Int_t numOfCluster=(*fresult).GetNumClust();
70a9314b 608 cout<<" progressive track number = "<<j<<"\n"; // stampa
04b80329 609 Long_t labITS=(*fresult).GetLabel();
70a9314b 610 cout << " ITS track label = " << labITS << "\n"; // stampa
13888096 611 int lab=track->GetLabel();
70a9314b 612 cout << " TPC track label = " << lab <<"\n"; // stampa
13888096 613
614
615//propagation to vertex
616
617 Double_t rbeam=3.;
618
04b80329 619 (*fresult).Propagation(rbeam);
13888096 620
04b80329 621 Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
622 (*fresult).GetCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44);
13888096 623
04b80329 624 Double_t pt=TMath::Abs((*fresult).GetPt());
625 Double_t dr=(*fresult).GetD();
626 Double_t z=(*fresult).GetZ();
627 Double_t tgl=(*fresult).GetTgl();
628 Double_t c=(*fresult).GetC();
629 Double_t cy=c/2.;
630 Double_t dz=z-(tgl/cy)*TMath::ASin((*fresult).Arga(rbeam));
631 dz-=vgeant(2);
13888096 632
04b80329 633 // cout<<" dr e dz alla fine = "<<dr<<" "<<dz<<"\n"; getchar();
634 Double_t phi=(*fresult).Getphi();
635 Double_t phivertex = phi - TMath::ASin((*fresult).ArgA(rbeam));
13888096 636 Double_t duepi=2.*TMath::Pi();
637 if(phivertex>duepi) phivertex-=duepi;
638 if(phivertex<0.) phivertex+=duepi;
04b80329 639 Double_t dtot=TMath::Sqrt(dr*dr+dz*dz);
13888096 640
641//////////////////////////////////////////////////////////////////////////////////////////
642
643 Int_t idmodule,idpoint;
04b80329 644 if(numOfCluster >=5) { // cinque - sei
645 //if(numOfCluster ==6) { // cinque - sei
13888096 646
647
04b80329 648 AliITSIOTrack outTrack;
13888096 649
04b80329 650 ioTrack=&outTrack;
13888096 651
04b80329 652 ioTrack->SetStatePhi(phi);
653 ioTrack->SetStateZ(z);
654 ioTrack->SetStateD(dr);
655 ioTrack->SetStateTgl(tgl);
656 ioTrack->SetStateC(c);
657 Double_t radius=(*fresult).Getrtrack();
658 ioTrack->SetRadius(radius);
13888096 659 Int_t charge;
04b80329 660 if(c>0.) charge=-1; else charge=1;
661 ioTrack->SetCharge(charge);
13888096 662
663
664
04b80329 665 ioTrack->SetCovMatrix(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44);
13888096 666
667 Double_t px=pt*TMath::Cos(phivertex);
668 Double_t py=pt*TMath::Sin(phivertex);
669 Double_t pz=pt*tgl;
670
04b80329 671 Double_t xtrack=dr*TMath::Sin(phivertex);
672 Double_t ytrack=dr*TMath::Cos(phivertex);
673 Double_t ztrack=dz+vgeant(2);
674
675
676 ioTrack->SetPx(px);
677 ioTrack->SetPy(py);
678 ioTrack->SetPz(pz);
679 ioTrack->SetX(xtrack);
680 ioTrack->SetY(ytrack);
681 ioTrack->SetZ(ztrack);
682 ioTrack->SetLabel(labITS);
13888096 683
684 Int_t il;
685 for(il=0;il<6; il++){
04b80329 686 ioTrack->SetIdPoint(il,(*fresult).GetIdPoint(il));
687 ioTrack->SetIdModule(il,(*fresult).GetIdModule(il));
13888096 688 }
689 tracktree1.Fill();
690
691 //cout<<" labITS = "<<labITS<<"\n";
04b80329 692 //cout<<" phi z dr tgl c = "<<phi<<" "<<z<<" "<<dr<<" "<<tgl<<" "<<c<<"\n"; getchar();
13888096 693
04b80329 694 dataOut(kkk) = ptg; kkk++; dataOut(kkk)=labITS; kkk++; dataOut(kkk)=lab; kkk++;
13888096 695
696 for (il=0;il<6;il++) {
04b80329 697 idpoint=(*fresult).GetIdPoint(il);
698 idmodule=(*fresult).GetIdModule(il);
699 *(fvettid[idmodule]+idpoint)=1;
700 ioTrack->SetIdPoint(il,idpoint);
701 ioTrack->SetIdModule(il,idmodule);
13888096 702 }
703
cbcfa38d 704 //cout<<" +++++++++++++ pt e ptg = "<<pt<<" "<<ptg<<" ++++++++++\n"; getchar();
04b80329 705
706 ///////////////////////////////
cbcfa38d 707 Double_t difpt= (pt-ptg)/ptg*100.;
708 //cout<<" difpt = "<<difpt<<"\n"; getchar();
04b80329 709 dataOut(kkk)=difpt; kkk++;
13888096 710 Double_t lambdag=TMath::ATan(pzg/ptg);
711 Double_t lam=TMath::ATan(tgl);
712 Double_t diflam = (lam - lambdag)*1000.;
04b80329 713 dataOut(kkk) = diflam; kkk++;
13888096 714 Double_t phig=TMath::ATan2(pyg,pxg); if(phig<0) phig=2.*TMath::Pi()+phig;
715 Double_t phi=phivertex;
716
717 Double_t difphi = (phi - phig)*1000.;
04b80329 718 dataOut(kkk)=difphi; kkk++;
719 dataOut(kkk)=dtot*1.e4; kkk++;
720 dataOut(kkk)=dr*1.e4; kkk++;
721 dataOut(kkk)=dz*1.e4; kkk++;
13888096 722 Int_t r;
04b80329 723 for (r=0; r<9; r++) { out<<dataOut(r)<<" ";}
13888096 724 out<<"\n";
725 kkk=0;
726
727
04b80329 728 } // end if on numOfCluster
13888096 729 //gObjectTable->Print(); // stampa memoria
04b80329 730 } // end for (int j=minTr; j<=maxTr; j++)
13888096 731
732 out.close();
04b80329 733
13888096 734
735 static Bool_t first=kTRUE;
736 static TFile *tfile;
737
738 if(first) {
739 tfile=new TFile("itstracks.root","RECREATE");
740 //cout<<"I have opened itstracks.root file "<<endl;
741 }
742 first=kFALSE;
743 tfile->cd();
744 tfile->ls();
745
746 char hname[30];
747 sprintf(hname,"TreeT%d",evNumber);
748
749 tracktree1.Write(hname);
750
751
752
753 TTree *fAli=gAlice->TreeK();
754 TFile *fileAli=0;
755
756 if (fAli) fileAli =fAli->GetCurrentFile();
757 fileAli->cd();
758
759 ////////////////////////////////////////////////////////////////////////////////////////////////
760
761 printf("delete vectors\n");
762 if(np) delete [] np;
04b80329 763 if(fvettid) delete [] fvettid;
764 if(fresult) delete fresult;
765
766}
767
768
769
770void AliITSTrackerV1::RecursiveTracking(TList *trackITSlist) {
771
772/////////////////////// This function perform the recursive tracking in ITS detectors /////////////////////
773/////////////////////// reference is a pointer to the final best track /////////////////////
774//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
775// The authors thank Mariana Bondila to have help them to resolve some problems. July-2000
776
777 //Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; Rlayer[3]=23.8; Rlayer[4]=39.1; Rlayer[5]=43.6; //vecchio
778
70a9314b 779 Float_t sigmaphil[6], sigmazl[6];
780 sigmaphil[0]=1.44e-6/(fAvrad[0]*fAvrad[0]);
781 sigmaphil[1]=1.44e-6/(fAvrad[1]*fAvrad[1]);
782 sigmaphil[2]=1.444e-5/(fAvrad[2]*fAvrad[2]);
783 sigmaphil[3]=1.444e-5/(fAvrad[3]*fAvrad[3]);
784 sigmaphil[4]=4e-6/(fAvrad[4]*fAvrad[4]);
785 sigmaphil[5]=4e-6/(fAvrad[5]*fAvrad[5]);
786
787 sigmazl[0]=1e-2;
788 sigmazl[1]=1e-2;
789 sigmazl[2]=7.84e-4;
790 sigmazl[3]=7.84e-4;
791 sigmazl[4]=0.6889;
792 sigmazl[5]=0.6889;
793
794 Int_t index;
795 AliITSgeom *g1 = fITS->GetITSgeom();
796 AliITSRecPoint *recp;
04b80329 797 for(index =0; index<trackITSlist->GetSize(); index++) {
798 AliITSTrackV1 *trackITS = (AliITSTrackV1 *) trackITSlist->At(index);
799
800 if((*trackITS).GetLayer()==7) fresult->SetChi2(10.223e140);
801 // cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n";
802 // cout<<"fvtrack =" <<"\n";
803 // cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" "<<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" "<<(*trackITS)(4)<<"\n";
804 // cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n";
805 // cout<< " Pt = "<<(*trackITS).GetPt()<<"\n";
806 // getchar();
807 Double_t chi2Now, chi2Ref;
70a9314b 808 Float_t numClustRef = fresult->GetNumClust();
04b80329 809 if((*trackITS).GetLayer()==1 ) {
810 chi2Now = trackITS->GetChi2();
811 Float_t numClustNow = trackITS->GetNumClust();
812 if(trackITS->GetNumClust()) chi2Now /= (Double_t )trackITS->GetNumClust();
70a9314b 813 chi2Ref = fresult->GetChi2();
04b80329 814 if(fresult->GetNumClust()) chi2Ref /= (Double_t )fresult->GetNumClust();
815 //cout<<" chi2Now and chi2Ref = "<<chi2Now<<" "<<chi2Ref<<"\n";
816 if( numClustNow > numClustRef ) {*fresult = *trackITS;}
817 if((numClustNow == numClustRef )&& (chi2Now < chi2Ref)) {*fresult = *trackITS;}
818 continue;
819 }
70a9314b 820 if(trackITS->Getfnoclust()>=2) continue;
04b80329 821 Float_t numClustNow = trackITS->GetNumClust();
822 if(numClustNow) {
823 chi2Now = trackITS->GetChi2();
70a9314b 824
825 if( numClustRef >3 && chi2Now>fresult->GetChi2()) continue;
04b80329 826 //cout<<" chi2Now = "<<chi2Now<<"\n";
70a9314b 827 // commentato il 30-7-2001
828 /*
829 chi2Now/=numClustNow; //commentato il 30-7-2001
830
04b80329 831 if(fPtref > 1.0 && chi2Now > 30.) continue;
832 if((fPtref >= 0.6 && fPtref<=1.0) && chi2Now > 40.) continue;
833 if((fPtref <= 0.6 && fPtref>0.2)&& chi2Now > 40.) continue;
70a9314b 834 if(fPtref <= 0.2 && chi2Now > 8.) continue;
835 */
836
837
04b80329 838 }
839
840 Int_t layerInit = (*trackITS).GetLayer();
841 Int_t layernew = layerInit - 2; // -1 for new layer, -1 for matrix index
04b80329 842
843 TList listoftrack;
844 Int_t ladp, ladm, detp,detm,ladinters,detinters;
845 Int_t layerfin=layerInit-1;
04b80329 846 // cout<<"Prima di intersection \n";
cbcfa38d 847
848 //if(!fTimerIntersection) fTimerIntersection = new TStopwatch(); // timer
849 //fTimerIntersection->Continue(); // timer
850 Int_t outinters=Intersection(*trackITS, layerfin, ladinters, detinters);
851 //fTimerIntersection->Stop(); // timer
04b80329 852
853 // cout<<" outinters = "<<outinters<<"\n";
854 // cout<<" Layer ladder detector intersection ="<<layerfin<<" "<<ladinters<<" "<<detinters<<"\n";
855 // cout << " phiinters zinters = "<<(*trackITS)(0) << " "<<(*trackITS)(1)<<"\n"; getchar();
856
857 if(outinters==-1) continue;
858
859 Int_t flaghit=0;
860 if(outinters==0){
861 TVector toucLad(9), toucDet(9);
862 Int_t lycur=layerfin;
863 ladp=ladinters+1;
864 ladm=ladinters-1;
865 if(ladm <= 0) ladm=fNlad[layerfin-1];
866 if(ladp > fNlad[layerfin-1]) ladp=1;
867 detp=detinters+1;
868 detm=detinters-1;
869 Int_t idetot=1;
70a9314b 870 /* commentato provvisoriamente il 1-8-2001
04b80329 871 toucLad(0)=ladinters; toucLad(1)=ladm; toucLad(2)=ladp;
872 toucLad(3)=ladinters; toucLad(4)=ladm; toucLad(5)=ladp;
873 toucLad(6)=ladinters; toucLad(7)=ladm; toucLad(8)=ladp;
874 toucDet(0)=detinters; toucDet(1)=detinters; toucDet(2)=detinters;
875 if(detm > 0 && detp <= fNdet[layerfin-1]) {
876 idetot=9;
877 toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm;
878 toucDet(6)=detp; toucDet(7)=detp; toucDet(8)=detp;
879 }
880
881 if(detm > 0 && detp > fNdet[layerfin-1]) {
882 idetot=6;
883 toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm;
884 }
885
886 if(detm <= 0 && detp <= fNdet[layerfin-1]) {
887 idetot=6;
888 toucDet(3)=detp; toucDet(4)=detp; toucDet(5)=detp;
889 }
70a9314b 890 */
891
892////////////////////////////////////////////////////////////////////////////
893//// nuova definizione idetot e toucLad e toucDet si puo' trasformare in un metodo
894 Float_t pigre=TMath::Pi();
b03a3bf1 895 Float_t rangephi=0.0, rangez=0.0;
70a9314b 896 if(layerfin==1 || layerfin ==2){
897 rangephi=30.*fepsphi*TMath::Sqrt(sigmaphil[layerfin-1]+(*trackITS).GetSigmaphi());
898 rangez = 30.*fepsz*TMath::Sqrt(sigmazl[layerfin-1]+(*trackITS).GetSigmaZ());
899 }
900 if(layerfin==3 || layerfin ==4){
901 rangephi=30.*fepsphi*TMath::Sqrt(sigmaphil[layerfin-1]+(*trackITS).GetSigmaphi());
902 rangez = 40.*fepsz*TMath::Sqrt(sigmazl[layerfin-1]+(*trackITS).GetSigmaZ());
903 }
904 if(layerfin==5 || layerfin ==6){
905 rangephi=20.*fepsphi*TMath::Sqrt(sigmaphil[layerfin-1]+(*trackITS).GetSigmaphi());
906 rangez =5.*fepsz*TMath::Sqrt(sigmazl[layerfin-1]+(*trackITS).GetSigmaZ());
907 }
908 Float_t phinters, zinters;
909 phinters=(*trackITS).Getphi();
910 zinters=(*trackITS).GetZ();
911 Float_t distz;
912 Float_t phicm, phicp, distphim, distphip;
913 phicm=phinters;
914 if(phinters>fphimax[layerfin-1][ladm]) phicm=phinters-2*pigre;
915 distphim=TMath::Abs(phicm-fphimax[layerfin-1][ladm]);
916 phicp=phinters;
917 if(phinters>fphimin[layerfin-1][ladp]) phicp=phinters-2.*pigre;
918 distphip=TMath::Abs(phicp-fphimin[layerfin-1][ladp]);
919 Int_t flagzmin=0;
920 Int_t flagzmax=0;
921 idetot=1;
922 toucLad(0)=ladinters; toucDet(0)=detinters;
923 if(detm>0) distz=TMath::Abs(zinters-fzmax[layerfin-1][detm-1]);
924 if(detm>0 && rangez>=distz){
925 flagzmin=1;
926 idetot++; toucLad(idetot-1)=ladinters; toucDet(idetot-1)=detm;
927 if(rangephi>=distphim){
928 idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detinters;
929 idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detm;
930 }
931 if(rangephi>=distphip){
932 idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detinters;
933 idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detm;
934 }
935 } //end detm>0....
936 if(detp<=fNdet[layerfin-1]) distz=TMath::Abs(zinters-fzmin[layerfin-1][detp-1]);
937 if(detp<=fNdet[layerfin-1] && rangez>=distz){
938 flagzmax=1;
939 idetot++; toucLad(idetot-1)=ladinters; toucDet(idetot-1)=detp;
940 if(rangephi>=distphim){
941 idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detp;
942 if(flagzmin == 0) {idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detinters;}
943 }
944 if(rangephi>=distphip){
945 idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detp;
946 if(flagzmin == 0) {idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detinters;}
947 }
948 } //end detm<fNdet[.......
949
950
951 if(flagzmin == 0 && flagzmax==0){
952 if(rangephi>=distphim){idetot++; toucLad(idetot-1)=ladm; toucDet(idetot-1)=detinters;}
953 if(rangephi>=distphip){idetot++; toucLad(idetot-1)=ladp; toucDet(idetot-1)=detinters;}
954 }
955
956///////////////////////////////////////////////////////////////////////////////////////////////////
957
958 Int_t iriv;
04b80329 959 for (iriv=0; iriv<idetot; iriv++) { //for on detectors
04b80329 960
961 ////////////////////////////////////////////////////////////////////////////////////////////////
962
963 /*** Rec points sorted by module *****/
964 /**************************************/
965
966 Int_t index;
70a9314b 967 // AliITSRecPoint *recp; //spostata il 2-8-2001
968 index = g1->GetModuleIndex(lycur,toucLad(iriv),toucDet(iriv));
969 fITS->ResetRecPoints();
970 gAlice->TreeR()->GetEvent(index);
04b80329 971 Int_t npoints=frecPoints->GetEntries();
70a9314b 972 /*
973 Int_t *indlist=new Int_t[npoints+1];
04b80329 974 Int_t counter=0;
975 Int_t ind;
976 for (ind=0; ind<=npoints; ind++) {
977 indlist[ind]=-1;
978 if (*(fvettid[index]+ind)==0) {
979 indlist[counter]=ind;
980 counter++;
981 }
982 }
983
984 ind=-1;
985
986 for(;;) {
987 ind++;
988 if(indlist[ind] < 0) recp=0;
989 else recp = (AliITSRecPoint*)frecPoints->UncheckedAt(indlist[ind]);
990
991 if((!recp) ) break;
70a9314b 992 */
993
994 Int_t indnuovo;
995 for(indnuovo=0; indnuovo<npoints; indnuovo++){
996 if (*(fvettid[index]+indnuovo)==0) recp = (AliITSRecPoint*)frecPoints->UncheckedAt(indnuovo);
997 else
998 continue;
999
04b80329 1000 TVector cluster(3),vecclust(9);
1001 vecclust(6)=vecclust(7)=vecclust(8)=-1.;
1002 Double_t sigma[2];
1003 // set veclust in global
1004 Float_t global[3], local[3];
1005 local[0]=recp->GetX();
1006 local[1]=0.;
1007 local[2]= recp->GetZ();
04b80329 1008 Int_t play = lycur;
1009 Int_t plad = TMath::Nint(toucLad(iriv));
70a9314b 1010 Int_t pdet = TMath::Nint(toucDet(iriv));
1011 // cout<<" lay lad det = "<<play<<" "<<plad<<" "<<pdet<<" x z = "<<local[0]<<" " <<local[2]<<"\n";
1012
04b80329 1013 g1->LtoG(play,plad,pdet,local,global);
1014
1015 vecclust(0)=global[0];
1016 vecclust(1)=global[1];
1017 vecclust(2)=global[2];
1018
1019
70a9314b 1020 vecclust(3) = (Float_t)recp->fTracks[0];
1021 vecclust(4) = (Float_t)indnuovo;
1022 // vecclust(4) = (float)indlist[ind];
1023 vecclust(5) = (Float_t)index;
1024 vecclust(6) = (Float_t)recp->fTracks[0];
1025 vecclust(7) = (Float_t)recp->fTracks[1];
1026 vecclust(8) = (Float_t)recp->fTracks[2];
04b80329 1027
1028 sigma[0] = (Double_t) recp->GetSigmaX2();
70a9314b 1029 sigma[1] = (Double_t) recp->GetSigmaZ2();
04b80329 1030
1031 //now we are in r,phi,z in global
1032 cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+vecclust(1)*vecclust(1));//r hit
70a9314b 1033 cluster(1) = TMath::ATan2(vecclust(1),vecclust(0)); if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi();
04b80329 1034 cluster(2) = vecclust(2); // z hit
1035
70a9314b 1036
04b80329 1037 Float_t sigmatotphi, sigmatotz;
1038
1039 //Float_t epsphi=3.2, epsz=3.;
70a9314b 1040 //Float_t epsphi=4.0, epsz=4.0; //commentato il 1-8-2001
1041 // if(fPtref<0.2) {epsphi=3.; epsz=3.;} //commentato il 1-8-2001
04b80329 1042
1043 Double_t rTrack=(*trackITS).Getrtrack();
70a9314b 1044 Double_t sigmaphi=sigma[0]/(cluster(0)*cluster(0));
1045 sigmatotphi=fepsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
04b80329 1046
70a9314b 1047 sigmatotz=fepsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
04b80329 1048 //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
1049 //if(vecclust(3)==481) getchar();
1050 if(cluster(1)<6. && (*trackITS).Getphi()>6.) cluster(1)=cluster(1)+(2.*TMath::Pi());
1051 if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi());
1052 if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue;
1053 // cout<<" supero sigmaphi \n";
1054 AliITSTrackV1 *newTrack = new AliITSTrackV1((*trackITS));
1055 (*newTrack).SetLayer((*trackITS).GetLayer()-1);
1056
1057 if (TMath::Abs(rTrack-cluster(0))/rTrack>1e-6)
1058 (*newTrack).Correct(Double_t(cluster(0)));
1059 //cout<<" cluster(2) e (*newTrack).GetZ() = "<<cluster(2)<<" "<< (*newTrack).GetZ()<<"\n";
1060 if(TMath::Abs(cluster(2)-(*newTrack).GetZ()) > sigmatotz){
1061 delete newTrack;
1062 continue;}
70a9314b 1063
1064 Double_t sigmanew[2];
1065 sigmanew[0]= sigmaphi;
1066 sigmanew[1]=sigma[1];
1067
1068 Double_t m[2];
1069 m[0]=cluster(1);
1070 m[1]=cluster(2);
1071 Double_t chi2pred=newTrack->GetPredChi2(m,sigmanew);
1072 // cout<<" chi2pred = "<<chi2pred<<"\n";
1073 /*
1074 if(fPtref>1.0 && chi2pred >30.) continue;
1075 if(fPtref<=1.0 && chi2pred >20.) continue;
1076 if(fPtref<0.4 && chi2pred >20.) continue;
1077 if(fPtref<0.2 && chi2pred >8.) continue;
1078 if(fPtref<0.1 && chi2pred >5.) continue;
1079 */
1080 if(chi2pred>fChi2max) continue; //aggiunto il 30-7-2001
1081
04b80329 1082 if(iriv == 0) flaghit=1;
1083
1084 (*newTrack).AddMS(frl); // add the multiple scattering matrix to the covariance matrix
70a9314b 1085 (*newTrack).AddEL(frl,1.,0);
04b80329 1086
70a9314b 1087
04b80329 1088
cbcfa38d 1089 //if(!fTimerKalman) fTimerKalman = new TStopwatch(); // timer
70a9314b 1090 //fTimerKalman->Continue(); // timer
1091
cbcfa38d 1092 if(fflagvert){
70a9314b 1093 // KalmanFilterVert(newTrack,cluster,sigmanew);
1094 KalmanFilterVert(newTrack,cluster,sigmanew,chi2pred); //modificata il 30-7-2001
cbcfa38d 1095 }
1096 else{
04b80329 1097 KalmanFilter(newTrack,cluster,sigmanew);
cbcfa38d 1098 }
1099 //fTimerKalman->Stop(); // timer
04b80329 1100
1101
1102 (*newTrack).PutCluster(layernew, vecclust);
1103 newTrack->AddClustInTrack();
1104
1105 listoftrack.AddLast(newTrack);
1106
70a9314b 1107 } // end of for(;;) on rec points
1108
1109 // delete [] indlist;
04b80329 1110
1111 } // end of for on detectors
1112
1113 }//end if(outinters==0)
1114
1115 if(flaghit==0 || outinters==-2) {
70a9314b 1116 AliITSTrackV1 *newTrack = new AliITSTrackV1(*trackITS);
1117 (*newTrack).Setfnoclust();
04b80329 1118 (*newTrack).SetLayer((*trackITS).GetLayer()-1);
1119 (*newTrack).AddMS(frl); // add the multiple scattering matrix to the covariance matrix
1120 (*newTrack).AddEL(frl,1.,0);
1121
1122 listoftrack.AddLast(newTrack);
1123 }
1124
1125
1126 //gObjectTable->Print(); // stampa memoria
1127
1128 RecursiveTracking(&listoftrack);
1129 listoftrack.Delete();
1130 } // end of for on tracks
1131
1132 //gObjectTable->Print(); // stampa memoria
1133
1134}
1135
cbcfa38d 1136Int_t AliITSTrackerV1::Intersection(AliITSTrackV1 &track, Int_t layer, Int_t &ladder, Int_t &detector) {
04b80329 1137//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
1138// Found the intersection and the detector
1139
cbcfa38d 1140 Double_t rk=fAvrad[layer-1];
04b80329 1141 if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;}
1142 track.Propagation(rk);
1143 Double_t zinters=track.GetZ();
1144 Double_t phinters=track.Getphi();
1145 //cout<<"zinters = "<<zinters<<" phinters = "<<phinters<<"\n";
04b80329 1146
1147 TVector det(9);
1148 TVector listDet(2);
1149 TVector distZCenter(2);
04b80329 1150
1151 Int_t iz=0;
04b80329 1152 Int_t iD;
cbcfa38d 1153
04b80329 1154 for(iD = 1; iD<= fNdet[layer-1]; iD++) {
cbcfa38d 1155 if(zinters > fzmin[layer-1][iD-1] && zinters <= fzmax[layer-1][iD-1]) {
1156 if(iz>1) {
1157 cout<< " Errore su iz in Intersection \n"; getchar();
1158 }
04b80329 1159 else {
1160 listDet(iz)= iD; distZCenter(iz)=TMath::Abs(zinters-det(2)); iz++;
1161 }
1162 }
1163 }
13888096 1164
04b80329 1165 if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;}
1166 detector=Int_t (listDet(0));
1167 if(iz>1 && (distZCenter(0)>distZCenter(1))) detector=Int_t (listDet(1));
cbcfa38d 1168
04b80329 1169 TVector listLad(2);
1170 TVector distPhiCenter(2);
1171 Int_t ip=0;
1172 Double_t pigre=TMath::Pi();
1173
1174 Int_t iLd;
cbcfa38d 1175 for(iLd = 1; iLd<= fNlad[layer-1]; iLd++) {
1176 Double_t phimin=fphimin[layer-1][iLd-1];
1177 Double_t phimax=fphimax[layer-1][iLd-1];
1178 Double_t phidet=fphidet[layer-1][iLd-1];
1179 Double_t phiconfr=phinters;
70a9314b 1180 if(phimin>phimax) {
1181 // if(phimin <5.5) {cout<<" Error in Intersection for phi \n"; getchar();}
1182 phimin-=(2.*pigre);
cbcfa38d 1183 if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
1184 if(phidet>(1.5*pigre)) phidet-=(2.*pigre);
04b80329 1185 }
cbcfa38d 1186 if(phiconfr>phimin && phiconfr<= phimax) {
1187 if(ip>1) {
1188 cout<< " Errore su ip in Intersection \n"; getchar();
1189 }
04b80329 1190 else {
1191 listLad(ip)= iLd; distPhiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++;
1192 }
1193 }
1194 }
1195 if(ip==0) { cout<< " No detector along phi \n"; getchar();}
1196 ladder=Int_t (listLad(0));
cbcfa38d 1197 if(ip>1 && (distPhiCenter(0)>distPhiCenter(1))) ladder=Int_t (listLad(1));
04b80329 1198
1199 return 0;
13888096 1200}
04b80329 1201
04b80329 1202void AliITSTrackerV1::KalmanFilter(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2]){
1203//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
1204// Kalman filter without vertex constraint
1205
1206
1207 ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
1208
1209 Double_t m[2];
1210 Double_t rk,phik,zk;
1211 rk=cluster(0); phik=cluster(1); zk=cluster(2);
1212 m[0]=phik; m[1]=zk;
1213
1214 ///////////////////////////////////// Evaluation of the error matrix V ///////////////////////////////
1215
1216 Double_t v00=sigma[0];
1217 Double_t v11=sigma[1];
1218
1219 ///////////////////////////////////////////////////////////////////////////////////////////
1220
1221
1222 Double_t cin00,cin10,cin20,cin30,cin40,cin11,cin21,cin31,cin41,cin22,cin32,cin42,cin33,cin43,cin44;
1223
1224 newTrack->GetCElements(cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40,
1225 cin41,cin42,cin43,cin44); //get C matrix
1226
1227 Double_t rold00=cin00+v00;
1228 Double_t rold10=cin10;
1229 Double_t rold11=cin11+v11;
1230
1231//////////////////////////////////// R matrix inversion ///////////////////////////////////////////////
1232
1233 Double_t det=rold00*rold11-rold10*rold10;
1234 Double_t r00=rold11/det;
1235 Double_t r10=-rold10/det;
1236 Double_t r11=rold00/det;
1237
1238////////////////////////////////////////////////////////////////////////////////////////////////////////
1239
1240 Double_t k00=cin00*r00+cin10*r10;
1241 Double_t k01=cin00*r10+cin10*r11;
1242 Double_t k10=cin10*r00+cin11*r10;
1243 Double_t k11=cin10*r10+cin11*r11;
1244 Double_t k20=cin20*r00+cin21*r10;
1245 Double_t k21=cin20*r10+cin21*r11;
1246 Double_t k30=cin30*r00+cin31*r10;
1247 Double_t k31=cin30*r10+cin31*r11;
1248 Double_t k40=cin40*r00+cin41*r10;
1249 Double_t k41=cin40*r10+cin41*r11;
1250
1251 Double_t x0,x1,x2,x3,x4;
1252 newTrack->GetXElements(x0,x1,x2,x3,x4); // get the state vector
1253
1254 Double_t savex0=x0, savex1=x1;
1255
1256 x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1);
1257 x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1);
1258 x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1);
1259 x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1);
1260 x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1);
1261
1262 Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44;
1263
1264 c00=cin00-k00*cin00-k01*cin10;
1265 c10=cin10-k00*cin10-k01*cin11;
1266 c20=cin20-k00*cin20-k01*cin21;
1267 c30=cin30-k00*cin30-k01*cin31;
1268 c40=cin40-k00*cin40-k01*cin41;
1269
1270 c11=cin11-k10*cin10-k11*cin11;
1271 c21=cin21-k10*cin20-k11*cin21;
1272 c31=cin31-k10*cin30-k11*cin31;
1273 c41=cin41-k10*cin40-k11*cin41;
1274
1275 c22=cin22-k20*cin20-k21*cin21;
1276 c32=cin32-k20*cin30-k21*cin31;
1277 c42=cin42-k20*cin40-k21*cin41;
1278
1279 c33=cin33-k30*cin30-k31*cin31;
1280 c43=cin43-k30*cin40-k31*cin41;
1281
1282 c44=cin44-k40*cin40-k41*cin41;
1283
1284 newTrack->PutXElements(x0,x1,x2,x3,x4); // put the new state vector
1285
1286 newTrack->PutCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44); // put in track the
1287 // new cov matrix
1288 Double_t vmcold00=v00-c00;
1289 Double_t vmcold10=-c10;
1290 Double_t vmcold11=v11-c11;
1291
1292///////////////////////////////////// Matrix vmc inversion ////////////////////////////////////////////////
1293
1294 det=vmcold00*vmcold11-vmcold10*vmcold10;
1295 Double_t vmc00=vmcold11/det;
1296 Double_t vmc10=-vmcold10/det;
1297 Double_t vmc11=vmcold00/det;
1298
1299////////////////////////////////////////////////////////////////////////////////////////////////////////////
1300
1301 Double_t chi2=(m[0]-x0)*( vmc00*(m[0]-x0) + 2.*vmc10*(m[1]-x1) ) +
1302 (m[1]-x1)*vmc11*(m[1]-x1);
1303
1304 newTrack->SetChi2(newTrack->GetChi2()+chi2);
1305
1306}
1307
1308
70a9314b 1309//void AliITSTrackerV1::KalmanFilterVert(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2]){
1310void AliITSTrackerV1::KalmanFilterVert(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2], Double_t chi2pred){
04b80329 1311//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
1312// Kalman filter with vertex constraint
1313
1314 ////////////////////////////// Evaluation of the measurement vector m ///////////////
1315
1316 Double_t m[4];
1317 Double_t rk,phik,zk;
1318 rk=cluster(0); phik=cluster(1); zk=cluster(2);
1319 m[0]=phik; m[1]=zk;
1320
1321 Double_t cc=(*newTrack).GetC();
1322 Double_t zv=(*newTrack).GetZv();
1323 Double_t dv=(*newTrack).GetDv();
1324 Double_t cy=cc/2.;
1325 Double_t tgl= (zk-zv)*cy/TMath::ASin(cy*rk);
1326 m[2]=dv; m[3]=tgl;
1327
1328 ///////////////////////////////////// Evaluation of the error matrix V //////////////
1329 Int_t layer=newTrack->GetLayer();
1330 Double_t v00=sigma[0];
1331 Double_t v11=sigma[1];
1332 Double_t v31=sigma[1]/rk;
1333 Double_t sigmaDv=newTrack->GetsigmaDv();
1334 Double_t v22=sigmaDv*sigmaDv + newTrack->Getd2(layer-1);
1335 Double_t v32=newTrack->Getdtgl(layer-1);
1336 Double_t sigmaZv=newTrack->GetsigmaZv();
1337 Double_t v33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(layer-1);
1338 ///////////////////////////////////////////////////////////////////////////////////////
1339
1340 Double_t cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40,cin41,cin42,cin43,cin44;
1341
1342 newTrack->GetCElements(cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40,
1343 cin41,cin42,cin43,cin44); //get C matrix
1344
1345 Double_t r[4][4];
1346 r[0][0]=cin00+v00;
1347 r[1][0]=cin10;
1348 r[2][0]=cin20;
1349 r[3][0]=cin30;
1350 r[1][1]=cin11+v11;
1351 r[2][1]=cin21;
1352 r[3][1]=cin31+sigma[1]/rk;
1353 r[2][2]=cin22+sigmaDv*sigmaDv+newTrack->Getd2(layer-1);
1354 r[3][2]=cin32+newTrack->Getdtgl(layer-1);
1355 r[3][3]=cin33+(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(layer-1);
1356
1357 r[0][1]=r[1][0]; r[0][2]=r[2][0]; r[0][3]=r[3][0]; r[1][2]=r[2][1]; r[1][3]=r[3][1];
1358 r[2][3]=r[3][2];
1359
1360///////////////////// Matrix R inversion ////////////////////////////////////////////
1361
1362 const Int_t kn=4;
1363 Double_t big, hold;
1364 Double_t d=1.;
1365 Int_t ll[kn],mm[kn];
1366
1367 Int_t i,j,k;
1368
1369 for(k=0; k<kn; k++) {
1370 ll[k]=k;
1371 mm[k]=k;
1372 big=r[k][k];
1373 for(j=k; j<kn ; j++) {
1374 for (i=j; i<kn; i++) {
1375 if(TMath::Abs(big) < TMath::Abs(r[i][j]) ) { big=r[i][j]; ll[k]=i; mm[k]=j; }
1376 }
1377 }
1378//
1379 j= ll[k];
1380 if(j > k) {
1381 for(i=0; i<kn; i++) { hold=-r[k][i]; r[k][i]=r[j][i]; r[j][i]=hold; }
1382
1383 }
1384//
1385 i=mm[k];
1386 if(i > k ) {
1387 for(j=0; j<kn; j++) { hold=-r[j][k]; r[j][k]=r[j][i]; r[j][i]=hold; }
1388 }
1389//
1390 if(!big) {
1391 d=0.;
1392 cout << "Singular matrix\n";
1393 }
1394 for(i=0; i<kn; i++) {
1395 if(i == k) { continue; }
1396 r[i][k]=r[i][k]/(-big);
1397 }
1398//
1399 for(i=0; i<kn; i++) {
1400 hold=r[i][k];
1401 for(j=0; j<kn; j++) {
1402 if(i == k || j == k) { continue; }
1403 r[i][j]=hold*r[k][j]+r[i][j];
1404 }
1405 }
1406//
1407 for(j=0; j<kn; j++) {
1408 if(j == k) { continue; }
1409 r[k][j]=r[k][j]/big;
1410 }
1411//
1412 d=d*big;
1413//
1414 r[k][k]=1./big;
1415 }
1416//
1417 for(k=kn-1; k>=0; k--) {
1418 i=ll[k];
1419 if(i > k) {
1420 for (j=0; j<kn; j++) {hold=r[j][k]; r[j][k]=-r[j][i]; r[j][i]=hold;}
1421 }
1422 j=mm[k];
1423 if(j > k) {
1424 for (i=0; i<kn; i++) {hold=r[k][i]; r[k][i]=-r[j][i]; r[j][i]=hold;}
1425 }
1426 }
1427//////////////////////////////////////////////////////////////////////////////////
1428
1429
1430 Double_t k00=cin00*r[0][0]+cin10*r[1][0]+cin20*r[2][0]+cin30*r[3][0];
1431 Double_t k01=cin00*r[1][0]+cin10*r[1][1]+cin20*r[2][1]+cin30*r[3][1];
1432 Double_t k02=cin00*r[2][0]+cin10*r[2][1]+cin20*r[2][2]+cin30*r[3][2];
1433 Double_t k03=cin00*r[3][0]+cin10*r[3][1]+cin20*r[3][2]+cin30*r[3][3];
1434 Double_t k10=cin10*r[0][0]+cin11*r[1][0]+cin21*r[2][0]+cin31*r[3][0];
1435 Double_t k11=cin10*r[1][0]+cin11*r[1][1]+cin21*r[2][1]+cin31*r[3][1];
1436 Double_t k12=cin10*r[2][0]+cin11*r[2][1]+cin21*r[2][2]+cin31*r[3][2];
1437 Double_t k13=cin10*r[3][0]+cin11*r[3][1]+cin21*r[3][2]+cin31*r[3][3];
1438 Double_t k20=cin20*r[0][0]+cin21*r[1][0]+cin22*r[2][0]+cin32*r[3][0];
1439 Double_t k21=cin20*r[1][0]+cin21*r[1][1]+cin22*r[2][1]+cin32*r[3][1];
1440 Double_t k22=cin20*r[2][0]+cin21*r[2][1]+cin22*r[2][2]+cin32*r[3][2];
1441 Double_t k23=cin20*r[3][0]+cin21*r[3][1]+cin22*r[3][2]+cin32*r[3][3];
1442 Double_t k30=cin30*r[0][0]+cin31*r[1][0]+cin32*r[2][0]+cin33*r[3][0];
1443 Double_t k31=cin30*r[1][0]+cin31*r[1][1]+cin32*r[2][1]+cin33*r[3][1];
1444 Double_t k32=cin30*r[2][0]+cin31*r[2][1]+cin32*r[2][2]+cin33*r[3][2];
1445 Double_t k33=cin30*r[3][0]+cin31*r[3][1]+cin32*r[3][2]+cin33*r[3][3];
1446 Double_t k40=cin40*r[0][0]+cin41*r[1][0]+cin42*r[2][0]+cin43*r[3][0];
1447 Double_t k41=cin40*r[1][0]+cin41*r[1][1]+cin42*r[2][1]+cin43*r[3][1];
1448 Double_t k42=cin40*r[2][0]+cin41*r[2][1]+cin42*r[2][2]+cin43*r[3][2];
1449 Double_t k43=cin40*r[3][0]+cin41*r[3][1]+cin42*r[3][2]+cin43*r[3][3];
1450
1451 Double_t x0,x1,x2,x3,x4;
1452 newTrack->GetXElements(x0,x1,x2,x3,x4); // get the state vector
1453
1454 Double_t savex0=x0, savex1=x1, savex2=x2, savex3=x3;
1455
1456 x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1)+k02*(m[2]-savex2)+
1457 k03*(m[3]-savex3);
1458 x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1)+k12*(m[2]-savex2)+
1459 k13*(m[3]-savex3);
1460 x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1)+k22*(m[2]-savex2)+
1461 k23*(m[3]-savex3);
1462 x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1)+k32*(m[2]-savex2)+
1463 k33*(m[3]-savex3);
1464 x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1)+k42*(m[2]-savex2)+
1465 k43*(m[3]-savex3);
1466
1467 Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44;
1468
1469 c00=cin00-k00*cin00-k01*cin10-k02*cin20-k03*cin30;
1470 c10=cin10-k00*cin10-k01*cin11-k02*cin21-k03*cin31;
1471 c20=cin20-k00*cin20-k01*cin21-k02*cin22-k03*cin32;
1472 c30=cin30-k00*cin30-k01*cin31-k02*cin32-k03*cin33;
1473 c40=cin40-k00*cin40-k01*cin41-k02*cin42-k03*cin43;
1474
1475 c11=cin11-k10*cin10-k11*cin11-k12*cin21-k13*cin31;
1476 c21=cin21-k10*cin20-k11*cin21-k12*cin22-k13*cin32;
1477 c31=cin31-k10*cin30-k11*cin31-k12*cin32-k13*cin33;
1478 c41=cin41-k10*cin40-k11*cin41-k12*cin42-k13*cin43;
1479
1480 c22=cin22-k20*cin20-k21*cin21-k22*cin22-k23*cin32;
1481 c32=cin32-k20*cin30-k21*cin31-k22*cin32-k23*cin33;
1482 c42=cin42-k20*cin40-k21*cin41-k22*cin42-k23*cin43;
1483
1484 c33=cin33-k30*cin30-k31*cin31-k32*cin32-k33*cin33;
1485 c43=cin43-k30*cin40-k31*cin41-k32*cin42-k33*cin43;
1486
1487 c44=cin44-k40*cin40-k41*cin41-k42*cin42-k43*cin43;
1488
1489 newTrack->PutXElements(x0,x1,x2,x3,x4); // put the new state vector
1490
1491 newTrack->PutCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44); // put in track the
1492 // new cov matrix
1493
1494 Double_t vmc[4][4];
1495
1496 vmc[0][0]=v00-c00; vmc[1][0]=-c10; vmc[2][0]=-c20; vmc[3][0]=-c30;
1497 vmc[1][1]=v11-c11; vmc[2][1]=-c21; vmc[3][1]=v31-c31;
1498 vmc[2][2]=v22-c22; vmc[3][2]=v32-c32;
1499 vmc[3][3]=v33-c33;
1500
1501 vmc[0][1]=vmc[1][0]; vmc[0][2]=vmc[2][0]; vmc[0][3]=vmc[3][0];
1502 vmc[1][2]=vmc[2][1]; vmc[1][3]=vmc[3][1];
1503 vmc[2][3]=vmc[3][2];
1504
1505
1506/////////////////////// vmc matrix inversion ///////////////////////////////////
1507
1508 d=1.;
1509
1510 for(k=0; k<kn; k++) {
1511 ll[k]=k;
1512 mm[k]=k;
1513 big=vmc[k][k];
1514 for(j=k; j<kn ; j++) {
1515 for (i=j; i<kn; i++) {
1516 if(TMath::Abs(big) < TMath::Abs(vmc[i][j]) ) { big=vmc[i][j]; ll[k]=i; mm[k]=j; }
1517 }
1518 }
1519//
1520 j= ll[k];
1521 if(j > k) {
1522 for(i=0; i<kn; i++) { hold=-vmc[k][i]; vmc[k][i]=vmc[j][i]; vmc[j][i]=hold; }
1523
1524 }
1525//
1526 i=mm[k];
1527 if(i > k ) {
1528 for(j=0; j<kn; j++) { hold=-vmc[j][k]; vmc[j][k]=vmc[j][i]; vmc[j][i]=hold; }
1529 }
1530//
1531 if(!big) {
1532 d=0.;
1533 cout << "Singular matrix\n";
1534 }
1535 for(i=0; i<kn; i++) {
1536 if(i == k) { continue; }
1537 vmc[i][k]=vmc[i][k]/(-big);
1538 }
1539//
1540 for(i=0; i<kn; i++) {
1541 hold=vmc[i][k];
1542 for(j=0; j<kn; j++) {
1543 if(i == k || j == k) { continue; }
1544 vmc[i][j]=hold*vmc[k][j]+vmc[i][j];
1545 }
1546 }
1547//
1548 for(j=0; j<kn; j++) {
1549 if(j == k) { continue; }
1550 vmc[k][j]=vmc[k][j]/big;
1551 }
1552//
1553 d=d*big;
1554//
1555 vmc[k][k]=1./big;
1556 }
1557//
1558 for(k=kn-1; k>=0; k--) {
1559 i=ll[k];
1560 if(i > k) {
1561 for (j=0; j<kn; j++) {hold=vmc[j][k]; vmc[j][k]=-vmc[j][i]; vmc[j][i]=hold;}
1562 }
1563 j=mm[k];
1564 if(j > k) {
1565 for (i=0; i<kn; i++) {hold=vmc[k][i]; vmc[k][i]=-vmc[j][i]; vmc[j][i]=hold;}
1566 }
1567 }
1568
1569
1570////////////////////////////////////////////////////////////////////////////////
1571
1572 Double_t chi2=(m[0]-x0)*( vmc[0][0]*(m[0]-x0) + 2.*vmc[1][0]*(m[1]-x1) +
1573 2.*vmc[2][0]*(m[2]-x2)+ 2.*vmc[3][0]*(m[3]-x3) ) +
1574 (m[1]-x1)* ( vmc[1][1]*(m[1]-x1) + 2.*vmc[2][1]*(m[2]-x2)+
1575 2.*vmc[3][1]*(m[3]-x3) ) +
1576 (m[2]-x2)* ( vmc[2][2]*(m[2]-x2)+ 2.*vmc[3][2]*(m[3]-x3) ) +
70a9314b 1577 (m[3]-x3)*vmc[3][3]*(m[3]-x3);
1578
1579 //cout<<" chi2 kalman = "<<chi2<<"\n"; getchar();
1580 newTrack->SetChi2(newTrack->GetChi2()+chi2); //commentata il 30-7-2001
1581 // newTrack->SetChi2(newTrack->GetChi2()+chi2pred); //aggiunta il 30-7-2001
04b80329 1582
1583}
1584