Line data Source code
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 : //
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
20 : //
21 : // grep "error" in output log to check abnormal termination
22 : //
23 : //
24 : // Xianguo Lu
25 : // lu@physi.uni-heidelberg.de
26 : // Xianguo.Lu@cern.ch
27 :
28 : //
29 :
30 : #include <TAxis.h>
31 : #include <TCanvas.h>
32 : #include <TGraph.h>
33 : #include <TTreeStream.h>
34 : #include <TVector3.h>
35 :
36 : #include "AliAnalysisManager.h"
37 : #include "AliESDEvent.h"
38 : #include "AliESDfriend.h"
39 : #include "AliESDtrack.h"
40 : #include "AliESDfriendTrack.h"
41 : #include "AliESDInputHandler.h"
42 : #include "AliTPCseed.h"
43 : #include "AliTrackerBase.h"
44 : #include "AliTrackPointArray.h"
45 :
46 : #include "AliTPCCosmicUtils.h"
47 :
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 0 : const Double_t dm = esd->GetMagneticField();
54 :
55 : Int_t mag = -999;
56 0 : if(dm>0)
57 0 : mag = (Int_t) (dm+0.5);
58 : else
59 0 : mag = (Int_t) (dm-0.5);
60 :
61 0 : if(TMath::Abs(mag)!=1 && TMath::Abs(mag)!=2 && TMath::Abs(mag)!=5){
62 0 : printf("error AliTPCCosmicUtils::GetB strange Bfield! %f %d\n", dm, mag); exit(1);
63 : }
64 :
65 0 : return mag;
66 : }
67 :
68 : Int_t AliTPCCosmicUtils::GetTrigger(const AliESDEvent *esd)
69 : {
70 : //
71 : //get cosmic trigger
72 : //
73 0 : const TString strig = esd->GetFiredTriggerClasses();
74 : Int_t btrig = 0;
75 :
76 0 : if( strig.Contains("C0OB0"))
77 0 : btrig += k0OB0;
78 0 : if( strig.Contains("C0OB1"))
79 0 : btrig += k0OB1;
80 0 : if( strig.Contains("C0HWU"))
81 0 : btrig += k0HWU;
82 0 : if( strig.Contains("CTRDCO2"))
83 0 : btrig += kTRDCO2;
84 0 : if( strig.Contains("AMU"))
85 0 : btrig += kAMU;
86 0 : if( strig.Contains("SCO"))
87 0 : btrig += kSCO;
88 :
89 : return btrig;
90 0 : }
91 :
92 : Bool_t AliTPCCosmicUtils::GetESD(AliESDEvent *& esdevent, AliESDfriend *& esdfriend)
93 : {
94 : //
95 : //get esdevent and esdfriend
96 : //
97 :
98 0 : esdevent = 0x0;
99 0 : esdfriend = 0x0;
100 :
101 0 : AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
102 0 : if(esdH){
103 0 : esdevent = (AliESDEvent*)esdH->GetEvent();
104 0 : }
105 :
106 0 : if(!esdevent){
107 0 : printf("AliAnalysisTaskdEdxCosmic::GetESD error No ESD Event\n");
108 0 : return kFALSE;
109 : }
110 :
111 0 : esdH->SetActiveBranches("ESDfriend*");
112 0 : esdfriend = (AliESDfriend *)esdevent->FindListObject("AliESDfriend");
113 :
114 0 : if(!esdfriend){
115 0 : printf("AliAnalysisTaskdEdxCosmic::GetESD error No ESD friend\n");
116 0 : return kFALSE;
117 : }
118 :
119 0 : esdevent->SetESDfriend(esdfriend);
120 :
121 0 : return kTRUE;
122 0 : }
123 :
124 :
125 : Double_t AliTPCCosmicUtils::GetMinPhi(const AliExternalTrackParam *params[])
126 : {
127 : //
128 : //minimum phi angle of the two tracks, 0: horizontal
129 : //
130 0 : const Double_t fsp[] = {TMath::Abs(TMath::Sin(params[0]->Phi())),
131 0 : TMath::Abs(TMath::Sin(params[1]->Phi()))};
132 0 : const Double_t minphi = TMath::ASin(TMath::Min(fsp[0], fsp[1])) * TMath::RadToDeg();
133 :
134 0 : return minphi;
135 : }
136 :
137 : Double_t AliTPCCosmicUtils::Point2LineDist(const TVector3 p0, const TVector3 l1, const TVector3 l2, TVector3 *vtx)
138 : {
139 : //
140 : //return distance of p0 to line (l2-l1)
141 : //
142 :
143 0 : const TVector3 va = p0 - l1;
144 0 : const TVector3 vb = l2 - l1;
145 :
146 : //calculate vtx position
147 0 : const Double_t ca = va.Dot(vb)/vb.Mag();
148 0 : const TVector3 imp = l1 + vb.Unit()*ca;
149 :
150 0 : if(vtx){
151 0 : (*vtx) = imp;
152 0 : }
153 :
154 : /* equivalent to imp.Mag()
155 : //calculate distance
156 : const TVector3 dd = va.Cross(vb);
157 : return dd.Mag()/vb.Mag();
158 : */
159 :
160 0 : return imp.Mag();
161 0 : }
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 1182 : const Double_t pmin = -TMath::Pi();
169 591 : const Double_t pmax = pmin+TMath::TwoPi();
170 :
171 1182 : while(phi<pmin)
172 0 : phi+= TMath::TwoPi();
173 :
174 1187 : while(phi>=pmax)
175 298 : phi-= TMath::TwoPi();
176 :
177 591 : return phi;
178 : }
179 :
180 : Bool_t AliTPCCosmicUtils::RotateSafe(AliExternalTrackParam *trackPar, const Double_t rawalpha)
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 :
197 : //the range-changing in AliExternalTrackParam::Rotate not precise enough
198 12 : const Double_t aa = AngleInRange(rawalpha);
199 :
200 6 : const Double_t a0 = trackPar->GetAlpha();
201 6 : const Double_t p2 = trackPar->GetParameter()[2];
202 :
203 : //copied from AliExternalTrackParam::Rotate
204 6 : const Double_t ca=TMath::Cos(aa-a0), sa=TMath::Sin(aa-a0);
205 6 : const Double_t sf=p2, cf=TMath::Sqrt((1.- p2)*(1.+p2));
206 :
207 6 : if((cf*ca+sf*sa)<0) {
208 0 : return kFALSE;
209 : }
210 :
211 6 : return trackPar->Rotate(aa);
212 :
213 : /*
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;
239 : */
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 6 : }
262 :
263 : void AliTPCCosmicUtils::PrintTrackParam(const Int_t id, const AliExternalTrackParam * trackpar, const char *tag)
264 : {
265 : //
266 : //print out TrackParam
267 : //
268 0 : Double_t xyz[3]={-999,-999,-999};
269 0 : trackpar->GetXYZ(xyz);
270 0 : Double_t pxpypz[3]={-999,-999,-999};
271 0 : trackpar->GetPxPyPz(pxpypz);
272 :
273 0 : 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 0 : }
275 :
276 : void AliTPCCosmicUtils::DrawTracks(AliESDtrack *esdtrks[], const TString tag, const TString outputformat)
277 : {
278 : //
279 : //draw esdtracks
280 : //
281 0 : const AliTPCseed *seeds[]={GetTPCseed(esdtrks[0]), GetTPCseed(esdtrks[1])};
282 0 : DrawSeeds(seeds, tag, outputformat);
283 0 : }
284 :
285 :
286 : void AliTPCCosmicUtils::DrawSeeds(const AliTPCseed * seeds[], const TString tag, const TString outputformat)
287 : {
288 : //
289 : //draw seed and output to file
290 : //
291 :
292 0 : if(!seeds[0] || !seeds[1])
293 : return;
294 :
295 0 : TGraph *grsyx[]={new TGraph, new TGraph};
296 0 : TGraph *grsyz[]={new TGraph, new TGraph};
297 :
298 0 : for(Int_t itrk=0; itrk<2; itrk++){
299 0 : if(!seeds[itrk])
300 : continue;
301 :
302 : Int_t ncl = 0;
303 0 : for(Int_t irow=0; irow<NRow(); irow++){
304 0 : const AliTPCclusterMI * cls = seeds[itrk]->GetClusterPointer(irow);
305 0 : if(!cls)
306 0 : continue;
307 :
308 0 : Float_t xyz[3]={-999,-999,-999};
309 0 : cls->GetGlobalXYZ(xyz);
310 : //printf("Test %f %f %f \n", xyz[0], xyz[1], xyz[2]);
311 0 : grsyx[itrk]->SetPoint(ncl, xyz[0], xyz[1]);
312 0 : grsyz[itrk]->SetPoint(ncl, xyz[2], xyz[1]);
313 0 : ncl++;
314 0 : }
315 0 : }
316 :
317 0 : grsyx[0]->SetTitle(tag);
318 0 : grsyx[0]->SetMaximum(250);
319 0 : grsyx[0]->SetMinimum(-250);
320 0 : grsyx[0]->GetXaxis()->SetLimits(-250,250);
321 0 : grsyx[0]->GetXaxis()->SetTitle("X (cm)");
322 0 : grsyx[0]->GetYaxis()->SetTitle("Y (cm)");
323 :
324 0 : grsyx[0]->SetMarkerStyle(20);
325 0 : grsyx[0]->SetMarkerColor(kRed);
326 0 : grsyx[1]->SetMarkerStyle(22);
327 0 : grsyx[1]->SetMarkerColor(kBlue);
328 :
329 0 : grsyz[0]->SetTitle(tag);
330 0 : grsyz[0]->SetMaximum(250);
331 0 : grsyz[0]->SetMinimum(-250);
332 0 : grsyz[0]->GetXaxis()->SetLimits(-250,250);
333 0 : grsyz[0]->GetXaxis()->SetTitle("Z (cm)");
334 0 : grsyz[0]->GetYaxis()->SetTitle("Y (cm)");
335 :
336 0 : grsyz[0]->SetMarkerStyle(20);
337 0 : grsyz[0]->SetMarkerColor(kRed);
338 0 : grsyz[1]->SetMarkerStyle(22);
339 0 : grsyz[1]->SetMarkerColor(kBlue);
340 :
341 0 : TCanvas *cc=new TCanvas("cc","",1200,600);
342 0 : cc->Divide(2);
343 0 : cc->cd(1);
344 0 : for(int itrk=0; itrk<2; itrk++){
345 0 : grsyx[itrk]->SetMarkerSize(1);
346 0 : grsyx[itrk]->Draw(itrk?"lp same":"alp");
347 : }
348 0 : cc->cd(2);
349 0 : for(int itrk=0; itrk<2; itrk++){
350 0 : grsyz[itrk]->SetMarkerSize(1);
351 0 : grsyz[itrk]->Draw(itrk?"lp same":"alp");
352 : }
353 :
354 0 : gErrorIgnoreLevel = 1001;
355 0 : cc->Print(Form("drawTrack%s.%s", tag.Data(), outputformat.Data()));
356 :
357 0 : for(Int_t ii=0; ii<2; ii++){
358 0 : delete grsyx[ii];
359 0 : delete grsyz[ii];
360 : }
361 0 : delete cc;
362 0 : }
363 :
364 : //================================================================================================
365 : AliTPCseed * AliTPCCosmicUtils::GetTPCseed(const AliESDtrack *esdtrack)
366 : {
367 : //
368 : //Get TPC seeds from ESDfriendTrack
369 : //
370 :
371 1404 : AliESDfriendTrack * friendtrk = (AliESDfriendTrack *) esdtrack->GetFriendTrack();
372 702 : if(!friendtrk){
373 : //printf("AliTPCCosmicUtils::GetTPCseed no friend track!\n"); exit(1);
374 0 : return 0x0;
375 : }
376 :
377 : TObject *calibObject=0x0;
378 : AliTPCseed *tseed = 0x0;
379 2864 : for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
380 3999 : if( (tseed=dynamic_cast<AliTPCseed*>(calibObject)) )
381 636 : break;
382 : }
383 : return tseed;
384 702 : }
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 :
393 8 : const Int_t rowstart = NRow()-1;
394 : const Int_t rowstop = 0;
395 : const Int_t drow = -1;
396 :
397 : //---
398 : const Int_t np = 3;
399 4 : AliTrackPoint *tpos[np];
400 32 : for(Int_t ii=0; ii<np; ii++)
401 12 : 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;
408 4 : const Int_t deltancls = NclsMin()/2-1;
409 4 : Int_t oldcl = -deltancls;
410 :
411 312 : for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){
412 156 : AliTPCclusterMI *cls = tseed->GetClusterPointer(irow);
413 156 : if(!cls) {
414 0 : continue;
415 : }
416 :
417 156 : if( icl == (oldcl+deltancls) ){
418 12 : Float_t txyz[3];
419 12 : cls->GetGlobalXYZ(txyz);
420 24 : 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 16 : if(npos==np) break;
425 20 : }
426 152 : icl++;
427 152 : }
428 4 : if(npos!=np){
429 0 : for(Int_t ii=0; ii<npos; ii++)
430 0 : delete tpos[ii];
431 :
432 0 : return 0x0;
433 : }
434 :
435 4 : AliExternalTrackParam * trackparam = AliTrackerBase::MakeSeed(*(tpos[0]), *(tpos[1]), *(tpos[2]));
436 8 : if(!trackparam || trackparam->Pt()==0){
437 0 : for(Int_t ii=0; ii<npos; ii++)
438 0 : delete tpos[ii];
439 0 : delete trackparam;
440 :
441 0 : return 0x0;
442 : }
443 :
444 : //------
445 :
446 4 : Double_t sxyz[3]={-999,-999,-999}, spxpypz[3]={-999,-999,-999};
447 4 : trackparam->GetXYZ(sxyz);
448 4 : trackparam->GetPxPyPz(spxpypz);
449 4 : Double_t scov[21];
450 4 : 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 176 : for(Int_t ii=0; ii<21; ii++) {
454 84 : scov[ii]=0;
455 : }
456 4 : trackparam->Set(sxyz, spxpypz, scov, sign);
457 :
458 32 : for(Int_t ii=0; ii<np; ii++)
459 24 : delete tpos[ii];
460 :
461 : return trackparam;
462 8 : }
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 16 : Double_t acov[16];
473 256 : for (Int_t i=0;i<15;i++)
474 120 : acov[i]=0;
475 :
476 8 : acov[0]=ksigma*ksigma;
477 8 : acov[2]=ksigma*ksigma;
478 8 : acov[5]=ksigma*ksigma;
479 8 : acov[9]=ksigma*ksigma;
480 8 : acov[14]=0.2*0.2;
481 :
482 8 : acov[5] = ksigma*ksigma/(ncl*ncl);
483 8 : acov[9] = ksigma*ksigma/(ncl*ncl);
484 :
485 : const Double_t resetcov = 4;
486 :
487 8 : trackPar->ResetCovariance(resetcov);
488 : //the following helps a lot!!
489 8 : trackPar->AddCovariance(acov);
490 8 : }
491 :
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)
493 : {
494 : //
495 : //fit single track
496 : //
497 : //kinward is the true geometry of the track. Incomming track: 1; outgoing track: 0
498 0 : const Double_t inde = kinward? -1 : 1;
499 0 : const Double_t outde = kinward? 1 : -1;
500 :
501 : //PrintTrackParam(9000, trackOutOld);
502 :
503 0 : AliExternalTrackParam trackOut = *trackOutOld;
504 0 : AliExternalTrackParam trackIn;
505 0 : 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
509 :
510 0 : Int_t rowouter = NRow()-1-rowstartshift;
511 :
512 : //so that when reversed, the same rows are read! important!!!
513 0 : Int_t rowinner = rowouter - rowouter/rowstep * rowstep;
514 :
515 0 : TVector3 gposStart;
516 0 : TVector3 gposStop;
517 0 : TVector3 dpos;
518 0 : lfit = 0;
519 :
520 0 : for(Int_t ii=0; ii<Niter(); ii++){
521 0 : nmiss = 0;
522 :
523 0 : gposStart.SetXYZ(-999,-999,-999);
524 0 : ksite = -999;
525 0 : trackIn = trackOut;
526 0 : FitKernel(&trackIn, tseed, rowouter, rowinner, -rowstep, xmin, xmax, inde , ksite, nfit, nmiss, pchi2, gposStart, gposStop, 0x0, kTRUE);
527 : //PrintTrackParam(9010+ii, &trackIn);
528 :
529 0 : dpos = gposStart-gposStop;
530 0 : lfit += dpos.Pt();
531 :
532 : //---------------------------
533 :
534 0 : nfit = 0;
535 0 : pchi2 = 0;
536 :
537 0 : gposStart.SetXYZ(-999,-999,-999);
538 0 : ksite = -999;
539 0 : trackOut = trackIn;
540 0 : FitKernel(&trackOut, tseed, rowinner, rowouter, rowstep, xmin, xmax, outde, ksite, nfit, nmiss, pchi2, gposStart, gposStop, (ii==Niter()-1 ? debugstreamer : 0x0), kTRUE);
541 : //PrintTrackParam(90020+ii, &trackOut);
542 :
543 0 : dpos = gposStart-gposStop;
544 0 : lfit += dpos.Pt();
545 : }
546 :
547 0 : lfit /= 2.*Niter();
548 :
549 0 : if(trackInOld)
550 0 : (*trackInOld) = trackIn;
551 :
552 0 : (*trackOutOld) = trackOut;
553 0 : }
554 :
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)
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 4 : AliExternalTrackParam trackLow = *(trackPars[1]);
570 2 : AliExternalTrackParam trackUp;
571 :
572 12 : for(Int_t ii=0; ii<Niter(); ii++){
573 4 : lfit = 0;
574 4 : vtxD = 0;
575 4 : vtxZ = 0;
576 4 : Double_t tmpl = -999, tmpd = -999, tmpz = -999;
577 :
578 : //lower->upper
579 4 : trackUp = trackLow;
580 4 : SubCombined(&trackUp, seeds, 1, 0, rowstartshift, rowstep, xmin, xmax, upde, nfit, nmiss, pchi2, tmpl, tmpd, tmpz);
581 4 : lfit += tmpl;
582 4 : vtxD += tmpd;
583 4 : vtxZ += tmpz;
584 :
585 : //upper->lower
586 4 : trackLow = trackUp;
587 4 : SubCombined(&trackLow, seeds, 0, 1, rowstartshift, rowstep, xmin, xmax, lowde, nfit, nmiss, pchi2, tmpl, tmpd, tmpz, (ii==Niter()-1? debugstreamer : 0x0));
588 4 : lfit += tmpl;
589 4 : vtxD += tmpd;
590 4 : vtxZ += tmpz;
591 4 : }
592 :
593 2 : *(trackPars[0]) = trackUp;
594 2 : *(trackPars[1]) = trackLow;
595 :
596 : //only last iteration used
597 2 : lfit /= 2;
598 2 : vtxD /= 2;
599 2 : vtxZ /= 2;
600 2 : }
601 :
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)
603 : {
604 : //
605 : //sub-routine for combined propagation
606 : //
607 :
608 16 : IniCov(trackPar, seeds[0]->GetNumberOfClusters()+seeds[1]->GetNumberOfClusters());
609 :
610 : Int_t dtk=1;
611 8 : if(tk0>tk1)
612 : dtk=-1;
613 :
614 : //reset counters
615 8 : Int_t ksite = -999;
616 8 : nfit = 0;
617 8 : nmiss = 0;
618 8 : pchi2 = 0;
619 8 : vtxD = -1e10;
620 8 : vtxZ = -1e10;
621 :
622 : //always nrow -> 1 -> nrow
623 8 : Int_t rowstart= NRow()-1-rowstartshift;
624 :
625 : //so that when reversed, the same rows are read! important!!!
626 8 : Int_t rowstop = rowstart - rowstart/rowstep * rowstep;
627 8 : Int_t drow = -rowstep;
628 :
629 8 : TVector3 gposStart(-999,-999,-999);
630 8 : TVector3 gposStop(-999,-999,-999);
631 :
632 48 : for(Int_t itrk=tk0; dtk*itrk<=dtk*tk1; itrk+=dtk){
633 16 : if(itrk==tk1){
634 : Int_t tmprow = rowstart;
635 : rowstart = rowstop;
636 : rowstop = tmprow;
637 8 : drow *= -1;
638 8 : }
639 :
640 16 : 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 16 : if(itrk==tk0){
644 8 : AliExternalTrackParam vertex(*trackPar);
645 : const Double_t maxStep = 1;
646 : const Bool_t rotateTo = kFALSE;
647 : const Double_t maxSnp = 0.8;
648 16 : if(AliTrackerBase::PropagateTrackToBxByBz(&vertex, 0, Mass(), maxStep, rotateTo, maxSnp, (Int_t)eloss)){
649 8 : vtxD = TMath::Abs(vertex.GetParameter()[0]);
650 8 : vtxZ = vertex.GetParameter()[1];
651 8 : }
652 8 : }
653 : }
654 :
655 8 : TVector3 dpos = gposStart-gposStop;
656 16 : lfit = dpos.Pt();
657 8 : }
658 :
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)
660 : {
661 : //
662 : //routine for propagation
663 : //
664 :
665 : //only for SingleFit-->
666 16 : if(kinicov)
667 0 : IniCov(trackPar, tseed->GetNumberOfClusters());
668 : //<--
669 :
670 : Int_t checkstop = -999;
671 :
672 5120 : for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){
673 : checkstop = irow;
674 : //printf("test irow %d\n", irow);
675 :
676 2544 : AliTPCclusterMI *cl=tseed->GetClusterPointer(irow);
677 2692 : if (!cl) continue;
678 2396 : if (cl->GetX()< XMin()) continue;
679 :
680 : //cut on cluster X (i.e. r) to specify leverarm
681 2396 : if(cl->GetX()< xmin){
682 0 : continue;
683 : }
684 2396 : if(cl->GetX()> xmax){
685 0 : continue;
686 : }
687 : //if propagation not successful, trackPar is not changed
688 2396 : AliExternalTrackParam tmppar(*trackPar);
689 :
690 : //--------------------------- to cure large chi2 in z, 2011-05-11. Thandks to M. Ivanov! ---------------------------
691 4792 : const Int_t tmpsite = ((cl->GetDetector()%36)<18);
692 2412 : if(tmpsite!=ksite && ksite!=-999){
693 8 : Double_t *clscov=(Double_t*)tmppar.GetCovariance();
694 8 : clscov[2]+=3*3;
695 8 : }
696 2396 : 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 2396 : Float_t gxyz[3];
701 2396 : cl->GetGlobalXYZ(gxyz);
702 :
703 2396 : const TVector3 gptmp(gxyz[0], gxyz[1], gxyz[2]);
704 :
705 2396 : TVector3 ptogo(gxyz);
706 :
707 : //printf("test %d : %f %f %f\n", irow, gxyz[0], gxyz[1], gxyz[2]);
708 :
709 2396 : const Double_t trkalpha = tmppar.GetAlpha();
710 2396 : ptogo.RotateZ(-trkalpha);
711 2396 : const Double_t xTogo = ptogo.X();
712 2396 : const Double_t yTogo = ptogo.Y();
713 2396 : 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!!!
722 4792 : const Bool_t kpro = AliTrackerBase::PropagateTrackToBxByBz(&tmppar, xTogo, Mass(), maxStep, rotateTo, maxSnp, (Int_t)eloss);
723 2396 : if(!kpro){
724 0 : nmiss++;
725 0 : continue;
726 : }
727 :
728 : //--------------------------- set error ---------------------------
729 : //necessary, help at mostly low p, and also high p
730 2396 : Double_t cov[3]={0.01, 0., 0.01};
731 2396 : AliTPCseed::GetError(cl, &tmppar, cov[0],cov[2]);
732 2396 : cov[0]*=cov[0];
733 2396 : cov[2]*=cov[2];
734 :
735 2396 : if(debugstreamer){
736 0 : TVector3 tmpgxyz(gxyz);
737 0 : (*debugstreamer)<<"ErrParam"<<
738 0 : "Cl.="<<cl<<
739 0 : "T.="<<&tmppar<<
740 0 : "covy="<<cov[0]<<
741 0 : "covz="<<cov[2]<<
742 0 : "clpos.="<<&ptogo<<
743 0 : "gpos.="<<&tmpgxyz<<
744 : "\n";
745 0 : }
746 :
747 : //--------------------------- get chi2 and THEN update ---------------------------
748 2396 : Double_t yz[2]={yTogo, zTogo};
749 :
750 4792 : pchi2 += tmppar.GetPredictedChi2(yz, cov);
751 :
752 2396 : tmppar.Update(yz,cov);
753 :
754 : //--------------------------- save change ---------------------------
755 2396 : (*trackPar) = tmppar;
756 :
757 2396 : nfit++;
758 :
759 2396 : gposStop = gptmp;
760 2396 : if(gposStart.X()<-998)
761 8 : gposStart = gptmp;
762 4792 : }
763 :
764 : //to make sure rowstart and rowstop are the actual ending
765 16 : if(checkstop != rowstop){
766 0 : printf("AliTPCCosmicUtils::FitKernel error wrong rowstart, stop, drow!! %d %d %d\n", rowstart, rowstop, drow); exit(1);
767 : }
768 16 : }
|
