]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCtrackerMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removed compiler warning
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
CommitLineData
1c53abe2 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
91162307 16/*
17$Log$
8c51a605 18Revision 1.18 2003/10/17 15:42:14 kowal2
19Back to the previous version. The warning was erronously generated
20by the compiler
21
b364ca79 22Revision 1.17 2003/10/17 12:28:02 kowal2
23Removed "always true" comparison
24
d58e9f0e 25Revision 1.16 2003/10/17 12:01:16 kowal2
26Removed compiler warning.
27
176aff27 28Revision 1.15 2003/09/29 11:56:58 kowal2
29bug fix2
30
3e5d0aa2 31Revision 1.14 2003/09/29 11:39:43 kowal2
32bug fix
33
732b9e78 34Revision 1.13 2003/09/29 11:28:19 kowal2
35completly rewritten
36
91162307 37Revision 1.9.4.3 2003/06/23 14:47:10 hristov
38Minor fix
39
40Revision 1.9.4.2 2003/06/23 10:06:13 hristov
41Updated information about the overlapping clusters (M.Ivanov)
42
43Revision 1.9.4.1 2003/06/19 06:59:58 hristov
44Updated version of parallel tracking (M.Ivanov)
45
46Revision 1.9 2003/03/19 17:14:11 hristov
47Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
48
49Revision 1.8 2003/03/05 11:16:15 kowal2
50Logs added
51
52*/
53
54
55
56
57
58
1c53abe2 59
60/*
61 AliTPC parallel tracker -
62 How to use? -
63 run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
64 run AliTPCFindTracksMI.C macro - to find tracks
65 tracks are written to AliTPCtracks.root file
66 for comparison also seeds are written to the same file - to special branch
67*/
68
69//-------------------------------------------------------
70// Implementation of the TPC tracker
71//
72// Origin: Marian Ivanov Marian.Ivanov@cern.ch
73//
74//-------------------------------------------------------
1c53abe2 75#include <TObjArray.h>
76#include <TFile.h>
77#include <TTree.h>
cc5e9db0 78#include "Riostream.h"
1c53abe2 79
80#include "AliTPCtrackerMI.h"
81#include "AliTPCclusterMI.h"
82#include "AliTPCParam.h"
83#include "AliTPCClustersRow.h"
84#include "AliComplexCluster.h"
1627d1c4 85#include "AliTPCpolyTrack.h"
1c53abe2 86#include "TStopwatch.h"
91162307 87#include "AliESD.h"
88#include "AliHelix.h"
89#include "TGraph.h"
90//
91#include "AliRunLoader.h"
1c53abe2 92
93
c9427e08 94
91162307 95ClassImp(AliTPCseed)
96ClassImp(AliTPCtrackerMI)
c9427e08 97
98
91162307 99class TPCFastMath {
100public:
101 TPCFastMath();
102 static Double_t FastAsin(Double_t x);
103 private:
104 static Double_t fgFastAsin[20000];
105};
c9427e08 106
91162307 107Double_t TPCFastMath::fgFastAsin[20000];
c9427e08 108
91162307 109TPCFastMath::TPCFastMath(){
110 for (Int_t i=0;i<10000;i++){
111 fgFastAsin[2*i] = TMath::ASin(i/10000.);
112 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
113 }
c9427e08 114}
115
91162307 116Double_t TPCFastMath::FastAsin(Double_t x){
117 if (x>0){
118 Int_t index = int(x*10000);
119 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
120 }
121 x*=-1;
122 Int_t index = int(x*10000);
123 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
1c53abe2 124}
125
91162307 126TPCFastMath gTPCMath;
1c53abe2 127
128
129
91162307 130Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
1c53abe2 131
91162307 132 AliTPCclusterMI* c =track->fCurrentCluster;
133 if (accept>0) track->fCurrentClusterIndex1 |=0x8000; //sign not accepted clusters
c9427e08 134
91162307 135 UInt_t i = track->fCurrentClusterIndex1;
1c53abe2 136
137 Int_t sec=(i&0xff000000)>>24;
91162307 138 //Int_t row = (i&0x00ff0000)>>16;
1c53abe2 139 track->fRow=(i&0x00ff0000)>>16;
140 track->fSector = sec;
141 // Int_t index = i&0xFFFF;
142 if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
91162307 143 track->SetClusterIndex2(track->fRow, i);
144
145 //track->fFirstPoint = row;
146 //if ( track->fLastPoint<row) track->fLastPoint =row;
147 // if (track->fRow<0 || track->fRow>160) {
148 // printf("problem\n");
149 //}
150 if (track->fFirstPoint>track->fRow)
151 track->fFirstPoint = track->fRow;
152 if (track->fLastPoint<track->fRow)
153 track->fLastPoint = track->fRow;
154
155
156 track->fClusterPointer[track->fRow] = c;
1c53abe2 157 //
158
1c53abe2 159 Float_t angle2 = track->GetSnp()*track->GetSnp();
160 angle2 = TMath::Sqrt(angle2/(1-angle2));
161 //
162 //SET NEW Track Point
163 //
91162307 164 // if (debug)
165 {
166 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->fRow));
1c53abe2 167 //
91162307 168 point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
169 point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
170 point.SetErrY(sqrt(track->fErrorY2));
171 point.SetErrZ(sqrt(track->fErrorZ2));
1c53abe2 172 //
91162307 173 point.SetX(track->GetX());
174 point.SetY(track->GetY());
175 point.SetZ(track->GetZ());
176 point.SetAngleY(angle2);
177 point.SetAngleZ(track->GetTgl());
178 if (point.fIsShared){
179 track->fErrorY2 *= 4;
180 track->fErrorZ2 *= 4;
181 }
182 }
183
184 Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
185 //
186 track->fErrorY2 *= 1.3;
187 track->fErrorY2 += 0.01;
188 track->fErrorZ2 *= 1.3;
189 track->fErrorZ2 += 0.005;
190 //}
191 if (accept>0) return 0;
192 if (track->GetNumberOfClusters()%20==0){
193 // if (track->fHelixIn){
194 // TClonesArray & larr = *(track->fHelixIn);
195 // Int_t ihelix = larr.GetEntriesFast();
196 // new(larr[ihelix]) AliHelix(*track) ;
197 //}
1c53abe2 198 }
91162307 199 track->fNoCluster =0;
200 return track->Update(c,chi2,i);
201}
202
203
204
205Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
206 Float_t cory, Float_t corz)
207{
1c53abe2 208 //
91162307 209 // decide according desired precision to accept given
210 // cluster for tracking
211 Double_t sy2=ErrY2(seed,cluster)*cory;
212 Double_t sz2=ErrZ2(seed,cluster)*corz;
213 //sy2=ErrY2(seed,cluster)*cory;
214 //sz2=ErrZ2(seed,cluster)*cory;
1c53abe2 215
91162307 216 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
217 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
218
219 Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
220 (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
221 Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
222 (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
223
224 Double_t rdistance2 = rdistancey2+rdistancez2;
225 //Int_t accept =0;
1c53abe2 226
91162307 227 if (rdistance2>16) return 3;
228
229
230 if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
231 return 2; //suspisiouce - will be changed
232
233 if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
234 // strict cut on overlaped cluster
235 return 2; //suspisiouce - will be changed
236
237 if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
238 && cluster->GetType()<0){
239 seed->fNFoundable--;
240 return 2;
1c53abe2 241 }
91162307 242 return 0;
243}
244
245
1c53abe2 246
1c53abe2 247
1c53abe2 248//_____________________________________________________________________________
f8aae377 249AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
1c53abe2 250AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
251{
252 //---------------------------------------------------------------------
253 // The main TPC tracker constructor
254 //---------------------------------------------------------------------
255 fInnerSec=new AliTPCSector[fkNIS];
256 fOuterSec=new AliTPCSector[fkNOS];
91162307 257
1c53abe2 258 Int_t i;
259 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
260 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
261
262 fN=0; fSectors=0;
263
1c53abe2 264 fSeeds=0;
265 fNtracks = 0;
91162307 266 fParam = par;
267 Int_t nrowlow = par->GetNRowLow();
268 Int_t nrowup = par->GetNRowUp();
269
270
271 for (Int_t i=0;i<nrowlow;i++){
272 fXRow[i] = par->GetPadRowRadiiLow(i);
273 fPadLength[i]= par->GetPadPitchLength(0,i);
274 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
275 }
276
277
278 for (Int_t i=0;i<nrowup;i++){
279 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
280 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
281 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
282 }
283 fSeeds=0;
284 //
285 fInput = 0;
286 fOutput = 0;
287 fSeedTree = 0;
288 fTreeDebug =0;
289 fNewIO =0;
290 fDebug =0;
291 fEvent =0;
1c53abe2 292}
293
294//_____________________________________________________________________________
295AliTPCtrackerMI::~AliTPCtrackerMI() {
296 //------------------------------------------------------------------
297 // TPC tracker destructor
298 //------------------------------------------------------------------
299 delete[] fInnerSec;
300 delete[] fOuterSec;
301 if (fSeeds) {
302 fSeeds->Delete();
303 delete fSeeds;
304 }
305}
306
91162307 307void AliTPCtrackerMI::SetIO()
308{
1c53abe2 309 //
91162307 310 fNewIO = kTRUE;
311 fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::fgkDefaultEventFolderName);
312 fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::fgkDefaultEventFolderName);
313 AliTPCtrack *iotrack= new AliTPCtrack;
314 // iotrack->fHelixIn = new TClonesArray("AliHelix");
315 //iotrack->fHelixOut = new TClonesArray("AliHelix");
316 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
317 delete iotrack;
318}
1c53abe2 319
91162307 320void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
321{
1c53abe2 322
91162307 323 // set input
324 fNewIO = kFALSE;
325 fInput = 0;
326 fOutput = 0;
327 fSeedTree = 0;
328 fTreeDebug =0;
329 fInput = input;
330 if (input==0){
331 return;
332 }
333 //set output
334 fOutput = output;
335 if (output){
336 AliTPCtrack *iotrack= new AliTPCtrack;
337 // iotrack->fHelixIn = new TClonesArray("AliHelix");
338 //iotrack->fHelixOut = new TClonesArray("AliHelix");
339 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
340 delete iotrack;
341 }
342 if (output && (fDebug&2)){
343 //write the full seed information if specified in debug mode
344 //
345 fSeedTree = new TTree("Seeds","Seeds");
346 AliTPCseed * vseed = new AliTPCseed;
347 //
348 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
349 arrtr->ExpandCreateFast(160);
350 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
351 //
352 vseed->fPoints = arrtr;
353 vseed->fEPoints = arre;
354 // vseed->fClusterPoints = arrcl;
355 fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
356 delete arrtr;
357 delete arre;
358 fTreeDebug = new TTree("trackDebug","trackDebug");
359 TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
360 fTreeDebug->Branch("debug",&arrd,32000,99);
1c53abe2 361 }
1c53abe2 362
1c53abe2 363
91162307 364 //set ESD event
365 fEvent = event;
366}
1c53abe2 367
91162307 368void AliTPCtrackerMI::WriteTracks()
369{
370 //
371 // write tracks to the given output tree -
372 // output specified with SetIO routine
373 if (!fSeeds) return;
374 if (fEvent){
375 // write tracks to the event
376 // store index of the track
377 Int_t nseed=fSeeds->GetEntriesFast();
378 for (Int_t i=0; i<nseed; i++) {
379 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
380 if (!pt) continue;
381 AliESDtrack iotrack;
382 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
383 //iotrack.SetTPCindex(i);
384 fEvent->AddTrack(&iotrack);
385 }
1c53abe2 386 }
387
1c53abe2 388
91162307 389 if (fOutput){
390 AliTPCtrack *iotrack= 0;
391 Int_t nseed=fSeeds->GetEntriesFast();
392 for (Int_t i=0; i<nseed; i++) {
393 iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
394 if (iotrack) break;
395 }
396
397 //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
398 TBranch * br = fOutput->GetBranch("tracks");
399 br->SetAddress(&iotrack);
400 //
401 for (Int_t i=0; i<nseed; i++) {
402 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
403 if (!pt) continue;
404 iotrack = pt;
405 pt->fLab2 =i;
406 // br->SetAddress(&iotrack);
407 fOutput->Fill();
408 iotrack =0;
409 }
410 }
411 // delete iotrack;
412 //
413 if (fSeedTree){
414 //write the full seed information if specified in debug mode
415
416 AliTPCseed * vseed = new AliTPCseed;
417 //
418 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
419 arrtr->ExpandCreateFast(160);
420 //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
421 //arrcl->ExpandCreateFast(160);
422 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
423 //
424 vseed->fPoints = arrtr;
425 vseed->fEPoints = arre;
426 // vseed->fClusterPoints = arrcl;
427 //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
428 TBranch * brseed = fSeedTree->GetBranch("seeds");
429
430 Int_t nseed=fSeeds->GetEntriesFast();
431
432 for (Int_t i=0; i<nseed; i++) {
433 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
434 if (!pt) continue;
435 pt->fPoints = arrtr;
436 // pt->fClusterPoints = arrcl;
437 pt->fEPoints = arre;
438 pt->RebuildSeed();
439 vseed = pt;
440 brseed->SetAddress(&vseed);
441 fSeedTree->Fill();
442 pt->fPoints = 0;
443 pt->fEPoints = 0;
444 // pt->fClusterPoints = 0;
445 }
446 fSeedTree->Write();
447 if (fTreeDebug) fTreeDebug->Write();
448 }
1c53abe2 449
91162307 450}
1c53abe2 451
1c53abe2 452
1c53abe2 453
454
91162307 455Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
456 //
457 //
458 //seed->SetErrorY2(0.1);
459 //return 0.1;
460 //calculate look-up table at the beginning
461 static Bool_t ginit = kFALSE;
462 static Float_t gnoise1,gnoise2,gnoise3;
463 static Float_t ggg1[10000];
464 static Float_t ggg2[10000];
465 static Float_t ggg3[10000];
466 static Float_t glandau1[10000];
467 static Float_t glandau2[10000];
468 static Float_t glandau3[10000];
469 //
470 static Float_t gcor01[500];
471 static Float_t gcor02[500];
472 static Float_t gcorp[500];
473 //
1c53abe2 474
91162307 475 //
476 if (ginit==kFALSE){
477 for (Int_t i=1;i<500;i++){
478 Float_t rsigma = float(i)/100.;
479 gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
480 gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
481 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
482 }
1c53abe2 483
91162307 484 //
485 for (Int_t i=3;i<10000;i++){
486 //
487 //
488 // inner sector
489 Float_t amp = float(i);
490 Float_t padlength =0.75;
491 gnoise1 = 0.0004/padlength;
492 Float_t nel = 0.268*amp;
493 Float_t nprim = 0.155*amp;
494 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
495 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
496 if (glandau1[i]>1) glandau1[i]=1;
497 glandau1[i]*=padlength*padlength/12.;
498 //
499 // outer short
500 padlength =1.;
501 gnoise2 = 0.0004/padlength;
502 nel = 0.3*amp;
503 nprim = 0.133*amp;
504 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
505 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
506 if (glandau2[i]>1) glandau2[i]=1;
507 glandau2[i]*=padlength*padlength/12.;
508 //
509 //
510 // outer long
511 padlength =1.5;
512 gnoise3 = 0.0004/padlength;
513 nel = 0.3*amp;
514 nprim = 0.133*amp;
515 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
516 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
517 if (glandau3[i]>1) glandau3[i]=1;
518 glandau3[i]*=padlength*padlength/12.;
519 //
520 }
521 ginit = kTRUE;
522 }
1c53abe2 523 //
524 //
91162307 525 //
526 Int_t amp = int(TMath::Abs(cl->GetQ()));
527 if (amp>9999) {
528 seed->SetErrorY2(1.);
529 return 1.;
530 }
1c53abe2 531 Float_t snoise2;
c9427e08 532 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
91162307 533 Int_t ctype = cl->GetType();
534 Float_t padlength= GetPadPitchLength(seed->fRow);
535 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
536 angle2 = angle2/(1-angle2);
537 //
538 //cluster "quality"
539 Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
540 Float_t res;
1c53abe2 541 //
1c53abe2 542 if (fSectors==fInnerSec){
91162307 543 snoise2 = gnoise1;
544 res = ggg1[amp]*z+glandau1[amp]*angle2;
545 if (ctype==0) res *= gcor01[rsigmay];
546 if ((ctype>0)){
547 res+=0.002;
548 res*= gcorp[rsigmay];
549 }
1c53abe2 550 }
551 else {
91162307 552 if (padlength<1.1){
553 snoise2 = gnoise2;
554 res = ggg2[amp]*z+glandau2[amp]*angle2;
555 if (ctype==0) res *= gcor02[rsigmay];
556 if ((ctype>0)){
557 res+=0.002;
558 res*= gcorp[rsigmay];
559 }
560 }
561 else{
562 snoise2 = gnoise3;
563 res = ggg3[amp]*z+glandau3[amp]*angle2;
564 if (ctype==0) res *= gcor02[rsigmay];
565 if ((ctype>0)){
566 res+=0.002;
567 res*= gcorp[rsigmay];
568 }
569 }
570 }
1c53abe2 571
91162307 572 if (ctype<0){
573 res+=0.005;
574 res*=2.4; // overestimate error 2 times
575 }
576 res+= snoise2;
577
1c53abe2 578 if (res<2*snoise2)
91162307 579 res = 2*snoise2;
580
581 seed->SetErrorY2(res);
1c53abe2 582 return res;
1c53abe2 583
584
91162307 585}
c9427e08 586
587
588
91162307 589Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
590 //
591 //
592 //seed->SetErrorY2(0.1);
593 //return 0.1;
594 //calculate look-up table at the beginning
595 static Bool_t ginit = kFALSE;
596 static Float_t gnoise1,gnoise2,gnoise3;
597 static Float_t ggg1[10000];
598 static Float_t ggg2[10000];
599 static Float_t ggg3[10000];
600 static Float_t glandau1[10000];
601 static Float_t glandau2[10000];
602 static Float_t glandau3[10000];
603 //
604 static Float_t gcor01[1000];
605 static Float_t gcor02[1000];
606 static Float_t gcorp[1000];
1627d1c4 607 //
1627d1c4 608
91162307 609 //
610 if (ginit==kFALSE){
611 for (Int_t i=1;i<1000;i++){
612 Float_t rsigma = float(i)/100.;
613 gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
614 gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
615 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
616 }
1c53abe2 617
91162307 618 //
619 for (Int_t i=3;i<10000;i++){
620 //
621 //
622 // inner sector
623 Float_t amp = float(i);
624 Float_t padlength =0.75;
625 gnoise1 = 0.0004/padlength;
626 Float_t nel = 0.268*amp;
627 Float_t nprim = 0.155*amp;
628 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
629 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
630 if (glandau1[i]>1) glandau1[i]=1;
631 glandau1[i]*=padlength*padlength/12.;
632 //
633 // outer short
634 padlength =1.;
635 gnoise2 = 0.0004/padlength;
636 nel = 0.3*amp;
637 nprim = 0.133*amp;
638 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
639 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
640 if (glandau2[i]>1) glandau2[i]=1;
641 glandau2[i]*=padlength*padlength/12.;
642 //
643 //
644 // outer long
645 padlength =1.5;
646 gnoise3 = 0.0004/padlength;
647 nel = 0.3*amp;
648 nprim = 0.133*amp;
649 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
650 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
651 if (glandau3[i]>1) glandau3[i]=1;
652 glandau3[i]*=padlength*padlength/12.;
653 //
654 }
655 ginit = kTRUE;
656 }
657 //
658 //
659 //
660 Int_t amp = int(TMath::Abs(cl->GetQ()));
661 if (amp>9999) {
662 seed->SetErrorY2(1.);
663 return 1.;
664 }
665 Float_t snoise2;
666 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
667 Int_t ctype = cl->GetType();
668 Float_t padlength= GetPadPitchLength(seed->fRow);
669 //
670 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
671 // if (angle2<0.6) angle2 = 0.6;
672 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
673 //
674 //cluster "quality"
675 Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
676 Float_t res;
677 //
678 if (fSectors==fInnerSec){
679 snoise2 = gnoise1;
680 res = ggg1[amp]*z+glandau1[amp]*angle2;
681 if (ctype==0) res *= gcor01[rsigmaz];
682 if ((ctype>0)){
683 res+=0.002;
684 res*= gcorp[rsigmaz];
685 }
686 }
687 else {
688 if (padlength<1.1){
689 snoise2 = gnoise2;
690 res = ggg2[amp]*z+glandau2[amp]*angle2;
691 if (ctype==0) res *= gcor02[rsigmaz];
692 if ((ctype>0)){
693 res+=0.002;
694 res*= gcorp[rsigmaz];
695 }
696 }
697 else{
698 snoise2 = gnoise3;
699 res = ggg3[amp]*z+glandau3[amp]*angle2;
700 if (ctype==0) res *= gcor02[rsigmaz];
701 if ((ctype>0)){
702 res+=0.002;
703 res*= gcorp[rsigmaz];
704 }
705 }
706 }
707
708 if (ctype<0){
709 res+=0.002;
710 res*=1.3;
711 }
712 if ((ctype<0) &&amp<70){
713 res+=0.002;
714 res*=1.3;
715 }
716 res += snoise2;
717 if (res<2*snoise2)
718 res = 2*snoise2;
719 if (res>3) res =3;
720 seed->SetErrorZ2(res);
721 return res;
722}
723
724
725
726/*
727Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
728 //
729 //
730 //seed->SetErrorZ2(0.1);
731 //return 0.1;
732
733 Float_t snoise2;
734 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
735 //
736 Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
737 Int_t ctype = cl->GetType();
738 Float_t amp = TMath::Abs(cl->GetQ());
739
740 Float_t nel;
741 Float_t nprim;
742 //
743 Float_t landau=2 ; //landau fluctuation part
744 Float_t gg=2; // gg fluctuation part
745 Float_t padlength= GetPadPitchLength(seed->GetX());
746
747 if (fSectors==fInnerSec){
748 snoise2 = 0.0004/padlength;
749 nel = 0.268*amp;
750 nprim = 0.155*amp;
751 gg = (2+0.001*nel/(padlength*padlength))/nel;
752 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
753 if (landau>1) landau=1;
754 }
755 else {
756 snoise2 = 0.0004/padlength;
757 nel = 0.3*amp;
758 nprim = 0.133*amp;
759 gg = (2+0.0008*nel/(padlength*padlength))/nel;
760 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
761 if (landau>1) landau=1;
762 }
763 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
764
765 //
766 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
767 angle2 = TMath::Sqrt((1-angle2));
768 if (angle2<0.6) angle2 = 0.6;
769 //angle2 = 1;
770
771 Float_t angle = seed->GetTgl()/angle2;
772 Float_t angular = landau*angle*angle*padlength*padlength/12.;
773 Float_t res = sdiff + angular;
774
775
776 if ((ctype==0) && (fSectors ==fOuterSec))
777 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
778
779 if ((ctype==0) && (fSectors ==fInnerSec))
780 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
781
782 if ((ctype>0)){
783 res+=0.005;
784 res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype;
785 }
786 if (ctype<0){
787 res+=0.002;
788 res*=1.3;
789 }
790 if ((ctype<0) &&amp<70){
791 res+=0.002;
792 res*=1.3;
793 }
794 res += snoise2;
795 if (res<2*snoise2)
796 res = 2*snoise2;
797
798 seed->SetErrorZ2(res);
799 return res;
800}
801*/
802
803
804
805void AliTPCseed::Reset(Bool_t all)
806{
807 //
808 //
809 SetNumberOfClusters(0);
810 fNFoundable = 0;
811 SetChi2(0);
812 ResetCovariance();
813 /*
814 if (fTrackPoints){
815 for (Int_t i=0;i<8;i++){
816 delete [] fTrackPoints[i];
817 }
818 delete fTrackPoints;
819 fTrackPoints =0;
820 }
821 */
822
823 if (all){
824 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
825 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
826 }
827
828}
829
830
831void AliTPCseed::Modify(Double_t factor)
832{
833
834 //------------------------------------------------------------------
835 //This function makes a track forget its history :)
836 //------------------------------------------------------------------
837 if (factor<=0) {
838 ResetCovariance();
839 return;
840 }
841 fC00*=factor;
842 fC10*=factor; fC11*=factor;
843 fC20*=factor; fC21*=factor; fC22*=factor;
844 fC30*=factor; fC31*=factor; fC32*=factor; fC33*=factor;
845 fC40*=factor; fC41*=factor; fC42*=factor; fC43*=factor; fC44*=factor;
846 SetNumberOfClusters(0);
847 fNFoundable =0;
848 SetChi2(0);
849 fRemoval = 0;
850 fCurrentSigmaY2 = 0.000005;
851 fCurrentSigmaZ2 = 0.000005;
852 fNoCluster = 0;
853 //fFirstPoint = 160;
854 //fLastPoint = 0;
855}
856
857
858
859
860Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
861{
862 //-----------------------------------------------------------------
863 // This function find proloncation of a track to a reference plane x=xk.
864 // doesn't change internal state of the track
865 //-----------------------------------------------------------------
866
1c53abe2 867 Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
91162307 868
869 if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
870 return 0;
871 }
872
1c53abe2 873 // Double_t y1=fP0, z1=fP1;
874 Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
875 Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
876
877 y = fP0;
878 z = fP1;
91162307 879 //y += dx*(c1+c2)/(r1+r2);
880 //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
881
882 Double_t dy = dx*(c1+c2)/(r1+r2);
883 Double_t dz = 0;
884 //
885 Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
886 /*
887 if (TMath::Abs(delta)>0.0001){
888 dz = fP3*TMath::ASin(delta)/fP4;
889 }else{
890 dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
891 }
892 */
893 dz = fP3*TPCFastMath::FastAsin(delta)/fP4;
894 //
895 y+=dy;
896 z+=dz;
897
898
899 return 1;
1c53abe2 900}
901
902
903//_____________________________________________________________________________
904Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
905{
906 //-----------------------------------------------------------------
907 // This function calculates a predicted chi2 increment.
908 //-----------------------------------------------------------------
909 //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
910 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
911 r00+=fC00; r01+=fC10; r11+=fC11;
912
913 Double_t det=r00*r11 - r01*r01;
914 if (TMath::Abs(det) < 1.e-10) {
915 Int_t n=GetNumberOfClusters();
916 if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
917 return 1e10;
918 }
919 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
920
921 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
922
923 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
924}
925
926
927//_________________________________________________________________________________________
928
929
930Int_t AliTPCseed::Compare(const TObject *o) const {
931 //-----------------------------------------------------------------
932 // This function compares tracks according to the sector - for given sector according z
933 //-----------------------------------------------------------------
934 AliTPCseed *t=(AliTPCseed*)o;
1c53abe2 935
91162307 936 if (fSort == 0){
937 if (t->fRelativeSector>fRelativeSector) return -1;
938 if (t->fRelativeSector<fRelativeSector) return 1;
939 Double_t z2 = t->GetZ();
940 Double_t z1 = GetZ();
941 if (z2>z1) return 1;
942 if (z2<z1) return -1;
943 return 0;
944 }
945 else {
946 Float_t f2 =1;
947 f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
948 if (t->fBConstrain) f2=1.2;
949
950 Float_t f1 =1;
951 f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
952
953 if (fBConstrain) f1=1.2;
954
955 if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
956 else return +1;
957 }
1c53abe2 958}
959
c9427e08 960void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
961{
962 //rotate to track "local coordinata
963 Float_t x = seed->GetX();
964 Float_t y = seed->GetY();
965 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
91162307 966
c9427e08 967 if (y > ymax) {
968 seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
969 if (!seed->Rotate(fSectors->GetAlpha()))
970 return;
971 } else if (y <-ymax) {
972 seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
973 if (!seed->Rotate(-fSectors->GetAlpha()))
974 return;
975 }
1c53abe2 976
c9427e08 977}
1c53abe2 978
979
980
981
982//_____________________________________________________________________________
176aff27 983Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
1c53abe2 984 //-----------------------------------------------------------------
985 // This function associates a cluster with this track.
986 //-----------------------------------------------------------------
1c53abe2 987 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
988
989 r00+=fC00; r01+=fC10; r11+=fC11;
990 Double_t det=r00*r11 - r01*r01;
991 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
992
993 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
994 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
995 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
996 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
997 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
998
999 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1000 Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1627d1c4 1001 if (TMath::Abs(cur*fX-eta) >= 0.9) {
1c53abe2 1002 return 0;
1003 }
1004
1005 fP0 += k00*dy + k01*dz;
1006 fP1 += k10*dy + k11*dz;
1007 fP2 = eta;
1008 fP3 += k30*dy + k31*dz;
1009 fP4 = cur;
1010
1011 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1012 Double_t c12=fC21, c13=fC31, c14=fC41;
1013
1014 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1015 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
1016 fC40-=k00*c04+k01*c14;
1017
1018 fC11-=k10*c01+k11*fC11;
1019 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
1020 fC41-=k10*c04+k11*c14;
1021
1022 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
1023 fC42-=k20*c04+k21*c14;
1024
1025 fC33-=k30*c03+k31*c13;
1026 fC43-=k40*c03+k41*c13;
1027
1028 fC44-=k40*c04+k41*c14;
1029
1030 Int_t n=GetNumberOfClusters();
91162307 1031 // fIndex[n]=index;
1c53abe2 1032 SetNumberOfClusters(n+1);
1033 SetChi2(GetChi2()+chisq);
1034
1035 return 1;
1036}
1037
1038
1039
1040//_____________________________________________________________________________
91162307 1041Double_t AliTPCtrackerMI::f1old(Double_t x1,Double_t y1,
1c53abe2 1042 Double_t x2,Double_t y2,
1043 Double_t x3,Double_t y3)
1044{
1045 //-----------------------------------------------------------------
1046 // Initial approximation of the track curvature
1047 //-----------------------------------------------------------------
1048 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1049 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1050 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1051 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1052 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1053
1054 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
91162307 1055 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1c53abe2 1056 return -xr*yr/sqrt(xr*xr+yr*yr);
1057}
1058
1059
91162307 1060
1c53abe2 1061//_____________________________________________________________________________
91162307 1062Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
1063 Double_t x2,Double_t y2,
1064 Double_t x3,Double_t y3)
1065{
1066 //-----------------------------------------------------------------
1067 // Initial approximation of the track curvature
1068 //-----------------------------------------------------------------
1069 x3 -=x1;
1070 x2 -=x1;
1071 y3 -=y1;
1072 y2 -=y1;
1073 //
1074 Double_t det = x3*y2-x2*y3;
1075 if (det==0) {
1076 return 100;
1077 }
1078 //
1079 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1080 Double_t x0 = x3*0.5-y3*u;
1081 Double_t y0 = y3*0.5+x3*u;
1082 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1083 if (det<0) c2*=-1;
1084 return c2;
1085}
1086
1087
1c53abe2 1088Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
1089 Double_t x2,Double_t y2,
1090 Double_t x3,Double_t y3)
91162307 1091{
1092 //-----------------------------------------------------------------
1093 // Initial approximation of the track curvature
1094 //-----------------------------------------------------------------
1095 x3 -=x1;
1096 x2 -=x1;
1097 y3 -=y1;
1098 y2 -=y1;
1099 //
1100 Double_t det = x3*y2-x2*y3;
1101 if (det==0) {
1102 return 100;
1103 }
1104 //
1105 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1106 Double_t x0 = x3*0.5-y3*u;
1107 Double_t y0 = y3*0.5+x3*u;
1108 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1109 if (det<0) c2*=-1;
1110 x0+=x1;
1111 x0*=c2;
1112 return x0;
1113}
1114
1115
1116
1117//_____________________________________________________________________________
1118Double_t AliTPCtrackerMI::f2old(Double_t x1,Double_t y1,
1119 Double_t x2,Double_t y2,
1120 Double_t x3,Double_t y3)
1c53abe2 1121{
1122 //-----------------------------------------------------------------
1123 // Initial approximation of the track curvature times center of curvature
1124 //-----------------------------------------------------------------
1125 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1126 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1127 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1128 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1129 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1130
1131 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1132
1133 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1134}
1135
1136//_____________________________________________________________________________
1137Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1,
1138 Double_t x2,Double_t y2,
1139 Double_t z1,Double_t z2)
1140{
1141 //-----------------------------------------------------------------
1142 // Initial approximation of the tangent of the track dip angle
1143 //-----------------------------------------------------------------
1144 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1145}
1146
1147
91162307 1148Double_t AliTPCtrackerMI::f3n(Double_t x1,Double_t y1,
1149 Double_t x2,Double_t y2,
1150 Double_t z1,Double_t z2, Double_t c)
1c53abe2 1151{
91162307 1152 //-----------------------------------------------------------------
1153 // Initial approximation of the tangent of the track dip angle
1154 //-----------------------------------------------------------------
1155
1156 // Double_t angle1;
1157
1158 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1c53abe2 1159 //
91162307 1160 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1161 if (TMath::Abs(d*c*0.5)>1) return 0;
1162 // Double_t angle2 = TMath::ASin(d*c*0.5);
1163 // Double_t angle2 = TPCFastMath::FastAsin(d*c*0.5);
1164 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): TPCFastMath::FastAsin(d*c*0.5);
1165
1166 angle2 = (z1-z2)*c/(angle2*2.);
1167 return angle2;
1168}
1169
1170Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1171{//-----------------------------------------------------------------
1172 // This function find proloncation of a track to a reference plane x=x2.
1173 //-----------------------------------------------------------------
1174
1175 Double_t dx=x2-x1;
1176
1177 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1178 return kFALSE;
1c53abe2 1179 }
f8aae377 1180
91162307 1181 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1182 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1183 y = x[0];
1184 z = x[1];
1185
1186 Double_t dy = dx*(c1+c2)/(r1+r2);
1187 Double_t dz = 0;
1188 //
1189 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1190
1191 if (TMath::Abs(delta)>0.01){
1192 dz = x[3]*TMath::ASin(delta)/x[4];
1193 }else{
1194 dz = x[3]*TPCFastMath::FastAsin(delta)/x[4];
1195 }
1196
1197 //dz = x[3]*TPCFastMath::FastAsin(delta)/x[4];
f8aae377 1198
91162307 1199 y+=dy;
1200 z+=dz;
1201
1202 return kTRUE;
1c53abe2 1203}
1204
91162307 1205
1206
1207Int_t AliTPCtrackerMI::LoadClusters()
1c53abe2 1208{
1209 //
1210 // load clusters to the memory
91162307 1211 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1212 clrow->SetClass("AliTPCclusterMI");
1213 clrow->SetArray(0);
1214 clrow->GetArray()->ExpandCreateFast(10000);
1215 //
1216 // TTree * tree = fClustersArray.GetTree();
1217
1218 TTree * tree = fInput;
1219 TBranch * br = tree->GetBranch("Segment");
1220 br->SetAddress(&clrow);
1221 //
1222 Int_t j=Int_t(tree->GetEntries());
1c53abe2 1223 for (Int_t i=0; i<j; i++) {
91162307 1224 br->GetEntry(i);
1225 //
1226 Int_t sec,row;
1227 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1228 //
1229 AliTPCRow * tpcrow=0;
1230 Int_t left=0;
1231 if (sec<fkNIS*2){
1232 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1233 left = sec/fkNIS;
1234 }
1235 else{
1236 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1237 left = (sec-fkNIS*2)/fkNOS;
1238 }
1239 if (left ==0){
1240 tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1241 tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1242 for (Int_t i=0;i<tpcrow->fN1;i++)
1243 tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1244 }
1245 if (left ==1){
1246 tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1247 tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1248 for (Int_t i=0;i<tpcrow->fN2;i++)
1249 tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1250 }
1c53abe2 1251 }
91162307 1252 //
1253 delete clrow;
1254 LoadOuterSectors();
1255 LoadInnerSectors();
1256 return 0;
1c53abe2 1257}
1258
1259
91162307 1260void AliTPCtrackerMI::UnloadClusters()
1261{
1262 //
1263 // unload clusters from the memory
1264 //
1265 Int_t nrows = fOuterSec->GetNRows();
1266 for (Int_t sec = 0;sec<fkNOS;sec++)
1267 for (Int_t row = 0;row<nrows;row++){
1268 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1269 if (tpcrow){
1270 if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1271 if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1272 }
1c53abe2 1273 }
91162307 1274 //
1275 nrows = fInnerSec->GetNRows();
1276 for (Int_t sec = 0;sec<fkNIS;sec++)
1277 for (Int_t row = 0;row<nrows;row++){
1278 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1279 if (tpcrow){
1280 if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1281 if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1282 }
1283 }
1284
1285 return ;
1c53abe2 1286}
1287
1288
1289//_____________________________________________________________________________
91162307 1290Int_t AliTPCtrackerMI::LoadOuterSectors() {
1c53abe2 1291 //-----------------------------------------------------------------
91162307 1292 // This function fills outer TPC sectors with clusters.
1c53abe2 1293 //-----------------------------------------------------------------
91162307 1294 Int_t nrows = fOuterSec->GetNRows();
1295 UInt_t index=0;
1296 for (Int_t sec = 0;sec<fkNOS;sec++)
1297 for (Int_t row = 0;row<nrows;row++){
1298 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1299 Int_t sec2 = sec+2*fkNIS;
1300 //left
1301 Int_t ncl = tpcrow->fN1;
1302 while (ncl--) {
1303 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1304 index=(((sec2<<8)+row)<<16)+ncl;
1305 tpcrow->InsertCluster(c,index);
1306 }
1307 //right
1308 ncl = tpcrow->fN2;
1309 while (ncl--) {
1310 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1311 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1312 tpcrow->InsertCluster(c,index);
1313 }
1314 //
1315 // write indexes for fast acces
1316 //
1317 for (Int_t i=0;i<510;i++)
1318 tpcrow->fFastCluster[i]=-1;
1319 for (Int_t i=0;i<tpcrow->GetN();i++){
1320 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1321 tpcrow->fFastCluster[zi]=i; // write index
1322 }
1323 Int_t last = 0;
1324 for (Int_t i=0;i<510;i++){
1325 if (tpcrow->fFastCluster[i]<0)
1326 tpcrow->fFastCluster[i] = last;
1327 else
1328 last = tpcrow->fFastCluster[i];
1329 }
1330 }
1331 fN=fkNOS;
1332 fSectors=fOuterSec;
1333 return 0;
1334}
1335
1336
1337//_____________________________________________________________________________
1338Int_t AliTPCtrackerMI::LoadInnerSectors() {
1339 //-----------------------------------------------------------------
1340 // This function fills inner TPC sectors with clusters.
1341 //-----------------------------------------------------------------
1342 Int_t nrows = fInnerSec->GetNRows();
1343 UInt_t index=0;
1344 for (Int_t sec = 0;sec<fkNIS;sec++)
1345 for (Int_t row = 0;row<nrows;row++){
1346 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1347 //
1348 //left
1349 Int_t ncl = tpcrow->fN1;
1350 while (ncl--) {
1351 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1352 index=(((sec<<8)+row)<<16)+ncl;
1353 tpcrow->InsertCluster(c,index);
1354 }
1355 //right
1356 ncl = tpcrow->fN2;
1357 while (ncl--) {
1358 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1359 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1360 tpcrow->InsertCluster(c,index);
1361 }
1362 //
1363 // write indexes for fast acces
1364 //
1365 for (Int_t i=0;i<510;i++)
1366 tpcrow->fFastCluster[i]=-1;
1367 for (Int_t i=0;i<tpcrow->GetN();i++){
1368 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1369 tpcrow->fFastCluster[zi]=i; // write index
1370 }
1371 Int_t last = 0;
1372 for (Int_t i=0;i<510;i++){
1373 if (tpcrow->fFastCluster[i]<0)
1374 tpcrow->fFastCluster[i] = last;
1375 else
1376 last = tpcrow->fFastCluster[i];
1377 }
1c53abe2 1378
91162307 1379 }
1380
1c53abe2 1381 fN=fkNIS;
1382 fSectors=fInnerSec;
91162307 1383 return 0;
1384}
1385
1386
1387
1388//_________________________________________________________________________
1389AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1390 //--------------------------------------------------------------------
1391 // Return pointer to a given cluster
1392 //--------------------------------------------------------------------
1393 Int_t sec=(index&0xff000000)>>24;
1394 Int_t row=(index&0x00ff0000)>>16;
1395 Int_t ncl=(index&0x0000ffff)>>00;
1396
1397 const AliTPCRow * tpcrow=0;
1398 AliTPCclusterMI * clrow =0;
1399 if (sec<fkNIS*2){
1400 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1401 if (sec<fkNIS)
1402 clrow = tpcrow->fClusters1;
1403 else
1404 clrow = tpcrow->fClusters2;
1405 }
1406 else{
1407 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1408 if (sec-2*fkNIS<fkNOS)
1409 clrow = tpcrow->fClusters1;
1410 else
1411 clrow = tpcrow->fClusters2;
1412 }
1413 if (tpcrow==0) return 0;
1414 if (tpcrow->GetN()<=ncl) return 0;
1415 // return (AliTPCclusterMI*)(*tpcrow)[ncl];
1416 return &(clrow[ncl]);
1417
1c53abe2 1418}
1419
91162307 1420
1421
1c53abe2 1422Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1423 //-----------------------------------------------------------------
1424 // This function tries to find a track prolongation to next pad row
1425 //-----------------------------------------------------------------
1c53abe2 1426 //
91162307 1427 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1627d1c4 1428
91162307 1429 // if (t.GetRadius()>x+10 ) return 0;
1430 // t.PropagateTo(x+0.02);
1431 //t.PropagateTo(x+0.01);
1627d1c4 1432 if (!t.PropagateTo(x)) {
91162307 1433 t.fRemoval = 10;
1627d1c4 1434 return 0;
1435 }
91162307 1436 //
1437 Double_t y=t.GetY(), z=t.GetZ();
1438 if (TMath::Abs(y)>ymax){
1439 if (y > ymax) {
1440 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1441 if (!t.Rotate(fSectors->GetAlpha()))
1442 return 0;
1443 } else if (y <-ymax) {
1444 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1445 if (!t.Rotate(-fSectors->GetAlpha()))
1446 return 0;
1447 }
1448 if (!t.PropagateTo(x)) {
1449 return 0;
1450 }
1451 y=t.GetY();
1452 }
1453 //
1454 // update current shape info every 3 pad-row
1455 if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1456 //t.fCurrentSigmaY = GetSigmaY(&t);
1457 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1458 GetShape(&t,nr);
1459 }
1c53abe2 1460 //
1461 AliTPCclusterMI *cl=0;
1462 UInt_t index=0;
91162307 1463
1464
1465 //Int_t nr2 = nr;
1466 if (t.GetClusterIndex2(nr)>0){
1467 //
1468 //cl = GetClusterMI(t.GetClusterIndex2(nr));
1469 index = t.GetClusterIndex2(nr);
1470 cl = t.fClusterPointer[nr];
1471 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1472 t.fCurrentClusterIndex1 = index;
1473 }
1474
1475 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1476 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1477 Double_t roady =1.;
1478 Double_t roadz = 1.;
1479 //
1c53abe2 1480 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1481 t.fInDead = kTRUE;
91162307 1482 t.SetClusterIndex2(nr,-1);
1c53abe2 1483 return 0;
1484 }
1485 else
1486 {
1627d1c4 1487 if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
1488 else
1489 return 0;
1c53abe2 1490 }
1491 //calculate
91162307 1492 if (cl){
c9427e08 1493 t.fCurrentCluster = cl;
91162307 1494 t.fRow = nr;
1495 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1496 if (accept<3) {
1497 //if founded cluster is acceptible
1498 UpdateTrack(&t,accept);
1499 return 1;
1500 }
1501 }
c9427e08 1502
91162307 1503 if (krow) {
1504 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1505 cl = krow.FindNearest2(y,z,roady,roadz,index);
1506 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1507 }
1508 // t.fNoCluster++;
c9427e08 1509
91162307 1510 if (cl) {
1511 t.fCurrentCluster = cl;
1512 t.fRow = nr;
1513 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
c9427e08 1514
91162307 1515 if (t.fCurrentCluster->IsUsed(10)){
1516 //
1517 //
c9427e08 1518
91162307 1519 t.fNShared++;
1520 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1521 t.fRemoval =10;
1522 return 0;
1523 }
1524 }
1525
1526 if (accept<3) UpdateTrack(&t,accept);
c9427e08 1527
91162307 1528 } else {
1529 if (t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1530
1531 }
1532 return 1;
1533}
c9427e08 1534
91162307 1535Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1536 //-----------------------------------------------------------------
1537 // This function tries to find a track prolongation to next pad row
1538 //-----------------------------------------------------------------
1539 //
1540 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1541 Double_t y,z;
1542 if (!t.GetProlongation(x,y,z)) {
1543 t.fRemoval = 10;
1544 return 0;
1545 }
1546 //
1547 //
1548 if (TMath::Abs(y)>ymax){
1549 return 0;
1550
1c53abe2 1551 if (y > ymax) {
1552 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1553 if (!t.Rotate(fSectors->GetAlpha()))
1554 return 0;
1555 } else if (y <-ymax) {
1556 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1557 if (!t.Rotate(-fSectors->GetAlpha()))
1558 return 0;
91162307 1559 }
1560 if (!t.PropagateTo(x)) {
1561 return 0;
1562 }
1563 t.GetProlongation(x,y,z);
1564 }
1565 //
1566 // update current shape info every 3 pad-row
1567 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1568 // t.fCurrentSigmaY = GetSigmaY(&t);
1569 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1570 GetShape(&t,nr);
1571 }
1572 //
1573 AliTPCclusterMI *cl=0;
1574 UInt_t index=0;
1575
1576
1577 //Int_t nr2 = nr;
1578 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1579 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1580 Double_t roady =1.;
1581 Double_t roadz = 1.;
1582 //
1583 Int_t row = nr;
1584 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1585 t.fInDead = kTRUE;
1586 t.SetClusterIndex2(row,-1);
1587 return 0;
1588 }
1589 else
1590 {
1591 if (TMath::Abs(z)>(1.05*x+10)) t.SetClusterIndex2(row,-1);
1c53abe2 1592 }
91162307 1593 //calculate
1594
1595 if ((cl==0)&&(krow)) {
1596 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1597 cl = krow.FindNearest2(y,z,roady,roadz,index);
1598
1599 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1600 }
1601
1602 if (cl) {
1603 t.fCurrentCluster = cl;
1604 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1605 //if (accept<3){
1606 t.SetClusterIndex2(row,index);
1607 t.fClusterPointer[row] = cl;
1608 //}
1c53abe2 1609 }
1610 return 1;
1611}
1612
1613
91162307 1614
1615Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1c53abe2 1616 //-----------------------------------------------------------------
1617 // This function tries to find a track prolongation to next pad row
1618 //-----------------------------------------------------------------
1619 t.fCurrentCluster = 0;
1620 t.fCurrentClusterIndex1 = 0;
1c53abe2 1621
1622 Double_t xt=t.GetX();
91162307 1623 Int_t row = GetRowNumber(xt)-1;
1624 Double_t ymax= GetMaxY(nr);
1625
1c53abe2 1626 if (row < nr) return 1; // don't prolongate if not information until now -
91162307 1627 if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1628 t.fRemoval =10;
1629 return 0; // not prolongate strongly inclined tracks
1630 }
1631 if (TMath::Abs(t.GetSnp())>0.95) {
1632 t.fRemoval =10;
1633 return 0; // not prolongate strongly inclined tracks
1634 }
1635
1636 Double_t x= GetXrow(nr);
1637 Double_t y,z;
1638 //t.PropagateTo(x+0.02);
1639 //t.PropagateTo(x+0.01);
1627d1c4 1640 if (!t.PropagateTo(x)){
1627d1c4 1641 return 0;
1642 }
1c53abe2 1643 //
91162307 1644 y=t.GetY();
1645 z=t.GetZ();
1c53abe2 1646
91162307 1647 if (TMath::Abs(y)>ymax){
1648 if (y > ymax) {
1649 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1650 if (!t.Rotate(fSectors->GetAlpha()))
1651 return 0;
1652 } else if (y <-ymax) {
1653 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1654 if (!t.Rotate(-fSectors->GetAlpha()))
1655 return 0;
1656 }
1657 if (!t.PropagateTo(x)){
1658 return 0;
1c53abe2 1659 }
91162307 1660 y = t.GetY();
1c53abe2 1661 }
91162307 1662 //
1c53abe2 1663
91162307 1664 AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1c53abe2 1665
1666 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1667 t.fInDead = kTRUE;
91162307 1668 t.SetClusterIndex2(nr,-1);
1c53abe2 1669 return 0;
1670 }
1671 else
1672 {
1627d1c4 1673 if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
1674 else
1675 return 0;
1c53abe2 1676 }
1c53abe2 1677
91162307 1678 // update current
1679 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1680 // t.fCurrentSigmaY = GetSigmaY(&t);
1681 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1682 GetShape(&t,nr);
1683 }
1c53abe2 1684
91162307 1685 AliTPCclusterMI *cl=0;
1686 UInt_t index=0;
1687 //
1688 Double_t roady = 1.;
1689 Double_t roadz = 1.;
1690 //
1691 if (krow) {
1692 //cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1693 cl = krow.FindNearest2(y,z,roady,roadz,index);
1694 }
1695 t.fCurrentCluster = cl;
1696 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1697 return 1;
1698}
1c53abe2 1699
1c53abe2 1700
91162307 1701Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1702 //-----------------------------------------------------------------
1703 // This function tries to find a track prolongation to next pad row
1704 //-----------------------------------------------------------------
1c53abe2 1705
91162307 1706 //update error according neighborhoud
1c53abe2 1707
91162307 1708 if (t.fCurrentCluster) {
1709 t.fRow = nr;
8c51a605 1710 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
91162307 1711
1712 if (t.fCurrentCluster->IsUsed(10)){
1713 //
1714 //
1715 // t.fErrorZ2*=2;
1716 // t.fErrorY2*=2;
1717 t.fNShared++;
1718 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1719 t.fRemoval =10;
1720 return 0;
1721 }
b364ca79 1722 }
1723
1724 if (accept<3) UpdateTrack(&t,accept);
1725
1c53abe2 1726 } else {
91162307 1727 if (fIteration==0){
1728 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1729 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1730
1731 if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1c53abe2 1732 }
1733 }
1734 return 1;
1735}
1736
1737
1738
91162307 1739//_____________________________________________________________________________
1740Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1c53abe2 1741 //-----------------------------------------------------------------
91162307 1742 // This function tries to find a track prolongation.
1c53abe2 1743 //-----------------------------------------------------------------
1744 Double_t xt=t.GetX();
1745 //
91162307 1746 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1c53abe2 1747 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1748 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1749 //
1750 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1c53abe2 1751
91162307 1752 Int_t first = GetRowNumber(xt)-1;
1753 for (Int_t nr= first; nr>=rf; nr-=step) {
1754 if (FollowToNext(t,nr)==0)
1755 if (!t.IsActive()) return 0;
1756
1757 }
1758 return 1;
1759}
1760
1c53abe2 1761
1762//_____________________________________________________________________________
91162307 1763Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1c53abe2 1764 //-----------------------------------------------------------------
1765 // This function tries to find a track prolongation.
1766 //-----------------------------------------------------------------
1767 Double_t xt=t.GetX();
1768 //
1769 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1770 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1771 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1772 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1773
1774 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1775
1776 if (FollowToNextFast(t,nr)==0)
1777 if (!t.IsActive()) return 0;
1c53abe2 1778
1c53abe2 1779 }
1780 return 1;
1781}
1782
1783
91162307 1784
1785
1786
1c53abe2 1787Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1788 //-----------------------------------------------------------------
1789 // This function tries to find a track prolongation.
1790 //-----------------------------------------------------------------
91162307 1791 // Double_t xt=t.GetX();
1c53abe2 1792 //
1793 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1794 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1795 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 1796 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1c53abe2 1797
91162307 1798 Int_t first = 0;
1799 first = t.fFirstPoint+3;
1800 //
1801 if (first<0) first=0;
1802 for (Int_t nr=first+1; nr<=rf; nr++) {
1803 //if ( (t.GetSnp()<0.9))
1804 FollowToNext(t,nr);
1c53abe2 1805 }
1806 return 1;
1807}
1808
1809
1810
1811
1812
1813Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1814{
1815 //
1816 //
1817 sum1=0;
1818 sum2=0;
1819 Int_t sum=0;
1c53abe2 1820 //
1821 Float_t dz2 =(s1->GetZ() - s2->GetZ());
c9427e08 1822 dz2*=dz2;
91162307 1823
c9427e08 1824 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1c53abe2 1825 dy2*=dy2;
1826 Float_t distance = TMath::Sqrt(dz2+dy2);
c9427e08 1827 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1c53abe2 1828
91162307 1829 // Int_t offset =0;
c9427e08 1830 Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
1831 Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
c9427e08 1832 if (lastpoint>160)
1833 lastpoint =160;
1834 if (firstpoint<0)
1835 firstpoint = 0;
91162307 1836 if (firstpoint>lastpoint) {
1837 firstpoint =lastpoint;
1838 // lastpoint =160;
c9427e08 1839 }
1840
1841
91162307 1842 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1843 if (s1->GetClusterIndex2(i)>0) sum1++;
1844 if (s2->GetClusterIndex2(i)>0) sum2++;
1845 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1c53abe2 1846 sum++;
1847 }
1848 }
91162307 1849 if (sum<5) return 0;
1850
1627d1c4 1851 Float_t summin = TMath::Min(sum1+1,sum2+1);
1852 Float_t ratio = (sum+1)/Float_t(summin);
1c53abe2 1853 return ratio;
1854}
1855
1856void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1857{
1858 //
1859 //
91162307 1860 if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
1861 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
1862
1c53abe2 1863 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1864 dz2*=dz2;
1865 Float_t dy2 =(s1->GetY() - s2->GetY());
1866 dy2*=dy2;
91162307 1867 Float_t distance = dz2+dy2;
1868 if (distance>325.) return ; // if there are far away - not overlap - to reduce combinatorics
1869
1c53abe2 1870 //
91162307 1871 Int_t sumshared=0;
1872 //
1873 Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
1874 Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
1875 //
1876 if (firstpoint>=lastpoint-5) return;;
1c53abe2 1877
91162307 1878 for (Int_t i=firstpoint;i<lastpoint;i++){
1879 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1880 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1881 sumshared++;
1882 }
1883 }
1884 if (sumshared>4){
1885 // sign clusters
1886 //
1887 for (Int_t i=firstpoint;i<lastpoint;i++){
1888 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1889 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1890 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
1891 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
1892 if (s1->IsActive()&&s2->IsActive()){
1893 p1->fIsShared = kTRUE;
1894 p2->fIsShared = kTRUE;
1895 }
1896 }
1897 }
1898 }
1899 //
1900 if (sumshared>10){
1901 for (Int_t i=0;i<4;i++){
1902 if (s1->fOverlapLabels[3*i]==0){
1903 s1->fOverlapLabels[3*i] = s2->GetLabel();
1904 s1->fOverlapLabels[3*i+1] = sumshared;
1905 s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
1906 break;
1907 }
1908 }
1909 for (Int_t i=0;i<4;i++){
1910 if (s2->fOverlapLabels[3*i]==0){
1911 s2->fOverlapLabels[3*i] = s1->GetLabel();
1912 s2->fOverlapLabels[3*i+1] = sumshared;
1913 s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
1914 break;
1915 }
1916 }
1917 }
1c53abe2 1918
91162307 1919}
1c53abe2 1920
91162307 1921void AliTPCtrackerMI::SignShared(TObjArray * arr)
1922{
1c53abe2 1923 //
91162307 1924 //sort trackss according sectors
1925 //
c9427e08 1926 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1927 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 1928 if (!pt) continue;
1929 //if (pt) RotateToLocal(pt);
1930 pt->fSort = 0;
c9427e08 1931 }
91162307 1932 arr->UnSort();
1c53abe2 1933 arr->Sort(); // sorting according z
1934 arr->Expand(arr->GetEntries());
91162307 1935 //
1936 //
1c53abe2 1937 Int_t nseed=arr->GetEntriesFast();
1c53abe2 1938 for (Int_t i=0; i<nseed; i++) {
1939 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 1940 if (!pt) continue;
1941 for (Int_t j=0;j<=12;j++){
1942 pt->fOverlapLabels[j] =0;
1c53abe2 1943 }
91162307 1944 }
1945 for (Int_t i=0; i<nseed; i++) {
1946 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1947 if (!pt) continue;
1948 if (pt->fRemoval>10) continue;
1c53abe2 1949 for (Int_t j=i+1; j<nseed; j++){
1950 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
91162307 1951 // if (pt2){
1952 if (pt2->fRemoval<=10) {
1953 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
1954 SignShared(pt,pt2);
c9427e08 1955 }
91162307 1956 }
1957 }
1958}
1959
1960void AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
1961{
1962 //
1963 //sort trackss according sectors
1964 //
1965 if (fDebug&1) {
1966 printf("Number of tracks before double removal- %d\n",arr->GetEntries());
1967 }
1968 //
1969 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1970 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1971 if (!pt) continue;
1972 pt->fSort = 0;
1973 }
1974 arr->UnSort();
1975 arr->Sort(); // sorting according z
1976 arr->Expand(arr->GetEntries());
1977 //
1978 //reset overlap labels
1979 //
1980 Int_t nseed=arr->GetEntriesFast();
1981 for (Int_t i=0; i<nseed; i++) {
1982 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1983 if (!pt) continue;
1984 pt->SetUniqueID(i);
1985 for (Int_t j=0;j<=12;j++){
1986 pt->fOverlapLabels[j] =0;
1c53abe2 1987 }
1988 }
91162307 1989 //
1990 //sign shared tracks
1c53abe2 1991 for (Int_t i=0; i<nseed; i++) {
1992 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 1993 if (!pt) continue;
1994 if (pt->fRemoval>10) continue;
1995 Float_t deltac = pt->GetC()*0.1;
1996 for (Int_t j=i+1; j<nseed; j++){
1997 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1998 // if (pt2){
1999 if (pt2->fRemoval<=10) {
2000 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2001 if (TMath::Abs(pt->GetC() -pt2->GetC())>deltac) continue;
2002 if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2003 //
2004 SignShared(pt,pt2);
2005 }
1c53abe2 2006 }
1c53abe2 2007 }
91162307 2008 //
2009 // remove highly shared tracks
2010 for (Int_t i=0; i<nseed; i++) {
2011 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2012 if (!pt) continue;
2013 if (pt->fRemoval>10) continue;
2014 //
2015 Int_t sumshared =0;
2016 for (Int_t j=0;j<4;j++){
2017 sumshared = pt->fOverlapLabels[j*3+1];
2018 }
2019 Float_t factor = factor1;
2020 if (pt->fRemoval>0) factor = factor2;
2021 if (sumshared/pt->GetNumberOfClusters()>factor){
2022 for (Int_t j=0;j<4;j++){
2023 if (pt->fOverlapLabels[3*j]==0) continue;
2024 if (pt->fOverlapLabels[3*j+1]<5) continue;
2025 if (pt->fRemoval==removalindex) continue;
2026 AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2027 if (!pt2) continue;
2028 if (pt2->GetSigma2C()<pt->GetSigma2C()){
2029 // pt->fRemoval = removalindex;
2030 delete arr->RemoveAt(i);
2031 break;
1c53abe2 2032 }
91162307 2033 }
1c53abe2 2034 }
91162307 2035 }
2036 arr->Compress();
2037 if (fDebug&1) {
2038 printf("Number of tracks after double removal- %d\n",arr->GetEntries());
2039 }
2040}
2041
2042
2043
1c53abe2 2044
2045
91162307 2046
2047void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode)
2048{
2049 //
2050 //sort tracks in array according mode criteria
2051 Int_t nseed = arr->GetEntriesFast();
2052 for (Int_t i=0; i<nseed; i++) {
2053 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2054 if (!pt) {
2055 continue;
2056 }
2057 pt->fSort = mode;
2058 }
2059 arr->UnSort();
2060 arr->Sort();
1c53abe2 2061}
c9427e08 2062
91162307 2063void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
c9427e08 2064{
2065
2066 //Loop over all tracks and remove "overlaps"
2067 //
2068 //
2069 Int_t nseed = arr->GetEntriesFast();
2070 Int_t good =0;
91162307 2071
c9427e08 2072 for (Int_t i=0; i<nseed; i++) {
2073 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2074 if (!pt) {
91162307 2075 delete arr->RemoveAt(i);
c9427e08 2076 }
91162307 2077 else{
2078 pt->fSort =1;
2079 pt->fBSigned = kFALSE;
c9427e08 2080 }
91162307 2081 }
2082 arr->Compress();
2083 nseed = arr->GetEntriesFast();
2084 arr->UnSort();
2085 arr->Sort();
2086 //
2087 //unsign used
2088 UnsignClusters();
2089 //
2090 for (Int_t i=0; i<nseed; i++) {
2091 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2092 if (!pt) {
2093 continue;
2094 }
2095 Int_t found,foundable,shared;
2096 if (pt->IsActive())
2097 pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2098 else
2099 pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2100 //
2101 Double_t factor = factor2;
2102 if (pt->fBConstrain) factor = factor1;
2103
2104 if ((Float_t(shared)/Float_t(found))>factor){
c9427e08 2105 pt->Desactivate(removalindex);
91162307 2106 continue;
2107 }
2108
2109 good++;
2110 for (Int_t i=0; i<160; i++) {
2111 Int_t index=pt->GetClusterIndex2(i);
2112 if (index<0 || index&0x8000 ) continue;
2113 AliTPCclusterMI *c= pt->fClusterPointer[i];
2114 if (!c) continue;
2115 // if (!c->IsUsed(10)) c->Use(10);
2116 //if (pt->IsActive())
2117 c->Use(10);
2118 //else
2119 // c->Use(5);
c9427e08 2120 }
91162307 2121
c9427e08 2122 }
2123 fNtracks = good;
c9427e08 2124
91162307 2125 printf("\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
c9427e08 2126}
2127
91162307 2128void AliTPCtrackerMI::UnsignClusters()
1c53abe2 2129{
91162307 2130 //
2131 // loop over all clusters and unsign them
2132 //
2133
2134 for (Int_t sec=0;sec<fkNIS;sec++){
2135 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2136 AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2137 for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2138 // if (cl[icl].IsUsed(10))
2139 cl[icl].Use(-1);
2140 cl = fInnerSec[sec][row].fClusters2;
2141 for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2142 //if (cl[icl].IsUsed(10))
2143 cl[icl].Use(-1);
2144 }
2145 }
2146
2147 for (Int_t sec=0;sec<fkNOS;sec++){
2148 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2149 AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2150 for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2151 //if (cl[icl].IsUsed(10))
2152 cl[icl].Use(-1);
2153 cl = fOuterSec[sec][row].fClusters2;
2154 for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2155 //if (cl[icl].IsUsed(10))
2156 cl[icl].Use(-1);
2157 }
2158 }
2159
1c53abe2 2160}
2161
91162307 2162
2163
732b9e78 2164void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
1c53abe2 2165{
2166 //
91162307 2167 //sign clusters to be "used"
2168 //
2169 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2170 // loop over "primaries"
2171
2172 Float_t sumdens=0;
2173 Float_t sumdens2=0;
2174 Float_t sumn =0;
2175 Float_t sumn2 =0;
2176 Float_t sumchi =0;
2177 Float_t sumchi2 =0;
2178
2179 Float_t sum =0;
2180
1c53abe2 2181 TStopwatch timer;
91162307 2182 timer.Start();
1c53abe2 2183
91162307 2184 Int_t nseed = arr->GetEntriesFast();
2185 for (Int_t i=0; i<nseed; i++) {
2186 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2187 if (!pt) {
2188 continue;
2189 }
2190 if (!(pt->IsActive())) continue;
2191 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2192 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2193 sumdens += dens;
2194 sumdens2+= dens*dens;
2195 sumn += pt->GetNumberOfClusters();
2196 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2197 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2198 if (chi2>5) chi2=5;
2199 sumchi +=chi2;
2200 sumchi2 +=chi2*chi2;
2201 sum++;
2202 }
1627d1c4 2203 }
91162307 2204
2205 Float_t mdensity = 0.9;
2206 Float_t meann = 130;
2207 Float_t meanchi = 1;
2208 Float_t sdensity = 0.1;
2209 Float_t smeann = 10;
2210 Float_t smeanchi =0.4;
1627d1c4 2211
91162307 2212
2213 if (sum>20){
2214 mdensity = sumdens/sum;
2215 meann = sumn/sum;
2216 meanchi = sumchi/sum;
2217 //
2218 sdensity = sumdens2/sum-mdensity*mdensity;
2219 sdensity = TMath::Sqrt(sdensity);
2220 //
2221 smeann = sumn2/sum-meann*meann;
2222 smeann = TMath::Sqrt(smeann);
2223 //
2224 smeanchi = sumchi2/sum - meanchi*meanchi;
2225 smeanchi = TMath::Sqrt(smeanchi);
2226 }
2227
2228
2229 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2230 //
2231 for (Int_t i=0; i<nseed; i++) {
2232 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2233 if (!pt) {
2234 continue;
1c53abe2 2235 }
91162307 2236 if (pt->fBSigned) continue;
2237 if (pt->fBConstrain) continue;
2238 //if (!(pt->IsActive())) continue;
2239 /*
2240 Int_t found,foundable,shared;
2241 pt->GetClusterStatistic(0,160,found, foundable,shared);
2242 if (shared/float(found)>0.3) {
2243 if (shared/float(found)>0.9 ){
2244 //delete arr->RemoveAt(i);
2245 }
2246 continue;
c9427e08 2247 }
91162307 2248 */
2249 Bool_t isok =kFALSE;
2250 if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2251 isok = kTRUE;
2252 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2253 isok =kTRUE;
2254 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2255 isok =kTRUE;
2256 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2257 isok =kTRUE;
2258
2259 if (isok)
2260 for (Int_t i=0; i<160; i++) {
2261 Int_t index=pt->GetClusterIndex2(i);
2262 if (index<0) continue;
2263 AliTPCclusterMI *c= pt->fClusterPointer[i];
2264 if (!c) continue;
2265 //if (!(c->IsUsed(10))) c->Use();
2266 c->Use(10);
2267 }
2268 }
2269
c9427e08 2270
1c53abe2 2271 //
91162307 2272 Double_t maxchi = meanchi+2.*smeanchi;
2273
2274 for (Int_t i=0; i<nseed; i++) {
2275 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2276 if (!pt) {
1c53abe2 2277 continue;
91162307 2278 }
2279 //if (!(pt->IsActive())) continue;
2280 if (pt->fBSigned) continue;
2281 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2282 if (chi>maxchi) continue;
2283
2284 Float_t bfactor=1;
2285 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2286
2287 //sign only tracks with enoug big density at the beginning
2288
2289 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2290
2291
3e5d0aa2 2292 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
91162307 2293 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2294
2295 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2296 if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2297 minn=0;
2298
2299 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2300 //Int_t noc=pt->GetNumberOfClusters();
2301 pt->fBSigned = kTRUE;
2302 for (Int_t i=0; i<160; i++) {
2303
2304 Int_t index=pt->GetClusterIndex2(i);
2305 if (index<0) continue;
2306 AliTPCclusterMI *c= pt->fClusterPointer[i];
2307 if (!c) continue;
2308 // if (!(c->IsUsed(10))) c->Use();
2309 c->Use(10);
2310 }
1c53abe2 2311 }
91162307 2312 }
2313 // gLastCheck = nseed;
2314 // arr->Compress();
2315 if (fDebug>0){
2316 timer.Print();
2317 }
1c53abe2 2318}
2319
2320
91162307 2321void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2)
2322{
2323 // stop not active tracks
2324 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2325 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2326 Int_t nseed = arr->GetEntriesFast();
2327 //
2328 for (Int_t i=0; i<nseed; i++) {
2329 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2330 if (!pt) {
2331 continue;
2332 }
2333 if (!(pt->IsActive())) continue;
2334 StopNotActive(pt,row0,th0, th1,th2);
2335 }
2336}
1c53abe2 2337
1c53abe2 2338
1c53abe2 2339
91162307 2340void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2341 Float_t th2)
2342{
2343 // stop not active tracks
2344 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2345 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2346 Int_t sumgood1 = 0;
2347 Int_t sumgood2 = 0;
2348 Int_t foundable = 0;
2349 Int_t maxindex = seed->fLastPoint; //last foundable row
2350 if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2351 seed->Desactivate(10) ;
2352 return;
2353 }
1c53abe2 2354
3e5d0aa2 2355 for (Int_t i=row0; i<maxindex; i++){
91162307 2356 Int_t index = seed->GetClusterIndex2(i);
2357 if (index!=-1) foundable++;
2358 //if (!c) continue;
2359 if (foundable<=30) sumgood1++;
2360 if (foundable<=50) {
2361 sumgood2++;
2362 }
2363 else{
2364 break;
2365 }
2366 }
2367 if (foundable>=30.){
2368 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2369 }
2370 if (foundable>=50)
2371 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2372}
1c53abe2 2373
1c53abe2 2374
1c53abe2 2375
91162307 2376Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2377{
2378 //
2379 // back propagation of ESD tracks
2380 //
1c53abe2 2381
91162307 2382 fEvent = event;
2383 ReadSeeds(event);
2384 PropagateBack(fSeeds);
2385 Int_t nseed = fSeeds->GetEntriesFast();
2386 for (Int_t i=0;i<nseed;i++){
2387 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2388 AliESDtrack *esd=event->GetTrack(i);
2389 seed->CookdEdx(0.02,0.06);
2390 CookLabel(seed,0.1); //For comparison only
2391 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2392 }
2393 fEvent =0;
2394 WriteTracks();
2395 return 0;
2396}
2397
2398
2399void AliTPCtrackerMI::DeleteSeeds()
2400{
2401 Int_t nseed = fSeeds->GetEntriesFast();
2402 for (Int_t i=0;i<nseed;i++){
2403 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2404 if (seed) delete fSeeds->RemoveAt(i);
2405 }
2406 delete fSeeds;
2407 fSeeds =0;
2408}
2409
2410void AliTPCtrackerMI::ReadSeeds(AliESD *event)
2411{
2412 //
2413 //read seeds from the event
2414
2415 Int_t nentr=event->GetNumberOfTracks();
2416 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
2417 if (fSeeds)
2418 DeleteSeeds();
2419 if (!fSeeds){
2420 fSeeds = new TObjArray;
2421 }
2422
2423 // Int_t ntrk=0;
2424 for (Int_t i=0; i<nentr; i++) {
2425 AliESDtrack *esd=event->GetTrack(i);
2426 ULong_t status=esd->GetStatus();
2427 const AliTPCtrack t(*esd);
2428 AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2429 if (status==AliESDtrack::kTPCin) seed->Modify(0.8);
2430 //
2431 //
2432 // rotate to the local coordinate system
2433
2434 fSectors=fInnerSec; fN=fkNIS;
2435
2436 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2437 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2438 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2439 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2440 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2441 alpha-=seed->GetAlpha();
2442 if (!seed->Rotate(alpha)) continue;
2443 //
2444 seed->PropagateTo(fSectors->GetX(0));
2445 //
2446 // Int_t index = esd->GetTPCindex();
2447 //AliTPCseed * seed2= (AliTPCseed*)fSeeds->At(index);
2448
2449 fSeeds->AddLast(seed);
2450 }
2451}
2452
2453
2454
2455//_____________________________________________________________________________
2456void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2457 Float_t deltay, Int_t ddsec) {
2458 //-----------------------------------------------------------------
2459 // This function creates track seeds.
2460 // SEEDING WITH VERTEX CONSTRAIN
2461 //-----------------------------------------------------------------
2462 // cuts[0] - fP4 cut
2463 // cuts[1] - tan(phi) cut
2464 // cuts[2] - zvertex cut
2465 // cuts[3] - fP3 cut
2466 Int_t nin0 = 0;
2467 Int_t nin1 = 0;
2468 Int_t nin2 = 0;
2469 Int_t nin = 0;
2470 Int_t nout1 = 0;
2471 Int_t nout2 = 0;
2472
2473 Double_t x[5], c[15];
2474 // Int_t di = i1-i2;
2475 //
2476 AliTPCseed * seed = new AliTPCseed;
2477 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2478 Double_t cs=cos(alpha), sn=sin(alpha);
2479 //
2480 // Double_t x1 =fOuterSec->GetX(i1);
2481 //Double_t xx2=fOuterSec->GetX(i2);
2482
2483 Double_t x1 =GetXrow(i1);
2484 Double_t xx2=GetXrow(i2);
2485
2486 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2487
2488 Int_t imiddle = (i2+i1)/2; //middle pad row index
2489 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2490 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2491 //
2492 Int_t ns =sec;
2493
2494 const AliTPCRow& kr1=GetRow(ns,i1);
2495 Double_t ymax = GetMaxY(i1)-kr1.fDeadZone-1.5;
2496 Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2497
2498 //
2499 // change cut on curvature if it can't reach this layer
2500 // maximal curvature set to reach it
2501 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2502 if (dvertexmax*0.5*cuts[0]>0.85){
2503 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2504 }
2505 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2506
2507 // Int_t ddsec = 1;
2508 if (deltay>0) ddsec = 0;
2509 // loop over clusters
2510 for (Int_t is=0; is < kr1; is++) {
2511 //
2512 if (kr1[is]->IsUsed(10)) continue;
2513 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2514 //if (TMath::Abs(y1)>ymax) continue;
2515
2516 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2517
2518 // find possible directions
2519 Float_t anglez = (z1-z3)/(x1-x3);
2520 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2521 //
2522 //
2523 //find rotation angles relative to line given by vertex and point 1
2524 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2525 Double_t dvertex = TMath::Sqrt(dvertex2);
2526 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2527 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2528
2529 //
2530 // loop over 2 sectors
2531 Int_t dsec1=-ddsec;
2532 Int_t dsec2= ddsec;
2533 if (y1<0) dsec2= 0;
2534 if (y1>0) dsec1= 0;
2535
2536 Double_t dddz1=0; // direction of delta inclination in z axis
2537 Double_t dddz2=0;
2538 if ( (z1-z3)>0)
2539 dddz1 =1;
2540 else
2541 dddz2 =1;
2542 //
2543 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2544 Int_t sec2 = sec + dsec;
2545 //
2546 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2547 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2548 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2549 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2550 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2551 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2552
2553 // rotation angles to p1-p3
2554 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2555 Double_t x2, y2, z2;
2556 //
2557 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2558
2559 //
2560 Double_t dxx0 = (xx2-x3)*cs13r;
2561 Double_t dyy0 = (xx2-x3)*sn13r;
2562 for (Int_t js=index1; js < index2; js++) {
2563 const AliTPCclusterMI *kcl = kr2[js];
2564 if (kcl->IsUsed(10)) continue;
2565 //
2566 //calcutate parameters
2567 //
2568 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2569 // stright track
2570 if (TMath::Abs(yy0)<0.000001) continue;
2571 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2572 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2573 Double_t r02 = (0.25+y0*y0)*dvertex2;
2574 //curvature (radius) cut
2575 if (r02<r2min) continue;
2576
2577 nin0++;
2578 //
2579 Double_t c0 = 1/TMath::Sqrt(r02);
2580 if (yy0>0) c0*=-1.;
2581
2582
2583 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2584 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2585 Double_t dfi0 = 2.*TPCFastMath::FastAsin(dvertex*c0*0.5);
2586 Double_t dfi1 = 2.*TPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2587 //
2588 //
2589 Double_t z0 = kcl->GetZ();
2590 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2591 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2592 nin1++;
2593 //
2594 Double_t dip = (z1-z0)*c0/dfi1;
2595 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2596 //
2597 y2 = kcl->GetY();
2598 if (dsec==0){
2599 x2 = xx2;
2600 z2 = kcl->GetZ();
2601 }
2602 else
2603 {
2604 // rotation
2605 z2 = kcl->GetZ();
2606 x2= xx2*cs-y2*sn*dsec;
2607 y2=+xx2*sn*dsec+y2*cs;
2608 }
2609
2610 x[0] = y1;
2611 x[1] = z1;
2612 x[2] = x0;
2613 x[3] = dip;
2614 x[4] = c0;
2615 //
2616 //
2617 // do we have cluster at the middle ?
2618 Double_t ym,zm;
2619 GetProlongation(x1,xm,x,ym,zm);
2620 UInt_t dummy;
2621 AliTPCclusterMI * cm=0;
2622 if (TMath::Abs(ym)-ymaxm<0){
2623 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2624 if ((!cm) || (cm->IsUsed(10))) {
2625 continue;
2626 }
2627 }
2628 else{
2629 // rotate y1 to system 0
2630 // get state vector in rotated system
2631 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2632 Double_t xr2 = x0*cs+yr1*sn*dsec;
2633 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2634 //
2635 GetProlongation(xx2,xm,xr,ym,zm);
2636 if (TMath::Abs(ym)-ymaxm<0){
2637 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2638 if ((!cm) || (cm->IsUsed(10))) {
2639 continue;
2640 }
2641 }
2642 }
2643
2644
2645 Double_t dym = 0;
2646 Double_t dzm = 0;
2647 if (cm){
2648 dym = ym - cm->GetY();
2649 dzm = zm - cm->GetZ();
2650 }
2651 nin2++;
2652
2653
2654 //
2655 //
2656 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2657 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2658 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2659 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2660 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2661
2662 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2663 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2664 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2665 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2666 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1c53abe2 2667 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
91162307 2668
1c53abe2 2669 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2670 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2671 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2672 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
91162307 2673
1c53abe2 2674 c[0]=sy1;
2675 c[1]=0.; c[2]=sz1;
2676 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2677 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2678 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2679 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2680 c[13]=f30*sy1*f40+f32*sy2*f42;
2681 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
91162307 2682
2683 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2684
1c53abe2 2685 UInt_t index=kr1.GetIndex(is);
91162307 2686 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
2687
1c53abe2 2688 track->fIsSeeding = kTRUE;
91162307 2689 track->fSeed1 = i1;
2690 track->fSeed2 = i2;
2691 track->fSeedType=3;
c9427e08 2692
91162307 2693
2694 //if (dsec==0) {
2695 FollowProlongation(*track, (i1+i2)/2,1);
2696 Int_t foundable,found,shared;
2697 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
2698 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
2699 seed->Reset();
2700 seed->~AliTPCseed();
2701 continue;
2702 }
2703 //}
2704
2705 nin++;
2706 FollowProlongation(*track, i2,1);
2707
2708
2709 //Int_t rc = 1;
2710 track->fBConstrain =1;
2711 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2712 track->fLastPoint = i1; // first cluster in track position
2713 track->fFirstPoint = track->fLastPoint;
2714
2715 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2716 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
2717 track->fNShared>0.4*track->GetNumberOfClusters() ) {
2718 seed->Reset();
2719 seed->~AliTPCseed();
c9427e08 2720 continue;
2721 }
91162307 2722 nout1++;
2723 // Z VERTEX CONDITION
2724 Double_t zv;
2725 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2726 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2727 if (TMath::Abs(zv-z3)>cuts[2]) {
2728 FollowProlongation(*track, TMath::Max(i2-20,0));
2729 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2730 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2731 if (TMath::Abs(zv-z3)>cuts[2]){
2732 FollowProlongation(*track, TMath::Max(i2-40,0));
2733 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2734 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2735 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
2736 // make seed without constrain
2737 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
2738 FollowProlongation(*track2, i2,1);
2739 track2->fBConstrain = kFALSE;
2740 track2->fSeedType = 1;
2741 arr->AddLast(track2);
2742 seed->Reset();
2743 seed->~AliTPCseed();
2744 continue;
2745 }
2746 else{
2747 seed->Reset();
2748 seed->~AliTPCseed();
2749 continue;
2750
2751 }
2752 }
c9427e08 2753 }
1c53abe2 2754
91162307 2755 track->fSeedType =0;
2756 arr->AddLast(track);
2757 seed = new AliTPCseed;
2758 nout2++;
2759 // don't consider other combinations
2760 if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
2761 break;
1c53abe2 2762 }
2763 }
2764 }
91162307 2765 if (fDebug>1){
2766 // printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
2767 }
2768 delete seed;
1c53abe2 2769}
2770
1627d1c4 2771
91162307 2772void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2773 Float_t deltay) {
2774
2775
2776
1627d1c4 2777 //-----------------------------------------------------------------
91162307 2778 // This function creates track seeds.
1627d1c4 2779 //-----------------------------------------------------------------
91162307 2780 // cuts[0] - fP4 cut
2781 // cuts[1] - tan(phi) cut
2782 // cuts[2] - zvertex cut
2783 // cuts[3] - fP3 cut
2784
2785
2786 Int_t nin0 = 0;
2787 Int_t nin1 = 0;
2788 Int_t nin2 = 0;
2789 Int_t nin = 0;
2790 Int_t nout1 = 0;
2791 Int_t nout2 = 0;
2792 Int_t nout3 =0;
2793 Double_t x[5], c[15];
2794 //
2795 // make temporary seed
2796 AliTPCseed * seed = new AliTPCseed;
1627d1c4 2797 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
2798 // Double_t cs=cos(alpha), sn=sin(alpha);
91162307 2799 //
2800 //
1627d1c4 2801
91162307 2802 // first 3 padrows
2803 Double_t x1 = GetXrow(i1-1);
2804 const AliTPCRow& kr1=GetRow(sec,i1-1);
2805 Double_t y1max = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
2806 //
2807 Double_t x1p = GetXrow(i1);
2808 const AliTPCRow& kr1p=GetRow(sec,i1);
2809 //
2810 Double_t x1m = GetXrow(i1-2);
2811 const AliTPCRow& kr1m=GetRow(sec,i1-2);
1627d1c4 2812
91162307 2813 //
2814 //last 3 padrow for seeding
2815 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
2816 Double_t x3 = GetXrow(i1-7);
2817 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
2818 //
2819 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
2820 Double_t x3p = GetXrow(i1-6);
2821 //
2822 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
2823 Double_t x3m = GetXrow(i1-8);
1627d1c4 2824
91162307 2825 //
2826 //
2827 // middle padrow
2828 Int_t im = i1-4; //middle pad row index
2829 Double_t xm = GetXrow(im); // radius of middle pad-row
2830 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
2831 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
2832 //
2833 //
2834 Double_t deltax = x1-x3;
2835 Double_t dymax = deltax*cuts[1];
2836 Double_t dzmax = deltax*cuts[3];
2837 //
2838 // loop over clusters
2839 for (Int_t is=0; is < kr1; is++) {
1627d1c4 2840 //
91162307 2841 if (kr1[is]->IsUsed(10)) continue;
2842 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1627d1c4 2843 //
91162307 2844 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2845 //
2846 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
2847 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
2848 //
2849 Double_t y3, z3;
1627d1c4 2850 //
1627d1c4 2851 //
91162307 2852 UInt_t index;
2853 for (Int_t js=index1; js < index2; js++) {
2854 const AliTPCclusterMI *kcl = kr3[js];
2855 if (kcl->IsUsed(10)) continue;
2856 y3 = kcl->GetY();
2857 // apply angular cuts
2858 if (TMath::Abs(y1-y3)>dymax) continue;
2859 x3 = x3;
2860 z3 = kcl->GetZ();
2861 if (TMath::Abs(z1-z3)>dzmax) continue;
2862 //
2863 Double_t angley = (y1-y3)/(x1-x3);
2864 Double_t anglez = (z1-z3)/(x1-x3);
2865 //
2866 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
2867 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
2868 //
2869 Double_t yyym = angley*(xm-x1)+y1;
2870 Double_t zzzm = anglez*(xm-x1)+z1;
2871
2872 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
2873 if (!kcm) continue;
2874 if (kcm->IsUsed(10)) continue;
2875
2876 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
2877 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
2878 //
2879 //
2880 //
2881 Int_t used =0;
2882 Int_t found =0;
2883 //
2884 // look around first
2885 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
2886 anglez*(x1m-x1)+z1,
2887 erry,errz,index);
2888 //
2889 if (kc1m){
2890 found++;
2891 if (kc1m->IsUsed(10)) used++;
1627d1c4 2892 }
91162307 2893 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
2894 anglez*(x1p-x1)+z1,
2895 erry,errz,index);
1627d1c4 2896 //
91162307 2897 if (kc1p){
2898 found++;
2899 if (kc1p->IsUsed(10)) used++;
1627d1c4 2900 }
91162307 2901 if (used>1) continue;
2902 if (found<1) continue;
1627d1c4 2903
91162307 2904 //
2905 // look around last
2906 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
2907 anglez*(x3m-x3)+z3,
2908 erry,errz,index);
2909 //
2910 if (kc3m){
2911 found++;
2912 if (kc3m->IsUsed(10)) used++;
2913 }
2914 else
2915 continue;
2916 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
2917 anglez*(x3p-x3)+z3,
2918 erry,errz,index);
2919 //
2920 if (kc3p){
2921 found++;
2922 if (kc3p->IsUsed(10)) used++;
2923 }
2924 else
2925 continue;
2926 if (used>1) continue;
2927 if (found<3) continue;
2928 //
2929 Double_t x2,y2,z2;
2930 x2 = xm;
2931 y2 = kcm->GetY();
2932 z2 = kcm->GetZ();
2933 //
2934
1627d1c4 2935 x[0]=y1;
2936 x[1]=z1;
2937 x[4]=f1(x1,y1,x2,y2,x3,y3);
91162307 2938 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
2939 nin0++;
2940 //
1627d1c4 2941 x[2]=f2(x1,y1,x2,y2,x3,y3);
91162307 2942 nin1++;
2943 //
2944 x[3]=f3n(x1,y1,x2,y2,z1,z2,x[4]);
2945 //if (TMath::Abs(x[3]) > cuts[3]) continue;
2946 nin2++;
2947 //
2948 //
2949 Double_t sy1=0.1, sz1=0.1;
2950 Double_t sy2=0.1, sz2=0.1;
2951 Double_t sy3=0.1, sy=0.1, sz=0.1;
1627d1c4 2952
2953 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2954 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2955 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2956 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2957 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2958 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
91162307 2959
1627d1c4 2960 Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2961 Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2962 Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2963 Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2964
2965 c[0]=sy1;
91162307 2966 c[1]=0.; c[2]=sz1;
1627d1c4 2967 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2968 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2969 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2970 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2971 c[13]=f30*sy1*f40+f32*sy2*f42;
2972 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2973
91162307 2974 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2975
2976 UInt_t index=kr1.GetIndex(is);
2977 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
2978
2979 track->fIsSeeding = kTRUE;
2980
2981 nin++;
2982 FollowProlongation(*track, i1-7,1);
2983 if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
2984 track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
2985 seed->Reset();
2986 seed->~AliTPCseed();
2987 continue;
2988 }
2989 nout1++;
2990 nout2++;
2991 //Int_t rc = 1;
2992 FollowProlongation(*track, i2,1);
2993 track->fBConstrain =0;
2994 track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2995 track->fFirstPoint = track->fLastPoint;
2996
2997 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2998 track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
2999 track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3000 seed->Reset();
3001 seed->~AliTPCseed();
3002 continue;
3003 }
3004
3005 {
3006 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3007 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3008 FollowProlongation(*track2, i2,1);
3009 track2->fBConstrain = kFALSE;
3010 track2->fSeedType = 4;
3011 arr->AddLast(track2);
3012 seed->Reset();
3013 seed->~AliTPCseed();
3014 }
3015
3016
3017 //arr->AddLast(track);
3018 //seed = new AliTPCseed;
3019 nout3++;
3020 }
3021 }
3022
3023 if (fDebug>1){
3024 // printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3025 }
3026 delete seed;
3027}
3028
3029
3030//_____________________________________________________________________________
176aff27 3031void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3032 Float_t deltay, Bool_t /*bconstrain*/) {
91162307 3033 //-----------------------------------------------------------------
3034 // This function creates track seeds - without vertex constraint
3035 //-----------------------------------------------------------------
3036 // cuts[0] - fP4 cut - not applied
3037 // cuts[1] - tan(phi) cut
3038 // cuts[2] - zvertex cut - not applied
3039 // cuts[3] - fP3 cut
3040 Int_t nin0=0;
3041 Int_t nin1=0;
3042 Int_t nin2=0;
3043 Int_t nin3=0;
3044 // Int_t nin4=0;
3045 //Int_t nin5=0;
3046
3047
3048
3049 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3050 // Double_t cs=cos(alpha), sn=sin(alpha);
3051 Int_t row0 = (i1+i2)/2;
3052 Int_t drow = (i1-i2)/2;
3053 const AliTPCRow& kr0=fSectors[sec][row0];
3054 AliTPCRow * kr=0;
3055
3056 AliTPCpolyTrack polytrack;
3057 Int_t nclusters=fSectors[sec][row0];
3058 AliTPCseed * seed = new AliTPCseed;
3059
3060 Int_t sumused=0;
3061 Int_t cused=0;
3062 Int_t cnused=0;
3063 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3064 Int_t nfound =0;
3065 Int_t nfoundable =0;
3066 for (Int_t iter =1; iter<2; iter++){ //iterations
3067 const AliTPCRow& krm=fSectors[sec][row0-iter];
3068 const AliTPCRow& krp=fSectors[sec][row0+iter];
3069 const AliTPCclusterMI * cl= kr0[is];
3070
3071 if (cl->IsUsed(10)) {
3072 cused++;
3073 }
3074 else{
3075 cnused++;
3076 }
3077 Double_t x = kr0.GetX();
3078 // Initialization of the polytrack
3079 nfound =0;
3080 nfoundable =0;
3081 polytrack.Reset();
3082 //
3083 Double_t y0= cl->GetY();
3084 Double_t z0= cl->GetZ();
3085 Float_t erry = 0;
3086 Float_t errz = 0;
3087
3088 Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3089 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3090
3091 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3092 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3093 polytrack.AddPoint(x,y0,z0,erry, errz);
3094
3095 sumused=0;
3096 if (cl->IsUsed(10)) sumused++;
3097
3098
3099 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3100 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3101 //
3102 x = krm.GetX();
3103 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3104 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3105 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3106 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3107 if (cl1->IsUsed(10)) sumused++;
3108 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3109 }
3110 //
3111 x = krp.GetX();
3112 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3113 if (cl2) {
3114 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3115 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3116 if (cl2->IsUsed(10)) sumused++;
3117 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3118 }
3119 //
3120 if (sumused>0) continue;
3121 nin0++;
3122 polytrack.UpdateParameters();
3123 // follow polytrack
3124 roadz = 1.2;
3125 roady = 1.2;
3126 //
3127 Double_t yn,zn;
3128 nfoundable = polytrack.GetN();
3129 nfound = nfoundable;
3130 //
3131 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3132 Float_t maxdist = 0.8*(1.+3./(ddrow));
3133 for (Int_t delta = -1;delta<=1;delta+=2){
3134 Int_t row = row0+ddrow*delta;
3135 kr = &(fSectors[sec][row]);
3136 Double_t xn = kr->GetX();
3137 Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3138 polytrack.GetFitPoint(xn,yn,zn);
3139 if (TMath::Abs(yn)>ymax) continue;
3140 nfoundable++;
3141 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3142 if (cln) {
3143 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3144 if (dist<maxdist){
3145 /*
3146 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3147 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3148 if (cln->IsUsed(10)) {
3149 // printf("used\n");
3150 sumused++;
3151 erry*=2;
3152 errz*=2;
3153 }
3154 */
3155 erry=0.1;
3156 errz=0.1;
3157 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3158 nfound++;
3159 }
3160 }
3161 }
3162 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3163 polytrack.UpdateParameters();
3164 }
3165 }
3166 if ( (sumused>3) || (sumused>0.5*nfound)) {
3167 //printf("sumused %d\n",sumused);
3168 continue;
3169 }
3170 nin1++;
3171 Double_t dy,dz;
3172 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3173 AliTPCpolyTrack track2;
3174
3175 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3176 if (track2.GetN()<0.5*nfoundable) continue;
3177 nin2++;
3178
3179 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3180 //
3181 // test seed with and without constrain
3182 for (Int_t constrain=0; constrain<=0;constrain++){
3183 // add polytrack candidate
3184
3185 Double_t x[5], c[15];
3186 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3187 track2.GetBoundaries(x3,x1);
3188 x2 = (x1+x3)/2.;
3189 track2.GetFitPoint(x1,y1,z1);
3190 track2.GetFitPoint(x2,y2,z2);
3191 track2.GetFitPoint(x3,y3,z3);
3192 //
3193 //is track pointing to the vertex ?
3194 Double_t x0,y0,z0;
3195 x0=0;
3196 polytrack.GetFitPoint(x0,y0,z0);
3197
3198 if (constrain) {
3199 x2 = x3;
3200 y2 = y3;
3201 z2 = z3;
3202
3203 x3 = 0;
3204 y3 = 0;
3205 z3 = 0;
3206 }
3207 x[0]=y1;
3208 x[1]=z1;
3209 x[4]=f1(x1,y1,x2,y2,x3,y3);
3210
3211 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3212 x[2]=f2(x1,y1,x2,y2,x3,y3);
3213
3214 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3215 //x[3]=f3(x1,y1,x2,y2,z1,z2);
3216 x[3]=f3n(x1,y1,x3,y3,z1,z3,x[4]);
3217 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3218
3219
3220 Double_t sy =0.1, sz =0.1;
3221 Double_t sy1=0.02, sz1=0.02;
3222 Double_t sy2=0.02, sz2=0.02;
3223 Double_t sy3=0.02;
3224
3225 if (constrain){
3226 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3227 }
3228
3229 Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3230 Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3231 Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3232 Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3233 Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3234 Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3235
3236 Double_t f30=(f3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3237 Double_t f31=(f3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3238 Double_t f32=(f3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3239 Double_t f34=(f3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3240
3241
3242 c[0]=sy1;
3243 c[1]=0.; c[2]=sz1;
3244 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3245 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3246 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3247 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3248 c[13]=f30*sy1*f40+f32*sy2*f42;
3249 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3250
3251 //Int_t row1 = fSectors->GetRowNumber(x1);
3252 Int_t row1 = GetRowNumber(x1);
3253
3254 UInt_t index=0;
3255 //kr0.GetIndex(is);
3256 AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3257 track->fIsSeeding = kTRUE;
3258 Int_t rc=FollowProlongation(*track, i2);
3259 if (constrain) track->fBConstrain =1;
3260 else
3261 track->fBConstrain =0;
3262 track->fLastPoint = row1+fInnerSec->GetNRows(); // first cluster in track position
3263 track->fFirstPoint = track->fLastPoint;
3264
3265 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3266 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3267 track->fNShared>0.4*track->GetNumberOfClusters()) {
3268 //delete track;
3269 seed->Reset();
3270 seed->~AliTPCseed();
3271 }
3272 else {
3273 arr->AddLast(track);
3274 seed = new AliTPCseed;
3275 }
3276 nin3++;
3277 }
3278 } // if accepted seed
3279 }
3280 if (fDebug>1){
3281 printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3282 }
3283 delete seed;
3284}
3285
3286
3287AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3288{
3289 //
3290 //
3291 //reseed
3292 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3293 Int_t p1 = int(r1*track->GetNumberOfClusters());
3294 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
176aff27 3295 Int_t pp2=0;
91162307 3296 Double_t x0[3],x1[3],x2[3];
3297 x0[0]=-1;
3298 x0[0]=-1;
3299 x0[0]=-1;
3300
3301 // find track position at given ratio of the length
3302 Int_t sec0, sec1, sec2;
176aff27 3303 sec0=0;
3304 sec1=0;
3305 sec2=0;
91162307 3306 Int_t index=-1;
3307 Int_t clindex;
3308 for (Int_t i=0;i<160;i++){
3309 if (track->fClusterPointer[i]){
3310 index++;
3311 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3312 if ( (index<p0) || x0[0]<0 ){
3313 if (trpoint->GetX()>1){
3314 clindex = track->GetClusterIndex2(i);
3315 if (clindex>0){
3316 x0[0] = trpoint->GetX();
3317 x0[1] = trpoint->GetY();
3318 x0[2] = trpoint->GetZ();
3319 sec0 = ((clindex&0xff000000)>>24)%18;
3320 }
3321 }
3322 }
3323
3324 if ( (index<p1) &&(trpoint->GetX()>1)){
3325 clindex = track->GetClusterIndex2(i);
3326 if (clindex>0){
3327 x1[0] = trpoint->GetX();
3328 x1[1] = trpoint->GetY();
3329 x1[2] = trpoint->GetZ();
3330 sec1 = ((clindex&0xff000000)>>24)%18;
3331 }
3332 }
3333 if ( (index<p2) &&(trpoint->GetX()>1)){
3334 clindex = track->GetClusterIndex2(i);
3335 if (clindex>0){
3336 x2[0] = trpoint->GetX();
3337 x2[1] = trpoint->GetY();
3338 x2[2] = trpoint->GetZ();
3339 sec2 = ((clindex&0xff000000)>>24)%18;
3340 pp2 = i;
3341 }
3342 }
3343 }
3344 }
3345
3346 Double_t alpha, cs,sn, xx2,yy2;
3347 //
3348 alpha = (sec1-sec2)*fSectors->GetAlpha();
3349 cs = TMath::Cos(alpha);
3350 sn = TMath::Sin(alpha);
3351 xx2= x1[0]*cs-x1[1]*sn;
3352 yy2= x1[0]*sn+x1[1]*cs;
3353 x1[0] = xx2;
3354 x1[1] = yy2;
3355 //
3356 alpha = (sec0-sec2)*fSectors->GetAlpha();
3357 cs = TMath::Cos(alpha);
3358 sn = TMath::Sin(alpha);
3359 xx2= x0[0]*cs-x0[1]*sn;
3360 yy2= x0[0]*sn+x0[1]*cs;
3361 x0[0] = xx2;
3362 x0[1] = yy2;
3363 //
3364 //
3365 //
3366 Double_t x[5],c[15];
3367 //
3368 x[0]=x2[1];
3369 x[1]=x2[2];
3370 x[4]=f1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3371 // if (x[4]>1) return 0;
3372 x[2]=f2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3373 x[3]=f3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3374 //if (TMath::Abs(x[3]) > 2.2) return 0;
3375 //if (TMath::Abs(x[2]) > 1.99) return 0;
3376 //
3377 Double_t sy =0.1, sz =0.1;
3378 //
3379 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3380 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3381 Double_t sy3=0.01+track->GetSigmaY2();
3382 //
3383 Double_t f40=(f1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3384 Double_t f42=(f1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3385 Double_t f43=(f1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3386 Double_t f20=(f2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3387 Double_t f22=(f2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3388 Double_t f23=(f2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3389 //
3390 Double_t f30=(f3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3391 Double_t f31=(f3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3392 Double_t f32=(f3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3393 Double_t f34=(f3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3394
3395
3396 c[0]=sy1;
3397 c[1]=0.; c[2]=sz1;
3398 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3399 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3400 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3401 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3402 c[13]=f30*sy1*f40+f32*sy2*f42;
3403 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3404
3405 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3406 AliTPCseed *seed=new AliTPCseed(0, x, c, x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3407 // Double_t y0,z0,y1,z1, y2,z2;
3408 //seed->GetProlongation(x0[0],y0,z0);
3409 // seed->GetProlongation(x1[0],y1,z1);
3410 //seed->GetProlongation(x2[0],y2,z2);
3411 // seed =0;
3412 seed->fLastPoint = pp2;
3413 seed->fFirstPoint = pp2;
3414
3415
3416 return seed;
3417}
3418
3419Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed, Float_t th)
3420{
3421 //
3422 //
3423 //
3424 for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3425 //
3426 if (TMath::Abs(seed->GetC())>0.01) return 0;
3427 //
3428
3429 Float_t x[160], y[160], erry[160], z[160], errz[160];
3430 Int_t sec[160];
3431 Float_t xt[160], yt[160], zt[160];
3432 Int_t i1 = 200;
3433 Int_t i2 = 0;
3434 Int_t secm = -1;
3435 Int_t padm = -1;
3436 Int_t middle = seed->GetNumberOfClusters()/2;
3437 //
3438 //
3439 // find central sector, get local cooordinates
3440 Int_t count = 0;
3441 for (Int_t i=seed->fFirstPoint;i<=seed->fLastPoint;i++) {
3442 sec[i]= seed->GetClusterSector(i)%18;
3443 x[i] = GetXrow(i);
3444 if (sec[i]>=0) {
3445 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3446 // if (cl==0) cl = GetClusterMI(seed->GetClusterIndex2(i));
3447 if (cl==0) {
3448 sec[i] = -1;
3449 continue;
3450 }
3451 //
3452 //
3453 if (i>i2) i2 = i; //last point with cluster
3454 if (i2<i1) i1 = i; //first point with cluster
3455 y[i] = cl->GetY();
3456 z[i] = cl->GetZ();
3457 AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3458 xt[i] = x[i];
3459 yt[i] = point->GetY();
3460 zt[i] = point->GetZ();
3461
3462 if (point->GetX()>0){
3463 erry[i] = point->GetErrY();
3464 errz[i] = point->GetErrZ();
3465 }
3466
3467 count++;
3468 if (count<middle) {
3469 secm = sec[i]; //central sector
3470 padm = i; //middle point with cluster
3471 }
3472 }
3473 }
3474 //
3475 // rotate position to global coordinate system connected to sector at last the point
3476 //
3477 for (Int_t i=i1;i<=i2;i++){
3478 //
3479 if (sec[i]<0) continue;
3480 Double_t alpha = (sec[i2]-sec[i])*fSectors->GetAlpha();
3481 Double_t cs = TMath::Cos(alpha);
3482 Double_t sn = TMath::Sin(alpha);
3483 Float_t xx2= x[i]*cs+y[i]*sn;
3484 Float_t yy2= -x[i]*sn+y[i]*cs;
3485 x[i] = xx2;
3486 y[i] = yy2;
3487 //
3488 xx2= xt[i]*cs+yt[i]*sn;
3489 yy2= -xt[i]*sn+yt[i]*cs;
3490 xt[i] = xx2;
3491 yt[i] = yy2;
3492
3493 }
3494 //get "state" vector
3495 Double_t xh[5],xm = x[padm];
3496 xh[0]=yt[i2];
3497 xh[1]=zt[i2];
3498 xh[4]=f1(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3499 xh[2]=f2(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3500 xh[3]=f3n(xt[i2],yt[i2],xt[i1],yt[i1],zt[i2],zt[i1],xh[4]);
3501 //
3502 //
3503 for (Int_t i=i1;i<=i2;i++){
3504 Double_t yy,zz;
3505 if (sec[i]<0) continue;
3506 GetProlongation(x[i2], x[i],xh,yy,zz);
3507 if (TMath::Abs(y[i]-yy)>4||TMath::Abs(z[i]-zz)>4){
3508 //Double_t xxh[5];
3509 //xxh[4]=f1old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3510 //xxh[2]=f2old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3511 printf("problem\n");
3512 }
3513 y[i] = y[i] - yy;
3514 z[i] = z[i] - zz;
3515 }
3516 Float_t dyup[160],dydown[160], dzup[160], dzdown[160];
3517 Float_t yup[160], ydown[160], zup[160], zdown[160];
3518
3519 AliTPCpolyTrack ptrack1,ptrack2;
3520 //
3521 // derivation up
3522 for (Int_t i=i1;i<=i2;i++){
3523 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3524 if (!cl) continue;
3525 if (cl->GetType()<0) continue;
3526 if (cl->GetType()>10) continue;
3527
3528 if (sec[i]>=0){
3529 ptrack1.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3530 }
3531 if (ptrack1.GetN()>4.){
3532 ptrack1.UpdateParameters();
3533 Double_t ddy,ddz;
3534 ptrack1.GetFitDerivation(x[i]-xm,ddy,ddz);
3535 Double_t yy,zz;
3536 ptrack1.GetFitPoint(x[i]-xm,yy,zz);
3537
3538 dyup[i] = ddy;
3539 dzup[i] = ddz;
3540 yup[i] = yy;
3541 zup[i] = zz;
3542
3543 }
3544 else{
3545 dyup[i]=0.; //not enough points
3546 }
3547 }
3548 //
3549 // derivation down
3550 for (Int_t i=i2;i>=i1;i--){
3551 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3552 if (!cl) continue;
3553 if (cl->GetType()<0) continue;
3554 if (cl->GetType()>10) continue;
3555 if (sec[i]>=0){
3556 ptrack2.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3557 }
3558 if (ptrack2.GetN()>4){
3559 ptrack2.UpdateParameters();
3560 Double_t ddy,ddz;
3561 ptrack2.GetFitDerivation(x[i]-xm,ddy,ddz);
3562 Double_t yy,zz;
3563 ptrack2.GetFitPoint(x[i]-xm,yy,zz);
3564
3565 dydown[i] = ddy;
3566 dzdown[i] = ddz;
3567 ydown[i] = yy;
3568 zdown[i] = zz;
3569 }
3570 else{
3571 dydown[i]=0.; //not enough points
3572 }
3573 }
3574 //
3575 //
3576 // find maximal difference of the derivation
3577 for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3578
3579
3580 for (Int_t i=i1+10;i<i2-10;i++){
3581 if ( (TMath::Abs(dydown[i])<0.00000001) || (TMath::Abs(dyup[i])<0.00000001) ||i<30)continue;
3582 // printf("%f\t%f\t%f\t%f\t%f\n",x[i],dydown[i],dyup[i],dzdown[i],dzup[i]);
3583 //
3584 Float_t ddy = TMath::Abs(dydown[i]-dyup[i]);
3585 Float_t ddz = TMath::Abs(dzdown[i]-dzup[i]);
3586 if ( (ddy+ddz)> th){
3587 seed->fKinkPoint[0] = i;
3588 seed->fKinkPoint[1] = ddy;
3589 seed->fKinkPoint[2] = ddz;
3590 th = ddy+ddz;
3591 }
3592 }
3593
3594 if (fTreeDebug){
3595 //
3596 //write information to the debug tree
3597 TBranch * br = fTreeDebug->GetBranch("debug");
3598 TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
3599 arr->ExpandCreateFast(i2-i1);
3600 br->SetAddress(&arr);
3601 //
3602 AliTPCclusterMI cldummy;
3603 cldummy.SetQ(0);
3604 AliTPCTrackPoint2 pdummy;
3605 pdummy.GetTPoint().fIsShared = 10;
3606 //
3607 Double_t alpha = sec[i2]*fSectors->GetAlpha();
3608 Double_t cs = TMath::Cos(alpha);
3609 Double_t sn = TMath::Sin(alpha);
3610
3611 for (Int_t i=i1;i<i2;i++){
3612 AliTPCTrackPoint2 *trpoint = (AliTPCTrackPoint2*)arr->UncheckedAt(i-i1);
3613 //cluster info
3614 AliTPCclusterMI * cl0 = seed->fClusterPointer[i];
3615 //
3616 AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3617
3618 if (cl0){
3619 Double_t x = GetXrow(i);
3620 trpoint->GetTPoint() = *point;
3621 trpoint->GetCPoint() = *cl0;
3622 trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
3623 trpoint->fID = seed->GetUniqueID();
3624 trpoint->fLab = seed->GetLabel();
3625 //
3626 trpoint->fGX = cs *x + sn*point->GetY();
3627 trpoint->fGY = -sn *x + cs*point->GetY() ;
3628 trpoint->fGZ = point->GetZ();
3629 //
3630 trpoint->fDY = y[i];
3631 trpoint->fDZ = z[i];
3632 //
3633 trpoint->fDYU = dyup[i];
3634 trpoint->fDZU = dzup[i];
3635 //
3636 trpoint->fDYD = dydown[i];
3637 trpoint->fDZD = dzdown[i];
3638 //
3639 if (TMath::Abs(dyup[i])>0.00000000001 &&TMath::Abs(dydown[i])>0.00000000001){
3640 trpoint->fDDY = dydown[i]-dyup[i];
3641 trpoint->fDDZ = dzdown[i]-dzup[i];
3642 }else{
3643 trpoint->fDDY = 0.;
3644 trpoint->fDDZ = 0.;
3645 }
3646 }
3647 else{
3648 *trpoint = pdummy;
3649 trpoint->GetCPoint()= cldummy;
3650 trpoint->fID = -1;
3651 }
3652 //
1627d1c4 3653 }
91162307 3654 fTreeDebug->Fill();
1627d1c4 3655 }
3656
3657
91162307 3658 return 0;
1627d1c4 3659
3660}
3661
3662
3663
3664
3665
91162307 3666AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
3667{
3668 //
3669 // reseed - refit - track
3670 //
3671 Int_t first = 0;
3672 // Int_t last = fSectors->GetNRows()-1;
3673 //
3674 if (fSectors == fOuterSec){
3675 first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows());
3676 //last =
3677 }
3678 else
3679 first = t->fFirstPoint;
3680 //
3681 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
3682 FollowBackProlongation(*t,fSectors->GetNRows()-1);
3683 t->Reset(kFALSE);
3684 FollowProlongation(*t,first);
3685 return seed;
3686}
3687
3688
3689
3690
3691
3692
3693
1c53abe2 3694//_____________________________________________________________________________
3695Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
3696 //-----------------------------------------------------------------
3697 // This function reades track seeds.
3698 //-----------------------------------------------------------------
3699 TDirectory *savedir=gDirectory;
3700
3701 TFile *in=(TFile*)inp;
3702 if (!in->IsOpen()) {
3703 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
3704 return 1;
3705 }
3706
3707 in->cd();
3708 TTree *seedTree=(TTree*)in->Get("Seeds");
3709 if (!seedTree) {
3710 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
3711 cerr<<"can't get a tree with track seeds !\n";
3712 return 2;
3713 }
3714 AliTPCtrack *seed=new AliTPCtrack;
3715 seedTree->SetBranchAddress("tracks",&seed);
3716
3717 if (fSeeds==0) fSeeds=new TObjArray(15000);
3718
3719 Int_t n=(Int_t)seedTree->GetEntries();
3720 for (Int_t i=0; i<n; i++) {
3721 seedTree->GetEvent(i);
3722 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
3723 }
3724
3725 delete seed;
3726 delete seedTree;
3727 savedir->cd();
3728 return 0;
3729}
3730
3731//_____________________________________________________________________________
f8aae377 3732Int_t AliTPCtrackerMI::Clusters2Tracks() {
1c53abe2 3733 //-----------------------------------------------------------------
3734 // This is a track finder.
3735 //-----------------------------------------------------------------
91162307 3736 TDirectory *savedir=gDirectory;
1c53abe2 3737 TStopwatch timer;
91162307 3738 //
3739 if (!fInput) SetIO(); //set default IO using loaders
3740 if (!fInput){
3741 cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
3742 return 1;
3743 }
1c53abe2 3744 LoadClusters();
91162307 3745 //
3746 fIteration = 0;
3747 fSeeds = Tracking();
1c53abe2 3748
1c53abe2 3749
91162307 3750 printf("Time for tracking: \t");timer.Print();timer.Start();
3751
3752 //activate again some tracks
3753 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
3754 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3755 if (!pt) continue;
3756 Int_t nc=t.GetNumberOfClusters();
3757 if (nc<20) {
3758 delete fSeeds->RemoveAt(i);
3759 continue;
3760 }
3761 if (pt->fRemoval==10) {
3762 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
3763 pt->Desactivate(10); // make track again active
3764 else{
3765 pt->Desactivate(20);
3766 delete fSeeds->RemoveAt(i);
3767 }
3768 }
3769 }
3770 RemoveDouble(fSeeds,0.2,0.6,11);
3771 //RemoveUsed(fSeeds,0.9,0.9,6);
3772 //RemoveUsed(fSeeds,0.8,0.8,6);
3773 //RemoveUsed(fSeeds,0.7,0.7,6);
3774 RemoveUsed(fSeeds,0.5,0.5,6);
c9427e08 3775
1c53abe2 3776 //
91162307 3777 Int_t nseed=fSeeds->GetEntriesFast();
3778 Int_t found = 0;
3779 for (Int_t i=0; i<nseed; i++) {
3780 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3781 if (!pt) continue;
3782 Int_t nc=t.GetNumberOfClusters();
3783 if (nc<15) {
3784 delete fSeeds->RemoveAt(i);
3785 continue;
3786 }
3787 CookLabel(pt,0.1); //For comparison only
3788 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3789 if ((pt->IsActive() || (pt->fRemoval==10) )){
3790 cerr<<found++<<'\r';
3791 }
3792 else
3793 delete fSeeds->RemoveAt(i);
3794 pt->fLab2 = i;
3795 }
3796
3797
3798 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
3799 SignShared(fSeeds);
3800 //RemoveUsed(fSeeds,0.9,0.9,6);
1c53abe2 3801 //
91162307 3802 nseed=fSeeds->GetEntriesFast();
3803 found = 0;
3804 for (Int_t i=0; i<nseed; i++) {
3805 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3806 if (!pt) continue;
3807 Int_t nc=t.GetNumberOfClusters();
3808 if (nc<15) {
3809 delete fSeeds->RemoveAt(i);
3810 continue;
3811 }
3812 t.SetUniqueID(i);
3813 t.CookdEdx(0.02,0.6);
3814 // CheckKinkPoint(&t,0.05);
3815 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3816 if ((pt->IsActive() || (pt->fRemoval==10) )){
3817 cerr<<found++<<'\r';
3818 }
3819 else
3820 delete fSeeds->RemoveAt(i);
3821 pt->fLab2 = i;
3822 }
3823
3824 SortTracks(fSeeds, 1);
1c53abe2 3825
91162307 3826 /*
3827 fIteration = 1;
3828 PrepareForBackProlongation(fSeeds,0.5);
3829 PropagateBack(fSeeds);
3830 printf("Time for back propagation: \t");timer.Print();timer.Start();
3831
3832 fIteration = 2;
1c53abe2 3833
91162307 3834 PrepareForProlongation(fSeeds,1.);
3835 PropagateForward();
3836
3837 fSectors = fOuterSec;
3838 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
3839 fSectors = fInnerSec;
3840 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
3841 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
3842 // RemoveUsed(fSeeds,0.7,0.7,6);
3843 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
3844
1c53abe2 3845 nseed=fSeeds->GetEntriesFast();
91162307 3846 found = 0;
3847 for (Int_t i=0; i<nseed; i++) {
1c53abe2 3848 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3849 if (!pt) continue;
3850 Int_t nc=t.GetNumberOfClusters();
91162307 3851 if (nc<15) {
3852 delete fSeeds->RemoveAt(i);
3853 continue;
3854 }
1c53abe2 3855 t.CookdEdx(0.02,0.6);
91162307 3856 // CookLabel(pt,0.1); //For comparison only
3857 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3858 if ((pt->IsActive() || (pt->fRemoval==10) )){
3859 cerr<<found++<<'\r';
3860 }
3861 else
3862 delete fSeeds->RemoveAt(i);
3863 pt->fLab2 = i;
1c53abe2 3864 }
91162307 3865 */
3866
c9427e08 3867 // fNTracks = found;
91162307 3868 printf("Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
3869 //
3870 if (fOutput) {
3871 WriteTracks();
3872 }
3873 if (!fNewIO) fOutput->Write();
3874 else
3875 AliRunLoader::GetDetectorLoader("TPC",AliConfig::fgkDefaultEventFolderName)->WriteTracks("OVERWRITE");
1c53abe2 3876
1c53abe2 3877
91162307 3878 cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
3879 savedir->cd();
3880 //if (seedtree) delete seedtree;
3881 // UnloadClusters();
3882 //printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1c53abe2 3883
1c53abe2 3884 return 0;
3885}
3886
91162307 3887void AliTPCtrackerMI::Tracking(TObjArray * arr)
3888{
3889 //
3890 // tracking of the seeds
3891 //
3892
3893 fSectors = fOuterSec;
3894 ParallelTracking(arr,150,63);
3895 fSectors = fOuterSec;
3896 ParallelTracking(arr,63,0);
3897}
3898
3899TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
3900{
3901 //
3902 //
3903 //tracking routine
3904 TObjArray * arr = new TObjArray;
3905 //
3906 fSectors = fOuterSec;
3907 TStopwatch timer;
3908 timer.Start();
3909 for (Int_t sec=0;sec<fkNOS;sec++){
3910 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
3911 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
3912 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
3913 }
3914 if (fDebug>0){
3915 printf("\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
3916 timer.Print();
3917 timer.Start();
3918 }
3919 Tracking(arr);
3920 if (fDebug>0){
3921 timer.Print();
3922 }
3923
3924 return arr;
3925}
3926
3927TObjArray * AliTPCtrackerMI::Tracking()
3928{
3929 //
3930 //
3931 TStopwatch timer;
3932 timer.Start();
3933 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
3934
3935 TObjArray * seeds = new TObjArray;
3936 TObjArray * arr=0;
3937
3938 Int_t gap =20;
3939 Float_t cuts[4];
3940 cuts[0] = 0.002;
3941 cuts[1] = 1.5;
3942 cuts[2] = 3.;
3943 cuts[3] = 3.;
3944 Float_t fnumber = 3.0;
3945 Float_t fdensity = 3.0;
3946
3947 //
3948 //find primaries
3949 cuts[0]=0.0066;
3950 for (Int_t delta = 0; delta<18; delta+=6){
3951 //
3952 cuts[0]=0.0070;
3953 cuts[1] = 1.5;
3954 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
3955 SumTracks(seeds,arr);
3956 SignClusters(seeds,fnumber,fdensity);
3957 //
3958 for (Int_t i=2;i<6;i+=2){
3959 // seed high pt tracks
3960 cuts[0]=0.0022;
3961 cuts[1]=0.3;
3962 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
3963 SumTracks(seeds,arr);
3964 SignClusters(seeds,fnumber,fdensity);
3965 }
3966 }
3967 fnumber = 4;
3968 fdensity = 4.;
3969 // RemoveUsed(seeds,0.9,0.9,1);
3970 // UnsignClusters();
3971 // SignClusters(seeds,fnumber,fdensity);
3972
3973 //find primaries
3974 cuts[0]=0.0077;
3975 for (Int_t delta = 20; delta<120; delta+=10){
3976 //
3977 // seed high pt tracks
3978 cuts[0]=0.0060;
3979 cuts[1]=0.3;
3980 cuts[2]=6.;
3981 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
3982 SumTracks(seeds,arr);
3983 SignClusters(seeds,fnumber,fdensity);
3984
3985 cuts[0]=0.003;
3986 cuts[1]=0.3;
3987 cuts[2]=6.;
3988 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
3989 SumTracks(seeds,arr);
3990 SignClusters(seeds,fnumber,fdensity);
3991 }
3992
3993 cuts[0] = 0.01;
3994 cuts[1] = 2.0;
3995 cuts[2] = 3.;
3996 cuts[3] = 2.0;
3997 fnumber = 2.;
3998 fdensity = 2.;
3999
4000 if (fDebug>0){
4001 printf("\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
4002 timer.Print();
4003 timer.Start();
4004 }
4005 // RemoveUsed(seeds,0.75,0.75,1);
4006 //UnsignClusters();
4007 //SignClusters(seeds,fnumber,fdensity);
4008
4009 // find secondaries
4010
4011 cuts[0] = 0.3;
4012 cuts[1] = 1.5;
4013 cuts[2] = 3.;
4014 cuts[3] = 1.5;
4015
4016 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
4017 SumTracks(seeds,arr);
4018 SignClusters(seeds,fnumber,fdensity);
4019 //
4020 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
4021 SumTracks(seeds,arr);
4022 SignClusters(seeds,fnumber,fdensity);
4023 //
4024 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
4025 SumTracks(seeds,arr);
4026 SignClusters(seeds,fnumber,fdensity);
4027 //
4028
4029
4030 for (Int_t delta = 3; delta<30; delta+=5){
4031 //
4032 cuts[0] = 0.3;
4033 cuts[1] = 1.5;
4034 cuts[2] = 3.;
4035 cuts[3] = 1.5;
4036 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4037 SumTracks(seeds,arr);
4038 SignClusters(seeds,fnumber,fdensity);
4039 //
4040 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
4041 SumTracks(seeds,arr);
4042 SignClusters(seeds,fnumber,fdensity);
4043 //
4044 }
4045 fnumber = 1;
4046 fdensity = 1;
4047 //
4048 // change cuts
4049 fnumber = 2.;
4050 fdensity = 2.;
4051 cuts[0]=0.0080;
4052
4053 // find secondaries
4054 for (Int_t delta = 30; delta<70; delta+=10){
4055 //
4056 cuts[0] = 0.3;
4057 cuts[1] = 1.5;
4058 cuts[2] = 3.;
4059 cuts[3] = 1.5;
4060 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4061 SumTracks(seeds,arr);
4062 SignClusters(seeds,fnumber,fdensity);
4063 //
4064 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
4065 SumTracks(seeds,arr);
4066 SignClusters(seeds,fnumber,fdensity);
4067 }
4068
4069 if (fDebug>0){
4070 printf("\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
4071 timer.Print();
4072 timer.Start();
4073 }
4074
4075 return seeds;
4076 //
4077
4078}
4079
4080
4081void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2)
4082{
4083 //
4084 //sum tracks to common container
4085 //remove suspicious tracks
4086 Int_t nseed = arr2->GetEntriesFast();
4087 for (Int_t i=0;i<nseed;i++){
4088 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
4089 if (pt){
4090
4091 // NORMAL ACTIVE TRACK
4092 if (pt->IsActive()){
4093 arr1->AddLast(arr2->RemoveAt(i));
4094 continue;
4095 }
4096 //remove not usable tracks
4097 if (pt->fRemoval!=10){
4098 delete arr2->RemoveAt(i);
4099 continue;
4100 }
4101 // REMOVE VERY SHORT TRACKS
4102 if (pt->GetNumberOfClusters()<20){
4103 delete arr2->RemoveAt(i);
4104 continue;
4105 }
4106 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
4107 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
4108 arr1->AddLast(arr2->RemoveAt(i));
4109 else{
4110 delete arr2->RemoveAt(i);
4111 }
4112 }
4113 }
4114 delete arr2;
4115}
4116
4117
1c53abe2 4118
91162307 4119void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
1c53abe2 4120{
4121 //
4122 // try to track in parralel
4123
91162307 4124 Int_t nseed=arr->GetEntriesFast();
1c53abe2 4125 //prepare seeds for tracking
4126 for (Int_t i=0; i<nseed; i++) {
91162307 4127 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1c53abe2 4128 if (!pt) continue;
4129 if (!t.IsActive()) continue;
4130 // follow prolongation to the first layer
91162307 4131 if ( (fSectors ==fInnerSec) || (t.fFirstPoint-fParam->GetNRowLow()>rfirst+1) )
c9427e08 4132 FollowProlongation(t, rfirst+1);
1c53abe2 4133 }
4134
4135
4136 //
4137 for (Int_t nr=rfirst; nr>=rlast; nr--){
4138 // make indexes with the cluster tracks for given
1c53abe2 4139
4140 // find nearest cluster
4141 for (Int_t i=0; i<nseed; i++) {
91162307 4142 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1c53abe2 4143 if (!pt) continue;
4144 if (!pt->IsActive()) continue;
91162307 4145 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
1627d1c4 4146 if (pt->fRelativeSector>17) {
4147 continue;
4148 }
91162307 4149 UpdateClusters(t,nr);
1c53abe2 4150 }
4151 // prolonagate to the nearest cluster - if founded
4152 for (Int_t i=0; i<nseed; i++) {
91162307 4153 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1c53abe2 4154 if (!pt) continue;
1627d1c4 4155 if (!pt->IsActive()) continue;
91162307 4156 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
1627d1c4 4157 if (pt->fRelativeSector>17) {
4158 continue;
4159 }
91162307 4160 FollowToNextCluster(*pt,nr);
1c53abe2 4161 }
91162307 4162 }
4163}
4164
4165void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac)
4166{
4167 //
4168 //
4169 // if we use TPC track itself we have to "update" covariance
4170 //
4171 Int_t nseed= arr->GetEntriesFast();
4172 for (Int_t i=0;i<nseed;i++){
4173 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4174 if (pt) {
4175 pt->Modify(fac);
4176 //
4177 //rotate to current local system at first accepted point
4178 Int_t index = pt->GetClusterIndex2(pt->fFirstPoint);
4179 Int_t sec = (index&0xff000000)>>24;
4180 sec = sec%18;
4181 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
4182 if (angle1>TMath::Pi())
4183 angle1-=2.*TMath::Pi();
4184 Float_t angle2 = pt->GetAlpha();
4185
4186 if (TMath::Abs(angle1-angle2)>0.001){
4187 pt->Rotate(angle1-angle2);
4188 //angle2 = pt->GetAlpha();
4189 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
4190 //if (pt->GetAlpha()<0)
4191 // pt->fRelativeSector+=18;
4192 //sec = pt->fRelativeSector;
4193 }
4194
4195 }
4196
4197 }
4198
4199
4200}
4201void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac)
4202{
4203 //
4204 //
4205 // if we use TPC track itself we have to "update" covariance
4206 //
4207 Int_t nseed= arr->GetEntriesFast();
4208 for (Int_t i=0;i<nseed;i++){
4209 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4210 if (pt) {
4211 pt->Modify(fac);
4212 pt->fFirstPoint = pt->fLastPoint;
4213 }
4214
4215 }
4216
4217
4218}
4219
4220Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
4221{
4222 //
4223 // make back propagation
4224 //
4225 Int_t nseed= arr->GetEntriesFast();
4226 for (Int_t i=0;i<nseed;i++){
4227 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4228 if (pt) {
4229 AliTPCseed *pt2 = new AliTPCseed(*pt);
4230 fSectors = fInnerSec;
4231 FollowBackProlongation(*pt,fSectors->GetNRows()-1);
4232 fSectors = fOuterSec;
4233 FollowBackProlongation(*pt,fSectors->GetNRows()-1);
4234 fSectors = fOuterSec;
4235 if (pt->GetNumberOfClusters()<35 && pt->GetLabel()>0 ){
4236 printf("\n%d",pt->GetLabel());
4237 fSectors = fInnerSec;
4238 FollowBackProlongation(*pt2,fSectors->GetNRows()-1);
4239 fSectors = fOuterSec;
4240 FollowBackProlongation(*pt2,fSectors->GetNRows()-1);
4241 fSectors = fOuterSec;
4242 }
4243 }
4244 }
4245 return 0;
4246}
4247
4248
4249Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
4250{
4251 //
4252 // make forward propagation
4253 //
4254 Int_t nseed= arr->GetEntriesFast();
4255 for (Int_t i=0;i<nseed;i++){
4256 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4257 if (pt) {
4258 AliTPCseed *pt2 = new AliTPCseed(*pt);
4259 fSectors = fOuterSec;
4260 FollowProlongation(*pt,0);
4261 fSectors = fOuterSec;
4262 FollowProlongation(*pt,0);
4263 fSectors = fInnerSec;
4264 if (pt->GetNumberOfClusters()<35 && pt->GetLabel()>0 ){
4265 printf("\n%d",pt->GetLabel());
4266 fSectors = fOuterSec;
4267 FollowProlongation(*pt2,0);
4268 fSectors = fOuterSec;
4269 FollowProlongation(*pt2,0);
4270 fSectors = fOuterSec;
4271 }
4272 }
4273 }
4274 return 0;
4275}
4276
4277
4278Int_t AliTPCtrackerMI::PropagateForward()
4279{
4280 fSectors = fOuterSec;
4281 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
4282 fSectors = fInnerSec;
4283 ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
4284 //WriteTracks();
4285 return 1;
4286}
4287
4288
4289
4290
4291
4292
4293Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
4294{
4295 //
4296 // make back propagation, in between row0 and row1
4297 //
1c53abe2 4298
91162307 4299 if (pt) {
4300 fSectors = fInnerSec;
4301 Int_t r1;
4302 //
4303 if (row1<fSectors->GetNRows())
4304 r1 = row1;
4305 else
4306 r1 = fSectors->GetNRows()-1;
4307
4308 if (row0<fSectors->GetNRows()&& r1>0 )
4309 FollowBackProlongation(*pt,r1);
4310 if (row1<=fSectors->GetNRows())
4311 return 0;
4312 //
4313 r1 = row1 - fSectors->GetNRows();
4314 if (r1<=0) return 0;
4315 if (r1>=fOuterSec->GetNRows()) return 0;
4316 fSectors = fOuterSec;
4317 return FollowBackProlongation(*pt,r1);
4318 }
4319 return 0;
4320}
4321
4322
4323
4324
4325void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
4326{
4327 //
4328 //
4329 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4330 // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4331 Float_t padlength = GetPadPitchLength(row);
4332 //
4333 Float_t sresy = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4334 Float_t angulary = seed->GetSnp();
4335 angulary = angulary*angulary/(1-angulary*angulary);
4336 seed->fCurrentSigmaY2 = sd2+padlength*padlength*angulary/12.+sresy*sresy;
4337 //
4338 Float_t sresz = fParam->GetZSigma();
4339 Float_t angularz = seed->GetTgl();
4340 seed->fCurrentSigmaZ2 = sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz;
4341 /*
4342 Float_t wy = GetSigmaY(seed);
4343 Float_t wz = GetSigmaZ(seed);
4344 wy*=wy;
4345 wz*=wz;
4346 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
4347 printf("problem\n");
4348 }
4349 */
1c53abe2 4350}
4351
91162307 4352
1c53abe2 4353Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
4354{
4355 //
4356 //
c9427e08 4357 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
1c53abe2 4358 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4359 Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4360 Float_t angular = seed->GetSnp();
4361 angular = angular*angular/(1-angular*angular);
4362 // angular*=angular;
4363 //angular = TMath::Sqrt(angular/(1-angular));
4364 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
4365 return res;
4366}
4367Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
4368{
4369 //
4370 //
c9427e08 4371 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
1c53abe2 4372 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4373 Float_t sres = fParam->GetZSigma();
4374 Float_t angular = seed->GetTgl();
4375 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
4376 return res;
4377}
4378
4379
4380
1c53abe2 4381
4382//__________________________________________________________________________
91162307 4383void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
1c53abe2 4384 //--------------------------------------------------------------------
4385 //This function "cooks" a track label. If label<0, this track is fake.
4386 //--------------------------------------------------------------------
4387 Int_t noc=t->GetNumberOfClusters();
91162307 4388 if (noc<10){
4389 printf("\nnot founded prolongation\n\n\n");
4390 t->Dump();
4391 return ;
4392 }
4393 Int_t lb[160];
4394 Int_t mx[160];
4395 AliTPCclusterMI *clusters[160];
4396 //
4397 for (Int_t i=0;i<160;i++) {
4398 clusters[i]=0;
4399 lb[i]=mx[i]=0;
4400 }
1c53abe2 4401
4402 Int_t i;
91162307 4403 Int_t current=0;
4404 for (i=0; i<160 && current<noc; i++) {
4405
4406 Int_t index=t->GetClusterIndex2(i);
4407 if (index<=0) continue;
4408 if (index&0x8000) continue;
4409 //
4410 //clusters[current]=GetClusterMI(index);
4411 if (t->fClusterPointer[i]){
4412 clusters[current]=t->fClusterPointer[i];
4413 current++;
4414 }
1c53abe2 4415 }
91162307 4416 noc = current;
1c53abe2 4417
4418 Int_t lab=123456789;
4419 for (i=0; i<noc; i++) {
4420 AliTPCclusterMI *c=clusters[i];
91162307 4421 if (!c) continue;
1c53abe2 4422 lab=TMath::Abs(c->GetLabel(0));
4423 Int_t j;
4424 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
4425 lb[j]=lab;
4426 (mx[j])++;
4427 }
4428
4429 Int_t max=0;
4430 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
4431
4432 for (i=0; i<noc; i++) {
9918f10a 4433 AliTPCclusterMI *c=clusters[i];
91162307 4434 if (!c) continue;
1c53abe2 4435 if (TMath::Abs(c->GetLabel(1)) == lab ||
4436 TMath::Abs(c->GetLabel(2)) == lab ) max++;
4437 }
4438
4439 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
4440
4441 else {
4442 Int_t tail=Int_t(0.10*noc);
4443 max=0;
91162307 4444 Int_t ind=0;
4445 for (i=1; i<=160&&ind<tail; i++) {
4446 // AliTPCclusterMI *c=clusters[noc-i];
4447 AliTPCclusterMI *c=clusters[i];
4448 if (!c) continue;
1c53abe2 4449 if (lab == TMath::Abs(c->GetLabel(0)) ||
4450 lab == TMath::Abs(c->GetLabel(1)) ||
4451 lab == TMath::Abs(c->GetLabel(2))) max++;
91162307 4452 ind++;
1c53abe2 4453 }
4454 if (max < Int_t(0.5*tail)) lab=-lab;
4455 }
4456
4457 t->SetLabel(lab);
4458
91162307 4459 // delete[] lb;
4460 //delete[] mx;
4461 //delete[] clusters;
1c53abe2 4462}
4463
4464//_________________________________________________________________________
4465void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
4466 //-----------------------------------------------------------------------
4467 // Setup inner sector
4468 //-----------------------------------------------------------------------
4469 if (f==0) {
4470 fAlpha=par->GetInnerAngle();
4471 fAlphaShift=par->GetInnerAngleShift();
4472 fPadPitchWidth=par->GetInnerPadPitchWidth();
4473 fPadPitchLength=par->GetInnerPadPitchLength();
4474 fN=par->GetNRowLow();
4475 fRow=new AliTPCRow[fN];
4476 for (Int_t i=0; i<fN; i++) {
4477 fRow[i].SetX(par->GetPadRowRadiiLow(i));
4478 fRow[i].fDeadZone =1.5; //1.5 cm of dead zone
4479 }
4480 } else {
4481 fAlpha=par->GetOuterAngle();
4482 fAlphaShift=par->GetOuterAngleShift();
4483 fPadPitchWidth = par->GetOuterPadPitchWidth();
4484 fPadPitchLength = par->GetOuter1PadPitchLength();
4485 f1PadPitchLength = par->GetOuter1PadPitchLength();
4486 f2PadPitchLength = par->GetOuter2PadPitchLength();
4487
4488 fN=par->GetNRowUp();
4489 fRow=new AliTPCRow[fN];
4490 for (Int_t i=0; i<fN; i++) {
4491 fRow[i].SetX(par->GetPadRowRadiiUp(i));
4492 fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone
4493 }
4494 }
4495}
4496
4497
4498AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
4499 //
1c53abe2 4500}
4501
4502
4503
1c53abe2 4504//_________________________________________________________________________
4505void
4506AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
4507 //-----------------------------------------------------------------------
4508 // Insert a cluster into this pad row in accordence with its y-coordinate
4509 //-----------------------------------------------------------------------
4510 if (fN==kMaxClusterPerRow) {
4511 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
4512 }
4513 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
4514 Int_t i=Find(c->GetZ());
4515 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
4516 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
4517 fIndex[i]=index; fClusters[i]=c; fN++;
4518}
4519
91162307 4520
1c53abe2 4521//___________________________________________________________________
4522Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
4523 //-----------------------------------------------------------------------
4524 // Return the index of the nearest cluster
4525 //-----------------------------------------------------------------------
4526 if (fN==0) return 0;
4527 if (z <= fClusters[0]->GetZ()) return 0;
4528 if (z > fClusters[fN-1]->GetZ()) return fN;
4529 Int_t b=0, e=fN-1, m=(b+e)/2;
4530 for (; b<e; m=(b+e)/2) {
4531 if (z > fClusters[m]->GetZ()) b=m+1;
4532 else e=m;
4533 }
4534 return m;
4535}
4536
4537
4538
1627d1c4 4539//___________________________________________________________________
4540AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
4541 //-----------------------------------------------------------------------
4542 // Return the index of the nearest cluster in z y
4543 //-----------------------------------------------------------------------
4544 Float_t maxdistance = roady*roady + roadz*roadz;
4545
4546 AliTPCclusterMI *cl =0;
4547 for (Int_t i=Find(z-roadz); i<fN; i++) {
4548 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4549 if (c->GetZ() > z+roadz) break;
4550 if ( (c->GetY()-y) > roady ) continue;
4551 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4552 if (maxdistance>distance) {
4553 maxdistance = distance;
4554 cl=c;
4555 }
4556 }
4557 return cl;
4558}
4559
91162307 4560AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4561{
4562 //-----------------------------------------------------------------------
4563 // Return the index of the nearest cluster in z y
4564 //-----------------------------------------------------------------------
4565 Float_t maxdistance = roady*roady + roadz*roadz;
4566 Int_t iz1 = TMath::Max(fFastCluster[Int_t(z-roadz+254.5)]-1,0);
4567 Int_t iz2 = TMath::Min(fFastCluster[Int_t(z+roadz+255.5)]+1,fN);
4568
4569 AliTPCclusterMI *cl =0;
4570 //FindNearest3(y,z,roady,roadz,index);
4571 // for (Int_t i=Find(z-roadz); i<fN; i++) {
4572 for (Int_t i=iz1; i<iz2; i++) {
4573 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4574 if (c->GetZ() > z+roadz) break;
4575 if ( c->GetY()-y > roady ) continue;
4576 if ( y-c->GetY() > roady ) continue;
4577 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4578 if (maxdistance>distance) {
4579 maxdistance = distance;
4580 cl=c;
4581 index =i;
4582 //roady = TMath::Sqrt(maxdistance);
4583 }
4584 }
4585 return cl;
4586}
4587
4588
4589
4590AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4591{
4592 //-----------------------------------------------------------------------
4593 // Return the index of the nearest cluster in z y
4594 //-----------------------------------------------------------------------
4595 Float_t maxdistance = roady*roady + roadz*roadz;
4596 // Int_t iz = Int_t(z+255.);
4597 AliTPCclusterMI *cl =0;
4598 for (Int_t i=Find(z-roadz); i<fN; i++) {
4599 //for (Int_t i=fFastCluster[iz-2]; i<fFastCluster[iz+2]; i++) {
4600 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4601 if (c->GetZ() > z+roadz) break;
4602 if ( c->GetY()-y > roady ) continue;
4603 if ( y-c->GetY() > roady ) continue;
4604 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4605 if (maxdistance>distance) {
4606 maxdistance = distance;
4607 cl=c;
4608 index =i;
4609 //roady = TMath::Sqrt(maxdistance);
4610 }
4611 }
4612 return cl;
4613}
1627d1c4 4614
4615
4616
4617
1c53abe2 4618AliTPCseed::AliTPCseed():AliTPCtrack(){
4619 //
4620 fRow=0;
4621 fRemoval =0;
91162307 4622 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
4623 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
4624
1c53abe2 4625 fPoints = 0;
4626 fEPoints = 0;
4627 fNFoundable =0;
4628 fNShared =0;
91162307 4629 // fTrackPoints =0;
1c53abe2 4630 fRemoval = 0;
91162307 4631 fSort =0;
4632 fFirstPoint =0;
4633 fNoCluster =0;
4634 fBSigned = kFALSE;
4635 fSeed1 =-1;
4636 fSeed2 =-1;
1c53abe2 4637}
4638
4639AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
4640 fPoints = 0;
4641 fEPoints = 0;
4642 fNShared =0;
91162307 4643 // fTrackPoints =0;
1c53abe2 4644 fRemoval =0;
91162307 4645 fSort =0;
4646 for (Int_t i=0;i<160;i++) {
4647 fClusterPointer[i] = 0;
4648 Int_t index = t.GetClusterIndex(i);
4649 if (index>0) {
4650 SetClusterIndex2(i,index);
4651 }
4652 else{
4653 SetClusterIndex2(i,-3);
4654 }
4655 }
4656 fFirstPoint =0;
4657 fNoCluster =0;
4658 fBSigned = kFALSE;
4659 fSeed1 =-1;
4660 fSeed2 =-1;
1c53abe2 4661}
4662
4663AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
4664 fRow=0;
91162307 4665 for (Int_t i=0;i<160;i++) {
4666 fClusterPointer[i] = 0;
4667 Int_t index = t.GetClusterIndex(i);
4668 SetClusterIndex2(i,index);
4669 }
4670
1c53abe2 4671 fPoints = 0;
4672 fEPoints = 0;
4673 fNFoundable =0;
4674 fNShared =0;
91162307 4675 // fTrackPoints =0;
1c53abe2 4676 fRemoval =0;
91162307 4677 fSort = 0;
4678 fFirstPoint =0;
4679 fNoCluster =0;
4680 fBSigned = kFALSE;
4681 fSeed1 =-1;
4682 fSeed2 =-1;
1c53abe2 4683}
4684
4685AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
4686 Double_t xr, Double_t alpha):
4687 AliTPCtrack(index, xx, cc, xr, alpha) {
4688 //
4689 //
4690 fRow =0;
91162307 4691 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
4692 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
1c53abe2 4693 fPoints = 0;
4694 fEPoints = 0;
4695 fNFoundable =0;
4696 fNShared = 0;
91162307 4697 // fTrackPoints =0;
1c53abe2 4698 fRemoval =0;
91162307 4699 fSort =0;
4700 fFirstPoint =0;
4701 // fHelixIn = new TClonesArray("AliHelix",0);
4702 //fHelixOut = new TClonesArray("AliHelix",0);
4703 fNoCluster =0;
4704 fBSigned = kFALSE;
4705 fSeed1 =-1;
4706 fSeed2 =-1;
1c53abe2 4707}
4708
4709AliTPCseed::~AliTPCseed(){
4710 if (fPoints) delete fPoints;
4711 fPoints =0;
91162307 4712 if (fEPoints) delete fEPoints;
1c53abe2 4713 fEPoints = 0;
91162307 4714 fNoCluster =0;
1c53abe2 4715}
4716
91162307 4717AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
1c53abe2 4718{
4719 //
4720 //
91162307 4721 return &fTrackPoints[i];
1c53abe2 4722}
4723
4724void AliTPCseed::RebuildSeed()
4725{
4726 //
4727 // rebuild seed to be ready for storing
91162307 4728 AliTPCclusterMI cldummy;
4729 cldummy.SetQ(0);
4730 AliTPCTrackPoint pdummy;
4731 pdummy.GetTPoint().fIsShared = 10;
1c53abe2 4732 for (Int_t i=0;i<160;i++){
91162307 4733 AliTPCclusterMI * cl0 = fClusterPointer[i];
4734 AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
4735 if (cl0){
4736 trpoint->GetTPoint() = *(GetTrackPoint(i));
4737 trpoint->GetCPoint() = *cl0;
4738 trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
4739 }
4740 else{
4741 *trpoint = pdummy;
4742 trpoint->GetCPoint()= cldummy;
4743 }
4744
4745 }
4746
4747}
4748
4749
4750Double_t AliTPCseed::GetDensityFirst(Int_t n)
4751{
4752 //
4753 //
4754 // return cluster for n rows bellow first point
4755 Int_t nfoundable = 1;
4756 Int_t nfound = 1;
4757 for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
4758 Int_t index = GetClusterIndex2(i);
4759 if (index!=-1) nfoundable++;
4760 if (index>0) nfound++;
1c53abe2 4761 }
91162307 4762 if (nfoundable<n) return 0;
4763 return Double_t(nfound)/Double_t(nfoundable);
4764
4765}
1c53abe2 4766
91162307 4767
4768void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
4769{
4770 // get cluster stat. on given region
4771 //
4772 found = 0;
4773 foundable = 0;
4774 shared =0;
4775 for (Int_t i=first;i<last; i++){
4776 Int_t index = GetClusterIndex2(i);
4777 if (index!=-1) foundable++;
4778 if (fClusterPointer[i]) {
4779 found++;
4780 }
4781 else
4782 continue;
4783
4784 if (fClusterPointer[i]->IsUsed(10)) {
4785 shared++;
4786 continue;
4787 }
4788 if (!plus2) continue; //take also neighborhoud
4789 //
4790 if ( (i>0) && fClusterPointer[i-1]){
4791 if (fClusterPointer[i-1]->IsUsed(10)) {
4792 shared++;
4793 continue;
4794 }
4795 }
4796 if ( fClusterPointer[i+1]){
4797 if (fClusterPointer[i+1]->IsUsed(10)) {
4798 shared++;
4799 continue;
4800 }
4801 }
4802
4803 }
4804 if (shared>found){
4805 printf("problem\n");
4806 }
1c53abe2 4807}
4808
4809//_____________________________________________________________________________
91162307 4810void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
1c53abe2 4811 //-----------------------------------------------------------------
4812 // This funtion calculates dE/dX within the "low" and "up" cuts.
4813 //-----------------------------------------------------------------
4814
4815 Float_t amp[200];
4816 Float_t angular[200];
4817 Float_t weight[200];
4818 Int_t index[200];
4819 //Int_t nc = 0;
4820 // TClonesArray & arr = *fPoints;
4821 Float_t meanlog = 100.;
4822
4823 Float_t mean[4] = {0,0,0,0};
4824 Float_t sigma[4] = {1000,1000,1000,1000};
4825 Int_t nc[4] = {0,0,0,0};
4826 Float_t norm[4] = {1000,1000,1000,1000};
4827 //
4828 //
4829 fNShared =0;
4830
4831 for (Int_t of =0; of<4; of++){
91162307 4832 for (Int_t i=of+i1;i<i2;i+=4)
1c53abe2 4833 {
91162307 4834 Int_t index = fIndex[i];
4835 if (index<0||index&0x8000) continue;
4836
1c53abe2 4837 //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
91162307 4838 AliTPCTrackerPoint * point = GetTrackPoint(i);
4839 //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
4840 //AliTPCTrackerPoint * pointp = 0;
4841 //if (i<159) pointp = GetTrackPoint(i+1);
4842
1c53abe2 4843 if (point==0) continue;
91162307 4844 AliTPCclusterMI * cl = fClusterPointer[i];
4845 if (cl==0) continue;
4846 if (onlyused && (!cl->IsUsed(10))) continue;
4847 if (cl->IsUsed(11)) {
1c53abe2 4848 fNShared++;
4849 continue;
4850 }
91162307 4851 Int_t type = cl->GetType();
4852 //if (point->fIsShared){
4853 // fNShared++;
4854 // continue;
4855 //}
4856 //if (pointm)
4857 // if (pointm->fIsShared) continue;
4858 //if (pointp)
4859 // if (pointp->fIsShared) continue;
4860
4861 if (type<0) continue;
4862 //if (type>10) continue;
4863 //if (point->GetErrY()==0) continue;
4864 //if (point->GetErrZ()==0) continue;
4865
4866 //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
4867 //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
4868 //if ((ddy*ddy+ddz*ddz)>10) continue;
4869
4870
4871 // if (point->GetCPoint().GetMax()<5) continue;
4872 if (cl->GetMax()<5) continue;
4873 Float_t angley = point->GetAngleY();
4874 Float_t anglez = point->GetAngleZ();
4875
4876 Float_t rsigmay2 = point->GetSigmaY();
4877 Float_t rsigmaz2 = point->GetSigmaZ();
4878 /*
4879 Float_t ns = 1.;
4880 if (pointm){
4881 rsigmay += pointm->GetTPoint().GetSigmaY();
4882 rsigmaz += pointm->GetTPoint().GetSigmaZ();
4883 ns+=1.;
4884 }
4885 if (pointp){
4886 rsigmay += pointp->GetTPoint().GetSigmaY();
4887 rsigmaz += pointp->GetTPoint().GetSigmaZ();
4888 ns+=1.;
4889 }
4890 rsigmay/=ns;
4891 rsigmaz/=ns;
4892 */
4893
4894 Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
1c53abe2 4895
4896 Float_t ampc = 0; // normalization to the number of electrons
4897 if (i>64){
91162307 4898 // ampc = 1.*point->GetCPoint().GetMax();
4899 ampc = 1.*cl->GetMax();
1c53abe2 4900 //ampc = 1.*point->GetCPoint().GetQ();
4901 // AliTPCClusterPoint & p = point->GetCPoint();
4902 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
4903 // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
4904 //Float_t dz =
4905 // TMath::Abs( Int_t(iz) - iz + 0.5);
4906 //ampc *= 1.15*(1-0.3*dy);
4907 //ampc *= 1.15*(1-0.3*dz);
1627d1c4 4908 // Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
4909 //ampc *=zfactor;
1c53abe2 4910 }
4911 else{
91162307 4912 //ampc = 1.0*point->GetCPoint().GetMax();
4913 ampc = 1.0*cl->GetMax();
1c53abe2 4914 //ampc = 1.0*point->GetCPoint().GetQ();
4915 //AliTPCClusterPoint & p = point->GetCPoint();
4916 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
4917 //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
4918 //Float_t dz =
4919 // TMath::Abs( Int_t(iz) - iz + 0.5);
4920
4921 //ampc *= 1.15*(1-0.3*dy);
4922 //ampc *= 1.15*(1-0.3*dz);
1627d1c4 4923 // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
4924 //ampc *=zfactor;
1c53abe2 4925
4926 }
4927 ampc *= 2.0; // put mean value to channel 50
4928 //ampc *= 0.58; // put mean value to channel 50
4929 Float_t w = 1.;
4930 // if (type>0) w = 1./(type/2.-0.5);
91162307 4931 // Float_t z = TMath::Abs(cl->GetZ());
1c53abe2 4932 if (i<64) {
4933 ampc /= 0.6;
1627d1c4 4934 //ampc /= (1+0.0008*z);
4935 } else
4936 if (i>128){
4937 ampc /=1.5;
4938 //ampc /= (1+0.0008*z);
4939 }else{
4940 //ampc /= (1+0.0008*z);
4941 }
4942
1c53abe2 4943 if (type<0) { //amp at the border - lower weight
4944 // w*= 2.;
1627d1c4 4945
1c53abe2 4946 continue;
4947 }
4948 if (rsigma>1.5) ampc/=1.3; // if big backround
4949 amp[nc[of]] = ampc;
4950 angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
4951 weight[nc[of]] = w;
4952 nc[of]++;
4953 }
4954
4955 TMath::Sort(nc[of],amp,index,kFALSE);
4956 Float_t sumamp=0;
4957 Float_t sumamp2=0;
4958 Float_t sumw=0;
4959 //meanlog = amp[index[Int_t(nc[of]*0.33)]];
91162307 4960 meanlog = 50;
1c53abe2 4961 for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
4962 Float_t ampl = amp[index[i]]/angular[index[i]];
4963 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
4964 //
4965 sumw += weight[index[i]];
4966 sumamp += weight[index[i]]*ampl;
4967 sumamp2 += weight[index[i]]*ampl*ampl;
4968 norm[of] += angular[index[i]]*weight[index[i]];
4969 }
4970 if (sumw<1){
4971 SetdEdx(0);
4972 }
4973 else {
4974 norm[of] /= sumw;
4975 mean[of] = sumamp/sumw;
4976 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
4977 if (sigma[of]>0.1)
4978 sigma[of] = TMath::Sqrt(sigma[of]);
4979 else
4980 sigma[of] = 1000;
4981
4982 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
4983 //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
4984 //mean *=(1-0.1*(norm-1.));
4985 }
4986 }
4987
4988 Float_t dedx =0;
4989 fSdEdx =0;
4990 fMAngular =0;
4991 // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
4992 // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
4993
4994
4995 // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
4996 // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
4997
4998 Int_t norm2 = 0;
4999 Int_t norm3 = 0;
5000 for (Int_t i =0;i<4;i++){
5001 if (nc[i]>2&&nc[i]<1000){
5002 dedx += mean[i] *nc[i];
5003 fSdEdx += sigma[i]*(nc[i]-2);
5004 fMAngular += norm[i] *nc[i];
5005 norm2 += nc[i];
5006 norm3 += nc[i]-2;
5007 }
5008 fDEDX[i] = mean[i];
5009 fSDEDX[i] = sigma[i];
5010 fNCDEDX[i]= nc[i];
5011 }
5012
5013 if (norm3>0){
5014 dedx /=norm2;
5015 fSdEdx /=norm3;
5016 fMAngular/=norm2;
5017 }
5018 else{
5019 SetdEdx(0);
5020 return;
5021 }
5022 // Float_t dedx1 =dedx;
91162307 5023 /*
1c53abe2 5024 dedx =0;
5025 for (Int_t i =0;i<4;i++){
5026 if (nc[i]>2&&nc[i]<1000){
5027 mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
5028 dedx += mean[i] *nc[i];
5029 }
5030 fDEDX[i] = mean[i];
5031 }
5032 dedx /= norm2;
91162307 5033 */
1c53abe2 5034
5035
5036 SetdEdx(dedx);
5037
5038 //mi deDX
5039
5040
5041
5042 //Very rough PID
5043 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
5044
5045 if (p<0.6) {
5046 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5047 if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
5048 SetMass(0.93827); return;
5049 }
5050
5051 if (p<1.2) {
5052 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5053 SetMass(0.93827); return;
5054 }
5055
5056 SetMass(0.13957); return;
5057
5058}
5059
5060
91162307 5061
5062/*
5063
5064
5065
5066void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
5067 //-----------------------------------------------------------------
5068 // This funtion calculates dE/dX within the "low" and "up" cuts.
5069 //-----------------------------------------------------------------
5070
5071 Float_t amp[200];
5072 Float_t angular[200];
5073 Float_t weight[200];
5074 Int_t index[200];
5075 Bool_t inlimit[200];
5076 for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
5077 for (Int_t i=0;i<200;i++) amp[i]=10000;
5078 for (Int_t i=0;i<200;i++) angular[i]= 1;;
5079
5080
5081 //
5082 Float_t meanlog = 100.;
5083 Int_t indexde[4]={0,64,128,160};
5084
5085 Float_t amean =0;
5086 Float_t asigma =0;
5087 Float_t anc =0;
5088 Float_t anorm =0;
5089
5090 Float_t mean[4] = {0,0,0,0};
5091 Float_t sigma[4] = {1000,1000,1000,1000};
5092 Int_t nc[4] = {0,0,0,0};
5093 Float_t norm[4] = {1000,1000,1000,1000};
5094 //
5095 //
5096 fNShared =0;
5097
5098 // for (Int_t of =0; of<3; of++){
5099 // for (Int_t i=indexde[of];i<indexde[of+1];i++)
5100 for (Int_t i =0; i<160;i++)
5101 {
5102 AliTPCTrackPoint * point = GetTrackPoint(i);
5103 if (point==0) continue;
5104 if (point->fIsShared){
5105 fNShared++;
5106 continue;
5107 }
5108 Int_t type = point->GetCPoint().GetType();
5109 if (type<0) continue;
5110 if (point->GetCPoint().GetMax()<5) continue;
5111 Float_t angley = point->GetTPoint().GetAngleY();
5112 Float_t anglez = point->GetTPoint().GetAngleZ();
5113 Float_t rsigmay = point->GetCPoint().GetSigmaY();
5114 Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
5115 Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
5116
5117 Float_t ampc = 0; // normalization to the number of electrons
5118 if (i>64){
5119 ampc = point->GetCPoint().GetMax();
5120 }
5121 else{
5122 ampc = point->GetCPoint().GetMax();
5123 }
5124 ampc *= 2.0; // put mean value to channel 50
5125 // ampc *= 0.565; // put mean value to channel 50
5126
5127 Float_t w = 1.;
5128 Float_t z = TMath::Abs(point->GetCPoint().GetZ());
5129 if (i<64) {
5130 ampc /= 0.63;
5131 } else
5132 if (i>128){
5133 ampc /=1.51;
5134 }
5135 if (type<0) { //amp at the border - lower weight
5136 continue;
5137 }
5138 if (rsigma>1.5) ampc/=1.3; // if big backround
5139 angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
5140 amp[i] = ampc/angular[i];
5141 weight[i] = w;
5142 anc++;
5143 }
5144
5145 TMath::Sort(159,amp,index,kFALSE);
5146 for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
5147 inlimit[index[i]] = kTRUE; // take all clusters
5148 }
5149
5150 // meanlog = amp[index[Int_t(anc*0.3)]];
5151 meanlog =10000.;
5152 for (Int_t of =0; of<3; of++){
5153 Float_t sumamp=0;
5154 Float_t sumamp2=0;
5155 Float_t sumw=0;
5156 for (Int_t i=indexde[of];i<indexde[of+1];i++)
5157 {
5158 if (inlimit[i]==kFALSE) continue;
5159 Float_t ampl = amp[i];
5160 ///angular[i];
5161 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
5162 //
5163 sumw += weight[i];
5164 sumamp += weight[i]*ampl;
5165 sumamp2 += weight[i]*ampl*ampl;
5166 norm[of] += angular[i]*weight[i];
5167 nc[of]++;
5168 }
5169 if (sumw<1){
5170 SetdEdx(0);
5171 }
5172 else {
5173 norm[of] /= sumw;
5174 mean[of] = sumamp/sumw;
5175 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
5176 if (sigma[of]>0.1)
5177 sigma[of] = TMath::Sqrt(sigma[of]);
5178 else
5179 sigma[of] = 1000;
5180 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
5181 }
5182 }
5183
5184 Float_t dedx =0;
5185 fSdEdx =0;
5186 fMAngular =0;
5187 //
5188 Int_t norm2 = 0;
5189 Int_t norm3 = 0;
5190 Float_t www[3] = {12.,14.,17.};
5191 //Float_t www[3] = {1.,1.,1.};
5192
5193 for (Int_t i =0;i<3;i++){
5194 if (nc[i]>2&&nc[i]<1000){
5195 dedx += mean[i] *nc[i]*www[i]/sigma[i];
5196 fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
5197 fMAngular += norm[i] *nc[i];
5198 norm2 += nc[i]*www[i]/sigma[i];
5199 norm3 += (nc[i]-2)*www[i]/sigma[i];
5200 }
5201 fDEDX[i] = mean[i];
5202 fSDEDX[i] = sigma[i];
5203 fNCDEDX[i]= nc[i];
5204 }
5205
5206 if (norm3>0){
5207 dedx /=norm2;
5208 fSdEdx /=norm3;
5209 fMAngular/=norm2;
5210 }
5211 else{
5212 SetdEdx(0);
5213 return;
5214 }
5215 // Float_t dedx1 =dedx;
5216
5217 dedx =0;
5218 Float_t norm4 = 0;
5219 for (Int_t i =0;i<3;i++){
5220 if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
5221 //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
5222 dedx += mean[i] *(nc[i])/(sigma[i]);
5223 norm4 += (nc[i])/(sigma[i]);
5224 }
5225 fDEDX[i] = mean[i];
5226 }
5227 if (norm4>0) dedx /= norm4;
5228
5229
5230
5231 SetdEdx(dedx);
5232
5233 //mi deDX
5234
5235}
5236
5237*/