]>
Commit | Line | Data |
---|---|---|
1f30d65c | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | // | |
15d37333 | 16 | // Static function member which can be used in standalone cases |
17 | // especially as utils for AliTPCCosmicTrackfit | |
18 | // | |
19 | // detailed description can be found inside individual function | |
1f30d65c | 20 | // |
3aa6a136 | 21 | // grep "error" in output log to check abnormal termination |
1f30d65c | 22 | // |
15d37333 | 23 | // |
24 | // Xianguo Lu | |
25 | // lu@physi.uni-heidelberg.de | |
26 | // Xianguo.Lu@cern.ch | |
27 | ||
1f30d65c | 28 | // |
29 | ||
30 | #include <TAxis.h> | |
31 | #include <TCanvas.h> | |
32 | #include <TGraph.h> | |
33 | #include <TTreeStream.h> | |
34 | #include <TVector3.h> | |
35 | ||
3aa6a136 | 36 | #include "AliAnalysisManager.h" |
37 | #include "AliESDEvent.h" | |
38 | #include "AliESDfriend.h" | |
1f30d65c | 39 | #include "AliESDtrack.h" |
40 | #include "AliESDfriendTrack.h" | |
3aa6a136 | 41 | #include "AliESDInputHandler.h" |
1f30d65c | 42 | #include "AliTPCseed.h" |
43 | #include "AliTrackerBase.h" | |
44 | #include "AliTrackPointArray.h" | |
45 | ||
46 | #include "AliTPCCosmicUtils.h" | |
47 | ||
3aa6a136 | 48 | Int_t AliTPCCosmicUtils::GetBField(const AliESDEvent *esd) |
49 | { | |
50 | // | |
51 | //get b-field in unit of kg, expect 1, 2, 5kg only | |
52 | // | |
53 | const Double_t dm = esd->GetMagneticField(); | |
54 | ||
55 | Int_t mag = -999; | |
56 | if(dm>0) | |
57 | mag = (Int_t) (dm+0.5); | |
58 | else | |
59 | mag = (Int_t) (dm-0.5); | |
60 | ||
61 | if(TMath::Abs(mag)!=1 && TMath::Abs(mag)!=2 && TMath::Abs(mag)!=5){ | |
62 | printf("error AliTPCCosmicUtils::GetB strange Bfield! %f %d\n", dm, mag); exit(1); | |
63 | } | |
64 | ||
65 | return mag; | |
66 | } | |
67 | ||
68 | Int_t AliTPCCosmicUtils::GetTrigger(const AliESDEvent *esd) | |
69 | { | |
70 | // | |
71 | //get cosmic trigger | |
72 | // | |
73 | const TString strig = esd->GetFiredTriggerClasses(); | |
74 | Int_t btrig = 0; | |
75 | ||
76 | if( strig.Contains("C0OB0")) | |
77 | btrig += k0OB0; | |
78 | if( strig.Contains("C0OB1")) | |
79 | btrig += k0OB1; | |
80 | if( strig.Contains("C0HWU")) | |
81 | btrig += k0HWU; | |
82 | if( strig.Contains("CTRDCO2")) | |
83 | btrig += kTRDCO2; | |
84 | if( strig.Contains("AMU")) | |
85 | btrig += kAMU; | |
86 | if( strig.Contains("SCO")) | |
87 | btrig += kSCO; | |
88 | ||
89 | return btrig; | |
90 | } | |
91 | ||
92 | Bool_t AliTPCCosmicUtils::GetESD(AliESDEvent *& esdevent, AliESDfriend *& esdfriend) | |
93 | { | |
94 | // | |
95 | //get esdevent and esdfriend | |
96 | // | |
97 | ||
98 | esdevent = 0x0; | |
99 | esdfriend = 0x0; | |
100 | ||
101 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
102 | if(esdH){ | |
103 | esdevent = esdH->GetEvent(); | |
104 | } | |
105 | ||
106 | if(!esdevent){ | |
107 | printf("AliAnalysisTaskdEdxCosmic::GetESD error No ESD Event\n"); | |
108 | return kFALSE; | |
109 | } | |
110 | ||
111 | esdH->SetActiveBranches("ESDfriend*"); | |
112 | esdfriend = (AliESDfriend *)esdevent->FindListObject("AliESDfriend"); | |
113 | ||
114 | if(!esdfriend){ | |
115 | printf("AliAnalysisTaskdEdxCosmic::GetESD error No ESD friend\n"); | |
116 | return kFALSE; | |
117 | } | |
118 | ||
119 | esdevent->SetESDfriend(esdfriend); | |
120 | ||
121 | return kTRUE; | |
122 | } | |
123 | ||
124 | ||
125 | Double_t AliTPCCosmicUtils::GetMinPhi(const AliExternalTrackParam *params[]) | |
126 | { | |
127 | // | |
128 | //minimum phi angle of the two tracks, 0: horizontal | |
129 | // | |
130 | const Double_t fsp[] = {TMath::Abs(TMath::Sin(params[0]->Phi())), | |
131 | TMath::Abs(TMath::Sin(params[1]->Phi()))}; | |
132 | const Double_t minphi = TMath::ASin(TMath::Min(fsp[0], fsp[1])) * TMath::RadToDeg(); | |
133 | ||
134 | return minphi; | |
135 | } | |
136 | ||
137 | Double_t AliTPCCosmicUtils::Point2LineDist(const TVector3 p0, const TVector3 l1, const TVector3 l2, TVector3 *vtx) | |
15d37333 | 138 | { |
139 | // | |
140 | //return distance of p0 to line (l2-l1) | |
141 | // | |
142 | ||
143 | const TVector3 va = p0 - l1; | |
144 | const TVector3 vb = l2 - l1; | |
145 | ||
3aa6a136 | 146 | //calculate vtx position |
147 | const Double_t ca = va.Dot(vb)/vb.Mag(); | |
148 | const TVector3 imp = l1 + vb.Unit()*ca; | |
15d37333 | 149 | |
3aa6a136 | 150 | if(vtx){ |
151 | (*vtx) = imp; | |
152 | } | |
153 | ||
154 | /* equivalent to imp.Mag() | |
155 | //calculate distance | |
156 | const TVector3 dd = va.Cross(vb); | |
15d37333 | 157 | return dd.Mag()/vb.Mag(); |
3aa6a136 | 158 | */ |
159 | ||
160 | return imp.Mag(); | |
15d37333 | 161 | } |
162 | ||
163 | Double_t AliTPCCosmicUtils::AngleInRange(Double_t phi) | |
164 | { | |
165 | // | |
166 | //get the phi value (in rad) in -pi ~ pi, so that fabs() works naiively | |
167 | // | |
168 | const Double_t pmin = -TMath::Pi(); | |
169 | const Double_t pmax = pmin+TMath::TwoPi(); | |
170 | ||
171 | while(phi<pmin) | |
172 | phi+= TMath::TwoPi(); | |
173 | ||
174 | while(phi>=pmax) | |
175 | phi-= TMath::TwoPi(); | |
176 | ||
177 | return phi; | |
178 | } | |
179 | ||
180 | Bool_t AliTPCCosmicUtils::RotateSafe(AliExternalTrackParam *trackPar, const Double_t rawalpha) | |
1f30d65c | 181 | { |
182 | //1. in AliExternalTrackParam::GetXYZ | |
183 | //r[0]=fX; r[1]=fP[0]; r[2]=fP[1]; | |
184 | //return Local2GlobalPosition(r,fAlpha); | |
185 | //in AliVParticle::Local2GlobalMomentum | |
186 | //Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha); | |
187 | //Double_t r=TMath::Sqrt((1. - p[1])*(1. + p[1])); | |
188 | //p[0]=pt*(r*cs - p[1]*sn); p[1]=pt*(p[1]*cs + r*sn); p[2]=pt*p[2]; | |
189 | //which is cos(phi_local=Snp+fAlpha) > 0 always. | |
190 | //2. also in Bool_t AliExternalTrackParam::Rotate(Double_t alpha) | |
191 | // Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); | |
192 | //in AliExternalTrackParam::Set | |
193 | //mom.RotateZ(-fAlpha); | |
194 | //fP[2] = TMath::Sin(mom.Phi()); | |
195 | //since only sin is used for mom.Phi(), that is assuming cos(mom.Phi())>0 | |
196 | ||
15d37333 | 197 | //the range-changing in AliExternalTrackParam::Rotate not precise enough |
198 | const Double_t aa = AngleInRange(rawalpha); | |
199 | ||
200 | const Double_t a0 = trackPar->GetAlpha(); | |
201 | const Double_t p2 = trackPar->GetParameter()[2]; | |
202 | ||
203 | //copied from AliExternalTrackParam::Rotate | |
204 | const Double_t ca=TMath::Cos(aa-a0), sa=TMath::Sin(aa-a0); | |
205 | const Double_t sf=p2, cf=TMath::Sqrt((1.- p2)*(1.+p2)); | |
206 | ||
207 | if((cf*ca+sf*sa)<0) { | |
208 | return kFALSE; | |
209 | } | |
210 | ||
211 | return trackPar->Rotate(aa); | |
212 | ||
213 | /* | |
1f30d65c | 214 | Double_t xyz[3], pxpypz[3]; |
215 | trackPar->GetXYZ(xyz); | |
216 | trackPar->GetPxPyPz(pxpypz); | |
217 | const TVector3 pos0(xyz); | |
218 | const TVector3 mom0(pxpypz); | |
219 | TVector3 pos1, mom1, dpos, dmom; | |
220 | const Double_t eps = 1e-6; | |
221 | AliExternalTrackParam tmppar = (*trackPar); | |
222 | ||
223 | if(!tmppar.Rotate(aa)){ | |
224 | return 0; | |
225 | } | |
226 | ||
227 | tmppar.GetXYZ(xyz); | |
228 | tmppar.GetPxPyPz(pxpypz); | |
229 | pos1.SetXYZ(xyz[0], xyz[1], xyz[2]); | |
230 | mom1.SetXYZ(pxpypz[0], pxpypz[1], pxpypz[2]); | |
231 | dpos = pos1-pos0; | |
232 | dmom = mom1-mom0; | |
233 | if(dpos.Mag()<eps && dmom.Mag()<eps){ | |
234 | (*trackPar)=tmppar; | |
235 | return 1; | |
236 | } | |
237 | else | |
238 | return 0; | |
15d37333 | 239 | */ |
1f30d65c | 240 | |
241 | /* | |
242 | //flip | |
243 | tmppar = (*trackPar); | |
244 | if(!tmppar.Rotate(aa+TMath::Pi())){ | |
245 | return 0; | |
246 | } | |
247 | ||
248 | tmppar.GetXYZ(xyz); | |
249 | tmppar.GetPxPyPz(pxpypz); | |
250 | pos1.SetXYZ(xyz[0], xyz[1], xyz[2]); | |
251 | mom1.SetXYZ(pxpypz[0], pxpypz[1], pxpypz[2]); | |
252 | dpos = pos1-pos0; | |
253 | dmom = mom1-mom0; | |
254 | if(dpos.Mag()<eps && dmom.Mag()<eps){ | |
255 | (*trackPar)=tmppar; | |
256 | return -1; | |
257 | } | |
258 | ||
259 | return 0; | |
260 | */ | |
261 | } | |
262 | ||
263 | void AliTPCCosmicUtils::PrintTrackParam(const Int_t id, const AliExternalTrackParam * trackpar, const char *tag) | |
264 | { | |
265 | // | |
266 | //print out TrackParam | |
267 | // | |
268 | Double_t xyz[3]={-999,-999,-999}; | |
269 | trackpar->GetXYZ(xyz); | |
270 | Double_t pxpypz[3]={-999,-999,-999}; | |
271 | trackpar->GetPxPyPz(pxpypz); | |
272 | ||
15d37333 | 273 | printf("PrintTrackPar %2s %3d : %6.1f %5.1f : %6.6f %6.6f %6.6f : %6.6f : %8.3f %8.3f %8.3f : (%2d) %8.3f %11.6f %.2f\n", tag, id, trackpar->GetX(), sqrt(pow(xyz[0],2)+pow(xyz[1],2)), xyz[0], xyz[1], xyz[2], trackpar->Phi(), pxpypz[0], pxpypz[1], pxpypz[2], trackpar->Charge(), trackpar->Pt(), trackpar->P(), trackpar->GetAlpha()*TMath::RadToDeg()); |
274 | } | |
275 | ||
276 | void AliTPCCosmicUtils::DrawTracks(AliESDtrack *esdtrks[], const TString tag, const TString outputformat) | |
277 | { | |
278 | // | |
279 | //draw esdtracks | |
280 | // | |
281 | const AliTPCseed *seeds[]={GetTPCseed(esdtrks[0]), GetTPCseed(esdtrks[1])}; | |
282 | DrawSeeds(seeds, tag, outputformat); | |
1f30d65c | 283 | } |
284 | ||
15d37333 | 285 | |
1f30d65c | 286 | void AliTPCCosmicUtils::DrawSeeds(const AliTPCseed * seeds[], const TString tag, const TString outputformat) |
287 | { | |
288 | // | |
289 | //draw seed and output to file | |
290 | // | |
15d37333 | 291 | |
292 | if(!seeds[0] || !seeds[1]) | |
293 | return; | |
294 | ||
1f30d65c | 295 | TGraph *grsyx[]={new TGraph, new TGraph}; |
296 | TGraph *grsyz[]={new TGraph, new TGraph}; | |
297 | ||
298 | for(Int_t itrk=0; itrk<2; itrk++){ | |
299 | if(!seeds[itrk]) | |
300 | continue; | |
301 | ||
302 | Int_t ncl = 0; | |
3aa6a136 | 303 | for(Int_t irow=0; irow<NRow(); irow++){ |
1f30d65c | 304 | const AliTPCclusterMI * cls = seeds[itrk]->GetClusterPointer(irow); |
305 | if(!cls) | |
306 | continue; | |
307 | ||
308 | Float_t xyz[3]={-999,-999,-999}; | |
309 | cls->GetGlobalXYZ(xyz); | |
310 | //printf("Test %f %f %f \n", xyz[0], xyz[1], xyz[2]); | |
311 | grsyx[itrk]->SetPoint(ncl, xyz[0], xyz[1]); | |
312 | grsyz[itrk]->SetPoint(ncl, xyz[2], xyz[1]); | |
313 | ncl++; | |
314 | } | |
315 | } | |
316 | ||
317 | grsyx[0]->SetTitle(tag); | |
318 | grsyx[0]->SetMaximum(250); | |
319 | grsyx[0]->SetMinimum(-250); | |
320 | grsyx[0]->GetXaxis()->SetLimits(-250,250); | |
321 | grsyx[0]->GetXaxis()->SetTitle("X (cm)"); | |
322 | grsyx[0]->GetYaxis()->SetTitle("Y (cm)"); | |
323 | ||
324 | grsyx[0]->SetMarkerStyle(20); | |
325 | grsyx[0]->SetMarkerColor(kRed); | |
326 | grsyx[1]->SetMarkerStyle(22); | |
327 | grsyx[1]->SetMarkerColor(kBlue); | |
328 | ||
329 | grsyz[0]->SetTitle(tag); | |
330 | grsyz[0]->SetMaximum(250); | |
331 | grsyz[0]->SetMinimum(-250); | |
332 | grsyz[0]->GetXaxis()->SetLimits(-250,250); | |
333 | grsyz[0]->GetXaxis()->SetTitle("Z (cm)"); | |
334 | grsyz[0]->GetYaxis()->SetTitle("Y (cm)"); | |
335 | ||
336 | grsyz[0]->SetMarkerStyle(20); | |
337 | grsyz[0]->SetMarkerColor(kRed); | |
338 | grsyz[1]->SetMarkerStyle(22); | |
339 | grsyz[1]->SetMarkerColor(kBlue); | |
340 | ||
341 | TCanvas *cc=new TCanvas("cc","",1200,600); | |
342 | cc->Divide(2); | |
343 | cc->cd(1); | |
344 | for(int itrk=0; itrk<2; itrk++){ | |
345 | grsyx[itrk]->SetMarkerSize(1); | |
346 | grsyx[itrk]->Draw(itrk?"lp same":"alp"); | |
347 | } | |
348 | cc->cd(2); | |
349 | for(int itrk=0; itrk<2; itrk++){ | |
350 | grsyz[itrk]->SetMarkerSize(1); | |
351 | grsyz[itrk]->Draw(itrk?"lp same":"alp"); | |
352 | } | |
15d37333 | 353 | |
354 | gErrorIgnoreLevel = 1001; | |
1f30d65c | 355 | cc->Print(Form("drawTrack%s.%s", tag.Data(), outputformat.Data())); |
356 | ||
357 | for(Int_t ii=0; ii<2; ii++){ | |
358 | delete grsyx[ii]; | |
359 | delete grsyz[ii]; | |
360 | } | |
361 | delete cc; | |
362 | } | |
363 | ||
364 | //================================================================================================ | |
365 | AliTPCseed * AliTPCCosmicUtils::GetTPCseed(const AliESDtrack *esdtrack) | |
366 | { | |
367 | // | |
368 | //Get TPC seeds from ESDfriendTrack | |
369 | // | |
370 | ||
371 | AliESDfriendTrack * friendtrk = (AliESDfriendTrack *) esdtrack->GetFriendTrack(); | |
372 | if(!friendtrk){ | |
3aa6a136 | 373 | //printf("AliTPCCosmicUtils::GetTPCseed no friend track!\n"); exit(1); |
374 | return 0x0; | |
1f30d65c | 375 | } |
3aa6a136 | 376 | |
1f30d65c | 377 | TObject *calibObject=0x0; |
378 | AliTPCseed *tseed = 0x0; | |
379 | for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) { | |
380 | if( (tseed=dynamic_cast<AliTPCseed*>(calibObject)) ) | |
381 | break; | |
382 | } | |
383 | return tseed; | |
384 | } | |
385 | ||
386 | //================================================================================================ | |
387 | AliExternalTrackParam *AliTPCCosmicUtils::MakeSeed(const AliTPCseed *tseed) | |
388 | { | |
389 | // | |
390 | //make seed for propagation of TrackParam, using np = 3 outer clusters (separated by deltancls clusters) in TPCseed | |
391 | // | |
392 | ||
3aa6a136 | 393 | const Int_t rowstart = NRow()-1; |
1f30d65c | 394 | const Int_t rowstop = 0; |
395 | const Int_t drow = -1; | |
396 | ||
397 | //--- | |
398 | const Int_t np = 3; | |
399 | AliTrackPoint *tpos[np]; | |
400 | for(Int_t ii=0; ii<np; ii++) | |
401 | tpos[ii] = 0x0; | |
402 | ||
403 | const Float_t cov[6]={0,0,0, 0.01*0.01 ,0, 0.01*0.01}; | |
404 | ||
405 | //--- | |
406 | Int_t npos = 0; | |
407 | Int_t icl = 0; | |
3aa6a136 | 408 | const Int_t deltancls = NclsMin()/2-1; |
1f30d65c | 409 | Int_t oldcl = -deltancls; |
410 | ||
411 | for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){ | |
412 | AliTPCclusterMI *cls = tseed->GetClusterPointer(irow); | |
413 | if(!cls) { | |
414 | continue; | |
415 | } | |
416 | ||
417 | if( icl == (oldcl+deltancls) ){ | |
418 | Float_t txyz[3]; | |
419 | cls->GetGlobalXYZ(txyz); | |
420 | tpos[npos++] = new AliTrackPoint(txyz, cov, 0); | |
421 | //printf("------ %d %f %f %f\n", npos, txyz[0], txyz[1], txyz[2]); | |
422 | ||
423 | oldcl = icl; | |
424 | if(npos==np) break; | |
425 | } | |
426 | icl++; | |
427 | } | |
428 | if(npos!=np){ | |
429 | for(Int_t ii=0; ii<npos; ii++) | |
430 | delete tpos[ii]; | |
431 | ||
432 | return 0x0; | |
433 | } | |
434 | ||
435 | AliExternalTrackParam * trackparam = AliTrackerBase::MakeSeed(*(tpos[0]), *(tpos[1]), *(tpos[2])); | |
436 | if(!trackparam || trackparam->Pt()==0){ | |
437 | for(Int_t ii=0; ii<npos; ii++) | |
438 | delete tpos[ii]; | |
439 | delete trackparam; | |
440 | ||
441 | return 0x0; | |
442 | } | |
443 | ||
444 | //------ | |
445 | ||
446 | Double_t sxyz[3]={-999,-999,-999}, spxpypz[3]={-999,-999,-999}; | |
447 | trackparam->GetXYZ(sxyz); | |
448 | trackparam->GetPxPyPz(spxpypz); | |
449 | Double_t scov[21]; | |
450 | Int_t sign = trackparam->Charge(); | |
451 | ||
452 | //reset covariance matrix -- necessary, otherwise bad fitting result: trackparam->GetCovariance()[ii] has problem: some are strange values | |
453 | for(Int_t ii=0; ii<21; ii++) { | |
454 | scov[ii]=0; | |
455 | } | |
456 | trackparam->Set(sxyz, spxpypz, scov, sign); | |
457 | ||
458 | for(Int_t ii=0; ii<np; ii++) | |
459 | delete tpos[ii]; | |
460 | ||
461 | return trackparam; | |
462 | } | |
463 | ||
464 | //================================================================================================ | |
465 | void AliTPCCosmicUtils::IniCov(AliExternalTrackParam *trackPar, const Double_t ncl) | |
466 | { | |
467 | // | |
468 | //initialize covariance matrix | |
469 | // | |
470 | ||
471 | const Double_t ksigma=5.; | |
472 | Double_t acov[16]; | |
473 | for (Int_t i=0;i<15;i++) | |
474 | acov[i]=0; | |
475 | ||
476 | acov[0]=ksigma*ksigma; | |
477 | acov[2]=ksigma*ksigma; | |
478 | acov[5]=ksigma*ksigma; | |
479 | acov[9]=ksigma*ksigma; | |
480 | acov[14]=0.2*0.2; | |
481 | ||
482 | acov[5] = ksigma*ksigma/(ncl*ncl); | |
483 | acov[9] = ksigma*ksigma/(ncl*ncl); | |
484 | ||
485 | const Double_t resetcov = 4; | |
486 | ||
487 | trackPar->ResetCovariance(resetcov); | |
488 | //the following helps a lot!! | |
489 | trackPar->AddCovariance(acov); | |
490 | } | |
491 | ||
15d37333 | 492 | void AliTPCCosmicUtils::SingleFit(AliExternalTrackParam * trackInOld, AliExternalTrackParam * trackOutOld, const AliTPCseed *tseed, const Bool_t kinward, const Int_t rowstartshift, const Int_t rowstep, const Double_t xmin, const Double_t xmax, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, Double_t &lfit, TTreeSRedirector *debugstreamer) |
1f30d65c | 493 | { |
494 | // | |
495 | //fit single track | |
496 | // | |
497 | //kinward is the true geometry of the track. Incomming track: 1; outgoing track: 0 | |
498 | const Double_t inde = kinward? -1 : 1; | |
499 | const Double_t outde = kinward? 1 : -1; | |
500 | ||
501 | //PrintTrackParam(9000, trackOutOld); | |
502 | ||
503 | AliExternalTrackParam trackOut = *trackOutOld; | |
504 | AliExternalTrackParam trackIn; | |
505 | Int_t ksite = -999; | |
506 | ||
507 | //nmiss is from the 2 FitKernel of the last iteration | |
508 | //nfit, pchi2 is from the last FitKernel of the last iteration | |
15d37333 | 509 | |
3aa6a136 | 510 | Int_t rowouter = NRow()-1-rowstartshift; |
15d37333 | 511 | |
512 | //so that when reversed, the same rows are read! important!!! | |
513 | Int_t rowinner = rowouter - rowouter/rowstep * rowstep; | |
514 | ||
515 | TVector3 gposStart; | |
516 | TVector3 gposStop; | |
517 | TVector3 dpos; | |
518 | lfit = 0; | |
519 | ||
3aa6a136 | 520 | for(Int_t ii=0; ii<Niter(); ii++){ |
1f30d65c | 521 | nmiss = 0; |
522 | ||
15d37333 | 523 | gposStart.SetXYZ(-999,-999,-999); |
1f30d65c | 524 | ksite = -999; |
525 | trackIn = trackOut; | |
15d37333 | 526 | FitKernel(&trackIn, tseed, rowouter, rowinner, -rowstep, xmin, xmax, inde , ksite, nfit, nmiss, pchi2, gposStart, gposStop, 0x0, kTRUE); |
527 | //PrintTrackParam(9010+ii, &trackIn); | |
528 | ||
529 | dpos = gposStart-gposStop; | |
530 | lfit += dpos.Pt(); | |
531 | ||
532 | //--------------------------- | |
1f30d65c | 533 | |
534 | nfit = 0; | |
535 | pchi2 = 0; | |
536 | ||
15d37333 | 537 | gposStart.SetXYZ(-999,-999,-999); |
1f30d65c | 538 | ksite = -999; |
539 | trackOut = trackIn; | |
3aa6a136 | 540 | FitKernel(&trackOut, tseed, rowinner, rowouter, rowstep, xmin, xmax, outde, ksite, nfit, nmiss, pchi2, gposStart, gposStop, (ii==Niter()-1 ? debugstreamer : 0x0), kTRUE); |
1f30d65c | 541 | //PrintTrackParam(90020+ii, &trackOut); |
15d37333 | 542 | |
543 | dpos = gposStart-gposStop; | |
544 | lfit += dpos.Pt(); | |
1f30d65c | 545 | } |
546 | ||
3aa6a136 | 547 | lfit /= 2.*Niter(); |
15d37333 | 548 | |
1f30d65c | 549 | if(trackInOld) |
550 | (*trackInOld) = trackIn; | |
551 | ||
552 | (*trackOutOld) = trackOut; | |
553 | } | |
554 | ||
15d37333 | 555 | void AliTPCCosmicUtils::CombinedFit(AliExternalTrackParam *trackPars[], const AliTPCseed *seeds[], const Int_t rowstartshift, const Int_t rowstep, const Double_t xmin, const Double_t xmax, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, Double_t &lfit, Double_t &vtxD, Double_t &vtxZ, TTreeSRedirector *debugstreamer) |
1f30d65c | 556 | { |
557 | // | |
558 | //combined propagation | |
559 | // | |
560 | ||
561 | //seen from 1/pt_fit - 1/pt_gen, for differnt de: | |
562 | //u+d+: combUp good, combLow bad | |
563 | //u+d-: combUp good, combLow good | |
564 | //u-d-: combUp bad, combLow good | |
565 | //u-d+: combUp bad, combLow bad | |
566 | const Double_t upde = 1; | |
567 | const Double_t lowde = -1; | |
568 | ||
569 | AliExternalTrackParam trackLow = *(trackPars[1]); | |
570 | AliExternalTrackParam trackUp; | |
571 | ||
3aa6a136 | 572 | for(Int_t ii=0; ii<Niter(); ii++){ |
15d37333 | 573 | lfit = 0; |
574 | vtxD = 0; | |
575 | vtxZ = 0; | |
576 | Double_t tmpl = -999, tmpd = -999, tmpz = -999; | |
1f30d65c | 577 | |
578 | //lower->upper | |
579 | trackUp = trackLow; | |
15d37333 | 580 | SubCombined(&trackUp, seeds, 1, 0, rowstartshift, rowstep, xmin, xmax, upde, nfit, nmiss, pchi2, tmpl, tmpd, tmpz); |
581 | lfit += tmpl; | |
582 | vtxD += tmpd; | |
583 | vtxZ += tmpz; | |
1f30d65c | 584 | |
1f30d65c | 585 | //upper->lower |
586 | trackLow = trackUp; | |
3aa6a136 | 587 | SubCombined(&trackLow, seeds, 0, 1, rowstartshift, rowstep, xmin, xmax, lowde, nfit, nmiss, pchi2, tmpl, tmpd, tmpz, (ii==Niter()-1? debugstreamer : 0x0)); |
15d37333 | 588 | lfit += tmpl; |
589 | vtxD += tmpd; | |
590 | vtxZ += tmpz; | |
1f30d65c | 591 | } |
592 | ||
593 | *(trackPars[0]) = trackUp; | |
594 | *(trackPars[1]) = trackLow; | |
15d37333 | 595 | |
596 | //only last iteration used | |
597 | lfit /= 2; | |
598 | vtxD /= 2; | |
599 | vtxZ /= 2; | |
1f30d65c | 600 | } |
601 | ||
15d37333 | 602 | void AliTPCCosmicUtils::SubCombined(AliExternalTrackParam *trackPar, const AliTPCseed *seeds[], const Int_t tk0, const Int_t tk1, const Int_t rowstartshift, const Int_t rowstep, const Double_t xmin, const Double_t xmax, const Double_t eloss, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, Double_t &lfit, Double_t &vtxD, Double_t &vtxZ, TTreeSRedirector *debugstreamer) |
1f30d65c | 603 | { |
604 | // | |
605 | //sub-routine for combined propagation | |
606 | // | |
607 | ||
608 | IniCov(trackPar, seeds[0]->GetNumberOfClusters()+seeds[1]->GetNumberOfClusters()); | |
609 | ||
610 | Int_t dtk=1; | |
611 | if(tk0>tk1) | |
612 | dtk=-1; | |
613 | ||
614 | //reset counters | |
615 | Int_t ksite = -999; | |
616 | nfit = 0; | |
617 | nmiss = 0; | |
618 | pchi2 = 0; | |
15d37333 | 619 | vtxD = -1e10; |
620 | vtxZ = -1e10; | |
1f30d65c | 621 | |
622 | //always nrow -> 1 -> nrow | |
3aa6a136 | 623 | Int_t rowstart= NRow()-1-rowstartshift; |
15d37333 | 624 | |
625 | //so that when reversed, the same rows are read! important!!! | |
626 | Int_t rowstop = rowstart - rowstart/rowstep * rowstep; | |
627 | Int_t drow = -rowstep; | |
628 | ||
629 | TVector3 gposStart(-999,-999,-999); | |
630 | TVector3 gposStop(-999,-999,-999); | |
631 | ||
1f30d65c | 632 | for(Int_t itrk=tk0; dtk*itrk<=dtk*tk1; itrk+=dtk){ |
633 | if(itrk==tk1){ | |
634 | Int_t tmprow = rowstart; | |
635 | rowstart = rowstop; | |
636 | rowstop = tmprow; | |
15d37333 | 637 | drow *= -1; |
638 | } | |
639 | ||
640 | FitKernel(trackPar, seeds[itrk], rowstart, rowstop, drow, xmin, xmax, eloss, ksite, nfit, nmiss, pchi2, gposStart, gposStop, debugstreamer, kFALSE); | |
641 | ||
642 | //get the impact parameters at the end of the first propagation (X=0) | |
643 | if(itrk==tk0){ | |
644 | AliExternalTrackParam vertex(*trackPar); | |
645 | const Double_t maxStep = 1; | |
646 | const Bool_t rotateTo = kFALSE; | |
647 | const Double_t maxSnp = 0.8; | |
3aa6a136 | 648 | if(AliTrackerBase::PropagateTrackToBxByBz(&vertex, 0, Mass(), maxStep, rotateTo, maxSnp, (Int_t)eloss)){ |
15d37333 | 649 | vtxD = TMath::Abs(vertex.GetParameter()[0]); |
650 | vtxZ = vertex.GetParameter()[1]; | |
651 | } | |
1f30d65c | 652 | } |
1f30d65c | 653 | } |
15d37333 | 654 | |
655 | TVector3 dpos = gposStart-gposStop; | |
656 | lfit = dpos.Pt(); | |
1f30d65c | 657 | } |
658 | ||
15d37333 | 659 | void AliTPCCosmicUtils::FitKernel(AliExternalTrackParam *trackPar, const AliTPCseed *tseed, const Int_t rowstart, const Int_t rowstop, const Int_t drow, const Double_t xmin, const Double_t xmax, const Double_t eloss, Int_t &ksite, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, TVector3 &gposStart, TVector3 &gposStop, TTreeSRedirector *debugstreamer, const Bool_t kinicov) |
1f30d65c | 660 | { |
661 | // | |
662 | //routine for propagation | |
663 | // | |
664 | ||
665 | //only for SingleFit--> | |
666 | if(kinicov) | |
667 | IniCov(trackPar, tseed->GetNumberOfClusters()); | |
668 | //<-- | |
669 | ||
15d37333 | 670 | Int_t checkstop = -999; |
1f30d65c | 671 | |
672 | for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){ | |
15d37333 | 673 | checkstop = irow; |
674 | //printf("test irow %d\n", irow); | |
1f30d65c | 675 | |
676 | AliTPCclusterMI *cl=tseed->GetClusterPointer(irow); | |
677 | if (!cl) continue; | |
3aa6a136 | 678 | if (cl->GetX()< XMin()) continue; |
1f30d65c | 679 | |
15d37333 | 680 | //cut on cluster X (i.e. r) to specify leverarm |
681 | if(cl->GetX()< xmin){ | |
682 | continue; | |
683 | } | |
684 | if(cl->GetX()> xmax){ | |
685 | continue; | |
686 | } | |
1f30d65c | 687 | //if propagation not successful, trackPar is not changed |
688 | AliExternalTrackParam tmppar(*trackPar); | |
689 | ||
690 | //--------------------------- to cure large chi2 in z, 2011-05-11. Thandks to M. Ivanov! --------------------------- | |
691 | const Int_t tmpsite = ((cl->GetDetector()%36)<18); | |
692 | if(tmpsite!=ksite && ksite!=-999){ | |
693 | Double_t *clscov=(Double_t*)tmppar.GetCovariance(); | |
694 | clscov[2]+=3*3; | |
695 | } | |
696 | ksite = tmpsite; | |
697 | ||
698 | //--------------------------- rotate cluster position to trackPar position --------------------------- | |
699 | //DO NOT rotate trackPar, there is "bug" in trackparam::rotate!! stay in the initial one defined as alpha = ATan2(posy0,posx0) | |
700 | Float_t gxyz[3]; | |
701 | cl->GetGlobalXYZ(gxyz); | |
15d37333 | 702 | |
703 | const TVector3 gptmp(gxyz[0], gxyz[1], gxyz[2]); | |
704 | ||
1f30d65c | 705 | TVector3 ptogo(gxyz); |
706 | ||
707 | //printf("test %d : %f %f %f\n", irow, gxyz[0], gxyz[1], gxyz[2]); | |
708 | ||
709 | const Double_t trkalpha = tmppar.GetAlpha(); | |
710 | ptogo.RotateZ(-trkalpha); | |
711 | const Double_t xTogo = ptogo.X(); | |
712 | const Double_t yTogo = ptogo.Y(); | |
713 | const Double_t zTogo = ptogo.Z(); | |
714 | ||
715 | //--------------------------- propagate --------------------------- | |
716 | //AliTrackerBase::PropagateTrackToBxByBz Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); | |
717 | const Double_t maxStep = 1; | |
718 | const Bool_t rotateTo = kFALSE; | |
719 | const Double_t maxSnp = 0.8; | |
720 | ||
721 | //eloss critical only for combine fit!!! | |
3aa6a136 | 722 | const Bool_t kpro = AliTrackerBase::PropagateTrackToBxByBz(&tmppar, xTogo, Mass(), maxStep, rotateTo, maxSnp, (Int_t)eloss); |
1f30d65c | 723 | if(!kpro){ |
724 | nmiss++; | |
725 | continue; | |
726 | } | |
727 | ||
728 | //--------------------------- set error --------------------------- | |
729 | //necessary, help at mostly low p, and also high p | |
730 | Double_t cov[3]={0.01, 0., 0.01}; | |
731 | AliTPCseed::GetError(cl, &tmppar, cov[0],cov[2]); | |
732 | cov[0]*=cov[0]; | |
733 | cov[2]*=cov[2]; | |
734 | ||
735 | if(debugstreamer){ | |
736 | TVector3 tmpgxyz(gxyz); | |
737 | (*debugstreamer)<<"ErrParam"<< | |
738 | "Cl.="<<cl<< | |
739 | "T.="<<&tmppar<< | |
740 | "covy="<<cov[0]<< | |
741 | "covz="<<cov[2]<< | |
742 | "clpos.="<<&ptogo<< | |
743 | "gpos.="<<&tmpgxyz<< | |
744 | "\n"; | |
745 | } | |
746 | ||
747 | //--------------------------- get chi2 and THEN update --------------------------- | |
748 | Double_t yz[2]={yTogo, zTogo}; | |
749 | ||
750 | pchi2 += tmppar.GetPredictedChi2(yz, cov); | |
751 | ||
752 | tmppar.Update(yz,cov); | |
753 | ||
754 | //--------------------------- save change --------------------------- | |
755 | (*trackPar) = tmppar; | |
756 | ||
757 | nfit++; | |
15d37333 | 758 | |
759 | gposStop = gptmp; | |
760 | if(gposStart.X()<-998) | |
761 | gposStart = gptmp; | |
762 | } | |
763 | ||
764 | //to make sure rowstart and rowstop are the actual ending | |
765 | if(checkstop != rowstop){ | |
3aa6a136 | 766 | printf("AliTPCCosmicUtils::FitKernel error wrong rowstart, stop, drow!! %d %d %d\n", rowstart, rowstop, drow); exit(1); |
1f30d65c | 767 | } |
768 | } |