]>
Commit | Line | Data |
---|---|---|
d953664a | 1 | #include <iostream.h> |
dce1757b | 2 | #include <TMath.h> |
3 | #include <TList.h> | |
4 | #include <TTree.h> | |
5 | #include <TVector.h> | |
6 | #include <TMatrix.h> | |
7 | #include <TObjectTable.h> | |
8 | ||
9 | #include "AliITStracking.h" | |
10 | #include "AliRun.h" | |
11 | #include "AliITS.h" | |
12 | #include "AliITSgeom.h" | |
13 | #include "AliITSRecPoint.h" | |
14 | #include "AliITStrack.h" | |
15 | ||
16 | ClassImp(AliITStracking) | |
17 | ||
18 | ||
19 | AliITStracking::AliITStracking(TList *trackITSlist, AliITStrack *reference, | |
8af13b4b | 20 | AliITS *aliITS, TObjArray *rpoints, Double_t Ptref, Int_t **vettid, Bool_t flagvert, AliITSRad *rl) { |
dce1757b | 21 | /////////////////////// This function perform the tracking in ITS detectors ///////////////////// |
22 | /////////////////////// reference is a pointer to the final best track ///////////////////// | |
23 | //Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it | |
24 | // The authors thank Mariana Bondila to have help them to resolve some problems. July-2000 | |
25 | ||
dce1757b | 26 | Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; Rlayer[3]=23.8; Rlayer[4]=39.1; Rlayer[5]=43.6; |
cc6e75dd | 27 | |
28 | Int_t index; | |
29 | for(index =0; index<trackITSlist->GetSize(); index++) { | |
dce1757b | 30 | AliITStrack *trackITS = (AliITStrack *) trackITSlist->At(index); |
cc6e75dd | 31 | |
dce1757b | 32 | if((*trackITS).GetLayer()==7) reference->SetChi2(10.223e140); |
cc6e75dd | 33 | // cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n"; |
34 | // cout<<"fvtrack =" <<"\n"; | |
35 | // cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" "<<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" "<<(*trackITS)(4)<<"\n"; | |
36 | // cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n"; | |
37 | // cout<< " Pt = "<<(*trackITS).GetPt()<<"\n"; | |
38 | // getchar(); | |
dce1757b | 39 | Double_t Chi2Now, Chi2Ref; |
40 | if((*trackITS).GetLayer()==1 ) { | |
41 | Chi2Now = trackITS->GetChi2(); | |
42 | Float_t NumClustNow = trackITS->GetNumClust(); | |
43 | if(trackITS->GetNumClust()) Chi2Now /= (Double_t )trackITS->GetNumClust(); | |
44 | Chi2Ref = reference->GetChi2(); | |
45 | Float_t NumClustRef = reference->GetNumClust(); | |
46 | if(reference->GetNumClust()) Chi2Ref /= (Double_t )reference->GetNumClust(); | |
47 | //cout<<" Chi2Now and Chi2Ref = "<<Chi2Now<<" "<<Chi2Ref<<"\n"; | |
48 | if( NumClustNow > NumClustRef ) {*reference = *trackITS;} | |
49 | if((NumClustNow == NumClustRef )&& (Chi2Now < Chi2Ref)) {*reference = *trackITS;} | |
50 | continue; | |
51 | } | |
dce1757b | 52 | Float_t NumClustNow = trackITS->GetNumClust(); |
53 | if(NumClustNow) { | |
54 | Chi2Now = trackITS->GetChi2(); | |
55 | Chi2Now/=NumClustNow; | |
56 | //cout<<" Chi2Now = "<<Chi2Now<<"\n"; | |
8af13b4b | 57 | /* |
dce1757b | 58 | // if(Ptref > 0.6 && Chi2Now > 20.) continue; |
59 | if(Ptref > 0.6 && Chi2Now > 30.) continue; | |
60 | if((Ptref <= 0.6 && Ptref>0.2)&& Chi2Now > 15.) continue; | |
61 | // if(Chi2Now>5.) continue; | |
62 | //if(Chi2Now>15.) continue; | |
63 | // if(Ptref <= 0.2 && Chi2Now > 10.) continue; | |
8af13b4b | 64 | if(Ptref <= 0.2 && Chi2Now > 8.) continue; |
65 | */ | |
66 | if(Ptref > 1.0 && Chi2Now > 30.) continue; | |
67 | if((Ptref >= 0.6 && Ptref<=1.0) && Chi2Now > 40.) continue; | |
68 | if((Ptref <= 0.6 && Ptref>0.2)&& Chi2Now > 40.) continue; | |
69 | if(Ptref <= 0.2 && Chi2Now > 8.) continue; | |
dce1757b | 70 | } |
71 | ||
72 | Int_t layerInit = (*trackITS).GetLayer(); | |
73 | Int_t layernew = layerInit - 2; // -1 for new layer, -1 for matrix index | |
74 | ||
75 | Int_t NLadder[]= {20, 40, 14, 22, 34, 38}; | |
8af13b4b | 76 | Int_t NDetector[]= {4, 4, 6, 8, 23, 26}; |
dce1757b | 77 | |
78 | TList listoftrack; | |
79 | Int_t ladp, ladm, detp,detm,ladinters,detinters; | |
80 | Int_t layerfin=layerInit-1; | |
81 | Double_t Rfin=Rlayer[layerfin-1]; | |
82 | // cout<<"Prima di intersection \n"; | |
cc6e75dd | 83 | |
dce1757b | 84 | Int_t outinters=NewIntersection(*trackITS, Rfin, layerfin, ladinters, detinters); |
85 | ||
86 | // cout<<" outinters = "<<outinters<<"\n"; | |
87 | // cout<<" Layer ladder detector intersection ="<<layerfin<<" "<<ladinters<<" "<<detinters<<"\n"; | |
88 | // cout << " phiinters zinters = "<<(*trackITS)(0) << " "<<(*trackITS)(1)<<"\n"; getchar(); | |
89 | ||
90 | if(outinters==-1) continue; | |
91 | ||
cc6e75dd | 92 | Int_t flaghit=0; |
dce1757b | 93 | if(outinters==0){ |
94 | TVector Touclad(9), Toucdet(9); | |
95 | Int_t lycur=layerfin; | |
96 | ladp=ladinters+1; | |
97 | ladm=ladinters-1; | |
98 | if(ladm <= 0) ladm=NLadder[layerfin-1]; | |
99 | if(ladp > NLadder[layerfin-1]) ladp=1; | |
100 | detp=detinters+1; | |
101 | detm=detinters-1; | |
102 | Int_t idetot=1; | |
103 | Touclad(0)=ladinters; Touclad(1)=ladm; Touclad(2)=ladp; | |
104 | Touclad(3)=ladinters; Touclad(4)=ladm; Touclad(5)=ladp; | |
105 | Touclad(6)=ladinters; Touclad(7)=ladm; Touclad(8)=ladp; | |
106 | Toucdet(0)=detinters; Toucdet(1)=detinters; Toucdet(2)=detinters; | |
107 | if(detm > 0 && detp <= NDetector[layerfin-1]) { | |
108 | idetot=9; | |
109 | Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm; | |
110 | Toucdet(6)=detp; Toucdet(7)=detp; Toucdet(8)=detp; | |
111 | } | |
112 | ||
113 | if(detm > 0 && detp > NDetector[layerfin-1]) { | |
114 | idetot=6; | |
115 | Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm; | |
116 | } | |
117 | ||
118 | if(detm <= 0 && detp <= NDetector[layerfin-1]) { | |
119 | idetot=6; | |
120 | Toucdet(3)=detp; Toucdet(4)=detp; Toucdet(5)=detp; | |
cc6e75dd | 121 | } |
122 | Int_t iriv; | |
123 | for (iriv=0; iriv<idetot; iriv++) { //for on detectors | |
dce1757b | 124 | AliITSgeom *g1 = aliITS->GetITSgeom(); |
125 | TVector CTF(9); | |
126 | g1->GetCenterThetaPhi(layerInit-1,(Int_t)Touclad(iriv),(Int_t)Toucdet(iriv),CTF); | |
127 | ||
128 | // cout<<" layer, ladder, det, xo, yo, zo = "<<layerInit-1<<" "<<(Int_t)Touclad(iriv)<< | |
129 | // " "<<(Int_t)Toucdet(iriv)<<" "<<CTF(0)<<" "<<CTF(1)<<" "<<CTF(2)<< " "<<CTF(6)<<"\n"; getchar(); | |
130 | ||
131 | //////////////////////////////////////////////////////////////////////////////////////////////// | |
132 | ||
133 | /*** Rec points sorted by module *****/ | |
134 | /**************************************/ | |
135 | ||
136 | Int_t index; | |
137 | AliITSRecPoint *recp; | |
138 | AliITSgeom *geom = aliITS->GetITSgeom(); | |
139 | index = geom->GetModuleIndex(lycur,Touclad(iriv),Toucdet(iriv)); | |
140 | Int_t lay,lad,det; | |
141 | geom->GetModuleId(index,lay,lad,det); | |
142 | aliITS->ResetRecPoints(); | |
340cd0d8 | 143 | //gAlice->TreeR()->GetEvent(index+1); //first entry in TreeR is empty |
144 | gAlice->TreeR()->GetEvent(index); //first entry in TreeR is empty | |
dce1757b | 145 | |
146 | Int_t npoints=rpoints->GetEntries(); | |
147 | Int_t *indlist=new Int_t[npoints+1]; | |
148 | Int_t counter=0; | |
149 | Int_t ind; | |
150 | for (ind=0; ind<=npoints; ind++) { | |
151 | indlist[ind]=-1; | |
152 | if (*(vettid[index]+ind)==0) { | |
153 | indlist[counter]=ind; | |
154 | counter++; | |
155 | } | |
156 | } | |
157 | ||
158 | ind=-1; | |
cc6e75dd | 159 | |
dce1757b | 160 | for(;;) { |
161 | ind++; | |
162 | if(indlist[ind] < 0) recp=0; | |
163 | else recp = (AliITSRecPoint*)rpoints->UncheckedAt(indlist[ind]); | |
164 | ||
165 | if((!recp) ) break; | |
166 | TVector cluster(3),vecclust(9); | |
167 | vecclust(6)=vecclust(7)=vecclust(8)=-1.; | |
53e7090a | 168 | Double_t sigma[2]; |
169 | //modificata 8-3-2001 | |
dce1757b | 170 | // set veclust in global |
171 | Float_t global[3], local[3]; | |
172 | local[0]=recp->GetX(); | |
173 | local[1]=0.; | |
174 | local[2]= recp->GetZ(); | |
175 | AliITSgeom *g1 = aliITS->GetITSgeom(); | |
176 | Int_t play = lycur; | |
177 | Int_t plad = TMath::Nint(Touclad(iriv)); | |
178 | Int_t pdet = TMath::Nint(Toucdet(iriv)); | |
179 | g1->LtoG(play,plad,pdet,local,global); | |
53e7090a | 180 | |
dce1757b | 181 | vecclust(0)=global[0]; |
182 | vecclust(1)=global[1]; | |
53e7090a | 183 | vecclust(2)=global[2]; |
184 | ||
185 | /* | |
186 | //// modificato 8-3-2001 | |
187 | vecclust(0)=recp->GetRhit();; | |
188 | vecclust(1)=recp->Getphi(); | |
189 | vecclust(2)=recp->GetZglobal(); | |
190 | */ | |
191 | ||
dce1757b | 192 | vecclust(3) = (float)recp->fTracks[0]; |
193 | vecclust(4) = (float)indlist[ind]; | |
194 | vecclust(5) = (float)index; | |
195 | vecclust(6) = (float)recp->fTracks[0]; | |
196 | vecclust(7) = (float)recp->fTracks[1]; | |
197 | vecclust(8) = (float)recp->fTracks[2]; | |
cc6e75dd | 198 | |
199 | sigma[0] = (Double_t) recp->GetSigmaX2(); | |
53e7090a | 200 | sigma[1] = (Double_t) recp->GetSigmaZ2(); |
201 | //commentato 8-3-2001 | |
dce1757b | 202 | //now we are in r,phi,z in global |
203 | cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+vecclust(1)*vecclust(1));//r hit | |
204 | cluster(1) = PhiDef(vecclust(0),vecclust(1)); // phi hit | |
53e7090a | 205 | cluster(2) = vecclust(2); // z hit |
206 | ||
207 | /* | |
208 | //modificato 8-3-2001 | |
209 | cluster(0) = vecclust(0);//r hit | |
210 | cluster(1) = vecclust(1); // phi hit | |
dce1757b | 211 | cluster(2) = vecclust(2); // z hit |
53e7090a | 212 | */ |
cc6e75dd | 213 | // cout<<" layer = "<<play<<"\n"; |
214 | // cout<<" cluster prima = "<<vecclust(0)<<" "<<vecclust(1)<<" " | |
215 | // <<vecclust(2)<<"\n"; getchar(); | |
dce1757b | 216 | //cluster(1)= cluster(1)-trackITS->Getalphaprov(); //provvisorio; |
217 | //if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi(); //provvisorio | |
218 | //cout<<" cluster(1) dopo = "<<cluster(1)<< " alphaprov = "<<trackITS->Getalphaprov()<<"\n"; | |
219 | Float_t sigmatotphi, sigmatotz; | |
220 | ||
221 | //Float_t epsphi=3.2, epsz=3.; | |
8af13b4b | 222 | Float_t epsphi=5.0, epsz=5.0; |
223 | if(Ptref<0.2) {epsphi=3.; epsz=3.;} | |
dce1757b | 224 | |
225 | Double_t Rtrack=(*trackITS).Getrtrack(); | |
226 | Double_t sigmaphi=sigma[0]/(Rtrack*Rtrack); | |
cc6e75dd | 227 | sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi()); |
228 | ||
dce1757b | 229 | sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ()); |
230 | //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n"; | |
cc6e75dd | 231 | //if(vecclust(3)==481) getchar(); |
dce1757b | 232 | if(cluster(1)<6. && (*trackITS).Getphi()>6.) cluster(1)=cluster(1)+(2.*TMath::Pi()); |
233 | if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi()); | |
234 | if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue; | |
235 | // cout<<" supero sigmaphi \n"; | |
236 | AliITStrack *newTrack = new AliITStrack((*trackITS)); | |
237 | (*newTrack).SetLayer((*trackITS).GetLayer()-1); | |
238 | ||
239 | if (TMath::Abs(Rtrack-cluster(0))/Rtrack>1e-6) | |
240 | (*newTrack).Correct(Double_t(cluster(0))); | |
241 | //cout<<" cluster(2) e (*newTrack).GetZ() = "<<cluster(2)<<" "<< (*newTrack).GetZ()<<"\n"; | |
242 | if(TMath::Abs(cluster(2)-(*newTrack).GetZ()) > sigmatotz){ | |
243 | delete newTrack; | |
244 | continue;} | |
245 | ||
246 | if(iriv == 0) flaghit=1; | |
cc6e75dd | 247 | |
8af13b4b | 248 | (*newTrack).AddMS(rl); // add the multiple scattering matrix to the covariance matrix |
249 | (*newTrack).AddEL(rl,1.,0); | |
cc6e75dd | 250 | |
dce1757b | 251 | Double_t sigmanew[2]; |
252 | sigmanew[0]= sigmaphi; | |
253 | sigmanew[1]=sigma[1]; | |
53e7090a | 254 | |
dce1757b | 255 | if(flagvert) |
256 | KalmanFilterVert(newTrack,cluster,sigmanew); | |
257 | else | |
258 | KalmanFilter(newTrack,cluster,sigmanew); | |
cc6e75dd | 259 | |
dce1757b | 260 | |
dce1757b | 261 | (*newTrack).PutCluster(layernew, vecclust); |
262 | newTrack->AddClustInTrack(); | |
263 | ||
264 | listoftrack.AddLast(newTrack); | |
265 | ||
266 | } // end of for(;;) on rec points | |
267 | ||
268 | delete [] indlist; | |
269 | ||
270 | } // end of for on detectors | |
cc6e75dd | 271 | |
dce1757b | 272 | }//end if(outinters==0) |
cc6e75dd | 273 | |
dce1757b | 274 | if(flaghit==0 || outinters==-2) { |
275 | AliITStrack *newTrack = new AliITStrack(*trackITS); | |
276 | (*newTrack).SetLayer((*trackITS).GetLayer()-1); | |
8af13b4b | 277 | (*newTrack).AddMS(rl); // add the multiple scattering matrix to the covariance matrix |
278 | (*newTrack).AddEL(rl,1.,0); | |
dce1757b | 279 | |
280 | listoftrack.AddLast(newTrack); | |
281 | } | |
282 | ||
283 | ||
284 | //gObjectTable->Print(); // stampa memoria | |
53e7090a | 285 | |
8af13b4b | 286 | AliITStracking(&listoftrack, reference, aliITS, rpoints,Ptref,vettid,flagvert,rl); |
dce1757b | 287 | listoftrack.Delete(); |
288 | } // end of for on tracks | |
cc6e75dd | 289 | |
dce1757b | 290 | //gObjectTable->Print(); // stampa memoria |
291 | ||
292 | } | |
293 | ||
294 | ||
295 | Int_t AliITStracking::NewIntersection(AliITStrack &track, Double_t rk,Int_t layer, Int_t &ladder, Int_t &detector) { | |
296 | //Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it | |
297 | // Found the intersection and the detector | |
298 | ||
299 | if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;} | |
300 | track.Propagation(rk); | |
301 | Double_t zinters=track.GetZ(); | |
302 | Double_t phinters=track.Getphi(); | |
303 | //phinters = phinters+track.Getalphaprov(); //provvisorio | |
304 | //if(phinters>2.*3.14) phinters=phinters-2.*3.14; //provvisorio | |
305 | //cout<<"zinters = "<<zinters<<" phinters = "<<phinters<<"\n"; | |
306 | ||
307 | ////////////////////////////////// limits for Geometry 5 ///////////////////////////// | |
308 | ||
309 | Int_t NLadder[]= {20, 40, 14, 22, 34, 38}; | |
8af13b4b | 310 | Int_t NDetector[]= {4, 4, 6, 8, 23, 26}; |
dce1757b | 311 | |
312 | Float_t Detx[]= {0.64, 0.64, 3.509, 3.509, 3.65, 3.65 }; | |
313 | Float_t Detz[]= {4.19, 4.19, 3.75 , 3.75 , 2 , 2 }; | |
314 | //////////////////////////////////////////////////////////////////////////////////////////////// | |
315 | ||
316 | TVector det(9); | |
317 | TVector ListDet(2); | |
318 | TVector DistZCenter(2); | |
319 | AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom(); | |
320 | Int_t iz=0; | |
321 | Double_t epsz=1.2; | |
322 | Double_t epszpixel=0.05; | |
323 | ||
cc6e75dd | 324 | Int_t iD; |
325 | for(iD = 1; iD<= NDetector[layer-1]; iD++) { | |
dce1757b | 326 | g1->GetCenterThetaPhi(layer,1,iD,det); |
327 | Double_t zmin=det(2)-Detz[layer-1]; | |
328 | if(iD==1) zmin=det(2)-(Detz[layer-1])*epsz; | |
329 | Double_t zmax=det(2)+Detz[layer-1]; | |
330 | if(iD==NDetector[layer-1]) zmax=det(2)+(Detz[layer-1])*epsz; | |
3d7c7ec8 | 331 | //added to take into account problem on drift |
99b69e4c | 332 | if(layer == 4 || layer==3) zmin=zmin-epszpixel; zmax=zmax+epszpixel; |
dce1757b | 333 | //cout<<"zmin zinters zmax det(2)= "<<zmin<<" "<<zinters<<" "<<zmax<<" "<<det(2)<<"\n"; |
334 | if(zinters > zmin && zinters <= zmax) { | |
335 | if(iz>1) {cout<< " Errore su iz in NewIntersection \n"; getchar();} | |
336 | else { | |
337 | ListDet(iz)= iD; DistZCenter(iz)=TMath::Abs(zinters-det(2)); iz++; | |
338 | } | |
339 | } | |
340 | } | |
341 | ||
342 | if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;} | |
343 | detector=Int_t (ListDet(0)); | |
344 | if(iz>1 && (DistZCenter(0)>DistZCenter(1))) detector=Int_t (ListDet(1)); | |
345 | ||
346 | AliITSgeom *g2 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom(); | |
347 | Float_t global[3]; | |
348 | Float_t local[3]; | |
349 | TVector ListLad(2); | |
350 | TVector DistphiCenter(2); | |
351 | Int_t ip=0; | |
352 | Double_t pigre=TMath::Pi(); | |
cc6e75dd | 353 | |
354 | Int_t iLd; | |
355 | for(iLd = 1; iLd<= NLadder[layer-1]; iLd++) { | |
dce1757b | 356 | g1->GetCenterThetaPhi(layer,iLd,detector,det); |
357 | Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1))); | |
dce1757b | 358 | // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar(); |
359 | Double_t xmin,ymin,xmax,ymax; | |
3d7c7ec8 | 360 | // Double_t phiconfr=0.0; |
361 | //cout<<" phiconfr inizio = "<<phiconfr <<"\n"; getchar(); | |
dce1757b | 362 | local[1]=local[2]=0.; |
363 | local[0]= -(Detx[layer-1]); | |
364 | if(layer==1) local[0]= (Detx[layer-1]); //take into account different reference system | |
365 | g2->LtoG(layer,iLd,detector,local,global); | |
366 | xmax=global[0]; ymax=global[1]; | |
367 | local[0]= (Detx[layer-1]); | |
368 | if(layer==1) local[0]= -(Detx[layer-1]); //take into account different reference system | |
369 | g2->LtoG(layer,iLd,detector,local,global); | |
370 | xmin=global[0]; ymin=global[1]; | |
371 | Double_t phimin=PhiDef(xmin,ymin); | |
372 | Double_t phimax=PhiDef(xmax,ymax); | |
373 | //cout<<" xmin ymin = "<<xmin<<" "<<ymin<<"\n"; | |
374 | // cout<<" xmax ymax = "<<xmax<<" "<<ymax<<"\n"; | |
3d7c7ec8 | 375 | // cout<<" iLd phimin phimax ="<<iLd<<" "<<phimin<<" "<<phimax<<"\n"; |
376 | ||
377 | Double_t phiconfr=phinters; | |
378 | if(phimin>phimax ){ | |
379 | if(phimin <5.5) {cout<<" Error in NewIntersection for phi \n"; getchar();} | |
dce1757b | 380 | phimin=phimin-(2.*pigre); |
3d7c7ec8 | 381 | if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre); |
dce1757b | 382 | if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre); |
3d7c7ec8 | 383 | } |
384 | // cout<<" phiconfr finale = "<<phiconfr<<"\n"; getchar(); | |
dce1757b | 385 | if(phiconfr>phimin && phiconfr<= phimax) { |
386 | if(ip>1) { | |
387 | cout<< " Errore su ip in NewIntersection \n"; getchar(); | |
388 | } | |
389 | else { | |
390 | ListLad(ip)= iLd; DistphiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++; | |
391 | } | |
392 | } | |
393 | } | |
394 | if(ip==0) { cout<< " No detector along phi \n"; getchar();} | |
395 | ladder=Int_t (ListLad(0)); | |
396 | if(ip>1 && (DistphiCenter(0)>DistphiCenter(1))) ladder=Int_t (ListLad(1)); | |
397 | ||
398 | return 0; | |
399 | } | |
400 | ||
401 | ||
402 | Double_t AliITStracking::PhiDef(Double_t x, Double_t y){ | |
403 | Double_t pigre= TMath::Pi(); | |
6204061f | 404 | Double_t phi=0.0; |
dce1757b | 405 | if(y == 0. || x == 0.) { |
406 | if(y == 0. && x == 0.) { | |
407 | cout << " Error in AliITStracking::PhiDef x=0 and y=0 \n"; getchar(); | |
408 | } | |
409 | if(y==0. && x>0.) phi=0.; | |
410 | if(y==0. && x<0.) phi=pigre; | |
411 | if(x==0 && y>0.) phi=pigre/2.; | |
412 | if(x==0 && y<0.) phi=1.5*pigre; | |
413 | } | |
414 | else { | |
415 | if (x>0. && y>0.) phi=TMath::ATan(y/x); | |
416 | if (x<0. && y>0.) phi=pigre+TMath::ATan(y/x); | |
417 | if (x<0. && y<0.) phi=pigre+TMath::ATan(y/x); | |
418 | if (x>0. && y<0.) phi=(2.*pigre)+TMath::ATan(y/x); | |
419 | } | |
420 | if(phi<0. || phi>(2*pigre)) { | |
421 | cout<<" Error on phi in AliITStracking::PhiDef \n"; getchar(); | |
422 | } | |
423 | return phi; | |
424 | } | |
425 | ||
426 | ||
427 | void AliITStracking::KalmanFilter(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){ | |
428 | //Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it | |
429 | // Kalman filter without vertex constraint | |
430 | ||
dce1757b | 431 | |
432 | ////////////////////////////// Evaluation of the measurement vector ///////////////////////////////////// | |
433 | ||
53e7090a | 434 | Double_t m[2]; |
dce1757b | 435 | Double_t rk,phik,zk; |
436 | rk=cluster(0); phik=cluster(1); zk=cluster(2); | |
53e7090a | 437 | m[0]=phik; m[1]=zk; |
438 | ||
dce1757b | 439 | ///////////////////////////////////// Evaluation of the error matrix V /////////////////////////////// |
440 | ||
53e7090a | 441 | Double_t V00=sigma[0]; |
442 | Double_t V11=sigma[1]; | |
443 | ||
dce1757b | 444 | /////////////////////////////////////////////////////////////////////////////////////////// |
445 | ||
53e7090a | 446 | |
447 | Double_t Cin00,Cin10,Cin20,Cin30,Cin40,Cin11,Cin21,Cin31,Cin41,Cin22,Cin32,Cin42,Cin33,Cin43,Cin44; | |
448 | ||
449 | newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40, | |
450 | Cin41,Cin42,Cin43,Cin44); //get C matrix | |
451 | ||
452 | Double_t Rold00=Cin00+V00; | |
453 | Double_t Rold10=Cin10; | |
454 | Double_t Rold11=Cin11+V11; | |
455 | ||
456 | //////////////////////////////////// R matrix inversion /////////////////////////////////////////////// | |
457 | ||
458 | Double_t det=Rold00*Rold11-Rold10*Rold10; | |
459 | Double_t R00=Rold11/det; | |
460 | Double_t R10=-Rold10/det; | |
461 | Double_t R11=Rold00/det; | |
462 | ||
463 | //////////////////////////////////////////////////////////////////////////////////////////////////////// | |
464 | ||
465 | Double_t K00=Cin00*R00+Cin10*R10; | |
466 | Double_t K01=Cin00*R10+Cin10*R11; | |
467 | Double_t K10=Cin10*R00+Cin11*R10; | |
468 | Double_t K11=Cin10*R10+Cin11*R11; | |
469 | Double_t K20=Cin20*R00+Cin21*R10; | |
470 | Double_t K21=Cin20*R10+Cin21*R11; | |
471 | Double_t K30=Cin30*R00+Cin31*R10; | |
472 | Double_t K31=Cin30*R10+Cin31*R11; | |
473 | Double_t K40=Cin40*R00+Cin41*R10; | |
474 | Double_t K41=Cin40*R10+Cin41*R11; | |
475 | ||
476 | Double_t X0,X1,X2,X3,X4; | |
477 | newTrack->GetXElements(X0,X1,X2,X3,X4); // get the state vector | |
478 | ||
479 | Double_t savex0=X0, savex1=X1; | |
480 | ||
481 | X0+=K00*(m[0]-savex0)+K01*(m[1]-savex1); | |
482 | X1+=K10*(m[0]-savex0)+K11*(m[1]-savex1); | |
483 | X2+=K20*(m[0]-savex0)+K21*(m[1]-savex1); | |
484 | X3+=K30*(m[0]-savex0)+K31*(m[1]-savex1); | |
485 | X4+=K40*(m[0]-savex0)+K41*(m[1]-savex1); | |
486 | ||
487 | Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44; | |
488 | ||
489 | C00=Cin00-K00*Cin00-K01*Cin10; | |
490 | C10=Cin10-K00*Cin10-K01*Cin11; | |
491 | C20=Cin20-K00*Cin20-K01*Cin21; | |
492 | C30=Cin30-K00*Cin30-K01*Cin31; | |
493 | C40=Cin40-K00*Cin40-K01*Cin41; | |
494 | ||
495 | C11=Cin11-K10*Cin10-K11*Cin11; | |
496 | C21=Cin21-K10*Cin20-K11*Cin21; | |
497 | C31=Cin31-K10*Cin30-K11*Cin31; | |
498 | C41=Cin41-K10*Cin40-K11*Cin41; | |
499 | ||
500 | C22=Cin22-K20*Cin20-K21*Cin21; | |
501 | C32=Cin32-K20*Cin30-K21*Cin31; | |
502 | C42=Cin42-K20*Cin40-K21*Cin41; | |
503 | ||
504 | C33=Cin33-K30*Cin30-K31*Cin31; | |
505 | C43=Cin43-K30*Cin40-K31*Cin41; | |
506 | ||
507 | C44=Cin44-K40*Cin40-K41*Cin41; | |
508 | ||
509 | newTrack->PutXElements(X0,X1,X2,X3,X4); // put the new state vector | |
510 | ||
511 | newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the | |
512 | // new cov matrix | |
513 | Double_t VMCold00=V00-C00; | |
514 | Double_t VMCold10=-C10; | |
515 | Double_t VMCold11=V11-C11; | |
516 | ||
517 | ///////////////////////////////////// Matrix VMC inversion //////////////////////////////////////////////// | |
518 | ||
519 | det=VMCold00*VMCold11-VMCold10*VMCold10; | |
520 | Double_t VMC00=VMCold11/det; | |
521 | Double_t VMC10=-VMCold10/det; | |
522 | Double_t VMC11=VMCold00/det; | |
523 | ||
524 | //////////////////////////////////////////////////////////////////////////////////////////////////////////// | |
525 | ||
526 | Double_t chi2=(m[0]-X0)*( VMC00*(m[0]-X0) + 2.*VMC10*(m[1]-X1) ) + | |
527 | (m[1]-X1)*VMC11*(m[1]-X1); | |
dce1757b | 528 | |
529 | newTrack->SetChi2(newTrack->GetChi2()+chi2); | |
530 | ||
531 | } | |
532 | ||
533 | ||
534 | void AliITStracking::KalmanFilterVert(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){ | |
535 | //Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it | |
536 | // Kalman filter with vertex constraint | |
dce1757b | 537 | |
53e7090a | 538 | ////////////////////////////// Evaluation of the measurement vector m /////////////// |
dce1757b | 539 | |
53e7090a | 540 | Double_t m[4]; |
dce1757b | 541 | Double_t rk,phik,zk; |
542 | rk=cluster(0); phik=cluster(1); zk=cluster(2); | |
53e7090a | 543 | m[0]=phik; m[1]=zk; |
dce1757b | 544 | |
545 | Double_t CC=(*newTrack).GetC(); | |
546 | Double_t Zv=(*newTrack).GetZv(); | |
547 | Double_t Dv=(*newTrack).GetDv(); | |
dce1757b | 548 | Double_t Cy=CC/2.; |
549 | Double_t tgl= (zk-Zv)*Cy/TMath::ASin(Cy*rk); | |
53e7090a | 550 | m[2]=Dv; m[3]=tgl; |
dce1757b | 551 | |
53e7090a | 552 | ///////////////////////////////////// Evaluation of the error matrix V ////////////// |
dce1757b | 553 | Int_t Layer=newTrack->GetLayer(); |
53e7090a | 554 | Double_t V00=sigma[0]; |
555 | Double_t V11=sigma[1]; | |
556 | Double_t V31=sigma[1]/rk; | |
dce1757b | 557 | Double_t sigmaDv=newTrack->GetsigmaDv(); |
53e7090a | 558 | Double_t V22=sigmaDv*sigmaDv + newTrack->Getd2(Layer-1); |
559 | Double_t V32=newTrack->Getdtgl(Layer-1); | |
dce1757b | 560 | Double_t sigmaZv=newTrack->GetsigmaZv(); |
53e7090a | 561 | Double_t V33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1); |
562 | /////////////////////////////////////////////////////////////////////////////////////// | |
563 | ||
564 | Double_t Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,Cin41,Cin42,Cin43,Cin44; | |
565 | ||
566 | newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40, | |
567 | Cin41,Cin42,Cin43,Cin44); //get C matrix | |
568 | ||
569 | Double_t R[4][4]; | |
570 | R[0][0]=Cin00+V00; | |
571 | R[1][0]=Cin10; | |
572 | R[2][0]=Cin20; | |
573 | R[3][0]=Cin30; | |
574 | R[1][1]=Cin11+V11; | |
575 | R[2][1]=Cin21; | |
576 | R[3][1]=Cin31+sigma[1]/rk; | |
577 | R[2][2]=Cin22+sigmaDv*sigmaDv+newTrack->Getd2(Layer-1); | |
578 | R[3][2]=Cin32+newTrack->Getdtgl(Layer-1); | |
579 | R[3][3]=Cin33+(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1); | |
580 | ||
581 | 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]; | |
582 | R[2][3]=R[3][2]; | |
583 | ||
584 | ///////////////////// Matrix R inversion //////////////////////////////////////////// | |
585 | ||
586 | Int_t n=4; | |
587 | Double_t big, hold; | |
588 | Double_t d=1.; | |
589 | Int_t ll[n],mm[n]; | |
590 | ||
591 | for(Int_t k=0; k<n; k++) { | |
592 | ll[k]=k; | |
593 | mm[k]=k; | |
594 | big=R[k][k]; | |
595 | for(Int_t j=k; j<n ; j++) { | |
596 | for (Int_t i=j; i<n; i++) { | |
597 | if(TMath::Abs(big) < TMath::Abs(R[i][j]) ) { big=R[i][j]; ll[k]=i; mm[k]=j; } | |
598 | } | |
599 | } | |
600 | // | |
601 | Int_t j= ll[k]; | |
602 | if(j > k) { | |
603 | for(Int_t i=0; i<n; i++) { hold=-R[k][i]; R[k][i]=R[j][i]; R[j][i]=hold; } | |
604 | ||
605 | } | |
606 | // | |
607 | Int_t i=mm[k]; | |
608 | if(i > k ) { | |
609 | for(int j=0; j<n; j++) { hold=-R[j][k]; R[j][k]=R[j][i]; R[j][i]=hold; } | |
610 | } | |
611 | // | |
612 | if(!big) { | |
613 | d=0.; | |
614 | cout << "Singular matrix\n"; | |
615 | } | |
616 | for(Int_t i=0; i<n; i++) { | |
617 | if(i == k) { continue; } | |
618 | R[i][k]=R[i][k]/(-big); | |
619 | } | |
620 | // | |
621 | for(Int_t i=0; i<n; i++) { | |
622 | hold=R[i][k]; | |
623 | for(Int_t j=0; j<n; j++) { | |
624 | if(i == k || j == k) { continue; } | |
625 | R[i][j]=hold*R[k][j]+R[i][j]; | |
626 | } | |
627 | } | |
628 | // | |
629 | for(Int_t j=0; j<n; j++) { | |
630 | if(j == k) { continue; } | |
631 | R[k][j]=R[k][j]/big; | |
632 | } | |
633 | // | |
634 | d=d*big; | |
635 | // | |
636 | R[k][k]=1./big; | |
637 | } | |
638 | // | |
639 | for(Int_t k=n-1; k>=0; k--) { | |
640 | Int_t i=ll[k]; | |
641 | if(i > k) { | |
642 | for (Int_t j=0; j<n; j++) {hold=R[j][k]; R[j][k]=-R[j][i]; R[j][i]=hold;} | |
643 | } | |
644 | Int_t j=mm[k]; | |
645 | if(j > k) { | |
646 | for (Int_t i=0; i<n; i++) {hold=R[k][i]; R[k][i]=-R[j][i]; R[j][i]=hold;} | |
647 | } | |
648 | } | |
649 | ////////////////////////////////////////////////////////////////////////////////// | |
650 | ||
651 | ||
652 | Double_t K00=Cin00*R[0][0]+Cin10*R[1][0]+Cin20*R[2][0]+Cin30*R[3][0]; | |
653 | Double_t K01=Cin00*R[1][0]+Cin10*R[1][1]+Cin20*R[2][1]+Cin30*R[3][1]; | |
654 | Double_t K02=Cin00*R[2][0]+Cin10*R[2][1]+Cin20*R[2][2]+Cin30*R[3][2]; | |
655 | Double_t K03=Cin00*R[3][0]+Cin10*R[3][1]+Cin20*R[3][2]+Cin30*R[3][3]; | |
656 | Double_t K10=Cin10*R[0][0]+Cin11*R[1][0]+Cin21*R[2][0]+Cin31*R[3][0]; | |
657 | Double_t K11=Cin10*R[1][0]+Cin11*R[1][1]+Cin21*R[2][1]+Cin31*R[3][1]; | |
658 | Double_t K12=Cin10*R[2][0]+Cin11*R[2][1]+Cin21*R[2][2]+Cin31*R[3][2]; | |
659 | Double_t K13=Cin10*R[3][0]+Cin11*R[3][1]+Cin21*R[3][2]+Cin31*R[3][3]; | |
660 | Double_t K20=Cin20*R[0][0]+Cin21*R[1][0]+Cin22*R[2][0]+Cin32*R[3][0]; | |
661 | Double_t K21=Cin20*R[1][0]+Cin21*R[1][1]+Cin22*R[2][1]+Cin32*R[3][1]; | |
662 | Double_t K22=Cin20*R[2][0]+Cin21*R[2][1]+Cin22*R[2][2]+Cin32*R[3][2]; | |
663 | Double_t K23=Cin20*R[3][0]+Cin21*R[3][1]+Cin22*R[3][2]+Cin32*R[3][3]; | |
664 | Double_t K30=Cin30*R[0][0]+Cin31*R[1][0]+Cin32*R[2][0]+Cin33*R[3][0]; | |
665 | Double_t K31=Cin30*R[1][0]+Cin31*R[1][1]+Cin32*R[2][1]+Cin33*R[3][1]; | |
666 | Double_t K32=Cin30*R[2][0]+Cin31*R[2][1]+Cin32*R[2][2]+Cin33*R[3][2]; | |
667 | Double_t K33=Cin30*R[3][0]+Cin31*R[3][1]+Cin32*R[3][2]+Cin33*R[3][3]; | |
668 | Double_t K40=Cin40*R[0][0]+Cin41*R[1][0]+Cin42*R[2][0]+Cin43*R[3][0]; | |
669 | Double_t K41=Cin40*R[1][0]+Cin41*R[1][1]+Cin42*R[2][1]+Cin43*R[3][1]; | |
670 | Double_t K42=Cin40*R[2][0]+Cin41*R[2][1]+Cin42*R[2][2]+Cin43*R[3][2]; | |
671 | Double_t K43=Cin40*R[3][0]+Cin41*R[3][1]+Cin42*R[3][2]+Cin43*R[3][3]; | |
672 | ||
673 | Double_t X0,X1,X2,X3,X4; | |
674 | newTrack->GetXElements(X0,X1,X2,X3,X4); // get the state vector | |
675 | ||
676 | Double_t savex0=X0, savex1=X1, savex2=X2, savex3=X3; | |
677 | ||
678 | X0+=K00*(m[0]-savex0)+K01*(m[1]-savex1)+K02*(m[2]-savex2)+ | |
679 | K03*(m[3]-savex3); | |
680 | X1+=K10*(m[0]-savex0)+K11*(m[1]-savex1)+K12*(m[2]-savex2)+ | |
681 | K13*(m[3]-savex3); | |
682 | X2+=K20*(m[0]-savex0)+K21*(m[1]-savex1)+K22*(m[2]-savex2)+ | |
683 | K23*(m[3]-savex3); | |
684 | X3+=K30*(m[0]-savex0)+K31*(m[1]-savex1)+K32*(m[2]-savex2)+ | |
685 | K33*(m[3]-savex3); | |
686 | X4+=K40*(m[0]-savex0)+K41*(m[1]-savex1)+K42*(m[2]-savex2)+ | |
687 | K43*(m[3]-savex3); | |
688 | ||
689 | Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44; | |
690 | ||
691 | C00=Cin00-K00*Cin00-K01*Cin10-K02*Cin20-K03*Cin30; | |
692 | C10=Cin10-K00*Cin10-K01*Cin11-K02*Cin21-K03*Cin31; | |
693 | C20=Cin20-K00*Cin20-K01*Cin21-K02*Cin22-K03*Cin32; | |
694 | C30=Cin30-K00*Cin30-K01*Cin31-K02*Cin32-K03*Cin33; | |
695 | C40=Cin40-K00*Cin40-K01*Cin41-K02*Cin42-K03*Cin43; | |
696 | ||
697 | C11=Cin11-K10*Cin10-K11*Cin11-K12*Cin21-K13*Cin31; | |
698 | C21=Cin21-K10*Cin20-K11*Cin21-K12*Cin22-K13*Cin32; | |
699 | C31=Cin31-K10*Cin30-K11*Cin31-K12*Cin32-K13*Cin33; | |
700 | C41=Cin41-K10*Cin40-K11*Cin41-K12*Cin42-K13*Cin43; | |
701 | ||
702 | C22=Cin22-K20*Cin20-K21*Cin21-K22*Cin22-K23*Cin32; | |
703 | C32=Cin32-K20*Cin30-K21*Cin31-K22*Cin32-K23*Cin33; | |
704 | C42=Cin42-K20*Cin40-K21*Cin41-K22*Cin42-K23*Cin43; | |
705 | ||
706 | C33=Cin33-K30*Cin30-K31*Cin31-K32*Cin32-K33*Cin33; | |
707 | C43=Cin43-K30*Cin40-K31*Cin41-K32*Cin42-K33*Cin43; | |
708 | ||
709 | C44=Cin44-K40*Cin40-K41*Cin41-K42*Cin42-K43*Cin43; | |
710 | ||
711 | newTrack->PutXElements(X0,X1,X2,X3,X4); // put the new state vector | |
712 | ||
713 | newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the | |
714 | // new cov matrix | |
715 | ||
716 | Double_t VMC[4][4]; | |
717 | ||
718 | VMC[0][0]=V00-C00; VMC[1][0]=-C10; VMC[2][0]=-C20; VMC[3][0]=-C30; | |
719 | VMC[1][1]=V11-C11; VMC[2][1]=-C21; VMC[3][1]=V31-C31; | |
720 | VMC[2][2]=V22-C22; VMC[3][2]=V32-C32; | |
721 | VMC[3][3]=V33-C33; | |
722 | ||
723 | VMC[0][1]=VMC[1][0]; VMC[0][2]=VMC[2][0]; VMC[0][3]=VMC[3][0]; | |
724 | VMC[1][2]=VMC[2][1]; VMC[1][3]=VMC[3][1]; | |
725 | VMC[2][3]=VMC[3][2]; | |
726 | ||
727 | ||
728 | /////////////////////// VMC matrix inversion /////////////////////////////////// | |
dce1757b | 729 | |
53e7090a | 730 | n=4; |
731 | d=1.; | |
732 | ||
733 | for(Int_t k=0; k<n; k++) { | |
734 | ll[k]=k; | |
735 | mm[k]=k; | |
736 | big=VMC[k][k]; | |
737 | for(Int_t j=k; j<n ; j++) { | |
738 | for (Int_t i=j; i<n; i++) { | |
739 | if(TMath::Abs(big) < TMath::Abs(VMC[i][j]) ) { big=VMC[i][j]; ll[k]=i; mm[k]=j; } | |
740 | } | |
741 | } | |
742 | // | |
743 | Int_t j= ll[k]; | |
744 | if(j > k) { | |
745 | for(Int_t i=0; i<n; i++) { hold=-VMC[k][i]; VMC[k][i]=VMC[j][i]; VMC[j][i]=hold; } | |
746 | ||
747 | } | |
748 | // | |
749 | Int_t i=mm[k]; | |
750 | if(i > k ) { | |
751 | for(int j=0; j<n; j++) { hold=-VMC[j][k]; VMC[j][k]=VMC[j][i]; VMC[j][i]=hold; } | |
752 | } | |
753 | // | |
754 | if(!big) { | |
755 | d=0.; | |
756 | cout << "Singular matrix\n"; | |
757 | } | |
758 | for(Int_t i=0; i<n; i++) { | |
759 | if(i == k) { continue; } | |
760 | VMC[i][k]=VMC[i][k]/(-big); | |
761 | } | |
762 | // | |
763 | for(Int_t i=0; i<n; i++) { | |
764 | hold=VMC[i][k]; | |
765 | for(Int_t j=0; j<n; j++) { | |
766 | if(i == k || j == k) { continue; } | |
767 | VMC[i][j]=hold*VMC[k][j]+VMC[i][j]; | |
768 | } | |
769 | } | |
770 | // | |
771 | for(Int_t j=0; j<n; j++) { | |
772 | if(j == k) { continue; } | |
773 | VMC[k][j]=VMC[k][j]/big; | |
774 | } | |
775 | // | |
776 | d=d*big; | |
777 | // | |
778 | VMC[k][k]=1./big; | |
779 | } | |
780 | // | |
781 | for(Int_t k=n-1; k>=0; k--) { | |
782 | Int_t i=ll[k]; | |
783 | if(i > k) { | |
784 | for (Int_t j=0; j<n; j++) {hold=VMC[j][k]; VMC[j][k]=-VMC[j][i]; VMC[j][i]=hold;} | |
785 | } | |
786 | Int_t j=mm[k]; | |
787 | if(j > k) { | |
788 | for (Int_t i=0; i<n; i++) {hold=VMC[k][i]; VMC[k][i]=-VMC[j][i]; VMC[j][i]=hold;} | |
789 | } | |
790 | } | |
791 | ||
792 | ||
793 | //////////////////////////////////////////////////////////////////////////////// | |
794 | ||
795 | Double_t chi2=(m[0]-X0)*( VMC[0][0]*(m[0]-X0) + 2.*VMC[1][0]*(m[1]-X1) + | |
796 | 2.*VMC[2][0]*(m[2]-X2)+ 2.*VMC[3][0]*(m[3]-X3) ) + | |
797 | (m[1]-X1)* ( VMC[1][1]*(m[1]-X1) + 2.*VMC[2][1]*(m[2]-X2)+ | |
798 | 2.*VMC[3][1]*(m[3]-X3) ) + | |
799 | (m[2]-X2)* ( VMC[2][2]*(m[2]-X2)+ 2.*VMC[3][2]*(m[3]-X3) ) + | |
800 | (m[3]-X3)*VMC[3][3]*(m[3]-X3); | |
dce1757b | 801 | |
802 | newTrack->SetChi2(newTrack->GetChi2()+chi2); | |
803 | ||
804 | } | |
805 |