Store seed in friends regardless on stream-level, since needed for calibration
[u/mrichter/AliRoot.git] / TPC / TPCrec / AliTPCtracker.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
1c53abe2 16
17//-------------------------------------------------------
18// Implementation of the TPC tracker
19//
20// Origin: Marian Ivanov Marian.Ivanov@cern.ch
21//
34acb742 22// AliTPC parallel tracker
6d493ea0 23//
dee67df8 24// The track fitting is based on Kalman filtering approach
6d493ea0 25// The track finding steps:
26// 1. Seeding - with and without vertex constraint
27// - seeding with vertex constain done at first n^2 proble
28// - seeding without vertex constraint n^3 problem
29// 2. Tracking - follow prolongation road - find cluster - update kalman track
6d493ea0 30// The seeding and tracking is repeated several times, in different seeding region.
31// This approach enables to find the track which cannot be seeded in some region of TPC
32// This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
33
34// With this approach we reach almost 100 % efficiency also for high occupancy events.
35// (If the seeding efficiency in a region is about 90 % than with logical or of several
36// regions we will reach 100% (in theory - supposing independence)
37
38// Repeating several seeding - tracking procedures some of the tracks can be find
39// several times.
40
41// The procedures to remove multi find tacks are impremented:
42// RemoveUsed2 - fast procedure n problem -
43// Algorithm - Sorting tracks according quality
44// remove tracks with some shared fraction
45// Sharing in respect to all tacks
46// Signing clusters in gold region
47// FindSplitted - slower algorithm n^2
48// Sort the tracks according quality
49// Loop over pair of tracks
50// If overlap with other track bigger than threshold - remove track
51//
52// FindCurling - Finds the pair of tracks which are curling
53// - About 10% of tracks can be find with this procedure
54// The combinatorial background is too big to be used in High
55// multiplicity environment
56// - n^2 problem - Slow procedure - currently it is disabled because of
57// low efficiency
58//
59// The number of splitted tracks can be reduced disabling the sharing of the cluster.
60// tpcRecoParam-> SetClusterSharing(kFALSE);
61// IT IS HIGHLY non recomended to use it in high flux enviroonment
62// Even using this switch some tracks can be found more than once
63// (because of multiple seeding and low quality tracks which will not cross full chamber)
64//
65//
66// The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
67// To enable storage of the TPC tracks in the ESD friend track
16299eac 68// use AliTPCReconstructor::SetStreamLevel(n);
6d493ea0 69//
65e67c9c 70// The debug level - different procedure produce tree for numerical debugging of code and data (see comments foEStreamFlags in AliTPCtracker.h )
6d493ea0 71//
92f513f5 72
73//
74// Adding systematic errors to the covariance:
75//
76// The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
77// of the tracks (not to the clusters as they are dependent):
78// The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
b87c2bbc 79// The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/GeV)
92f513f5 80// The default values are 0.
81//
65e67c9c 82// The systematic errors are added to the covariance matrix in following places:
92f513f5 83//
829455ad 84// 1. During fisrt itteration - AliTPCtracker::FillESD
92f513f5 85// 2. Second iteration -
829455ad 86// 2.a ITS->TPC - AliTPCtracker::ReadSeeds
87// 2.b TPC->TRD - AliTPCtracker::PropagateBack
92f513f5 88// 3. Third iteration -
829455ad 89// 3.a TRD->TPC - AliTPCtracker::ReadSeeds
90// 3.b TPC->ITS - AliTPCtracker::RefitInward
92f513f5 91//
47966a6d 92
93/* $Id$ */
94
cc5e9db0 95#include "Riostream.h"
6d171107 96#include <TClonesArray.h>
97#include <TFile.h>
98#include <TObjArray.h>
99#include <TTree.h>
8cc1870a 100#include <TMatrixD.h>
6d64657a 101#include <TGraphErrors.h>
661f340b 102#include <TTimeStamp.h>
a3232aae 103#include "AliLog.h"
47966a6d 104#include "AliComplexCluster.h"
af885e0f 105#include "AliESDEvent.h"
aad72f45 106#include "AliESDtrack.h"
107#include "AliESDVertex.h"
6c94f330 108#include "AliKink.h"
109#include "AliV0.h"
91162307 110#include "AliHelix.h"
91162307 111#include "AliRunLoader.h"
6d171107 112#include "AliTPCClustersRow.h"
113#include "AliTPCParam.h"
9996a03b 114#include "AliTPCReconstructor.h"
6d171107 115#include "AliTPCpolyTrack.h"
81e97e0d 116#include "AliTPCreco.h"
9350f379 117#include "AliTPCseed.h"
118
119#include "AliTPCtrackerSector.h"
829455ad 120#include "AliTPCtracker.h"
6d171107 121#include "TStopwatch.h"
81e97e0d 122#include "AliTPCReconstructor.h"
5d837844 123#include "AliAlignObj.h"
124#include "AliTrackPointArray.h"
6d493ea0 125#include "TRandom.h"
24db6af7 126#include "AliTPCcalibDB.h"
6d64657a 127#include "AliTPCcalibDButil.h"
24db6af7 128#include "AliTPCTransform.h"
fd065ea2 129#include "AliTPCClusterParam.h"
64bf5ca0 130#include "AliTPCdEdxInfo.h"
ec7e4ad6 131#include "AliDCSSensorArray.h"
132#include "AliDCSSensor.h"
133#include "AliDAQ.h"
3aa6a136 134#include "AliCosmicTracker.h"
44ca7282 135#include "AliTPCROC.h"
e731a3f8 136#include "AliMathBase.h"
6d171107 137//
c9427e08 138
a11596ad 139using std::cerr;
140using std::endl;
829455ad 141ClassImp(AliTPCtracker)
c9427e08 142
143
002af263 144
b67e07dc 145class AliTPCFastMath {
91162307 146public:
b67e07dc 147 AliTPCFastMath();
91162307 148 static Double_t FastAsin(Double_t x);
149 private:
b67e07dc 150 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
91162307 151};
c9427e08 152
b67e07dc 153Double_t AliTPCFastMath::fgFastAsin[20000];
2274b54b 154AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
c9427e08 155
b67e07dc 156AliTPCFastMath::AliTPCFastMath(){
157 //
158 // initialized lookup table;
91162307 159 for (Int_t i=0;i<10000;i++){
160 fgFastAsin[2*i] = TMath::ASin(i/10000.);
161 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
162 }
c9427e08 163}
164
b67e07dc 165Double_t AliTPCFastMath::FastAsin(Double_t x){
166 //
167 // return asin using lookup table
91162307 168 if (x>0){
169 Int_t index = int(x*10000);
170 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
171 }
172 x*=-1;
173 Int_t index = int(x*10000);
174 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
1c53abe2 175}
e046d791 176//__________________________________________________________________
829455ad 177AliTPCtracker::AliTPCtracker()
e046d791 178 :AliTracker(),
179 fkNIS(0),
180 fInnerSec(0),
181 fkNOS(0),
182 fOuterSec(0),
183 fN(0),
184 fSectors(0),
185 fInput(0),
186 fOutput(0),
187 fSeedTree(0),
188 fTreeDebug(0),
189 fEvent(0),
72e25240 190 fEventHLT(0),
e046d791 191 fDebug(0),
192 fNewIO(kFALSE),
193 fNtracks(0),
194 fSeeds(0),
195 fIteration(0),
47af7ca4 196 fkParam(0),
1598ba75 197 fDebugStreamer(0),
f06a1ff6 198 fUseHLTClusters(4),
8cc1870a 199 fCrossTalkSignalArray(0),
f06a1ff6 200 fSeedsPool(0),
ddfbc51a 201 fFreeSeedsID(500),
202 fNFreeSeeds(0),
203 fLastSeedID(-1)
e046d791 204{
205 //
206 // default constructor
207 //
0f923679 208 for (Int_t irow=0; irow<200; irow++){
209 fXRow[irow]=0;
210 fYMax[irow]=0;
211 fPadLength[irow]=0;
212 }
e046d791 213}
214//_____________________________________________________________________
1c53abe2 215
216
217
829455ad 218Int_t AliTPCtracker::UpdateTrack(AliTPCseed * track, Int_t accept){
b67e07dc 219 //
220 //update track information using current cluster - track->fCurrentCluster
221
1c53abe2 222
b9671574 223 AliTPCclusterMI* c =track->GetCurrentCluster();
f124f8bf 224 if (accept > 0) //sign not accepted clusters
225 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
226 else // unsign accpeted clusters
227 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
b9671574 228 UInt_t i = track->GetCurrentClusterIndex1();
1c53abe2 229
230 Int_t sec=(i&0xff000000)>>24;
91162307 231 //Int_t row = (i&0x00ff0000)>>16;
b9671574 232 track->SetRow((i&0x00ff0000)>>16);
233 track->SetSector(sec);
1c53abe2 234 // Int_t index = i&0xFFFF;
47af7ca4 235 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
b9671574 236 track->SetClusterIndex2(track->GetRow(), i);
91162307 237 //track->fFirstPoint = row;
238 //if ( track->fLastPoint<row) track->fLastPoint =row;
239 // if (track->fRow<0 || track->fRow>160) {
240 // printf("problem\n");
241 //}
b9671574 242 if (track->GetFirstPoint()>track->GetRow())
243 track->SetFirstPoint(track->GetRow());
244 if (track->GetLastPoint()<track->GetRow())
245 track->SetLastPoint(track->GetRow());
91162307 246
247
b9671574 248 track->SetClusterPointer(track->GetRow(),c);
1c53abe2 249 //
250
3f82c4f2 251 Double_t angle2 = track->GetSnp()*track->GetSnp();
1c53abe2 252 //
253 //SET NEW Track Point
254 //
85e2b57d 255 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
91162307 256 {
85e2b57d 257 angle2 = TMath::Sqrt(angle2/(1-angle2));
b9671574 258 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
1c53abe2 259 //
b9671574 260 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
261 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
262 point.SetErrY(sqrt(track->GetErrorY2()));
263 point.SetErrZ(sqrt(track->GetErrorZ2()));
1c53abe2 264 //
91162307 265 point.SetX(track->GetX());
266 point.SetY(track->GetY());
267 point.SetZ(track->GetZ());
268 point.SetAngleY(angle2);
269 point.SetAngleZ(track->GetTgl());
b9671574 270 if (point.IsShared()){
271 track->SetErrorY2(track->GetErrorY2()*4);
272 track->SetErrorZ2(track->GetErrorZ2()*4);
91162307 273 }
274 }
275
b9671574 276 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
91162307 277 //
1d91d749 278// track->SetErrorY2(track->GetErrorY2()*1.3);
279// track->SetErrorY2(track->GetErrorY2()+0.01);
280// track->SetErrorZ2(track->GetErrorZ2()*1.3);
281// track->SetErrorZ2(track->GetErrorZ2()+0.005);
91162307 282 //}
283 if (accept>0) return 0;
284 if (track->GetNumberOfClusters()%20==0){
285 // if (track->fHelixIn){
286 // TClonesArray & larr = *(track->fHelixIn);
287 // Int_t ihelix = larr.GetEntriesFast();
288 // new(larr[ihelix]) AliHelix(*track) ;
289 //}
1c53abe2 290 }
65e67c9c 291 if (AliTPCReconstructor::StreamLevel()&kStreamUpdateTrack) {
bda5ad6d 292 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
293 AliExternalTrackParam param(*track);
294 TTreeSRedirector &cstream = *fDebugStreamer;
65e67c9c 295 cstream<<"UpdateTrack"<<
bda5ad6d 296 "cl.="<<c<<
265b5e37 297 "event="<<event<<
bda5ad6d 298 "track.="<<&param<<
299 "\n";
300 }
b9671574 301 track->SetNoCluster(0);
91162307 302 return track->Update(c,chi2,i);
303}
304
305
306
829455ad 307Int_t AliTPCtracker::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
91162307 308{
1c53abe2 309 //
91162307 310 // decide according desired precision to accept given
311 // cluster for tracking
6fbe1e5c 312 Double_t yt=0,zt=0;
313 seed->GetProlongation(cluster->GetX(),yt,zt);
00055a22 314 Double_t sy2=ErrY2(seed,cluster);
315 Double_t sz2=ErrZ2(seed,cluster);
1c53abe2 316
91162307 317 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
318 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
6fbe1e5c 319 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
f47588e0 320 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
6fbe1e5c 321 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
322 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
323 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
324 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
91162307 325
326 Double_t rdistance2 = rdistancey2+rdistancez2;
327 //Int_t accept =0;
328
f5a565c3 329 if (AliTPCReconstructor::StreamLevel()>2 && ( (fIteration>0)|| (seed->GetNumberOfClusters()>20))) {
a9ff2f09 330 // if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
fd065ea2 331 Float_t rmsy2 = seed->GetCurrentSigmaY2();
332 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
e0e13b88 333 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
334 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
42aec1f1 335 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
336 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
6fbe1e5c 337 AliExternalTrackParam param(*seed);
eea6b724 338 static TVectorD gcl(3),gtr(3);
339 Float_t gclf[3];
340 param.GetXYZ(gcl.GetMatrixArray());
341 cluster->GetGlobalXYZ(gclf);
342 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
a9ff2f09 343 Int_t nclSeed=seed->GetNumberOfClusters();
b194b32c 344
65e67c9c 345 if (AliTPCReconstructor::StreamLevel()&kStreamErrParam) { // flag:stream in debug mode cluster and track extrapolation at given row together with error nad shape estimate
b9d88e03 346 Int_t eventNr = fEvent->GetEventNumberInFile();
347
fd065ea2 348 (*fDebugStreamer)<<"ErrParam"<<
498b9441 349 "iter="<<fIteration<<
b9d88e03 350 "eventNr="<<eventNr<<
e0e13b88 351 "Cl.="<<cluster<<
a9ff2f09 352 "nclSeed="<<nclSeed<<
e0e13b88 353 "T.="<<&param<<
6fbe1e5c 354 "dy="<<dy<<
355 "dz="<<dz<<
356 "yt="<<yt<<
357 "zt="<<zt<<
eea6b724 358 "gcl.="<<&gcl<<
359 "gtr.="<<&gtr<<
e0e13b88 360 "erry2="<<sy2<<
361 "errz2="<<sz2<<
362 "rmsy2="<<rmsy2<<
363 "rmsz2="<<rmsz2<<
364 "rmsy2p30="<<rmsy2p30<<
365 "rmsz2p30="<<rmsz2p30<<
42aec1f1 366 "rmsy2p30R="<<rmsy2p30R<<
367 "rmsz2p30R="<<rmsz2p30R<<
f47588e0 368 // normalize distance -
369 "rdisty="<<rdistancey2<<
370 "rdistz="<<rdistancez2<<
371 "rdist="<<rdistance2<< //
e0e13b88 372 "\n";
b194b32c 373 }
fd065ea2 374 }
6fbe1e5c 375 //return 0; // temporary
376 if (rdistance2>32) return 3;
91162307 377
378
00055a22 379 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
91162307 380 return 2; //suspisiouce - will be changed
381
00055a22 382 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
91162307 383 // strict cut on overlaped cluster
384 return 2; //suspisiouce - will be changed
1c53abe2 385
00055a22 386 if ( (rdistancey2>1. || rdistancez2>6.25 )
91162307 387 && cluster->GetType()<0){
b9671574 388 seed->SetNFoundable(seed->GetNFoundable()-1);
91162307 389 return 2;
1c53abe2 390 }
19dcf504 391
1598ba75 392 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
36b89541 393 if (fIteration==2){
394 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
395 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
396 return 2;
397 }
1598ba75 398 }
19dcf504 399 }
1598ba75 400
91162307 401 return 0;
402}
403
404
1c53abe2 405
1c53abe2 406
1af5da7e 407
1c53abe2 408//_____________________________________________________________________________
829455ad 409AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
e046d791 410AliTracker(),
411 fkNIS(par->GetNInnerSector()/2),
412 fInnerSec(0),
413 fkNOS(par->GetNOuterSector()/2),
414 fOuterSec(0),
415 fN(0),
416 fSectors(0),
417 fInput(0),
418 fOutput(0),
419 fSeedTree(0),
420 fTreeDebug(0),
421 fEvent(0),
72e25240 422 fEventHLT(0),
e046d791 423 fDebug(0),
424 fNewIO(0),
425 fNtracks(0),
426 fSeeds(0),
427 fIteration(0),
47af7ca4 428 fkParam(0),
8cc1870a 429 fDebugStreamer(0),
430 fUseHLTClusters(4),
431 fCrossTalkSignalArray(0),
8cc1870a 432 fSeedsPool(0),
ddfbc51a 433 fFreeSeedsID(500),
434 fNFreeSeeds(0),
435 fLastSeedID(-1)
1c53abe2 436{
437 //---------------------------------------------------------------------
438 // The main TPC tracker constructor
439 //---------------------------------------------------------------------
9350f379 440 fInnerSec=new AliTPCtrackerSector[fkNIS];
441 fOuterSec=new AliTPCtrackerSector[fkNOS];
91162307 442
1c53abe2 443 Int_t i;
444 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
445 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
446
47af7ca4 447 fkParam = par;
91162307 448 Int_t nrowlow = par->GetNRowLow();
449 Int_t nrowup = par->GetNRowUp();
450
451
77f88633 452 for (i=0;i<nrowlow;i++){
91162307 453 fXRow[i] = par->GetPadRowRadiiLow(i);
454 fPadLength[i]= par->GetPadPitchLength(0,i);
455 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
456 }
457
458
77f88633 459 for (i=0;i<nrowup;i++){
91162307 460 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
461 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
462 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
463 }
e046d791 464
b194b32c 465 if (AliTPCReconstructor::StreamLevel()>0) {
3d674021 466 fDebugStreamer = new TTreeSRedirector("TPCdebug.root","recreate");
b194b32c 467 }
f06a1ff6 468 //
ddfbc51a 469 fSeedsPool = new TClonesArray("AliTPCseed",1000);
265b5e37 470
471 // crosstalk array and matrix initialization
472 Int_t nROCs = 72;
473 Int_t nTimeBinsAll = par->GetMaxTBin();
474 Int_t nWireSegments = 11;
3c8c25a3 475 fCrossTalkSignalArray = new TObjArray(nROCs*4); //
265b5e37 476 fCrossTalkSignalArray->SetOwner(kTRUE);
3c8c25a3 477 for (Int_t isector=0; isector<4*nROCs; isector++){
39240815 478 TMatrixD * crossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll);
265b5e37 479 for (Int_t imatrix = 0; imatrix<11; imatrix++)
480 for (Int_t jmatrix = 0; jmatrix<nTimeBinsAll; jmatrix++){
39240815 481 (*crossTalkSignal)[imatrix][jmatrix]=0.;
265b5e37 482 }
39240815 483 fCrossTalkSignalArray->AddAt(crossTalkSignal,isector);
265b5e37 484 }
485
1c53abe2 486}
2fc0c115 487//________________________________________________________________________
829455ad 488AliTPCtracker::AliTPCtracker(const AliTPCtracker &t):
58251ea0 489 AliTracker(t),
e046d791 490 fkNIS(t.fkNIS),
491 fInnerSec(0),
492 fkNOS(t.fkNOS),
493 fOuterSec(0),
494 fN(0),
495 fSectors(0),
496 fInput(0),
497 fOutput(0),
498 fSeedTree(0),
499 fTreeDebug(0),
500 fEvent(0),
72e25240 501 fEventHLT(0),
e046d791 502 fDebug(0),
503 fNewIO(kFALSE),
504 fNtracks(0),
505 fSeeds(0),
506 fIteration(0),
47af7ca4 507 fkParam(0),
8cc1870a 508 fDebugStreamer(0),
509 fUseHLTClusters(4),
510 fCrossTalkSignalArray(0),
8cc1870a 511 fSeedsPool(0),
ddfbc51a 512 fFreeSeedsID(500),
513 fNFreeSeeds(0),
514 fLastSeedID(-1)
58251ea0 515{
2fc0c115 516 //------------------------------------
517 // dummy copy constructor
518 //------------------------------------------------------------------
e046d791 519 fOutput=t.fOutput;
17abbffb 520 for (Int_t irow=0; irow<200; irow++){
521 fXRow[irow]=0;
522 fYMax[irow]=0;
523 fPadLength[irow]=0;
524 }
525
2fc0c115 526}
829455ad 527AliTPCtracker & AliTPCtracker::operator=(const AliTPCtracker& /*r*/)
47af7ca4 528{
2fc0c115 529 //------------------------------
530 // dummy
531 //--------------------------------------------------------------
532 return *this;
533}
1c53abe2 534//_____________________________________________________________________________
829455ad 535AliTPCtracker::~AliTPCtracker() {
1c53abe2 536 //------------------------------------------------------------------
537 // TPC tracker destructor
538 //------------------------------------------------------------------
539 delete[] fInnerSec;
540 delete[] fOuterSec;
541 if (fSeeds) {
f06a1ff6 542 fSeeds->Clear();
1c53abe2 543 delete fSeeds;
544 }
8cc1870a 545 if (fCrossTalkSignalArray) delete fCrossTalkSignalArray;
81e97e0d 546 if (fDebugStreamer) delete fDebugStreamer;
ddfbc51a 547 if (fSeedsPool) delete fSeedsPool;
1c53abe2 548}
549
1c53abe2 550
829455ad 551void AliTPCtracker::FillESD(const TObjArray* arr)
91162307 552{
47966a6d 553 //
554 //
555 //fill esds using updated tracks
ada522d7 556
ec26e231 557 if (!fEvent) return;
ada522d7 558
559 AliESDtrack iotrack;
ec26e231 560
91162307 561 // write tracks to the event
562 // store index of the track
d26d9159 563 Int_t nseed=arr->GetEntriesFast();
51ad6848 564 //FindKinks(arr,fEvent);
91162307 565 for (Int_t i=0; i<nseed; i++) {
d26d9159 566 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 567 if (!pt) continue;
51ad6848 568 pt->UpdatePoints();
92f513f5 569 AddCovariance(pt);
65e67c9c 570 if (AliTPCReconstructor::StreamLevel()&kStreamFillESD) {
571 (*fDebugStreamer)<<"FillESD"<< // flag: stream track information in FillESD function (after track Iteration 0)
8cecaa87 572 "Tr0.="<<pt<<
573 "\n";
574 }
47af7ca4 575 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
eea478d3 576 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
47af7ca4 577 pt->PropagateTo(fkParam->GetInnerRadiusLow());
eea478d3 578 }
51ad6848 579
580 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
ada522d7 581 iotrack.~AliESDtrack();
582 new(&iotrack) AliESDtrack;
51ad6848 583 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
a2e215dd 584 iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
51ad6848 585 iotrack.SetTPCPoints(pt->GetPoints());
586 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 587 iotrack.SetV0Indexes(pt->GetV0Indexes());
588 // iotrack.SetTPCpid(pt->fTPCr);
f5cbf2ef 589 //iotrack.SetTPCindex(i);
590 MakeESDBitmaps(pt, &iotrack);
51ad6848 591 fEvent->AddTrack(&iotrack);
592 continue;
593 }
594
b9671574 595 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
ada522d7 596 iotrack.~AliESDtrack();
597 new(&iotrack) AliESDtrack;
51ad6848 598 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
a2e215dd 599 iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
51ad6848 600 iotrack.SetTPCPoints(pt->GetPoints());
d26d9159 601 //iotrack.SetTPCindex(i);
51ad6848 602 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 603 iotrack.SetV0Indexes(pt->GetV0Indexes());
f5cbf2ef 604 MakeESDBitmaps(pt, &iotrack);
81e97e0d 605 // iotrack.SetTPCpid(pt->fTPCr);
d26d9159 606 fEvent->AddTrack(&iotrack);
a42a6bae 607 continue;
608 }
51ad6848 609 //
610 // short tracks - maybe decays
611
b9671574 612 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
a42a6bae 613 Int_t found,foundable,shared;
614 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
b9671574 615 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
ada522d7 616 iotrack.~AliESDtrack();
617 new(&iotrack) AliESDtrack;
a42a6bae 618 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
a2e215dd 619 iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
a42a6bae 620 //iotrack.SetTPCindex(i);
51ad6848 621 iotrack.SetTPCPoints(pt->GetPoints());
622 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 623 iotrack.SetV0Indexes(pt->GetV0Indexes());
f5cbf2ef 624 MakeESDBitmaps(pt, &iotrack);
81e97e0d 625 //iotrack.SetTPCpid(pt->fTPCr);
a42a6bae 626 fEvent->AddTrack(&iotrack);
627 continue;
628 }
629 }
630
b9671574 631 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
a42a6bae 632 Int_t found,foundable,shared;
633 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
634 if (found<20) continue;
b9671574 635 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
a42a6bae 636 //
ada522d7 637 iotrack.~AliESDtrack();
638 new(&iotrack) AliESDtrack;
a42a6bae 639 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
a2e215dd 640 iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
51ad6848 641 iotrack.SetTPCPoints(pt->GetPoints());
642 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 643 iotrack.SetV0Indexes(pt->GetV0Indexes());
f5cbf2ef 644 MakeESDBitmaps(pt, &iotrack);
81e97e0d 645 //iotrack.SetTPCpid(pt->fTPCr);
51ad6848 646 //iotrack.SetTPCindex(i);
647 fEvent->AddTrack(&iotrack);
648 continue;
649 }
650 // short tracks - secondaties
651 //
652 if ( (pt->GetNumberOfClusters()>30) ) {
653 Int_t found,foundable,shared;
654 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
b9671574 655 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
ada522d7 656 iotrack.~AliESDtrack();
657 new(&iotrack) AliESDtrack;
51ad6848 658 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
a2e215dd 659 iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
51ad6848 660 iotrack.SetTPCPoints(pt->GetPoints());
661 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 662 iotrack.SetV0Indexes(pt->GetV0Indexes());
f5cbf2ef 663 MakeESDBitmaps(pt, &iotrack);
81e97e0d 664 //iotrack.SetTPCpid(pt->fTPCr);
51ad6848 665 //iotrack.SetTPCindex(i);
666 fEvent->AddTrack(&iotrack);
667 continue;
668 }
669 }
670
671 if ( (pt->GetNumberOfClusters()>15)) {
672 Int_t found,foundable,shared;
673 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
674 if (found<15) continue;
e7eb17e4 675 if (foundable<=0) continue;
b9671574 676 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
51ad6848 677 if (float(found)/float(foundable)<0.8) continue;
678 //
ada522d7 679 iotrack.~AliESDtrack();
680 new(&iotrack) AliESDtrack;
51ad6848 681 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
a2e215dd 682 iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
51ad6848 683 iotrack.SetTPCPoints(pt->GetPoints());
684 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
81e97e0d 685 iotrack.SetV0Indexes(pt->GetV0Indexes());
f5cbf2ef 686 MakeESDBitmaps(pt, &iotrack);
81e97e0d 687 // iotrack.SetTPCpid(pt->fTPCr);
a42a6bae 688 //iotrack.SetTPCindex(i);
689 fEvent->AddTrack(&iotrack);
690 continue;
691 }
91162307 692 }
6fbe1e5c 693 // >> account for suppressed tracks in the kink indices (RS)
694 int nESDtracks = fEvent->GetNumberOfTracks();
695 for (int it=nESDtracks;it--;) {
696 AliESDtrack* esdTr = fEvent->GetTrack(it);
697 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
698 for (int ik=0;ik<3;ik++) {
699 int knkId=0;
700 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
701 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
702 if (!kink) {
703 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
704 continue;
705 }
706 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
707 }
708 }
ada522d7 709
ec26e231 710 // << account for suppressed tracks in the kink indices (RS)
711 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
712
d26d9159 713}
714
d26d9159 715
1c53abe2 716
1c53abe2 717
718
829455ad 719Double_t AliTPCtracker::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
91162307 720 //
721 //
fd065ea2 722 // Use calibrated cluster error from OCDB
91162307 723 //
fd065ea2 724 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
91162307 725 //
47af7ca4 726 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
91162307 727 Int_t ctype = cl->GetType();
fd065ea2 728 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
729 Double_t angle = seed->GetSnp()*seed->GetSnp();
e0e13b88 730 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
fd065ea2 731 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
732 if (ctype<0) {
733 erry2+=0.5; // edge cluster
734 }
735 erry2*=erry2;
0b2c141b 736 Double_t addErr=0;
737 const Double_t *errInner = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorClusterInner();
f858edae 738 addErr=errInner[0]*TMath::Exp(-TMath::Abs((cl->GetX()-85.)/errInner[1]));
0b2c141b 739 erry2+=addErr*addErr;
44ca7282 740 const Double_t *errCluster = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorCluster();
741 erry2+=errCluster[0]*errCluster[0];
fd065ea2 742 seed->SetErrorY2(erry2);
743 //
744 return erry2;
745
746//calculate look-up table at the beginning
747// static Bool_t ginit = kFALSE;
748// static Float_t gnoise1,gnoise2,gnoise3;
749// static Float_t ggg1[10000];
750// static Float_t ggg2[10000];
751// static Float_t ggg3[10000];
752// static Float_t glandau1[10000];
753// static Float_t glandau2[10000];
754// static Float_t glandau3[10000];
755// //
756// static Float_t gcor01[500];
757// static Float_t gcor02[500];
758// static Float_t gcorp[500];
759// //
91162307 760
fd065ea2 761// //
762// if (ginit==kFALSE){
763// for (Int_t i=1;i<500;i++){
764// Float_t rsigma = float(i)/100.;
765// gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
766// gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
767// gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
768// }
769
770// //
771// for (Int_t i=3;i<10000;i++){
772// //
773// //
774// // inner sector
775// Float_t amp = float(i);
776// Float_t padlength =0.75;
777// gnoise1 = 0.0004/padlength;
778// Float_t nel = 0.268*amp;
779// Float_t nprim = 0.155*amp;
47af7ca4 780// ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
fd065ea2 781// glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
782// if (glandau1[i]>1) glandau1[i]=1;
783// glandau1[i]*=padlength*padlength/12.;
784// //
785// // outer short
786// padlength =1.;
787// gnoise2 = 0.0004/padlength;
788// nel = 0.3*amp;
789// nprim = 0.133*amp;
47af7ca4 790// ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
fd065ea2 791// glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
792// if (glandau2[i]>1) glandau2[i]=1;
793// glandau2[i]*=padlength*padlength/12.;
794// //
795// //
796// // outer long
797// padlength =1.5;
798// gnoise3 = 0.0004/padlength;
799// nel = 0.3*amp;
800// nprim = 0.133*amp;
47af7ca4 801// ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
fd065ea2 802// glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
803// if (glandau3[i]>1) glandau3[i]=1;
804// glandau3[i]*=padlength*padlength/12.;
805// //
806// }
807// ginit = kTRUE;
808// }
809// //
810// //
811// //
812// Int_t amp = int(TMath::Abs(cl->GetQ()));
813// if (amp>9999) {
814// seed->SetErrorY2(1.);
815// return 1.;
816// }
817// Float_t snoise2;
47af7ca4 818// Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
fd065ea2 819// Int_t ctype = cl->GetType();
820// Float_t padlength= GetPadPitchLength(seed->GetRow());
821// Double_t angle2 = seed->GetSnp()*seed->GetSnp();
822// angle2 = angle2/(1-angle2);
823// //
824// //cluster "quality"
825// Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
826// Float_t res;
827// //
828// if (fSectors==fInnerSec){
829// snoise2 = gnoise1;
830// res = ggg1[amp]*z+glandau1[amp]*angle2;
831// if (ctype==0) res *= gcor01[rsigmay];
832// if ((ctype>0)){
833// res+=0.002;
834// res*= gcorp[rsigmay];
835// }
836// }
837// else {
838// if (padlength<1.1){
839// snoise2 = gnoise2;
840// res = ggg2[amp]*z+glandau2[amp]*angle2;
841// if (ctype==0) res *= gcor02[rsigmay];
842// if ((ctype>0)){
843// res+=0.002;
844// res*= gcorp[rsigmay];
845// }
846// }
847// else{
848// snoise2 = gnoise3;
849// res = ggg3[amp]*z+glandau3[amp]*angle2;
850// if (ctype==0) res *= gcor02[rsigmay];
851// if ((ctype>0)){
852// res+=0.002;
853// res*= gcorp[rsigmay];
854// }
855// }
856// }
857
858// if (ctype<0){
859// res+=0.005;
860// res*=2.4; // overestimate error 2 times
861// }
862// res+= snoise2;
1c53abe2 863
fd065ea2 864// if (res<2*snoise2)
865// res = 2*snoise2;
91162307 866
fd065ea2 867// seed->SetErrorY2(res);
868// return res;
1c53abe2 869
870
91162307 871}
872
1c53abe2 873
874
829455ad 875Double_t AliTPCtracker::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
1c53abe2 876 //
877 //
fd065ea2 878 // Use calibrated cluster error from OCDB
91162307 879 //
fd065ea2 880 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
91162307 881 //
47af7ca4 882 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
91162307 883 Int_t ctype = cl->GetType();
fd065ea2 884 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
91162307 885 //
3f82c4f2 886 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
91162307 887 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
e0e13b88 888 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
fd065ea2 889 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
890 if (ctype<0) {
891 errz2+=0.5; // edge cluster
91162307 892 }
1af5da7e 893 errz2*=errz2;
0b2c141b 894 Double_t addErr=0;
895 const Double_t *errInner = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorClusterInner();
f858edae 896 addErr=errInner[0]*TMath::Exp(-TMath::Abs((cl->GetX()-85.)/errInner[1]));
0b2c141b 897 errz2+=addErr*addErr;
44ca7282 898 const Double_t *errCluster = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorCluster();
899 errz2+=errCluster[1]*errCluster[1];
fd065ea2 900 seed->SetErrorZ2(errz2);
901 //
902 return errz2;
1c53abe2 903
fd065ea2 904
905
906// //seed->SetErrorY2(0.1);
907// //return 0.1;
908// //calculate look-up table at the beginning
909// static Bool_t ginit = kFALSE;
910// static Float_t gnoise1,gnoise2,gnoise3;
911// static Float_t ggg1[10000];
912// static Float_t ggg2[10000];
913// static Float_t ggg3[10000];
914// static Float_t glandau1[10000];
915// static Float_t glandau2[10000];
916// static Float_t glandau3[10000];
917// //
918// static Float_t gcor01[1000];
919// static Float_t gcor02[1000];
920// static Float_t gcorp[1000];
921// //
922
923// //
924// if (ginit==kFALSE){
925// for (Int_t i=1;i<1000;i++){
926// Float_t rsigma = float(i)/100.;
927// gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
928// gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
929// gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
930// }
931
932// //
933// for (Int_t i=3;i<10000;i++){
934// //
935// //
936// // inner sector
937// Float_t amp = float(i);
938// Float_t padlength =0.75;
939// gnoise1 = 0.0004/padlength;
940// Float_t nel = 0.268*amp;
941// Float_t nprim = 0.155*amp;
47af7ca4 942// ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
fd065ea2 943// glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
944// if (glandau1[i]>1) glandau1[i]=1;
945// glandau1[i]*=padlength*padlength/12.;
946// //
947// // outer short
948// padlength =1.;
949// gnoise2 = 0.0004/padlength;
950// nel = 0.3*amp;
951// nprim = 0.133*amp;
47af7ca4 952// ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
fd065ea2 953// glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
954// if (glandau2[i]>1) glandau2[i]=1;
955// glandau2[i]*=padlength*padlength/12.;
956// //
957// //
958// // outer long
959// padlength =1.5;
960// gnoise3 = 0.0004/padlength;
961// nel = 0.3*amp;
962// nprim = 0.133*amp;
47af7ca4 963// ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
fd065ea2 964// glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
965// if (glandau3[i]>1) glandau3[i]=1;
966// glandau3[i]*=padlength*padlength/12.;
967// //
968// }
969// ginit = kTRUE;
970// }
971// //
972// //
973// //
974// Int_t amp = int(TMath::Abs(cl->GetQ()));
975// if (amp>9999) {
976// seed->SetErrorY2(1.);
977// return 1.;
978// }
979// Float_t snoise2;
47af7ca4 980// Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
fd065ea2 981// Int_t ctype = cl->GetType();
982// Float_t padlength= GetPadPitchLength(seed->GetRow());
983// //
984// Double_t angle2 = seed->GetSnp()*seed->GetSnp();
985// // if (angle2<0.6) angle2 = 0.6;
986// angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
987// //
988// //cluster "quality"
989// Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
990// Float_t res;
991// //
992// if (fSectors==fInnerSec){
993// snoise2 = gnoise1;
994// res = ggg1[amp]*z+glandau1[amp]*angle2;
995// if (ctype==0) res *= gcor01[rsigmaz];
996// if ((ctype>0)){
997// res+=0.002;
998// res*= gcorp[rsigmaz];
999// }
1000// }
1001// else {
1002// if (padlength<1.1){
1003// snoise2 = gnoise2;
1004// res = ggg2[amp]*z+glandau2[amp]*angle2;
1005// if (ctype==0) res *= gcor02[rsigmaz];
1006// if ((ctype>0)){
1007// res+=0.002;
1008// res*= gcorp[rsigmaz];
1009// }
1010// }
1011// else{
1012// snoise2 = gnoise3;
1013// res = ggg3[amp]*z+glandau3[amp]*angle2;
1014// if (ctype==0) res *= gcor02[rsigmaz];
1015// if ((ctype>0)){
1016// res+=0.002;
1017// res*= gcorp[rsigmaz];
1018// }
1019// }
1020// }
1021
1022// if (ctype<0){
1023// res+=0.002;
1024// res*=1.3;
1025// }
1026// if ((ctype<0) &&amp<70){
1027// res+=0.002;
1028// res*=1.3;
1029// }
1030// res += snoise2;
1031// if (res<2*snoise2)
1032// res = 2*snoise2;
1033// if (res>3) res =3;
1034// seed->SetErrorZ2(res);
1035// return res;
91162307 1036}
1037
1038
1039
c9427e08 1040
1c53abe2 1041
829455ad 1042void AliTPCtracker::RotateToLocal(AliTPCseed *seed)
c9427e08 1043{
1044 //rotate to track "local coordinata
1045 Float_t x = seed->GetX();
1046 Float_t y = seed->GetY();
1047 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
91162307 1048
c9427e08 1049 if (y > ymax) {
b9671574 1050 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
c9427e08 1051 if (!seed->Rotate(fSectors->GetAlpha()))
1052 return;
1053 } else if (y <-ymax) {
b9671574 1054 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
c9427e08 1055 if (!seed->Rotate(-fSectors->GetAlpha()))
1056 return;
1057 }
1c53abe2 1058
c9427e08 1059}
1c53abe2 1060
1061
1062
1c53abe2 1063//_____________________________________________________________________________
829455ad 1064Double_t AliTPCtracker::F1old(Double_t x1,Double_t y1,
1c53abe2 1065 Double_t x2,Double_t y2,
47af7ca4 1066 Double_t x3,Double_t y3) const
1c53abe2 1067{
1068 //-----------------------------------------------------------------
1069 // Initial approximation of the track curvature
1070 //-----------------------------------------------------------------
1071 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1072 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1073 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1074 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1075 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1076
1077 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
91162307 1078 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1c53abe2 1079 return -xr*yr/sqrt(xr*xr+yr*yr);
1080}
1081
1082
91162307 1083
1c53abe2 1084//_____________________________________________________________________________
829455ad 1085Double_t AliTPCtracker::F1(Double_t x1,Double_t y1,
91162307 1086 Double_t x2,Double_t y2,
47af7ca4 1087 Double_t x3,Double_t y3) const
91162307 1088{
1089 //-----------------------------------------------------------------
1090 // Initial approximation of the track curvature
1091 //-----------------------------------------------------------------
1092 x3 -=x1;
1093 x2 -=x1;
1094 y3 -=y1;
1095 y2 -=y1;
1096 //
1097 Double_t det = x3*y2-x2*y3;
6e23caff 1098 if (TMath::Abs(det)<1e-10){
91162307 1099 return 100;
1100 }
1101 //
1102 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1103 Double_t x0 = x3*0.5-y3*u;
1104 Double_t y0 = y3*0.5+x3*u;
1105 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1106 if (det<0) c2*=-1;
1107 return c2;
1108}
1109
1110
829455ad 1111Double_t AliTPCtracker::F2(Double_t x1,Double_t y1,
1c53abe2 1112 Double_t x2,Double_t y2,
47af7ca4 1113 Double_t x3,Double_t y3) const
1c53abe2 1114{
1115 //-----------------------------------------------------------------
91162307 1116 // Initial approximation of the track curvature
1117 //-----------------------------------------------------------------
1118 x3 -=x1;
1119 x2 -=x1;
1120 y3 -=y1;
1121 y2 -=y1;
1122 //
1123 Double_t det = x3*y2-x2*y3;
6e23caff 1124 if (TMath::Abs(det)<1e-10) {
91162307 1125 return 100;
1126 }
1127 //
1128 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1129 Double_t x0 = x3*0.5-y3*u;
1130 Double_t y0 = y3*0.5+x3*u;
1131 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1132 if (det<0) c2*=-1;
1133 x0+=x1;
1134 x0*=c2;
1135 return x0;
1136}
1137
1138
1139
1140//_____________________________________________________________________________
829455ad 1141Double_t AliTPCtracker::F2old(Double_t x1,Double_t y1,
91162307 1142 Double_t x2,Double_t y2,
47af7ca4 1143 Double_t x3,Double_t y3) const
91162307 1144{
1145 //-----------------------------------------------------------------
1c53abe2 1146 // Initial approximation of the track curvature times center of curvature
1147 //-----------------------------------------------------------------
1148 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1149 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1150 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1151 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1152 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1153
1154 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1155
1156 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1157}
1158
1159//_____________________________________________________________________________
829455ad 1160Double_t AliTPCtracker::F3(Double_t x1,Double_t y1,
1c53abe2 1161 Double_t x2,Double_t y2,
47af7ca4 1162 Double_t z1,Double_t z2) const
1c53abe2 1163{
1164 //-----------------------------------------------------------------
1165 // Initial approximation of the tangent of the track dip angle
1166 //-----------------------------------------------------------------
1167 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1168}
1169
1170
829455ad 1171Double_t AliTPCtracker::F3n(Double_t x1,Double_t y1,
91162307 1172 Double_t x2,Double_t y2,
47af7ca4 1173 Double_t z1,Double_t z2, Double_t c) const
1c53abe2 1174{
91162307 1175 //-----------------------------------------------------------------
1176 // Initial approximation of the tangent of the track dip angle
1177 //-----------------------------------------------------------------
1178
1179 // Double_t angle1;
1180
1181 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1c53abe2 1182 //
91162307 1183 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1184 if (TMath::Abs(d*c*0.5)>1) return 0;
1185 // Double_t angle2 = TMath::ASin(d*c*0.5);
b67e07dc 1186 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1187 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
91162307 1188
1189 angle2 = (z1-z2)*c/(angle2*2.);
1190 return angle2;
1191}
1192
829455ad 1193Bool_t AliTPCtracker::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
91162307 1194{//-----------------------------------------------------------------
1195 // This function find proloncation of a track to a reference plane x=x2.
1196 //-----------------------------------------------------------------
1197
1198 Double_t dx=x2-x1;
1199
1200 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1201 return kFALSE;
1c53abe2 1202 }
f8aae377 1203
bfd20868 1204 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1205 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
91162307 1206 y = x[0];
1207 z = x[1];
1208
1209 Double_t dy = dx*(c1+c2)/(r1+r2);
1210 Double_t dz = 0;
1211 //
1212 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1213
1214 if (TMath::Abs(delta)>0.01){
1215 dz = x[3]*TMath::ASin(delta)/x[4];
1216 }else{
b67e07dc 1217 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
91162307 1218 }
1219
b67e07dc 1220 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
f8aae377 1221
91162307 1222 y+=dy;
1223 z+=dz;
1224
1225 return kTRUE;
1c53abe2 1226}
1227
829455ad 1228Int_t AliTPCtracker::LoadClusters (TTree *const tree)
d26d9159 1229{
544c295f 1230 // load clusters
d26d9159 1231 //
1232 fInput = tree;
1233 return LoadClusters();
1234}
91162307 1235
af86c1fd 1236
829455ad 1237Int_t AliTPCtracker::LoadClusters(const TObjArray *arr)
af86c1fd 1238{
1239 //
1240 // load clusters to the memory
e7de6656 1241 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
af86c1fd 1242 Int_t lower = arr->LowerBound();
1243 Int_t entries = arr->GetEntriesFast();
77f88633 1244
af86c1fd 1245 for (Int_t i=lower; i<entries; i++) {
1246 clrow = (AliTPCClustersRow*) arr->At(i);
77f88633 1247 if(!clrow) continue;
1248 if(!clrow->GetArray()) continue;
1249
af86c1fd 1250 //
1251 Int_t sec,row;
47af7ca4 1252 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
77f88633 1253
af86c1fd 1254 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1255 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1256 }
1257 //
1258 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
bd26fa83 1259 AliTPCtrackerRow * tpcrow=0;
af86c1fd 1260 Int_t left=0;
1261 if (sec<fkNIS*2){
1262 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1263 left = sec/fkNIS;
1264 }
1265 else{
1266 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1267 left = (sec-fkNIS*2)/fkNOS;
1268 }
1269 if (left ==0){
1270 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
77f88633 1271 for (Int_t j=0;j<tpcrow->GetN1();++j)
1272 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
af86c1fd 1273 }
1274 if (left ==1){
1275 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
77f88633 1276 for (Int_t j=0;j<tpcrow->GetN2();++j)
bfa00fba 1277 tpcrow->SetCluster2(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
af86c1fd 1278 }
e7de6656 1279 clrow->GetArray()->Clear("C");
af86c1fd 1280 }
1281 //
1282 delete clrow;
1283 LoadOuterSectors();
1284 LoadInnerSectors();
1285 return 0;
1286}
1287
829455ad 1288Int_t AliTPCtracker::LoadClusters(const TClonesArray *arr)
aa7f1a5a 1289{
1290 //
1291 // load clusters to the memory from one
1292 // TClonesArray
1293 //
1294 AliTPCclusterMI *clust=0;
98ee6d31 1295 Int_t count[72][96] = { {0} , {0} };
aa7f1a5a 1296
1297 // loop over clusters
1298 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1299 clust = (AliTPCclusterMI*)arr->At(icl);
1300 if(!clust) continue;
1301 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1302
1303 // transform clusters
1304 Transform(clust);
1305
1306 // count clusters per pad row
1307 count[clust->GetDetector()][clust->GetRow()]++;
1308 }
1309
1310 // insert clusters to sectors
1311 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1312 clust = (AliTPCclusterMI*)arr->At(icl);
1313 if(!clust) continue;
1314
1315 Int_t sec = clust->GetDetector();
1316 Int_t row = clust->GetRow();
98ee6d31 1317
1318 // filter overlapping pad rows needed by HLT
1319 if(sec<fkNIS*2) { //IROCs
1320 if(row == 30) continue;
1321 }
1322 else { // OROCs
1323 if(row == 27 || row == 76) continue;
1324 }
1325
2a97785a 1326 // Int_t left=0;
aa7f1a5a 1327 if (sec<fkNIS*2){
2a97785a 1328 // left = sec/fkNIS;
47af7ca4 1329 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
aa7f1a5a 1330 }
1331 else{
2a97785a 1332 // left = (sec-fkNIS*2)/fkNOS;
47af7ca4 1333 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
aa7f1a5a 1334 }
1335 }
1336
98ee6d31 1337 // Load functions must be called behind LoadCluster(TClonesArray*)
1338 // needed by HLT
1339 //LoadOuterSectors();
1340 //LoadInnerSectors();
aa7f1a5a 1341
1342 return 0;
1343}
1344
af86c1fd 1345
829455ad 1346Int_t AliTPCtracker::LoadClusters()
1c53abe2 1347{
3c8c25a3 1348 //
1349 // load clusters to the memory
bfa00fba 1350 static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
91162307 1351 //
1352 // TTree * tree = fClustersArray.GetTree();
9a836cc2 1353 AliInfo("LoadClusters()\n");
91162307 1354
1355 TTree * tree = fInput;
1356 TBranch * br = tree->GetBranch("Segment");
1357 br->SetAddress(&clrow);
bad6eb00 1358
1359 // Conversion of pad, row coordinates in local tracking coords.
1360 // Could be skipped here; is already done in clusterfinder
1361
91162307 1362 Int_t j=Int_t(tree->GetEntries());
1c53abe2 1363 for (Int_t i=0; i<j; i++) {
91162307 1364 br->GetEntry(i);
1365 //
1366 Int_t sec,row;
47af7ca4 1367 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
0201b65c 1368 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
3c8c25a3 1369 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
0201b65c 1370 }
91162307 1371 //
bd26fa83 1372 AliTPCtrackerRow * tpcrow=0;
91162307 1373 Int_t left=0;
1374 if (sec<fkNIS*2){
1375 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1376 left = sec/fkNIS;
1377 }
1378 else{
1379 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1380 left = (sec-fkNIS*2)/fkNOS;
1381 }
1382 if (left ==0){
b9671574 1383 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
77f88633 1384 for (Int_t k=0;k<tpcrow->GetN1();++k)
1385 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
91162307 1386 }
1387 if (left ==1){
b9671574 1388 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
77f88633 1389 for (Int_t k=0;k<tpcrow->GetN2();++k)
bfa00fba 1390 tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
91162307 1391 }
1c53abe2 1392 }
3c8c25a3 1393 //
1394 clrow->Clear("C");
1395 LoadOuterSectors();
1396 LoadInnerSectors();
1397
1398 cout << " =================================================================================================================================== " << endl;
1399 cout << " AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() << endl;
1400 cout << " AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() << endl;
1401 cout << " =================================================================================================================================== " << endl;
1402
1403 if (AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection()) ApplyTailCancellation();
1404 if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) CalculateXtalkCorrection();
1405 if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) ApplyXtalkCorrection();
1406 //if (AliTPCReconstructor::GetRecoParam()->GetUseOulierClusterFilter()) FilterOutlierClusters();
1407 return 0;
1408}
1409
1410
1411void AliTPCtracker::CalculateXtalkCorrection(){
1412 //
6712c26f 1413 // Calculate crosstalk estimate
3c8c25a3 1414 //
1415 const Int_t nROCs = 72;
f5d2d8e0 1416 const Int_t nIterations=3; //
3c8c25a3 1417 // 0.) reset crosstalk matrix
1418 //
1419 for (Int_t isector=0; isector<nROCs*4; isector++){ //set all ellemts of crosstalk matrix to 0
1420 TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1421 if (crossTalkMatrix)(*crossTalkMatrix)*=0;
1422 }
1423
1424 //
1425 // 1.) Filling part -- loop over clusters
1426 //
1427 Double_t missingChargeFactor= AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrectionMissingCharge();
abf8ed8a 1428 for (Int_t iter=0; iter<nIterations;iter++){
3c8c25a3 1429 for (Int_t isector=0; isector<36; isector++){ // loop over sectors
1430 for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
1431 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1432 Int_t nrows = sector.GetNRows();
1433 Int_t sec=0;
1434 if (isector<18) sec=isector+18*iside;
abf8ed8a 1435 if (isector>=18) sec=18+isector+18*iside;
3c8c25a3 1436 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1437 //
1438 //
1439 Int_t wireSegmentID = fkParam->GetWireSegment(sec,row);
1440 Float_t nPadsPerSegment = (Float_t)(fkParam->GetNPadsPerSegment(wireSegmentID));
1441 TMatrixD &crossTalkSignal = *((TMatrixD*)fCrossTalkSignalArray->At(sec));
1442 TMatrixD &crossTalkSignalCache = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs*2)); // this is the cache value of the crosstalk from previous iteration
1443 TMatrixD &crossTalkSignalBelow = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs));
1444 Int_t nCols=crossTalkSignal.GetNcols();
1445 //
1446 AliTPCtrackerRow& tpcrow = sector[row];
1447 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1448 if (iside>0) ncl=tpcrow.GetN2();
1449 for (Int_t i=0;i<ncl;i++) { // loop over clusters
1450 AliTPCclusterMI *clXtalk= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1451
1452 Int_t timeBinXtalk = clXtalk->GetTimeBin();
1453 Double_t rmsPadMin2=0.5*0.5+(fkParam->GetDiffT()*fkParam->GetDiffT())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec)); // minimal PRF width - 0.5 is the PRF in the pad units - we should et it from fkparam getters
1454 Double_t rmsTimeMin2=1+(fkParam->GetDiffL()*fkParam->GetDiffL())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetZWidth()*fkParam->GetZWidth()); // minimal PRF width - 1 is the TRF in the time bin units - we should et it from fkParam getters
1455 Double_t rmsTime2 = clXtalk->GetSigmaZ2()/(fkParam->GetZWidth()*fkParam->GetZWidth());
1456 Double_t rmsPad2 = clXtalk->GetSigmaY2()/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec));
1457 if (rmsPadMin2>rmsPad2){
1458 rmsPad2=rmsPadMin2;
1459 }
1460 if (rmsTimeMin2>rmsTime2){
1461 rmsTime2=rmsTimeMin2;
1462 }
1463
1464 Double_t norm= 2.*TMath::Exp(-1.0/(2.*rmsTime2))+2.*TMath::Exp(-4.0/(2.*rmsTime2))+1.;
1465 Double_t qTotXtalk = 0.;
1466 Double_t qTotXtalkMissing = 0.;
1467 for (Int_t itb=timeBinXtalk-2, idelta=-2; itb<=timeBinXtalk+2; itb++,idelta++) {
1468 if (itb<0 || itb>=nCols) continue;
1469 Double_t missingCharge=0;
1470 Double_t trf= TMath::Exp(-idelta*idelta/(2.*rmsTime2));
1471 if (missingChargeFactor>0) {
1472 for (Int_t dpad=-2; dpad<=2; dpad++){
1473 Double_t qPad = clXtalk->GetMax()*TMath::Exp(-dpad*dpad/(2.*rmsPad2))*trf;
1474 if (TMath::Nint(qPad-crossTalkSignalCache[wireSegmentID][itb])<=fkParam->GetZeroSup()){
1475 missingCharge+=qPad+crossTalkSignalCache[wireSegmentID][itb];
1476 }else{
1477 missingCharge+=crossTalkSignalCache[wireSegmentID][itb];
1478 }
1479 }
1480 }
1481 qTotXtalk = clXtalk->GetQ()*trf/norm+missingCharge*missingChargeFactor;
1482 qTotXtalkMissing = missingCharge;
1483 crossTalkSignal[wireSegmentID][itb]+= qTotXtalk/nPadsPerSegment;
1484 crossTalkSignalBelow[wireSegmentID][itb]+= qTotXtalkMissing/nPadsPerSegment;
1485 } // end of time bin loop
1486 } // end of cluster loop
1487 } // end of rows loop
1488 } // end of side loop
1489 } // end of sector loop
1490 //
1491 // copy crosstalk matrix to cached used for next itteration
434cdbe1 1492 //
6712c26f 1493 //
1494 // 2.) dump the crosstalk matrices to tree for further investigation
1495 // a.) to estimate fluctuation of pedestal in indiviula wire segments
1496 // b.) to check correlation between regions
1497 // c.) to check relative conribution of signal below threshold to crosstalk
1498
1499 if (AliTPCReconstructor::StreamLevel()&kStreamCrosstalkMatrix) {
1500 for (Int_t isector=0; isector<nROCs; isector++){ //set all ellemts of crosstalk matrix to 0
1501 TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1502 TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
1503 TMatrixD * crossTalkMatrixCache = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs*2);
1504 TVectorD vecAll(crossTalkMatrix->GetNrows());
1505 TVectorD vecBelow(crossTalkMatrix->GetNrows());
1506 TVectorD vecCache(crossTalkMatrixCache->GetNrows());
1507 //
1508 for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
1509 for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
1510 vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
1511 vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
1512 vecCache[iwire]=(*crossTalkMatrixCache)(iwire,itime);
1513 }
1514 (*fDebugStreamer)<<"crosstalkMatrix"<<
1515 "iter="<<iter<< //iteration
1516 "sector="<<isector<< // sector
1517 "itime="<<itime<< // time bin index
1518 "vecAll.="<<&vecAll<< // crosstalk charge + charge below threshold
1519 "vecCache.="<<&vecCache<< // crosstalk charge + charge below threshold
1520 "vecBelow.="<<&vecBelow<< // crosstalk contribution from signal below threshold
1521 "\n";
434cdbe1 1522 }
434cdbe1 1523 }
1524 }
f5d2d8e0 1525 if (iter<nIterations-1) for (Int_t isector=0; isector<nROCs*2; isector++){ //set all ellemts of crosstalk matrix to 0
1526 TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1527 TMatrixD * crossTalkMatrixCache = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs*2);
1528 if (crossTalkMatrix){
1529 (*crossTalkMatrixCache)*=0;
1530 (*crossTalkMatrixCache)+=(*crossTalkMatrix);
1531 (*crossTalkMatrix)*=0;
1532 }
1533 }
434cdbe1 1534 }
1535
1536
1c53abe2 1537}
1538
3c8c25a3 1539
1540
1541
eea44233 1542void AliTPCtracker::FilterOutlierClusters(){
1543 //
1544 // filter outlier clusters
1545 //
1546 /*
1547 1.)..... booking part
1548 nSectors=72;
1549 nTimeBins=fParam->Get....
1550 TH2F hisTime("","", sector,0,sector, nTimeBins,0,nTimeBins);
1551 TH2F hisPadRow("","", sector,0,sector, nPadRows,0,nPadRows);
1552 2.) .... filling part
1553 .... cluster loop { hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin()); }
1554
1555 3.) ...filtering part
1556 sector loop { calculate median,mean80 and rms80 of the nclusters per time bin; calculate median,mean80 and rms80 of the nclusters per par row; .... export values to the debug streamers - to decide which threshold to be used... }
1557
1558 sector loop
1559 { disable clusters in time bins > mean+ n rms80+ offsetTime disable clusters in padRow > mean+ n rms80+ offsetPadRow // how to dislable clusters? - new bit to introduce }
1560 //
1561 4. Disabling clusters
1562
1563 */
1564
1565 //
1566 // 1.) booking part
e731a3f8 1567 //
3e1f1ce7 1568 // AliTPCcalibDB *db=AliTPCcalibDB::Instance();
e731a3f8 1569 Int_t nSectors=AliTPCROC::Instance()->GetNSectors();
3e1f1ce7 1570 Int_t nTimeBins= 1100; // *Bug here - we should get NTimeBins from ALTRO - Parameters not relyable
1571 Int_t nPadRows=(AliTPCROC::Instance()->GetNRows(0) + AliTPCROC::Instance()->GetNRows(36));
e731a3f8 1572 // parameters for filtering
1573 const Double_t nSigmaCut=9.; // should be in recoParam ?
1574 const Double_t offsetTime=100; // should be in RecoParam ? -
1575 const Double_t offsetPadRow=300; // should be in RecoParam ?
1576 const Double_t offsetTimeAccept=8; // should be in RecoParam ? - obtained as mean +1 rms in high IR pp
eea44233 1577 TH2F hisTime("hisSectorTime","hisSectorTime", nSectors,0,nSectors, nTimeBins,0,nTimeBins);
1578 TH2F hisPadRow("hisSectorRow","hisSectorRow", nSectors,0,nSectors, nPadRows,0,nPadRows);
1579 //
1580 // 2.) Filling part -- loop over clusters
1581 //
44ca7282 1582 for (Int_t isector=0; isector<36; isector++){ // loop over sectors
1583 for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
1584 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1585 Int_t nrows = sector.GetNRows();
1586 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1587 AliTPCtrackerRow& tpcrow = sector[row];
1588 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1589 if (iside>0) ncl=tpcrow.GetN2();
1590 for (Int_t i=0;i<ncl;i++) { // loop over clusters
1591 AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1592 hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin());
1593 hisPadRow.Fill(cluster->GetDetector(),cluster->GetRow());
1594 }
1595 }
1596 }
1597 }
44ca7282 1598
eea44233 1599 //
1600 // 3. Filtering part
1601 //
1602 TVectorD vecTime(nTimeBins);
44ca7282 1603 TVectorD vecPadRow(nPadRows);
c5dce091 1604 TVectorD vecMedianSectorTime(nSectors);
44ca7282 1605 TVectorD vecRMSSectorTime(nSectors);
c5dce091 1606 TVectorD vecMedianSectorTimeOut6(nSectors);
1607 TVectorD vecMedianSectorTimeOut9(nSectors);//
3e1f1ce7 1608 TVectorD vecMedianSectorTimeOut(nSectors);//
c5dce091 1609 TVectorD vecMedianSectorPadRow(nSectors);
44ca7282 1610 TVectorD vecRMSSectorPadRow(nSectors);
c5dce091 1611 TVectorD vecMedianSectorPadRowOut6(nSectors);
1612 TVectorD vecMedianSectorPadRowOut9(nSectors);
3e1f1ce7 1613 TVectorD vecMedianSectorPadRowOut(nSectors);
e731a3f8 1614 TVectorD vecSectorOut6(nSectors);
1615 TVectorD vecSectorOut9(nSectors);
3e1f1ce7 1616 TMatrixD matSectorCluster(nSectors,2);
e731a3f8 1617 //
1618 // 3.a) median, rms calculations for hisTime
1619 //
eea44233 1620 for (Int_t isec=0; isec<nSectors; isec++){
c5dce091 1621 vecMedianSectorTimeOut6[isec]=0;
1622 vecMedianSectorTimeOut9[isec]=0;
eea44233 1623 for (Int_t itime=0; itime<nTimeBins; itime++){
1624 vecTime[itime]=hisTime.GetBinContent(isec+1, itime+1);
1625 }
3e1f1ce7 1626 Double_t median= TMath::Mean(nTimeBins,vecTime.GetMatrixArray());
c5dce091 1627 Double_t rms= TMath::RMS(nTimeBins,vecTime.GetMatrixArray());
1628 vecMedianSectorTime[isec]=median;
44ca7282 1629 vecRMSSectorTime[isec]=rms;
3e1f1ce7 1630 if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector TimeStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
eea44233 1631 //
1632 // declare outliers
1633 for (Int_t itime=0; itime<nTimeBins; itime++){
1634 Double_t entries= hisTime.GetBinContent(isec+1, itime+1);
44ca7282 1635 if (entries>median+6.*rms+offsetTime) {
c5dce091 1636 vecMedianSectorTimeOut6[isec]+=1;
1637 }
44ca7282 1638 if (entries>median+9.*rms+offsetTime) {
c5dce091 1639 vecMedianSectorTimeOut9[isec]+=1;
eea44233 1640 }
1641 }
c5dce091 1642 }
e731a3f8 1643 //
1644 // 3.b) median, rms calculations for hisPadRow
1645 //
44ca7282 1646 for (Int_t isec=0; isec<nSectors; isec++){
1647 vecMedianSectorPadRowOut6[isec]=0;
1648 vecMedianSectorPadRowOut9[isec]=0;
1649 for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1650 vecPadRow[ipadrow]=hisPadRow.GetBinContent(isec+1, ipadrow+1);
1651 }
3e1f1ce7 1652 Int_t nPadRowsSector= AliTPCROC::Instance()->GetNRows(isec);
1653 Double_t median= TMath::Mean(nPadRowsSector,vecPadRow.GetMatrixArray());
1654 Double_t rms= TMath::RMS(nPadRowsSector,vecPadRow.GetMatrixArray());
44ca7282 1655 vecMedianSectorPadRow[isec]=median;
1656 vecRMSSectorPadRow[isec]=rms;
3e1f1ce7 1657 if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector PadRowStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
44ca7282 1658 //
1659 // declare outliers
1660 for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1661 Double_t entries= hisPadRow.GetBinContent(isec+1, ipadrow+1);
1662 if (entries>median+6.*rms+offsetPadRow) {
1663 vecMedianSectorPadRowOut6[isec]+=1;
1664 }
1665 if (entries>median+9.*rms+offsetPadRow) {
1666 vecMedianSectorPadRowOut9[isec]+=1;
1667 }
1668 }
1669 }
1670 //
e731a3f8 1671 // 3.c) filter outlier sectors
1672 //
1673 Double_t medianSectorTime = TMath::Median(nSectors, vecTime.GetMatrixArray());
1674 Double_t mean69SectorTime, rms69SectorTime=0;
1675 AliMathBase::EvaluateUni(nSectors, vecTime.GetMatrixArray(), mean69SectorTime,rms69SectorTime,69);
1676 for (Int_t isec=0; isec<nSectors; isec++){
1677 vecSectorOut6[isec]=0;
1678 vecSectorOut9[isec]=0;
3e1f1ce7 1679 matSectorCluster(isec,0)=0;
1680 matSectorCluster(isec,1)=0;
65e67c9c 1681 if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+6.*(rms69SectorTime+ offsetTimeAccept))) {
1682 vecSectorOut6[isec]=1;
1683 }
1684 if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+9.*(rms69SectorTime+ offsetTimeAccept))){
1685 vecSectorOut9[isec]=1;
1686 }
e731a3f8 1687 }
1688 // light version of export variable
1689 Int_t filteredSector= vecSectorOut9.Sum(); // light version of export variable
1690 Int_t filteredSectorTime= vecMedianSectorTimeOut9.Sum();
1691 Int_t filteredSectorPadRow= vecMedianSectorPadRowOut9.Sum();
b9d88e03 1692 if (fEvent) if (fEvent->GetHeader()){
1693 fEvent->GetHeader()->SetTPCNoiseFilterCounter(0,TMath::Min(filteredSector,255));
1694 fEvent->GetHeader()->SetTPCNoiseFilterCounter(1,TMath::Min(filteredSectorTime,255));
1695 fEvent->GetHeader()->SetTPCNoiseFilterCounter(2,TMath::Min(filteredSectorPadRow,255));
eea44233 1696 }
65e67c9c 1697
44ca7282 1698 //
1699 // 4. Disabling clusters in outlier layers
1700 //
65e67c9c 1701 Int_t counterAll=0;
1702 Int_t counterOut=0;
44ca7282 1703 for (Int_t isector=0; isector<36; isector++){ // loop over sectors
1704 for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
1705 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1706 Int_t nrows = sector.GetNRows();
1707 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1708 AliTPCtrackerRow& tpcrow = sector[row];
1709 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1710 if (iside>0) ncl=tpcrow.GetN2();
1711 for (Int_t i=0;i<ncl;i++) { // loop over clusters
1712 AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1713 Double_t medianTime=vecMedianSectorTime[cluster->GetDetector()];
1714 Double_t medianPadRow=vecMedianSectorPadRow[cluster->GetDetector()];
1715 Double_t rmsTime=vecRMSSectorTime[cluster->GetDetector()];
1716 Double_t rmsPadRow=vecRMSSectorPadRow[cluster->GetDetector()];
1717 Int_t entriesPadRow=hisPadRow.GetBinContent(cluster->GetDetector()+1, cluster->GetRow()+1);
1718 Int_t entriesTime=hisTime.GetBinContent(cluster->GetDetector()+1, cluster->GetTimeBin()+1);
1719 Bool_t isOut=kFALSE;
65e67c9c 1720 if (vecSectorOut9[cluster->GetDetector()]>0.5) {
1721 isOut=kTRUE;
1722 }
e731a3f8 1723
3e1f1ce7 1724 if (entriesTime>medianTime+nSigmaCut*rmsTime+offsetTime) {
1725 isOut=kTRUE;
1726 vecMedianSectorTimeOut[cluster->GetDetector()]++;
1727 }
1728 if (entriesPadRow>medianPadRow+nSigmaCut*rmsPadRow+offsetPadRow) {
1729 isOut=kTRUE;
1730 vecMedianSectorPadRowOut[cluster->GetDetector()]++;
1731 }
65e67c9c 1732 counterAll++;
3e1f1ce7 1733 matSectorCluster(cluster->GetDetector(),0)+=1;
44ca7282 1734 if (isOut){
1735 cluster->Disable();
65e67c9c 1736 counterOut++;
3e1f1ce7 1737 matSectorCluster(cluster->GetDetector(),1)+=1;
44ca7282 1738 }
1739 }
1740 }
1741 }
1742 }
3e1f1ce7 1743 for (Int_t isec=0; isec<nSectors; isec++){
1744 if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector Stat: %d\t%8.0f\t%8.0f",isec,matSectorCluster(isec,1),matSectorCluster(isec,0)).Data());
1745 }
65e67c9c 1746 //
1747 // dump info to streamer - for later tuning of cuts
1748 //
1749 if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) { // stream TPC data ouliers filtering infomation
3e1f1ce7 1750 AliLog::Flush();
1751 AliInfo(TString::Format("Cluster counter: (%d/%d) (Filtered/All)",counterOut,counterAll).Data());
1752 for (Int_t iSec=0; iSec<nSectors; iSec++){
1753 if (vecSectorOut9[iSec]>0 || matSectorCluster(iSec,1)>0) {
1754 AliInfo(TString::Format("Filtered sector\t%d",iSec).Data());
1755 Double_t vecMedTime =TMath::Median(72,vecMedianSectorTime.GetMatrixArray());
1756 Double_t vecMedPadRow =TMath::Median(72,vecMedianSectorPadRow.GetMatrixArray());
1757 Double_t vecMedCluster=(counterAll-counterOut)/72;
1758 AliInfo(TString::Format("VecMedianSectorTime\t(%4.4f/%4.4f/%4.4f)", vecMedianSectorTimeOut[iSec],vecMedianSectorTime[iSec],vecMedTime).Data());
1759 AliInfo(TString::Format("VecMedianSectorPadRow\t(%4.4f/%4.4f/%4.4f)", vecMedianSectorPadRowOut[iSec],vecMedianSectorPadRow[iSec],vecMedPadRow).Data());
1760 AliInfo(TString::Format("MatSectorCluster\t(%4.4f/%4.4f/%4.4f)\n", matSectorCluster(iSec,1), matSectorCluster(iSec,0), vecMedCluster).Data());
1761 AliLog::Flush();
1762 }
1763 }
1764 AliLog::Flush();
1765 Int_t eventNr = fEvent->GetEventNumberInFile();
65e67c9c 1766 (*fDebugStreamer)<<"filterClusterInfo"<<
1767 // minimal set variables for the ESDevent
3e1f1ce7 1768 "eventNr="<<eventNr<<
65e67c9c 1769 "counterAll="<<counterAll<<
1770 "counterOut="<<counterOut<<
3c8c25a3 1771 "matSectotCluster.="<<&matSectorCluster<< //
65e67c9c 1772 //
1773 "filteredSector="<<filteredSector<< // counter filtered sectors
1774 "filteredSectorTime="<<filteredSectorTime<< // counter filtered time bins
1775 "filteredSectorPadRow="<<filteredSectorPadRow<< // counter filtered pad-rows
1776 // per sector outlier information
1777 "medianSectorTime="<<medianSectorTime<< // median number of clusters per sector/timebin
1778 "mean69SectorTime="<<mean69SectorTime<< // LTM statistic mean of clusters per sector/timebin
1779 "rms69SectorTime="<<rms69SectorTime<< // LTM statistic RMS of clusters per sector/timebin
1780 "vecSectorOut6.="<<&vecSectorOut6<< // flag array sector - 6 sigma +accept margin outlier
1781 "vecSectorOut9.="<<&vecSectorOut9<< // flag array sector - 9 sigma + accept margin outlier
1782 // per sector/timebin outlier detection
1783 "vecMedianSectorTime.="<<&vecMedianSectorTime<<
1784 "vecRMSSectorTime.="<<&vecRMSSectorTime<<
1785 "vecMedianSectorTimeOut6.="<<&vecMedianSectorTimeOut6<<
1786 "vecMedianSectorTimeOut9.="<<&vecMedianSectorTimeOut9<<
3e1f1ce7 1787 "vecMedianSectorTimeOut0.="<<&vecMedianSectorTimeOut<<
65e67c9c 1788 // per sector/pad-row outlier detection
1789 "vecMedianSectorPadRow.="<<&vecMedianSectorPadRow<<
1790 "vecRMSSectorPadRow.="<<&vecRMSSectorPadRow<<
1791 "vecMedianSectorPadRowOut6.="<<&vecMedianSectorPadRowOut6<<
1792 "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut9<<
3e1f1ce7 1793 "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut<<
65e67c9c 1794 "\n";
1795 ((*fDebugStreamer)<<"filterClusterInfo").GetTree()->Write();
1796 fDebugStreamer->GetFile()->Flush();
1797 }
eea44233 1798}
1799
829455ad 1800void AliTPCtracker::UnloadClusters()
91162307 1801{
1802 //
1803 // unload clusters from the memory
1804 //
1805 Int_t nrows = fOuterSec->GetNRows();
1806 for (Int_t sec = 0;sec<fkNOS;sec++)
1807 for (Int_t row = 0;row<nrows;row++){
bd26fa83 1808 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
982aff31 1809 // if (tpcrow){
1810 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1811 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1812 //}
1813 tpcrow->ResetClusters();
91162307 1814 }
1815 //
1816 nrows = fInnerSec->GetNRows();
1817 for (Int_t sec = 0;sec<fkNIS;sec++)
1818 for (Int_t row = 0;row<nrows;row++){
bd26fa83 1819 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
982aff31 1820 //if (tpcrow){
1821 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1822 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1823 //}
1824 tpcrow->ResetClusters();
91162307 1825 }
1826
1827 return ;
1828}
1829
829455ad 1830void AliTPCtracker::FillClusterArray(TObjArray* array) const{
002af263 1831 //
1832 // Filling cluster to the array - For visualization purposes
1833 //
1834 Int_t nrows=0;
1835 nrows = fOuterSec->GetNRows();
1836 for (Int_t sec = 0;sec<fkNOS;sec++)
1837 for (Int_t row = 0;row<nrows;row++){
1838 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1839 if (!tpcrow) continue;
1840 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1841 array->AddLast((TObject*)((*tpcrow)[icl]));
1842 }
1843 }
1844 nrows = fInnerSec->GetNRows();
1845 for (Int_t sec = 0;sec<fkNIS;sec++)
1846 for (Int_t row = 0;row<nrows;row++){
1847 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1848 if (!tpcrow) continue;
1849 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1850 array->AddLast((TObject*)(*tpcrow)[icl]);
1851 }
1852 }
1853}
1854
1855
829455ad 1856void AliTPCtracker::Transform(AliTPCclusterMI * cluster){
893468d9 1857 //
544c295f 1858 // transformation
893468d9 1859 //
3f3549a3 1860 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1861 AliTPCTransform *transform = calibDB->GetTransform() ;
24db6af7 1862 if (!transform) {
1863 AliFatal("Tranformations not in calibDB");
ec26e231 1864 return;
24db6af7 1865 }
239f39b1 1866 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2942f542 1867 Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),static_cast<Double_t>(cluster->GetPad()),static_cast<Double_t>(cluster->GetTimeBin())};
24db6af7 1868 Int_t i[1]={cluster->GetDetector()};
022ee144 1869 transform->Transform(x,i,0,1);
afbaa016 1870 // if (cluster->GetDetector()%36>17){
1871 // x[1]*=-1;
1872 //}
022ee144 1873
24db6af7 1874 //
1875 // in debug mode check the transformation
1876 //
65e67c9c 1877 if ((AliTPCReconstructor::StreamLevel()&kStreamTransform)>0) {
af86c1fd 1878 Float_t gx[3];
1879 cluster->GetGlobalXYZ(gx);
002af263 1880 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
24db6af7 1881 TTreeSRedirector &cstream = *fDebugStreamer;
65e67c9c 1882 cstream<<"Transform"<< // needed for debugging of the cluster transformation, resp. used for later visualization
002af263 1883 "event="<<event<<
24db6af7 1884 "x0="<<x[0]<<
1885 "x1="<<x[1]<<
1886 "x2="<<x[2]<<
af86c1fd 1887 "gx0="<<gx[0]<<
1888 "gx1="<<gx[1]<<
1889 "gx2="<<gx[2]<<
24db6af7 1890 "Cl.="<<cluster<<
1891 "\n";
1892 }
1893 cluster->SetX(x[0]);
1894 cluster->SetY(x[1]);
1895 cluster->SetZ(x[2]);
19b00333 1896 // The old stuff:
0201b65c 1897 //
1898 //
1899 //
47af7ca4 1900 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
3f3549a3 1901 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1902 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1903 //TGeoHMatrix mat;
1904 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1905 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1906 if (mat) mat->LocalToMaster(pos,posC);
1907 else{
1908 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1909 }
1910 cluster->SetX(posC[0]);
1911 cluster->SetY(posC[1]);
1912 cluster->SetZ(posC[2]);
a7760332 1913 }
0201b65c 1914}
1c53abe2 1915
265b5e37 1916void AliTPCtracker::ApplyXtalkCorrection(){
1917 //
1918 // ApplyXtalk correction
1919 // Loop over all clusters
1920 // add to each cluster signal corresponding to common Xtalk mode for given time bin at given wire segment
265b5e37 1921 // cluster loop
1922 for (Int_t isector=0; isector<36; isector++){ //loop tracking sectors
1923 for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
1924 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1925 Int_t nrows = sector.GetNRows();
1926 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1927 AliTPCtrackerRow& tpcrow = sector[row]; // row object
1928 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1929 if (iside>0) ncl=tpcrow.GetN2();
1930 Int_t xSector=0; // sector number in the TPC convention 0-72
1931 if (isector<18){ //if IROC
1932 xSector=isector+(iside>0)*18;
1933 }else{
1934 xSector=isector+18; // isector -18 +36
1935 if (iside>0) xSector+=18;
1936 }
1937 TMatrixD &crossTalkMatrix= *((TMatrixD*)fCrossTalkSignalArray->At(xSector));
1938 Int_t wireSegmentID = fkParam->GetWireSegment(xSector,row);
1939 for (Int_t i=0;i<ncl;i++) {
1940 AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1941 Int_t iTimeBin=TMath::Nint(cluster->GetTimeBin());
1942 Double_t xTalk= crossTalkMatrix[wireSegmentID][iTimeBin];
1943 cluster->SetMax(cluster->GetMax()+xTalk);
6f8ff889 1944 const Double_t kDummy=4;
265b5e37 1945 Double_t sumxTalk=xTalk*kDummy; // should be calculated via time response function
1946 cluster->SetQ(cluster->GetQ()+sumxTalk);
6f8ff889 1947
1948
65e67c9c 1949 if ((AliTPCReconstructor::StreamLevel()&kStreamXtalk)>0) { // flag: stream crosstalk correctio as applied to cluster
6f8ff889 1950 TTreeSRedirector &cstream = *fDebugStreamer;
1951 if (gRandom->Rndm() > 0.){
1952 cstream<<"Xtalk"<<
1953 "isector=" << isector << // sector [0,36]
1954 "iside=" << iside << // side A or C
1955 "row=" << row << // padrow
1956 "i=" << i << // index of the cluster
1957 "xSector=" << xSector << // sector [0,72]
1958 "wireSegmentID=" << wireSegmentID << // anode wire segment id [0,10]
1959 "iTimeBin=" << iTimeBin << // timebin of the corrected cluster
1960 "xTalk=" << xTalk << // Xtalk contribution added to Qmax
1961 "sumxTalk=" << sumxTalk << // Xtalk contribution added to Qtot (roughly 3*Xtalk)
1962 "cluster.=" << cluster << // corrected cluster object
1963 "\n";
1964 }
1965 }// dump the results to the debug streamer if in debug mode
1966
1967
1968
1969
1970
265b5e37 1971 }
1972 }
1973 }
1974 }
1975}
1976
9a836cc2 1977void AliTPCtracker::ApplyTailCancellation(){
7bc5f28b 1978 //
9a836cc2 1979 // Correct the cluster charge for the ion tail effect
7bc5f28b 1980 // The TimeResponse function accessed via AliTPCcalibDB (TPC/Calib/IonTail)
1981 //
7bc5f28b 1982
9a836cc2 1983 // Retrieve
1984 TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
1985 if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
1986 TObject *rocFactorIROC = ionTailArr->FindObject("factorIROC");
1987 TObject *rocFactorOROC = ionTailArr->FindObject("factorOROC");
1988 Float_t factorIROC = (atof(rocFactorIROC->GetTitle()));
1989 Float_t factorOROC = (atof(rocFactorOROC->GetTitle()));
1990
1991 // find the number of clusters for the whole TPC (nclALL)
1992 Int_t nclALL=0;
1993 for (Int_t isector=0; isector<36; isector++){
1994 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1995 nclALL += sector.GetNClInSector(0);
1996 nclALL += sector.GetNClInSector(1);
1997 }
1998
1999 // start looping over all clusters
2000 for (Int_t iside=0; iside<2; iside++){ // loop over sides
7bc5f28b 2001 //
2002 //
9a836cc2 2003 for (Int_t secType=0; secType<2; secType++){ //loop over inner or outer sector
2004 // cache experimantal tuning factor for the different chamber type
2005 const Float_t ampfactor = (secType==0)?factorIROC:factorOROC;
2006 std::cout << " ampfactor = " << ampfactor << std::endl;
7bc5f28b 2007 //
9a836cc2 2008 for (Int_t sec = 0;sec<fkNOS;sec++){ //loop overs sectors
2009 //
2010 //
2011 // Cache time response functions and their positons to COG of the cluster
2012 TGraphErrors ** graphRes = new TGraphErrors *[20];
2013 Float_t * indexAmpGraphs = new Float_t[20];
2014 for (Int_t icache=0; icache<20; icache++)
2015 {
2016 graphRes[icache] = NULL;
2017 indexAmpGraphs[icache] = 0;
2018 }
2019 ///////////////////////////// --> position fo sie loop
d4559772 2020 if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(sec+36*secType+18*iside,graphRes,indexAmpGraphs))
2021 {
2022 continue;
2023 }
9a836cc2 2024
2025 AliTPCtrackerSector &sector= (secType==0)?fInnerSec[sec]:fOuterSec[sec];
2026 Int_t nrows = sector.GetNRows(); // number of rows
2027 Int_t nclSector = sector.GetNClInSector(iside); // ncl per sector to be used for debugging
2028
2029 for (Int_t row = 0;row<nrows;row++){ // loop over rows
2030
2031 AliTPCtrackerRow& tpcrow = sector[row]; // row object
2032 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
2033 if (iside>0) ncl=tpcrow.GetN2();
2034
2035 // Order clusters in time for the proper correction of ion tail
2036 Float_t qTotArray[ncl]; // arrays to be filled with modified Qtot and Qmax values in order to avoid float->int conversion
2037 Float_t qMaxArray[ncl];
2038 Int_t sortedClusterIndex[ncl];
2039 Float_t sortedClusterTimeBin[ncl];
2040 TObjArray *rowClusterArray = new TObjArray(ncl); // cache clusters for each row
2041 for (Int_t i=0;i<ncl;i++)
2042 {
2043 qTotArray[i]=0;
2044 qMaxArray[i]=0;
2045 sortedClusterIndex[i]=i;
2046 AliTPCclusterMI *rowcl= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
2047 if (rowcl) {
2048 rowClusterArray->AddAt(rowcl,i);
2049 } else {
2050 rowClusterArray->RemoveAt(i);
2051 }
2052 // Fill the timebin info to the array in order to sort wrt tb
2053 if (!rowcl) {
2054 sortedClusterTimeBin[i]=0.0;
2055 } else {
2056 sortedClusterTimeBin[i] = rowcl->GetTimeBin();
2057 }
2058
2059 }
2060 TMath::Sort(ncl,sortedClusterTimeBin,sortedClusterIndex,kFALSE); // sort clusters in time
2061
2062 // Main cluster correction loops over clusters
2063 for (Int_t icl0=0; icl0<ncl;icl0++){ // first loop over clusters
2064
2065 AliTPCclusterMI *cl0= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl0]));
2066
2067 if (!cl0) continue;
2068 Int_t nclPad=0;
265b5e37 2069 for (Int_t icl1=0; icl1<ncl;icl1++){ // second loop over clusters
9a836cc2 2070 AliTPCclusterMI *cl1= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl1]));
2071 if (!cl1) continue;
2072 if (TMath::Abs(cl0->GetPad()-cl1->GetPad())>4) continue; // no contribution if far away in pad direction
2073 if (cl0->GetTimeBin()<= cl1->GetTimeBin()) continue; // no contibution to the tail if later
2074 if (TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin())>600) continue; // out of the range of response function
2075
2076 if (TMath::Abs(cl0->GetPad()-cl1->GetPad())<4) nclPad++; // count ncl for every pad for debugging
2077
2078 // Get the correction values for Qmax and Qtot and find total correction for a given cluster
2079 Double_t ionTailMax=0.;
2080 Double_t ionTailTotal=0.;
2081 GetTailValue(ampfactor,ionTailMax,ionTailTotal,graphRes,indexAmpGraphs,cl0,cl1);
2082 ionTailMax=TMath::Abs(ionTailMax);
2083 ionTailTotal=TMath::Abs(ionTailTotal);
2084 qTotArray[icl0]+=ionTailTotal;
2085 qMaxArray[icl0]+=ionTailMax;
2086
2087 // Dump some info for debugging while clusters are being corrected
65e67c9c 2088 if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) { // flag: stream ion tail correction as applied to cluster
9a836cc2 2089 TTreeSRedirector &cstream = *fDebugStreamer;
2090 if (gRandom->Rndm() > 0.999){
2091 cstream<<"IonTail"<<
2092 "cl0.=" <<cl0 << // cluster 0 (to be corrected)
2093 "cl1.=" <<cl1 << // cluster 1 (previous cluster)
2094 "ionTailTotal=" <<ionTailTotal << // ion Tail from cluster 1 contribution to cluster0
2095 "ionTailMax=" <<ionTailMax << // ion Tail from cluster 1 contribution to cluster0
2096 "\n";
2097 }
2098 }// dump the results to the debug streamer if in debug mode
2099
2100 }//end of second loop over clusters
2101
2102 // Set corrected values of the corrected cluster
2103 cl0->SetQ(TMath::Nint(Float_t(cl0->GetQ())+Float_t(qTotArray[icl0])));
2104 cl0->SetMax(TMath::Nint(Float_t(cl0->GetMax())+qMaxArray[icl0]));
2105
2106 // Dump some info for debugging after clusters are corrected
65e67c9c 2107 if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {
9a836cc2 2108 TTreeSRedirector &cstream = *fDebugStreamer;
2109 if (gRandom->Rndm() > 0.999){
2110 cstream<<"IonTailCorrected"<<
2111 "cl0.=" << cl0 << // cluster 0 with huge Qmax
2112 "ionTailTotalPerCluster=" << qTotArray[icl0] <<
2113 "ionTailMaxPerCluster=" << qMaxArray[icl0] <<
2114 "nclALL=" << nclALL <<
2115 "nclSector=" << nclSector <<
2116 "nclRow=" << ncl <<
2117 "nclPad=" << nclPad <<
2118 "row=" << row <<
2119 "sector=" << sec <<
2120 "icl0=" << icl0 <<
2121 "\n";
2122 }
2123 }// dump the results to the debug streamer if in debug mode
2124
2125 }//end of first loop over cluster
2126 delete rowClusterArray;
2127 }//end of loop over rows
2128 for (int i=0; i<20; i++) delete graphRes[i];
2129 delete [] graphRes;
2130 delete [] indexAmpGraphs;
2131
2132 }//end of loop over sectors
2133 }//end of loop over IROC/OROC
2134 }// end of side loop
7bc5f28b 2135}
9a836cc2 2136//_____________________________________________________________________________
5a516e0a 2137void AliTPCtracker::GetTailValue(Float_t ampfactor,Double_t &ionTailMax, Double_t &ionTailTotal,TGraphErrors **graphRes,Float_t *indexAmpGraphs,AliTPCclusterMI *cl0,AliTPCclusterMI *cl1){
9a836cc2 2138
2139 //
2140 // Function in order to calculate the amount of the correction to be added for a given cluster, return values are ionTailTaoltal and ionTailMax
265b5e37 2141 // Parameters:
2142 // cl0 - cluster to be modified
2143 // cl1 - source cluster ion tail of this cluster will be added to the cl0 (accroding time and pad response function)
9a836cc2 2144 //
9a836cc2 2145 const Double_t kMinPRF = 0.5; // minimal PRF width
2146 ionTailTotal = 0.; // correction value to be added to Qtot of cl0
2147 ionTailMax = 0.; // correction value to be added to Qmax of cl0
2148
2149 Float_t qTot0 = cl0->GetQ(); // cl0 Qtot info
2150 Float_t qTot1 = cl1->GetQ(); // cl1 Qtot info
2151 Int_t sectorPad = cl1->GetDetector(); // sector number
2152 Int_t padcl0 = TMath::Nint(cl0->GetPad()); // pad0
2153 Int_t padcl1 = TMath::Nint(cl1->GetPad()); // pad1
2154 Float_t padWidth = (sectorPad < 36)?0.4:0.6; // pad width in cm
2155 const Int_t deltaTimebin = TMath::Nint(TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin()))+12; //distance between pads of cl1 and cl0 increased by 12 bins
2156 Double_t rmsPad1 = (cl1->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl1->GetSigmaY2())/padWidth);
2157 Double_t rmsPad0 = (cl0->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl0->GetSigmaY2())/padWidth);
265b5e37 2158
2159
a3986da2 2160
9a836cc2 2161 Double_t sumAmp1=0.;
2162 for (Int_t idelta =-2; idelta<=2;idelta++){
2163 sumAmp1+=TMath::Exp(-idelta*idelta/(2*rmsPad1));
2164 }
7bc5f28b 2165
9a836cc2 2166 Double_t sumAmp0=0.;
2167 for (Int_t idelta =-2; idelta<=2;idelta++){
2168 sumAmp0+=TMath::Exp(-idelta*idelta/(2*rmsPad0));
2169 }
7bc5f28b 2170
9a836cc2 2171 // Apply the correction --> cl1 corrects cl0 (loop over cl1's pads and find which pads of cl0 are going to be corrected)
2172 Int_t padScan=2; // +-2 pad-timebin window will be scanned
2173 for (Int_t ipad1=padcl1-padScan; ipad1<=padcl1+padScan; ipad1++) {
2174 //
2175 //
2176 Float_t deltaPad1 = TMath::Abs(cl1->GetPad()-(Float_t)ipad1);
2177 Double_t amp1 = (TMath::Exp(-(deltaPad1*deltaPad1)/(2*rmsPad1)))/sumAmp1; // normalized pad response function
2178 Float_t qTotPad1 = amp1*qTot1; // used as a factor to multipliy the response function
2179
2180 // find closest value of cl1 to COG (among the time response functions' amplitude array --> to select proper t.r.f.)
2181 Int_t ampIndex = 0;
2182 Float_t diffAmp = TMath::Abs(deltaPad1-indexAmpGraphs[0]);
2183 for (Int_t j=0;j<20;j++) {
2184 if (diffAmp > TMath::Abs(deltaPad1-indexAmpGraphs[j]) && indexAmpGraphs[j]!=0)
2185 {
2186 diffAmp = TMath::Abs(deltaPad1-indexAmpGraphs[j]);
2187 ampIndex = j;
2188 }
2189 }
2190 if (!graphRes[ampIndex]) continue;
b4742ddd 2191 if (deltaTimebin+2 >= graphRes[ampIndex]->GetN()) continue;
9a836cc2 2192 if (graphRes[ampIndex]->GetY()[deltaTimebin+2]>=0) continue;
2193
2194 for (Int_t ipad0=padcl0-padScan; ipad0<=padcl0+padScan; ipad0++) {
2195 //
2196 //
2197 if (ipad1!=ipad0) continue; // check if ipad1 channel sees ipad0 channel, if not no correction to be applied.
2198
2199 Float_t deltaPad0 = TMath::Abs(cl0->GetPad()-(Float_t)ipad0);
2200 Double_t amp0 = (TMath::Exp(-(deltaPad0*deltaPad0)/(2*rmsPad0)))/sumAmp0; // normalized pad resp function
2201 Float_t qMaxPad0 = amp0*qTot0;
2202
2203 // Add 5 timebin range contribution around the max peak (-+2 tb window)
2204 for (Int_t itb=deltaTimebin-2; itb<=deltaTimebin+2; itb++) {
2205
2206 if (itb<0) continue;
2207 if (itb>=graphRes[ampIndex]->GetN()) continue;
2208
2209 // calculate contribution to qTot
2210 Float_t tailCorr = TMath::Abs((qTotPad1*ampfactor)*(graphRes[ampIndex])->GetY()[itb]);
2211 if (ipad1!=padcl0) {
2212 ionTailTotal += TMath::Min(qMaxPad0,tailCorr); // for side pad
2213 } else {
2214 ionTailTotal += tailCorr; // for center pad
2215 }
2216 // calculate contribution to qMax
2217 if (itb == deltaTimebin && ipad1 == padcl0) ionTailMax += tailCorr;
2218
2219 } // end of tb correction loop which is applied over 5 tb range
2220
2221 } // end of cl0 loop
2222 } // end of cl1 loop
2223
2224}
7bc5f28b 2225
1c53abe2 2226//_____________________________________________________________________________
829455ad 2227Int_t AliTPCtracker::LoadOuterSectors() {
1c53abe2 2228 //-----------------------------------------------------------------
2229 // This function fills outer TPC sectors with clusters.
2230 //-----------------------------------------------------------------
91162307 2231 Int_t nrows = fOuterSec->GetNRows();
2232 UInt_t index=0;
2233 for (Int_t sec = 0;sec<fkNOS;sec++)
2234 for (Int_t row = 0;row<nrows;row++){
bd26fa83 2235 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
91162307 2236 Int_t sec2 = sec+2*fkNIS;
2237 //left
b9671574 2238 Int_t ncl = tpcrow->GetN1();
91162307 2239 while (ncl--) {
b9671574 2240 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
91162307 2241 index=(((sec2<<8)+row)<<16)+ncl;
2242 tpcrow->InsertCluster(c,index);
2243 }
2244 //right
b9671574 2245 ncl = tpcrow->GetN2();
91162307 2246 while (ncl--) {
b9671574 2247 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
91162307 2248 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
2249 tpcrow->InsertCluster(c,index);
2250 }
2251 //
2252 // write indexes for fast acces
2253 //
2254 for (Int_t i=0;i<510;i++)
b9671574 2255 tpcrow->SetFastCluster(i,-1);
91162307 2256 for (Int_t i=0;i<tpcrow->GetN();i++){
2257 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
b9671574 2258 tpcrow->SetFastCluster(zi,i); // write index
91162307 2259 }
2260 Int_t last = 0;
2261 for (Int_t i=0;i<510;i++){
b9671574 2262 if (tpcrow->GetFastCluster(i)<0)
2263 tpcrow->SetFastCluster(i,last);
91162307 2264 else
b9671574 2265 last = tpcrow->GetFastCluster(i);
91162307 2266 }
2267 }
1c53abe2 2268 fN=fkNOS;
2269 fSectors=fOuterSec;
91162307 2270 return 0;
1c53abe2 2271}
2272
2273
2274//_____________________________________________________________________________
829455ad 2275Int_t AliTPCtracker::LoadInnerSectors() {
1c53abe2 2276 //-----------------------------------------------------------------
2277 // This function fills inner TPC sectors with clusters.
2278 //-----------------------------------------------------------------
91162307 2279 Int_t nrows = fInnerSec->GetNRows();
2280 UInt_t index=0;
2281 for (Int_t sec = 0;sec<fkNIS;sec++)
2282 for (Int_t row = 0;row<nrows;row++){
bd26fa83 2283 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
91162307 2284 //
2285 //left
b9671574 2286 Int_t ncl = tpcrow->GetN1();
91162307 2287 while (ncl--) {
b9671574 2288 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
91162307 2289 index=(((sec<<8)+row)<<16)+ncl;
2290 tpcrow->InsertCluster(c,index);
2291 }
2292 //right
b9671574 2293 ncl = tpcrow->GetN2();
91162307 2294 while (ncl--) {
b9671574 2295 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
91162307 2296 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
2297 tpcrow->InsertCluster(c,index);
2298 }
2299 //
2300 // write indexes for fast acces
2301 //
2302 for (Int_t i=0;i<510;i++)
b9671574 2303 tpcrow->SetFastCluster(i,-1);
91162307 2304 for (Int_t i=0;i<tpcrow->GetN();i++){
2305 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
b9671574 2306 tpcrow->SetFastCluster(zi,i); // write index
91162307 2307 }
2308 Int_t last = 0;
2309 for (Int_t i=0;i<510;i++){
b9671574 2310 if (tpcrow->GetFastCluster(i)<0)
2311 tpcrow->SetFastCluster(i,last);
91162307 2312 else
b9671574 2313 last = tpcrow->GetFastCluster(i);
91162307 2314 }
1c53abe2 2315
91162307 2316 }
2317
1c53abe2 2318 fN=fkNIS;
2319 fSectors=fInnerSec;
91162307 2320 return 0;
2321}
2322
2323
2324
2325//_________________________________________________________________________
829455ad 2326AliTPCclusterMI *AliTPCtracker::GetClusterMI(Int_t index) const {
91162307 2327 //--------------------------------------------------------------------
2328 // Return pointer to a given cluster
2329 //--------------------------------------------------------------------
e7eb17e4 2330 if (index<0) return 0; // no cluster
91162307 2331 Int_t sec=(index&0xff000000)>>24;
2332 Int_t row=(index&0x00ff0000)>>16;
d26d9159 2333 Int_t ncl=(index&0x00007fff)>>00;
91162307 2334
bd26fa83 2335 const AliTPCtrackerRow * tpcrow=0;
bfa00fba 2336 TClonesArray * clrow =0;
3da363a1 2337
ad23441a 2338 if (sec<0 || sec>=fkNIS*4) {
2339 AliWarning(Form("Wrong sector %d",sec));
2340 return 0x0;
2341 }
2342
91162307 2343 if (sec<fkNIS*2){
462fb6e0 2344 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
2345 if (tracksec.GetNRows()<=row) return 0;
2346 tpcrow = &(tracksec[row]);
3da363a1 2347 if (tpcrow==0) return 0;
2348
2349 if (sec<fkNIS) {
b9671574 2350 if (tpcrow->GetN1()<=ncl) return 0;
2351 clrow = tpcrow->GetClusters1();
3da363a1 2352 }
2353 else {
b9671574 2354 if (tpcrow->GetN2()<=ncl) return 0;
2355 clrow = tpcrow->GetClusters2();
3da363a1 2356 }
91162307 2357 }
3da363a1 2358 else {
462fb6e0 2359 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
2360 if (tracksec.GetNRows()<=row) return 0;
2361 tpcrow = &(tracksec[row]);
3da363a1 2362 if (tpcrow==0) return 0;
2363
2364 if (sec-2*fkNIS<fkNOS) {
b9671574 2365 if (tpcrow->GetN1()<=ncl) return 0;
2366 clrow = tpcrow->GetClusters1();
3da363a1 2367 }
2368 else {
b9671574 2369 if (tpcrow->GetN2()<=ncl) return 0;
2370 clrow = tpcrow->GetClusters2();
3da363a1 2371 }
91162307 2372 }
3da363a1 2373
bfa00fba 2374 return (AliTPCclusterMI*)clrow->At(ncl);
91162307 2375
1c53abe2 2376}
2377
91162307 2378
2379
829455ad 2380Int_t AliTPCtracker::FollowToNext(AliTPCseed& t, Int_t nr) {
1c53abe2 2381 //-----------------------------------------------------------------
2382 // This function tries to find a track prolongation to next pad row
2383 //-----------------------------------------------------------------
1c53abe2 2384 //
91162307 2385 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
76d56fd6 2386 //
2387 //
4d158c36 2388 AliTPCclusterMI *cl=0;
2389 Int_t tpcindex= t.GetClusterIndex2(nr);
2390 //
2391 // update current shape info every 5 pad-row
2392 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
2393 GetShape(&t,nr);
2394 //}
2395 //
2396 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
2397 //
2398 if (tpcindex==-1) return 0; //track in dead zone
f124f8bf 2399 if (tpcindex >= 0){ //
b9671574 2400 cl = t.GetClusterPointer(nr);
97d77e7a 2401 //if (cl==0) cl = GetClusterMI(tpcindex);
2402 if (!cl) cl = GetClusterMI(tpcindex);
b9671574 2403 t.SetCurrentClusterIndex1(tpcindex);
4d158c36 2404 }
2405 if (cl){
2406 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
2407 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
2408 //
2409 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
2410 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
2411
2412 if (TMath::Abs(angle-t.GetAlpha())>0.001){
2413 Double_t rotation = angle-t.GetAlpha();
b9671574 2414 t.SetRelativeSector(relativesector);
1791d824 2415 if (!t.Rotate(rotation)) {
2416 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2417 return 0;
2418 }
2419 }
2420 if (!t.PropagateTo(x)) {
2421 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2422 return 0;
4d158c36 2423 }
4d158c36 2424 //
b9671574 2425 t.SetCurrentCluster(cl);
2426 t.SetRow(nr);
00055a22 2427 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
4d158c36 2428 if ((tpcindex&0x8000)==0) accept =0;
2429 if (accept<3) {
2430 //if founded cluster is acceptible
2431 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
b9671574 2432 t.SetErrorY2(t.GetErrorY2()+0.03);
2433 t.SetErrorZ2(t.GetErrorZ2()+0.03);
2434 t.SetErrorY2(t.GetErrorY2()*3);
2435 t.SetErrorZ2(t.GetErrorZ2()*3);
4d158c36 2436 }
b9671574 2437 t.SetNFoundable(t.GetNFoundable()+1);
4d158c36 2438 UpdateTrack(&t,accept);
2439 return 1;
f124f8bf 2440 }
2441 else { // Remove old cluster from track
2442 t.SetClusterIndex(nr, -3);
2443 t.SetClusterPointer(nr, 0);
2444 }
4d158c36 2445 }
1627d1c4 2446 }
3f82c4f2 2447 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
76d56fd6 2448 if (fIteration>1 && IsFindable(t)){
3f82c4f2 2449 // not look for new cluster during refitting
b9671574 2450 t.SetNFoundable(t.GetNFoundable()+1);
3f82c4f2 2451 return 0;
2452 }
91162307 2453 //
4d158c36 2454 UInt_t index=0;
ca142b1f 2455 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
f124f8bf 2456 if (!t.PropagateTo(x)) {
2457 if (fIteration==0) t.SetRemoval(10);
2458 return 0;
2459 }
2460 Double_t y = t.GetY();
91162307 2461 if (TMath::Abs(y)>ymax){
2462 if (y > ymax) {
b9671574 2463 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
91162307 2464 if (!t.Rotate(fSectors->GetAlpha()))
2465 return 0;
2466 } else if (y <-ymax) {
b9671574 2467 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
91162307 2468 if (!t.Rotate(-fSectors->GetAlpha()))
2469 return 0;
2470 }
f124f8bf 2471 if (!t.PropagateTo(x)) {
2472 if (fIteration==0) t.SetRemoval(10);
2473 return 0;
2474 }
2475 y = t.GetY();
91162307 2476 }
2477 //
4d158c36 2478 Double_t z=t.GetZ();
2479 //
a3232aae 2480
b9671574 2481 if (!IsActive(t.GetRelativeSector(),nr)) {
2482 t.SetInDead(kTRUE);
a3232aae 2483 t.SetClusterIndex2(nr,-1);
2484 return 0;
2485 }
2486 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
b9671574 2487 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
2488 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
a3232aae 2489
2490 if (!isActive || !isActive2) return 0;
2491
bd26fa83 2492 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
91162307 2493 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
2494 Double_t roady =1.;
2495 Double_t roadz = 1.;
2496 //
b9671574 2497 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2498 t.SetInDead(kTRUE);
91162307 2499 t.SetClusterIndex2(nr,-1);
1c53abe2 2500 return 0;
2501 }
2502 else
2503 {
76d56fd6 2504 if (IsFindable(t))
47af7ca4 2505 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
b9671574 2506 t.SetNFoundable(t.GetNFoundable()+1);
1627d1c4 2507 else
2508 return 0;
1c53abe2 2509 }
2510 //calculate
91162307 2511 if (krow) {
2512 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
2513 cl = krow.FindNearest2(y,z,roady,roadz,index);
b9671574 2514 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
91162307 2515 }
91162307 2516 if (cl) {
b9671574 2517 t.SetCurrentCluster(cl);
2518 t.SetRow(nr);
4d158c36 2519 if (fIteration==2&&cl->IsUsed(10)) return 0;
00055a22 2520 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
4d158c36 2521 if (fIteration==2&&cl->IsUsed(11)) {
b9671574 2522 t.SetErrorY2(t.GetErrorY2()+0.03);
2523 t.SetErrorZ2(t.GetErrorZ2()+0.03);
2524 t.SetErrorY2(t.GetErrorY2()*3);
2525 t.SetErrorZ2(t.GetErrorZ2()*3);
4d158c36 2526 }
d26d9159 2527 /*
91162307 2528 if (t.fCurrentCluster->IsUsed(10)){
2529 //
2530 //
c9427e08 2531
91162307 2532 t.fNShared++;
2533 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
2534 t.fRemoval =10;
2535 return 0;
2536 }
2537 }
d26d9159 2538 */
91162307 2539 if (accept<3) UpdateTrack(&t,accept);
c9427e08 2540
91162307 2541 } else {
b9671574 2542 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
91162307 2543
2544 }
2545 return 1;
2546}
c9427e08 2547
1c53abe2 2548
91162307 2549
5d837844 2550//_________________________________________________________________________
829455ad 2551Bool_t AliTPCtracker::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
5d837844 2552{
2553 // Get track space point by index
2554 // return false in case the cluster doesn't exist
2555 AliTPCclusterMI *cl = GetClusterMI(index);
2556 if (!cl) return kFALSE;
2557 Int_t sector = (index&0xff000000)>>24;
0201b65c 2558 // Int_t row = (index&0x00ff0000)>>16;
5d837844 2559 Float_t xyz[3];
47af7ca4 2560 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
0201b65c 2561 xyz[0] = cl->GetX();
5d837844 2562 xyz[1] = cl->GetY();
2563 xyz[2] = cl->GetZ();
2564 Float_t sin,cos;
47af7ca4 2565 fkParam->AdjustCosSin(sector,cos,sin);
5d837844 2566 Float_t x = cos*xyz[0]-sin*xyz[1];
2567 Float_t y = cos*xyz[1]+sin*xyz[0];
2568 Float_t cov[6];
2569 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
47af7ca4 2570 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
5d837844 2571 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
47af7ca4 2572 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
5d837844 2573 cov[0] = sin*sin*sigmaY2;
2574 cov[1] = -sin*cos*sigmaY2;
2575 cov[2] = 0.;
2576 cov[3] = cos*cos*sigmaY2;
2577 cov[4] = 0.;
2578 cov[5] = sigmaZ2;
2579 p.SetXYZ(x,y,xyz[2],cov);
ae079791 2580 AliGeomManager::ELayerID iLayer;
5d837844 2581 Int_t idet;
47af7ca4 2582 if (sector < fkParam->GetNInnerSector()) {
ae079791 2583 iLayer = AliGeomManager::kTPC1;
5d837844 2584 idet = sector;
2585 }
2586 else {
ae079791 2587 iLayer = AliGeomManager::kTPC2;
47af7ca4 2588 idet = sector - fkParam->GetNInnerSector();
5d837844 2589 }
ae079791 2590 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
5d837844 2591 p.SetVolumeID(volid);
2592 return kTRUE;
2593}
2594
2595
2596
829455ad 2597Int_t AliTPCtracker::UpdateClusters(AliTPCseed& t, Int_t nr) {
1c53abe2 2598 //-----------------------------------------------------------------
2599 // This function tries to find a track prolongation to next pad row
2600 //-----------------------------------------------------------------
b9671574 2601 t.SetCurrentCluster(0);
f124f8bf 2602 t.SetCurrentClusterIndex1(-3);
1c53abe2 2603
2604 Double_t xt=t.GetX();
91162307 2605 Int_t row = GetRowNumber(xt)-1;
2606 Double_t ymax= GetMaxY(nr);
2607
1c53abe2 2608 if (row < nr) return 1; // don't prolongate if not information until now -
ca142b1f 2609// if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
2610// t.fRemoval =10;
2611// return 0; // not prolongate strongly inclined tracks
2612// }
2613// if (TMath::Abs(t.GetSnp())>0.95) {
2614// t.fRemoval =10;
2615// return 0; // not prolongate strongly inclined tracks
2616// }// patch 28 fev 06
91162307 2617
2618 Double_t x= GetXrow(nr);
2619 Double_t y,z;
2620 //t.PropagateTo(x+0.02);
2621 //t.PropagateTo(x+0.01);
1627d1c4 2622 if (!t.PropagateTo(x)){
1627d1c4 2623 return 0;
2624 }
1c53abe2 2625 //
91162307 2626 y=t.GetY();
2627 z=t.GetZ();
1c53abe2 2628
91162307 2629 if (TMath::Abs(y)>ymax){
2630 if (y > ymax) {
b9671574 2631 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
91162307 2632 if (!t.Rotate(fSectors->GetAlpha()))
2633 return 0;
2634 } else if (y <-ymax) {
b9671574 2635 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
91162307 2636 if (!t.Rotate(-fSectors->GetAlpha()))
2637 return 0;
2638 }
982aff31 2639 // if (!t.PropagateTo(x)){
2640 // return 0;
2641 //}
2642 return 1;
2643 //y = t.GetY();
1c53abe2 2644 }
91162307 2645 //
3f82c4f2 2646 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
a3232aae 2647
b9671574 2648 if (!IsActive(t.GetRelativeSector(),nr)) {
2649 t.SetInDead(kTRUE);
a3232aae 2650 t.SetClusterIndex2(nr,-1);
2651 return 0;
2652 }
2653 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
2654
bd26fa83 2655 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1c53abe2 2656
b9671574 2657 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2658 t.SetInDead(kTRUE);
91162307 2659 t.SetClusterIndex2(nr,-1);
1c53abe2 2660 return 0;
2661 }
2662 else
2663 {
76d56fd6 2664
2665 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
2666 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1627d1c4 2667 else
2668 return 0;
1c53abe2 2669 }
1c53abe2 2670
91162307 2671 // update current
42aec1f1 2672 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
91162307 2673 // t.fCurrentSigmaY = GetSigmaY(&t);
2674 //t.fCurrentSigmaZ = GetSigmaZ(&t);
2675 GetShape(&t,nr);
2676 }
1c53abe2 2677
91162307 2678 AliTPCclusterMI *cl=0;
d184f295 2679 Int_t index=0;
91162307 2680 //
2681 Double_t roady = 1.;
2682 Double_t roadz = 1.;
2683 //
d26d9159 2684
2685 if (!cl){
2686 index = t.GetClusterIndex2(nr);
bad6eb00 2687 if ( (index >= 0) && (index&0x8000)==0){
b9671574 2688 cl = t.GetClusterPointer(nr);
bad6eb00 2689 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
b9671574 2690 t.SetCurrentClusterIndex1(index);
d26d9159 2691 if (cl) {
b9671574 2692 t.SetCurrentCluster(cl);
d26d9159 2693 return 1;
2694 }
2695 }
2696 }
2697
3d172d79 2698 // if (index<0) return 0;
2699 UInt_t uindex = TMath::Abs(index);
d184f295 2700
91162307 2701 if (krow) {
d184f295 2702 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
2703 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
91162307 2704 }
d26d9159 2705
b9671574 2706 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
2707 t.SetCurrentCluster(cl);
d26d9159 2708
91162307 2709 return 1;
2710}
1c53abe2 2711
1c53abe2 2712
829455ad 2713Int_t AliTPCtracker::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
91162307 2714 //-----------------------------------------------------------------
2715 // This function tries to find a track prolongation to next pad row
2716 //-----------------------------------------------------------------
1c53abe2 2717
91162307 2718 //update error according neighborhoud
1c53abe2 2719
b9671574 2720 if (t.GetCurrentCluster()) {
2721 t.SetRow(nr);
00055a22 2722 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
91162307 2723
b9671574 2724 if (t.GetCurrentCluster()->IsUsed(10)){
91162307 2725 //
2726 //
2727 // t.fErrorZ2*=2;
2728 // t.fErrorY2*=2;
b9671574 2729 t.SetNShared(t.GetNShared()+1);
2730 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
2731 t.SetRemoval(10);
91162307 2732 return 0;
2733 }
b364ca79 2734 }
d26d9159 2735 if (fIteration>0) accept = 0;
b364ca79 2736 if (accept<3) UpdateTrack(&t,accept);
2737
1c53abe2 2738 } else {
91162307 2739 if (fIteration==0){
f124f8bf 2740 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
2741 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
91162307 2742
b9671574 2743 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1c53abe2 2744 }
2745 }
2746 return 1;
2747}
2748
2749
2750
91162307 2751//_____________________________________________________________________________
829455ad 2752Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1c53abe2 2753 //-----------------------------------------------------------------
91162307 2754 // This function tries to find a track prolongation.
1c53abe2 2755 //-----------------------------------------------------------------
2756 Double_t xt=t.GetX();
2757 //
f124f8bf 2758 Double_t alpha=t.GetAlpha();
1c53abe2 2759 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2760 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 2761 //
a3f36f42 2762 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1c53abe2 2763
f124f8bf 2764 Int_t first = GetRowNumber(xt);
2765 if (!fromSeeds)
2766 first -= step;
a3f36f42 2767 if (first < 0)
2768 first = 0;
51ad6848 2769 for (Int_t nr= first; nr>=rf; nr-=step) {
2770 // update kink info
2771 if (t.GetKinkIndexes()[0]>0){
2772 for (Int_t i=0;i<3;i++){
2773 Int_t index = t.GetKinkIndexes()[i];
2774 if (index==0) break;
2775 if (index<0) continue;
2776 //
6c94f330 2777 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
51ad6848 2778 if (!kink){
2779 printf("PROBLEM\n");
2780 }
2781 else{
eea478d3 2782 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
51ad6848 2783 if (kinkrow==nr){
2784 AliExternalTrackParam paramd(t);
2785 kink->SetDaughter(paramd);
eea478d3 2786 kink->SetStatus(2,5);
51ad6848 2787 kink->Update();
51ad6848 2788 }
2789 }
2790 }
2791 }
2792
2793 if (nr==80) t.UpdateReference();
982aff31 2794 if (nr<fInnerSec->GetNRows())
2795 fSectors = fInnerSec;
2796 else
2797 fSectors = fOuterSec;
91162307 2798 if (FollowToNext(t,nr)==0)
4d158c36 2799 if (!t.IsActive())
2800 return 0;
91162307 2801
2802 }
2803 return 1;
2804}
2805
1c53abe2 2806
1c53abe2 2807
91162307 2808
2809
2810
829455ad 2811Int_t AliTPCtracker::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
1c53abe2 2812 //-----------------------------------------------------------------
2813 // This function tries to find a track prolongation.
2814 //-----------------------------------------------------------------
1c53abe2 2815 //
eea478d3 2816 Double_t xt=t.GetX();
f124f8bf 2817 Double_t alpha=t.GetAlpha();
1c53abe2 2818 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2819 if (alpha < 0. ) alpha += 2.*TMath::Pi();
bad6eb00 2820 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1c53abe2 2821
b9671574 2822 Int_t first = t.GetFirstPoint();
f124f8bf 2823 Int_t ri = GetRowNumber(xt);
2824 if (!fromSeeds)
2825 ri += 1;
2826
2827 if (first<ri) first = ri;
91162307 2828 //
2829 if (first<0) first=0;
4d158c36 2830 for (Int_t nr=first; nr<=rf; nr++) {
ca142b1f 2831 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
51ad6848 2832 if (t.GetKinkIndexes()[0]<0){
2833 for (Int_t i=0;i<3;i++){
2834 Int_t index = t.GetKinkIndexes()[i];
2835 if (index==0) break;
2836 if (index>0) continue;
2837 index = TMath::Abs(index);
6c94f330 2838 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
51ad6848 2839 if (!kink){
2840 printf("PROBLEM\n");
2841 }
2842 else{
eea478d3 2843 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
51ad6848 2844 if (kinkrow==nr){
2845 AliExternalTrackParam paramm(t);
2846 kink->SetMother(paramm);
eea478d3 2847 kink->SetStatus(2,1);
51ad6848 2848 kink->Update();
51ad6848 2849 }
2850 }
eea478d3 2851 }
51ad6848 2852 }
eea478d3 2853 //
d26d9159 2854 if (nr<fInnerSec->GetNRows())
2855 fSectors = fInnerSec;
2856 else
2857 fSectors = fOuterSec;
c274e255 2858
d26d9159 2859 FollowToNext(t,nr);
1c53abe2 2860 }
2861 return 1;
2862}
2863
2864
2865
2866
2867
829455ad 2868Float_t AliTPCtracker::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1c53abe2 2869{
544c295f 2870 // overlapping factor
1c53abe2 2871 //
2872 sum1=0;
2873 sum2=0;
2874 Int_t sum=0;
1c53abe2 2875 //
2876 Float_t dz2 =(s1->GetZ() - s2->GetZ());
c9427e08 2877 dz2*=dz2;
91162307 2878
c9427e08 2879 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1c53abe2 2880 dy2*=dy2;
2881 Float_t distance = TMath::Sqrt(dz2+dy2);
c9427e08 2882 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1c53abe2 2883
91162307 2884 // Int_t offset =0;
b9671574 2885 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2886 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
c9427e08 2887 if (lastpoint>160)
2888 lastpoint =160;
2889 if (firstpoint<0)
2890 firstpoint = 0;
91162307 2891 if (firstpoint>lastpoint) {
2892 firstpoint =lastpoint;
2893 // lastpoint =160;
c9427e08 2894 }
2895
2896
91162307 2897 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2898 if (s1->GetClusterIndex2(i)>0) sum1++;
2899 if (s2->GetClusterIndex2(i)>0) sum2++;
2900 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1c53abe2 2901 sum++;
2902 }
2903 }
91162307 2904 if (sum<5) return 0;
2905
1627d1c4 2906 Float_t summin = TMath::Min(sum1+1,sum2+1);
2907 Float_t ratio = (sum+1)/Float_t(summin);
1c53abe2 2908 return ratio;
2909}
2910
829455ad 2911void AliTPCtracker::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1c53abe2 2912{
544c295f 2913 // shared clusters
1c53abe2 2914 //
a0f4d6a6 2915 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2916 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2917 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2918 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
91162307 2919
1c53abe2 2920 //
91162307 2921 Int_t sumshared=0;
2922 //
a0f4d6a6 2923 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2924 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2925 Int_t firstpoint = 0;
2926 Int_t lastpoint = 160;
91162307 2927 //
a0f4d6a6 2928 // if (firstpoint>=lastpoint-5) return;;
1af5da7e 2929
91162307 2930 for (Int_t i=firstpoint;i<lastpoint;i++){
2931 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2932 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2933 sumshared++;
2934 }
2935 }
a0f4d6a6 2936 if (sumshared>cutN0){
91162307 2937 // sign clusters
2938 //
2939 for (Int_t i=firstpoint;i<lastpoint;i++){
2940 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2941 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2942 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2943 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2944 if (s1->IsActive()&&s2->IsActive()){
b9671574 2945 p1->SetShared(kTRUE);
2946 p2->SetShared(kTRUE);
91162307 2947 }
2948 }
2949 }
2950 }
2951 //
a0f4d6a6 2952 if (sumshared>cutN0){
91162307 2953 for (Int_t i=0;i<4;i++){
b9671574 2954 if (s1->GetOverlapLabel(3*i)==0){
2955 s1->SetOverlapLabel(3*i, s2->GetLabel());
2956 s1->SetOverlapLabel(3*i+1,sumshared);
2957 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
91162307 2958 break;
2959 }
2960 }
2961 for (Int_t i=0;i<4;i++){
b9671574 2962 if (s2->GetOverlapLabel(3*i)==0){
2963 s2->SetOverlapLabel(3*i, s1->GetLabel());
2964 s2->SetOverlapLabel(3*i+1,sumshared);
2965 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
91162307 2966 break;
2967 }
2968 }
2969 }
91162307 2970}
1c53abe2 2971
829455ad 2972void AliTPCtracker::SignShared(TObjArray * arr)
91162307 2973{
1c53abe2 2974 //
91162307 2975 //sort trackss according sectors
2976 //
c9427e08 2977 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2978 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 2979 if (!pt) continue;
2980 //if (pt) RotateToLocal(pt);
b9671574 2981 pt->SetSort(0);
c9427e08 2982 }
91162307 2983 arr->UnSort();
6d493ea0 2984 arr->Sort(); // sorting according relative sectors
1c53abe2 2985 arr->Expand(arr->GetEntries());
91162307 2986 //
2987 //
1c53abe2 2988 Int_t nseed=arr->GetEntriesFast();
1c53abe2 2989 for (Int_t i=0; i<nseed; i++) {
2990 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 2991 if (!pt) continue;
ec26e231 2992 for (Int_t j=0;j<12;j++){
b9671574 2993 pt->SetOverlapLabel(j,0);
1c53abe2 2994 }
91162307 2995 }
2996 for (Int_t i=0; i<nseed; i++) {
2997 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2998 if (!pt) continue;
b9671574 2999 if (pt->GetRemoval()>10) continue;
1c53abe2 3000 for (Int_t j=i+1; j<nseed; j++){
3001 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1af5da7e 3002 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
91162307 3003 // if (pt2){
b9671574 3004 if (pt2->GetRemoval()<=10) {
6d493ea0 3005 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
91162307 3006 SignShared(pt,pt2);
c9427e08 3007 }
91162307 3008 }
3009 }
3010}
3011
91162307 3012
829455ad 3013void AliTPCtracker::SortTracks(TObjArray * arr, Int_t mode) const
91162307 3014{
3015 //
3016 //sort tracks in array according mode criteria
3017 Int_t nseed = arr->GetEntriesFast();
3018 for (Int_t i=0; i<nseed; i++) {
3019 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
3020 if (!pt) {
3021 continue;
3022 }
b9671574 3023 pt->SetSort(mode);
91162307 3024 }
3025 arr->UnSort();
3026 arr->Sort();
1c53abe2 3027}
c9427e08 3028
51ad6848 3029
829455ad 3030void AliTPCtracker::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
51ad6848 3031{
51ad6848 3032 //
6d493ea0 3033 // Loop over all tracks and remove overlaped tracks (with lower quality)
3034 // Algorithm:
3035 // 1. Unsign clusters
3036 // 2. Sort tracks according quality
3037 // Quality is defined by the number of cluster between first and last points
3038 //
3039 // 3. Loop over tracks - decreasing quality order
3040 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
3041 // b.) remove - If the minimal number of clusters less than minimal and not ITS
30