PWGPP-71 Move TPC noise filtering information form the AliESDevent to the AliESDHeader
[u/mrichter/AliRoot.git] / TPC / Rec / 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;
434cdbe1 475 fCrossTalkSignalArray = new TObjArray(nROCs*2); //
265b5e37 476 fCrossTalkSignalArray->SetOwner(kTRUE);
434cdbe1 477 for (Int_t isector=0; isector<2*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{
1348 //
434cdbe1 1349 // load clusters to the memory
1350 Int_t nROCs = 72;
bfa00fba 1351 static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
91162307 1352 //
1353 // TTree * tree = fClustersArray.GetTree();
9a836cc2 1354 AliInfo("LoadClusters()\n");
39240815 1355 // reset crosstalk matrix
1356 //
434cdbe1 1357 for (Int_t isector=0; isector<nROCs*2; isector++){ //set all ellemts of crosstalk matrix to 0
39240815 1358 TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1359 if (crossTalkMatrix)(*crossTalkMatrix)*=0;
1360 }
91162307 1361
1362 TTree * tree = fInput;
1363 TBranch * br = tree->GetBranch("Segment");
1364 br->SetAddress(&clrow);
bad6eb00 1365
1366 // Conversion of pad, row coordinates in local tracking coords.
1367 // Could be skipped here; is already done in clusterfinder
1368
6f8ff889 1369
91162307 1370 Int_t j=Int_t(tree->GetEntries());
1c53abe2 1371 for (Int_t i=0; i<j; i++) {
91162307 1372 br->GetEntry(i);
1373 //
1374 Int_t sec,row;
47af7ca4 1375 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
8cc1870a 1376
1377 // wire segmentID and nPadsPerSegment to be used for Xtalk calculation
1378 Int_t wireSegmentID = fkParam->GetWireSegment(sec,row);
1379 Float_t nPadsPerSegment = (Float_t)(fkParam->GetNPadsPerSegment(wireSegmentID));
1380 TMatrixD &crossTalkSignal = *((TMatrixD*)fCrossTalkSignalArray->At(sec));
434cdbe1 1381 TMatrixD &crossTalkSignalBelow = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs));
265b5e37 1382 Int_t nCols=crossTalkSignal.GetNcols();
3449a97c 1383 Double_t missingChargeFactor= AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrectionMissingCharge();
0201b65c 1384 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
8cc1870a 1385 AliTPCclusterMI *clXtalk= static_cast<AliTPCclusterMI*>(clrow->GetArray()->At(icl));
8cc1870a 1386 Transform((AliTPCclusterMI*)(clXtalk));
434cdbe1 1387 Int_t timeBinXtalk = clXtalk->GetTimeBin();
1388 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
1389 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
6f8ff889 1390 Double_t rmsTime2 = clXtalk->GetSigmaZ2()/(fkParam->GetZWidth()*fkParam->GetZWidth());
3449a97c 1391 Double_t rmsPad2 = clXtalk->GetSigmaY2()/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec));
434cdbe1 1392 if (rmsPadMin2>rmsPad2){
1393 rmsPad2=rmsPadMin2;
1394 }
1395 if (rmsTimeMin2>rmsTime2){
1396 rmsTime2=rmsTimeMin2;
1397 }
3449a97c 1398
6f8ff889 1399 Double_t norm= 2.*TMath::Exp(-1.0/(2.*rmsTime2))+2.*TMath::Exp(-4.0/(2.*rmsTime2))+1.;
434cdbe1 1400 Double_t qTotXtalk = 0.;
1401 Double_t qTotXtalkMissing = 0.;
6f8ff889 1402 for (Int_t itb=timeBinXtalk-2, idelta=-2; itb<=timeBinXtalk+2; itb++,idelta++) {
265b5e37 1403 if (itb<0 || itb>=nCols) continue;
3449a97c 1404 Double_t missingCharge=0;
1405 Double_t trf= TMath::Exp(-idelta*idelta/(2.*rmsTime2));
434cdbe1 1406 if (missingChargeFactor>0) {
1407 for (Int_t dpad=-2; dpad<=2; dpad++){
1408 Double_t qPad = clXtalk->GetMax()*TMath::Exp(-dpad*dpad/(2.*rmsPad2))*trf;
1409 qPad-=2.*crossTalkSignal[wireSegmentID][itb]; // missing part due crosttalk - feedback part - factor 2 used
1410 if (TMath::Nint(qPad)<=fkParam->GetZeroSup()){
1411 missingCharge+=qPad+2.*crossTalkSignal[wireSegmentID][itb];
1412 }else{
1413 missingCharge+=2.*crossTalkSignal[wireSegmentID][itb];
1414 }
3449a97c 1415 }
1416 }
434cdbe1 1417 qTotXtalk = clXtalk->GetQ()*trf/norm+missingCharge*missingChargeFactor;
1418 qTotXtalkMissing = missingCharge;
1419 crossTalkSignal[wireSegmentID][itb]+= qTotXtalk/nPadsPerSegment;
1420 crossTalkSignalBelow[wireSegmentID][itb]+= qTotXtalkMissing/nPadsPerSegment;
1421 }
0201b65c 1422 }
434cdbe1 1423
91162307 1424 //
bd26fa83 1425 AliTPCtrackerRow * tpcrow=0;
91162307 1426 Int_t left=0;
1427 if (sec<fkNIS*2){
1428 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1429 left = sec/fkNIS;
1430 }
1431 else{
1432 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1433 left = (sec-fkNIS*2)/fkNOS;
1434 }
1435 if (left ==0){
b9671574 1436 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
77f88633 1437 for (Int_t k=0;k<tpcrow->GetN1();++k)
1438 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
91162307 1439 }
1440 if (left ==1){
b9671574 1441 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
77f88633 1442 for (Int_t k=0;k<tpcrow->GetN2();++k)
bfa00fba 1443 tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
91162307 1444 }
1c53abe2 1445 }
65e67c9c 1446 if (AliTPCReconstructor::StreamLevel()&kStreamCrosstalkMatrix) {
434cdbe1 1447 //
1448 // dump the crosstalk matrices to tree for further investigation
1449 // a.) to estimate fluctuation of pedestal in indiviula wire segments
1450 // b.) to check correlation between regions
1451 // c.) to check relative conribution of signal below threshold to crosstalk
1452 for (Int_t isector=0; isector<nROCs; isector++){ //set all ellemts of crosstalk matrix to 0
1453 TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
1454 TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
1455 TVectorD vecAll(crossTalkMatrix->GetNrows());
1456 TVectorD vecBelow(crossTalkMatrix->GetNrows());
1457 //
1458 for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
1459 for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
1460 vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
1461 vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
1462 }
1463 (*fDebugStreamer)<<"crosstalkMatrix"<<
1464 "sector="<<isector<<
1465 "itime="<<itime<<
1466 "vecAll.="<<&vecAll<< // crosstalk charge + charge below threshold
1467 "vecBelow.="<<&vecBelow<< // crosstalk contribution from signal below threshold
1468 "\n";
1469 }
1470 }
1471 }
1472
1473
91162307 1474 //
bfa00fba 1475 clrow->Clear("C");
91162307 1476 LoadOuterSectors();
1477 LoadInnerSectors();
6f8ff889 1478
1479 cout << " =================================================================================================================================== " << endl;
1480 cout << " AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() << endl;
1481 cout << " AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() << endl;
1482 cout << " =================================================================================================================================== " << endl;
1483
832977ba 1484 if (AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection()) ApplyTailCancellation();
265b5e37 1485 if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) ApplyXtalkCorrection();
65e67c9c 1486 //if (AliTPCReconstructor::GetRecoParam()->GetUseOulierClusterFilter()) FilterOutlierClusters();
91162307 1487 return 0;
1c53abe2 1488}
1489
eea44233 1490void AliTPCtracker::FilterOutlierClusters(){
1491 //
1492 // filter outlier clusters
1493 //
1494 /*
1495 1.)..... booking part
1496 nSectors=72;
1497 nTimeBins=fParam->Get....
1498 TH2F hisTime("","", sector,0,sector, nTimeBins,0,nTimeBins);
1499 TH2F hisPadRow("","", sector,0,sector, nPadRows,0,nPadRows);
1500 2.) .... filling part
1501 .... cluster loop { hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin()); }
1502
1503 3.) ...filtering part
1504 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... }
1505
1506 sector loop
1507 { 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 }
1508 //
1509 4. Disabling clusters
1510
1511 */
1512
1513 //
1514 // 1.) booking part
e731a3f8 1515 //
3e1f1ce7 1516 // AliTPCcalibDB *db=AliTPCcalibDB::Instance();
e731a3f8 1517 Int_t nSectors=AliTPCROC::Instance()->GetNSectors();
3e1f1ce7 1518 Int_t nTimeBins= 1100; // *Bug here - we should get NTimeBins from ALTRO - Parameters not relyable
1519 Int_t nPadRows=(AliTPCROC::Instance()->GetNRows(0) + AliTPCROC::Instance()->GetNRows(36));
e731a3f8 1520 // parameters for filtering
1521 const Double_t nSigmaCut=9.; // should be in recoParam ?
1522 const Double_t offsetTime=100; // should be in RecoParam ? -
1523 const Double_t offsetPadRow=300; // should be in RecoParam ?
1524 const Double_t offsetTimeAccept=8; // should be in RecoParam ? - obtained as mean +1 rms in high IR pp
eea44233 1525 TH2F hisTime("hisSectorTime","hisSectorTime", nSectors,0,nSectors, nTimeBins,0,nTimeBins);
1526 TH2F hisPadRow("hisSectorRow","hisSectorRow", nSectors,0,nSectors, nPadRows,0,nPadRows);
1527 //
1528 // 2.) Filling part -- loop over clusters
1529 //
44ca7282 1530 for (Int_t isector=0; isector<36; isector++){ // loop over sectors
1531 for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
1532 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1533 Int_t nrows = sector.GetNRows();
1534 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1535 AliTPCtrackerRow& tpcrow = sector[row];
1536 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1537 if (iside>0) ncl=tpcrow.GetN2();
1538 for (Int_t i=0;i<ncl;i++) { // loop over clusters
1539 AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1540 hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin());
1541 hisPadRow.Fill(cluster->GetDetector(),cluster->GetRow());
1542 }
1543 }
1544 }
1545 }
44ca7282 1546
eea44233 1547 //
1548 // 3. Filtering part
1549 //
1550 TVectorD vecTime(nTimeBins);
44ca7282 1551 TVectorD vecPadRow(nPadRows);
c5dce091 1552 TVectorD vecMedianSectorTime(nSectors);
44ca7282 1553 TVectorD vecRMSSectorTime(nSectors);
c5dce091 1554 TVectorD vecMedianSectorTimeOut6(nSectors);
1555 TVectorD vecMedianSectorTimeOut9(nSectors);//
3e1f1ce7 1556 TVectorD vecMedianSectorTimeOut(nSectors);//
c5dce091 1557 TVectorD vecMedianSectorPadRow(nSectors);
44ca7282 1558 TVectorD vecRMSSectorPadRow(nSectors);
c5dce091 1559 TVectorD vecMedianSectorPadRowOut6(nSectors);
1560 TVectorD vecMedianSectorPadRowOut9(nSectors);
3e1f1ce7 1561 TVectorD vecMedianSectorPadRowOut(nSectors);
e731a3f8 1562 TVectorD vecSectorOut6(nSectors);
1563 TVectorD vecSectorOut9(nSectors);
3e1f1ce7 1564 TMatrixD matSectorCluster(nSectors,2);
e731a3f8 1565 //
1566 // 3.a) median, rms calculations for hisTime
1567 //
eea44233 1568 for (Int_t isec=0; isec<nSectors; isec++){
c5dce091 1569 vecMedianSectorTimeOut6[isec]=0;
1570 vecMedianSectorTimeOut9[isec]=0;
eea44233 1571 for (Int_t itime=0; itime<nTimeBins; itime++){
1572 vecTime[itime]=hisTime.GetBinContent(isec+1, itime+1);
1573 }
3e1f1ce7 1574 Double_t median= TMath::Mean(nTimeBins,vecTime.GetMatrixArray());
c5dce091 1575 Double_t rms= TMath::RMS(nTimeBins,vecTime.GetMatrixArray());
1576 vecMedianSectorTime[isec]=median;
44ca7282 1577 vecRMSSectorTime[isec]=rms;
3e1f1ce7 1578 if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector TimeStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
eea44233 1579 //
1580 // declare outliers
1581 for (Int_t itime=0; itime<nTimeBins; itime++){
1582 Double_t entries= hisTime.GetBinContent(isec+1, itime+1);
44ca7282 1583 if (entries>median+6.*rms+offsetTime) {
c5dce091 1584 vecMedianSectorTimeOut6[isec]+=1;
1585 }
44ca7282 1586 if (entries>median+9.*rms+offsetTime) {
c5dce091 1587 vecMedianSectorTimeOut9[isec]+=1;
eea44233 1588 }
1589 }
c5dce091 1590 }
e731a3f8 1591 //
1592 // 3.b) median, rms calculations for hisPadRow
1593 //
44ca7282 1594 for (Int_t isec=0; isec<nSectors; isec++){
1595 vecMedianSectorPadRowOut6[isec]=0;
1596 vecMedianSectorPadRowOut9[isec]=0;
1597 for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1598 vecPadRow[ipadrow]=hisPadRow.GetBinContent(isec+1, ipadrow+1);
1599 }
3e1f1ce7 1600 Int_t nPadRowsSector= AliTPCROC::Instance()->GetNRows(isec);
1601 Double_t median= TMath::Mean(nPadRowsSector,vecPadRow.GetMatrixArray());
1602 Double_t rms= TMath::RMS(nPadRowsSector,vecPadRow.GetMatrixArray());
44ca7282 1603 vecMedianSectorPadRow[isec]=median;
1604 vecRMSSectorPadRow[isec]=rms;
3e1f1ce7 1605 if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector PadRowStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
44ca7282 1606 //
1607 // declare outliers
1608 for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
1609 Double_t entries= hisPadRow.GetBinContent(isec+1, ipadrow+1);
1610 if (entries>median+6.*rms+offsetPadRow) {
1611 vecMedianSectorPadRowOut6[isec]+=1;
1612 }
1613 if (entries>median+9.*rms+offsetPadRow) {
1614 vecMedianSectorPadRowOut9[isec]+=1;
1615 }
1616 }
1617 }
1618 //
e731a3f8 1619 // 3.c) filter outlier sectors
1620 //
1621 Double_t medianSectorTime = TMath::Median(nSectors, vecTime.GetMatrixArray());
1622 Double_t mean69SectorTime, rms69SectorTime=0;
1623 AliMathBase::EvaluateUni(nSectors, vecTime.GetMatrixArray(), mean69SectorTime,rms69SectorTime,69);
1624 for (Int_t isec=0; isec<nSectors; isec++){
1625 vecSectorOut6[isec]=0;
1626 vecSectorOut9[isec]=0;
3e1f1ce7 1627 matSectorCluster(isec,0)=0;
1628 matSectorCluster(isec,1)=0;
65e67c9c 1629 if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+6.*(rms69SectorTime+ offsetTimeAccept))) {
1630 vecSectorOut6[isec]=1;
1631 }
1632 if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+9.*(rms69SectorTime+ offsetTimeAccept))){
1633 vecSectorOut9[isec]=1;
1634 }
e731a3f8 1635 }
1636 // light version of export variable
1637 Int_t filteredSector= vecSectorOut9.Sum(); // light version of export variable
1638 Int_t filteredSectorTime= vecMedianSectorTimeOut9.Sum();
1639 Int_t filteredSectorPadRow= vecMedianSectorPadRowOut9.Sum();
b9d88e03 1640 if (fEvent) if (fEvent->GetHeader()){
1641 fEvent->GetHeader()->SetTPCNoiseFilterCounter(0,TMath::Min(filteredSector,255));
1642 fEvent->GetHeader()->SetTPCNoiseFilterCounter(1,TMath::Min(filteredSectorTime,255));
1643 fEvent->GetHeader()->SetTPCNoiseFilterCounter(2,TMath::Min(filteredSectorPadRow,255));
eea44233 1644 }
65e67c9c 1645
44ca7282 1646 //
1647 // 4. Disabling clusters in outlier layers
1648 //
65e67c9c 1649 Int_t counterAll=0;
1650 Int_t counterOut=0;
44ca7282 1651 for (Int_t isector=0; isector<36; isector++){ // loop over sectors
1652 for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
1653 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1654 Int_t nrows = sector.GetNRows();
1655 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1656 AliTPCtrackerRow& tpcrow = sector[row];
1657 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1658 if (iside>0) ncl=tpcrow.GetN2();
1659 for (Int_t i=0;i<ncl;i++) { // loop over clusters
1660 AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1661 Double_t medianTime=vecMedianSectorTime[cluster->GetDetector()];
1662 Double_t medianPadRow=vecMedianSectorPadRow[cluster->GetDetector()];
1663 Double_t rmsTime=vecRMSSectorTime[cluster->GetDetector()];
1664 Double_t rmsPadRow=vecRMSSectorPadRow[cluster->GetDetector()];
1665 Int_t entriesPadRow=hisPadRow.GetBinContent(cluster->GetDetector()+1, cluster->GetRow()+1);
1666 Int_t entriesTime=hisTime.GetBinContent(cluster->GetDetector()+1, cluster->GetTimeBin()+1);
1667 Bool_t isOut=kFALSE;
65e67c9c 1668 if (vecSectorOut9[cluster->GetDetector()]>0.5) {
1669 isOut=kTRUE;
1670 }
e731a3f8 1671
3e1f1ce7 1672 if (entriesTime>medianTime+nSigmaCut*rmsTime+offsetTime) {
1673 isOut=kTRUE;
1674 vecMedianSectorTimeOut[cluster->GetDetector()]++;
1675 }
1676 if (entriesPadRow>medianPadRow+nSigmaCut*rmsPadRow+offsetPadRow) {
1677 isOut=kTRUE;
1678 vecMedianSectorPadRowOut[cluster->GetDetector()]++;
1679 }
65e67c9c 1680 counterAll++;
3e1f1ce7 1681 matSectorCluster(cluster->GetDetector(),0)+=1;
44ca7282 1682 if (isOut){
1683 cluster->Disable();
65e67c9c 1684 counterOut++;
3e1f1ce7 1685 matSectorCluster(cluster->GetDetector(),1)+=1;
44ca7282 1686 }
1687 }
1688 }
1689 }
1690 }
3e1f1ce7 1691 for (Int_t isec=0; isec<nSectors; isec++){
1692 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());
1693 }
65e67c9c 1694 //
1695 // dump info to streamer - for later tuning of cuts
1696 //
1697 if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) { // stream TPC data ouliers filtering infomation
3e1f1ce7 1698 AliLog::Flush();
1699 AliInfo(TString::Format("Cluster counter: (%d/%d) (Filtered/All)",counterOut,counterAll).Data());
1700 for (Int_t iSec=0; iSec<nSectors; iSec++){
1701 if (vecSectorOut9[iSec]>0 || matSectorCluster(iSec,1)>0) {
1702 AliInfo(TString::Format("Filtered sector\t%d",iSec).Data());
1703 Double_t vecMedTime =TMath::Median(72,vecMedianSectorTime.GetMatrixArray());
1704 Double_t vecMedPadRow =TMath::Median(72,vecMedianSectorPadRow.GetMatrixArray());
1705 Double_t vecMedCluster=(counterAll-counterOut)/72;
1706 AliInfo(TString::Format("VecMedianSectorTime\t(%4.4f/%4.4f/%4.4f)", vecMedianSectorTimeOut[iSec],vecMedianSectorTime[iSec],vecMedTime).Data());
1707 AliInfo(TString::Format("VecMedianSectorPadRow\t(%4.4f/%4.4f/%4.4f)", vecMedianSectorPadRowOut[iSec],vecMedianSectorPadRow[iSec],vecMedPadRow).Data());
1708 AliInfo(TString::Format("MatSectorCluster\t(%4.4f/%4.4f/%4.4f)\n", matSectorCluster(iSec,1), matSectorCluster(iSec,0), vecMedCluster).Data());
1709 AliLog::Flush();
1710 }
1711 }
1712 AliLog::Flush();
1713 Int_t eventNr = fEvent->GetEventNumberInFile();
65e67c9c 1714 (*fDebugStreamer)<<"filterClusterInfo"<<
1715 // minimal set variables for the ESDevent
3e1f1ce7 1716 "eventNr="<<eventNr<<
65e67c9c 1717 "counterAll="<<counterAll<<
1718 "counterOut="<<counterOut<<
1719 //
1720 "filteredSector="<<filteredSector<< // counter filtered sectors
1721 "filteredSectorTime="<<filteredSectorTime<< // counter filtered time bins
1722 "filteredSectorPadRow="<<filteredSectorPadRow<< // counter filtered pad-rows
1723 // per sector outlier information
1724 "medianSectorTime="<<medianSectorTime<< // median number of clusters per sector/timebin
1725 "mean69SectorTime="<<mean69SectorTime<< // LTM statistic mean of clusters per sector/timebin
1726 "rms69SectorTime="<<rms69SectorTime<< // LTM statistic RMS of clusters per sector/timebin
1727 "vecSectorOut6.="<<&vecSectorOut6<< // flag array sector - 6 sigma +accept margin outlier
1728 "vecSectorOut9.="<<&vecSectorOut9<< // flag array sector - 9 sigma + accept margin outlier
1729 // per sector/timebin outlier detection
1730 "vecMedianSectorTime.="<<&vecMedianSectorTime<<
1731 "vecRMSSectorTime.="<<&vecRMSSectorTime<<
1732 "vecMedianSectorTimeOut6.="<<&vecMedianSectorTimeOut6<<
1733 "vecMedianSectorTimeOut9.="<<&vecMedianSectorTimeOut9<<
3e1f1ce7 1734 "vecMedianSectorTimeOut0.="<<&vecMedianSectorTimeOut<<
65e67c9c 1735 // per sector/pad-row outlier detection
1736 "vecMedianSectorPadRow.="<<&vecMedianSectorPadRow<<
1737 "vecRMSSectorPadRow.="<<&vecRMSSectorPadRow<<
1738 "vecMedianSectorPadRowOut6.="<<&vecMedianSectorPadRowOut6<<
1739 "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut9<<
3e1f1ce7 1740 "vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut<<
65e67c9c 1741 "\n";
1742 ((*fDebugStreamer)<<"filterClusterInfo").GetTree()->Write();
1743 fDebugStreamer->GetFile()->Flush();
1744 }
eea44233 1745}
1746
829455ad 1747void AliTPCtracker::UnloadClusters()
91162307 1748{
1749 //
1750 // unload clusters from the memory
1751 //
1752 Int_t nrows = fOuterSec->GetNRows();
1753 for (Int_t sec = 0;sec<fkNOS;sec++)
1754 for (Int_t row = 0;row<nrows;row++){
bd26fa83 1755 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
982aff31 1756 // if (tpcrow){
1757 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1758 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1759 //}
1760 tpcrow->ResetClusters();
91162307 1761 }
1762 //
1763 nrows = fInnerSec->GetNRows();
1764 for (Int_t sec = 0;sec<fkNIS;sec++)
1765 for (Int_t row = 0;row<nrows;row++){
bd26fa83 1766 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
982aff31 1767 //if (tpcrow){
1768 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1769 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1770 //}
1771 tpcrow->ResetClusters();
91162307 1772 }
1773
1774 return ;
1775}
1776
829455ad 1777void AliTPCtracker::FillClusterArray(TObjArray* array) const{
002af263 1778 //
1779 // Filling cluster to the array - For visualization purposes
1780 //
1781 Int_t nrows=0;
1782 nrows = fOuterSec->GetNRows();
1783 for (Int_t sec = 0;sec<fkNOS;sec++)
1784 for (Int_t row = 0;row<nrows;row++){
1785 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1786 if (!tpcrow) continue;
1787 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1788 array->AddLast((TObject*)((*tpcrow)[icl]));
1789 }
1790 }
1791 nrows = fInnerSec->GetNRows();
1792 for (Int_t sec = 0;sec<fkNIS;sec++)
1793 for (Int_t row = 0;row<nrows;row++){
1794 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1795 if (!tpcrow) continue;
1796 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1797 array->AddLast((TObject*)(*tpcrow)[icl]);
1798 }
1799 }
1800}
1801
1802
829455ad 1803void AliTPCtracker::Transform(AliTPCclusterMI * cluster){
893468d9 1804 //
544c295f 1805 // transformation
893468d9 1806 //
3f3549a3 1807 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1808 AliTPCTransform *transform = calibDB->GetTransform() ;
24db6af7 1809 if (!transform) {
1810 AliFatal("Tranformations not in calibDB");
ec26e231 1811 return;
24db6af7 1812 }
239f39b1 1813 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2942f542 1814 Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),static_cast<Double_t>(cluster->GetPad()),static_cast<Double_t>(cluster->GetTimeBin())};
24db6af7 1815 Int_t i[1]={cluster->GetDetector()};
022ee144 1816 transform->Transform(x,i,0,1);
afbaa016 1817 // if (cluster->GetDetector()%36>17){
1818 // x[1]*=-1;
1819 //}
022ee144 1820
24db6af7 1821 //
1822 // in debug mode check the transformation
1823 //
65e67c9c 1824 if ((AliTPCReconstructor::StreamLevel()&kStreamTransform)>0) {
af86c1fd 1825 Float_t gx[3];
1826 cluster->GetGlobalXYZ(gx);
002af263 1827 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
24db6af7 1828 TTreeSRedirector &cstream = *fDebugStreamer;
65e67c9c 1829 cstream<<"Transform"<< // needed for debugging of the cluster transformation, resp. used for later visualization
002af263 1830 "event="<<event<<
24db6af7 1831 "x0="<<x[0]<<
1832 "x1="<<x[1]<<
1833 "x2="<<x[2]<<
af86c1fd 1834 "gx0="<<gx[0]<<
1835 "gx1="<<gx[1]<<
1836 "gx2="<<gx[2]<<
24db6af7 1837 "Cl.="<<cluster<<
1838 "\n";
1839 }
1840 cluster->SetX(x[0]);
1841 cluster->SetY(x[1]);
1842 cluster->SetZ(x[2]);
19b00333 1843 // The old stuff:
0201b65c 1844 //
1845 //
1846 //
47af7ca4 1847 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
3f3549a3 1848 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1849 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1850 //TGeoHMatrix mat;
1851 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1852 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1853 if (mat) mat->LocalToMaster(pos,posC);
1854 else{
1855 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1856 }
1857 cluster->SetX(posC[0]);
1858 cluster->SetY(posC[1]);
1859 cluster->SetZ(posC[2]);
a7760332 1860 }
0201b65c 1861}
1c53abe2 1862
265b5e37 1863void AliTPCtracker::ApplyXtalkCorrection(){
1864 //
1865 // ApplyXtalk correction
1866 // Loop over all clusters
1867 // add to each cluster signal corresponding to common Xtalk mode for given time bin at given wire segment
265b5e37 1868 // cluster loop
1869 for (Int_t isector=0; isector<36; isector++){ //loop tracking sectors
1870 for (Int_t iside=0; iside<2; iside++){ // loop over sides A/C
1871 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1872 Int_t nrows = sector.GetNRows();
1873 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1874 AliTPCtrackerRow& tpcrow = sector[row]; // row object
1875 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1876 if (iside>0) ncl=tpcrow.GetN2();
1877 Int_t xSector=0; // sector number in the TPC convention 0-72
1878 if (isector<18){ //if IROC
1879 xSector=isector+(iside>0)*18;
1880 }else{
1881 xSector=isector+18; // isector -18 +36
1882 if (iside>0) xSector+=18;
1883 }
1884 TMatrixD &crossTalkMatrix= *((TMatrixD*)fCrossTalkSignalArray->At(xSector));
1885 Int_t wireSegmentID = fkParam->GetWireSegment(xSector,row);
1886 for (Int_t i=0;i<ncl;i++) {
1887 AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1888 Int_t iTimeBin=TMath::Nint(cluster->GetTimeBin());
1889 Double_t xTalk= crossTalkMatrix[wireSegmentID][iTimeBin];
1890 cluster->SetMax(cluster->GetMax()+xTalk);
6f8ff889 1891 const Double_t kDummy=4;
265b5e37 1892 Double_t sumxTalk=xTalk*kDummy; // should be calculated via time response function
1893 cluster->SetQ(cluster->GetQ()+sumxTalk);
6f8ff889 1894
1895
65e67c9c 1896 if ((AliTPCReconstructor::StreamLevel()&kStreamXtalk)>0) { // flag: stream crosstalk correctio as applied to cluster
6f8ff889 1897 TTreeSRedirector &cstream = *fDebugStreamer;
1898 if (gRandom->Rndm() > 0.){
1899 cstream<<"Xtalk"<<
1900 "isector=" << isector << // sector [0,36]
1901 "iside=" << iside << // side A or C
1902 "row=" << row << // padrow
1903 "i=" << i << // index of the cluster
1904 "xSector=" << xSector << // sector [0,72]
1905 "wireSegmentID=" << wireSegmentID << // anode wire segment id [0,10]
1906 "iTimeBin=" << iTimeBin << // timebin of the corrected cluster
1907 "xTalk=" << xTalk << // Xtalk contribution added to Qmax
1908 "sumxTalk=" << sumxTalk << // Xtalk contribution added to Qtot (roughly 3*Xtalk)
1909 "cluster.=" << cluster << // corrected cluster object
1910 "\n";
1911 }
1912 }// dump the results to the debug streamer if in debug mode
1913
1914
1915
1916
1917
265b5e37 1918 }
1919 }
1920 }
1921 }
1922}
1923
9a836cc2 1924void AliTPCtracker::ApplyTailCancellation(){
7bc5f28b 1925 //
9a836cc2 1926 // Correct the cluster charge for the ion tail effect
7bc5f28b 1927 // The TimeResponse function accessed via AliTPCcalibDB (TPC/Calib/IonTail)
1928 //
7bc5f28b 1929
9a836cc2 1930 // Retrieve
1931 TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
1932 if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
1933 TObject *rocFactorIROC = ionTailArr->FindObject("factorIROC");
1934 TObject *rocFactorOROC = ionTailArr->FindObject("factorOROC");
1935 Float_t factorIROC = (atof(rocFactorIROC->GetTitle()));
1936 Float_t factorOROC = (atof(rocFactorOROC->GetTitle()));
1937
1938 // find the number of clusters for the whole TPC (nclALL)
1939 Int_t nclALL=0;
1940 for (Int_t isector=0; isector<36; isector++){
1941 AliTPCtrackerSector &sector= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
1942 nclALL += sector.GetNClInSector(0);
1943 nclALL += sector.GetNClInSector(1);
1944 }
1945
1946 // start looping over all clusters
1947 for (Int_t iside=0; iside<2; iside++){ // loop over sides
7bc5f28b 1948 //
1949 //
9a836cc2 1950 for (Int_t secType=0; secType<2; secType++){ //loop over inner or outer sector
1951 // cache experimantal tuning factor for the different chamber type
1952 const Float_t ampfactor = (secType==0)?factorIROC:factorOROC;
1953 std::cout << " ampfactor = " << ampfactor << std::endl;
7bc5f28b 1954 //
9a836cc2 1955 for (Int_t sec = 0;sec<fkNOS;sec++){ //loop overs sectors
1956 //
1957 //
1958 // Cache time response functions and their positons to COG of the cluster
1959 TGraphErrors ** graphRes = new TGraphErrors *[20];
1960 Float_t * indexAmpGraphs = new Float_t[20];
1961 for (Int_t icache=0; icache<20; icache++)
1962 {
1963 graphRes[icache] = NULL;
1964 indexAmpGraphs[icache] = 0;
1965 }
1966 ///////////////////////////// --> position fo sie loop
d4559772 1967 if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(sec+36*secType+18*iside,graphRes,indexAmpGraphs))
1968 {
1969 continue;
1970 }
9a836cc2 1971
1972 AliTPCtrackerSector &sector= (secType==0)?fInnerSec[sec]:fOuterSec[sec];
1973 Int_t nrows = sector.GetNRows(); // number of rows
1974 Int_t nclSector = sector.GetNClInSector(iside); // ncl per sector to be used for debugging
1975
1976 for (Int_t row = 0;row<nrows;row++){ // loop over rows
1977
1978 AliTPCtrackerRow& tpcrow = sector[row]; // row object
1979 Int_t ncl = tpcrow.GetN1(); // number of clusters in the row
1980 if (iside>0) ncl=tpcrow.GetN2();
1981
1982 // Order clusters in time for the proper correction of ion tail
1983 Float_t qTotArray[ncl]; // arrays to be filled with modified Qtot and Qmax values in order to avoid float->int conversion
1984 Float_t qMaxArray[ncl];
1985 Int_t sortedClusterIndex[ncl];
1986 Float_t sortedClusterTimeBin[ncl];
1987 TObjArray *rowClusterArray = new TObjArray(ncl); // cache clusters for each row
1988 for (Int_t i=0;i<ncl;i++)
1989 {
1990 qTotArray[i]=0;
1991 qMaxArray[i]=0;
1992 sortedClusterIndex[i]=i;
1993 AliTPCclusterMI *rowcl= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
1994 if (rowcl) {
1995 rowClusterArray->AddAt(rowcl,i);
1996 } else {
1997 rowClusterArray->RemoveAt(i);
1998 }
1999 // Fill the timebin info to the array in order to sort wrt tb
2000 if (!rowcl) {
2001 sortedClusterTimeBin[i]=0.0;
2002 } else {
2003 sortedClusterTimeBin[i] = rowcl->GetTimeBin();
2004 }
2005
2006 }
2007 TMath::Sort(ncl,sortedClusterTimeBin,sortedClusterIndex,kFALSE); // sort clusters in time
2008
2009 // Main cluster correction loops over clusters
2010 for (Int_t icl0=0; icl0<ncl;icl0++){ // first loop over clusters
2011
2012 AliTPCclusterMI *cl0= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl0]));
2013
2014 if (!cl0) continue;
2015 Int_t nclPad=0;
265b5e37 2016 for (Int_t icl1=0; icl1<ncl;icl1++){ // second loop over clusters
9a836cc2 2017 AliTPCclusterMI *cl1= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl1]));
2018 if (!cl1) continue;
2019 if (TMath::Abs(cl0->GetPad()-cl1->GetPad())>4) continue; // no contribution if far away in pad direction
2020 if (cl0->GetTimeBin()<= cl1->GetTimeBin()) continue; // no contibution to the tail if later
2021 if (TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin())>600) continue; // out of the range of response function
2022
2023 if (TMath::Abs(cl0->GetPad()-cl1->GetPad())<4) nclPad++; // count ncl for every pad for debugging
2024
2025 // Get the correction values for Qmax and Qtot and find total correction for a given cluster
2026 Double_t ionTailMax=0.;
2027 Double_t ionTailTotal=0.;
2028 GetTailValue(ampfactor,ionTailMax,ionTailTotal,graphRes,indexAmpGraphs,cl0,cl1);
2029 ionTailMax=TMath::Abs(ionTailMax);
2030 ionTailTotal=TMath::Abs(ionTailTotal);
2031 qTotArray[icl0]+=ionTailTotal;
2032 qMaxArray[icl0]+=ionTailMax;
2033
2034 // Dump some info for debugging while clusters are being corrected
65e67c9c 2035 if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) { // flag: stream ion tail correction as applied to cluster
9a836cc2 2036 TTreeSRedirector &cstream = *fDebugStreamer;
2037 if (gRandom->Rndm() > 0.999){
2038 cstream<<"IonTail"<<
2039 "cl0.=" <<cl0 << // cluster 0 (to be corrected)
2040 "cl1.=" <<cl1 << // cluster 1 (previous cluster)
2041 "ionTailTotal=" <<ionTailTotal << // ion Tail from cluster 1 contribution to cluster0
2042 "ionTailMax=" <<ionTailMax << // ion Tail from cluster 1 contribution to cluster0
2043 "\n";
2044 }
2045 }// dump the results to the debug streamer if in debug mode
2046
2047 }//end of second loop over clusters
2048
2049 // Set corrected values of the corrected cluster
2050 cl0->SetQ(TMath::Nint(Float_t(cl0->GetQ())+Float_t(qTotArray[icl0])));
2051 cl0->SetMax(TMath::Nint(Float_t(cl0->GetMax())+qMaxArray[icl0]));
2052
2053 // Dump some info for debugging after clusters are corrected
65e67c9c 2054 if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {
9a836cc2 2055 TTreeSRedirector &cstream = *fDebugStreamer;
2056 if (gRandom->Rndm() > 0.999){
2057 cstream<<"IonTailCorrected"<<
2058 "cl0.=" << cl0 << // cluster 0 with huge Qmax
2059 "ionTailTotalPerCluster=" << qTotArray[icl0] <<
2060 "ionTailMaxPerCluster=" << qMaxArray[icl0] <<
2061 "nclALL=" << nclALL <<
2062 "nclSector=" << nclSector <<
2063 "nclRow=" << ncl <<
2064 "nclPad=" << nclPad <<
2065 "row=" << row <<
2066 "sector=" << sec <<
2067 "icl0=" << icl0 <<
2068 "\n";
2069 }
2070 }// dump the results to the debug streamer if in debug mode
2071
2072 }//end of first loop over cluster
2073 delete rowClusterArray;
2074 }//end of loop over rows
2075 for (int i=0; i<20; i++) delete graphRes[i];
2076 delete [] graphRes;
2077 delete [] indexAmpGraphs;
2078
2079 }//end of loop over sectors
2080 }//end of loop over IROC/OROC
2081 }// end of side loop
7bc5f28b 2082}
9a836cc2 2083//_____________________________________________________________________________
5a516e0a 2084void AliTPCtracker::GetTailValue(Float_t ampfactor,Double_t &ionTailMax, Double_t &ionTailTotal,TGraphErrors **graphRes,Float_t *indexAmpGraphs,AliTPCclusterMI *cl0,AliTPCclusterMI *cl1){
9a836cc2 2085
2086 //
2087 // Function in order to calculate the amount of the correction to be added for a given cluster, return values are ionTailTaoltal and ionTailMax
265b5e37 2088 // Parameters:
2089 // cl0 - cluster to be modified
2090 // cl1 - source cluster ion tail of this cluster will be added to the cl0 (accroding time and pad response function)
9a836cc2 2091 //
9a836cc2 2092 const Double_t kMinPRF = 0.5; // minimal PRF width
2093 ionTailTotal = 0.; // correction value to be added to Qtot of cl0
2094 ionTailMax = 0.; // correction value to be added to Qmax of cl0
2095
2096 Float_t qTot0 = cl0->GetQ(); // cl0 Qtot info
2097 Float_t qTot1 = cl1->GetQ(); // cl1 Qtot info
2098 Int_t sectorPad = cl1->GetDetector(); // sector number
2099 Int_t padcl0 = TMath::Nint(cl0->GetPad()); // pad0
2100 Int_t padcl1 = TMath::Nint(cl1->GetPad()); // pad1
2101 Float_t padWidth = (sectorPad < 36)?0.4:0.6; // pad width in cm
2102 const Int_t deltaTimebin = TMath::Nint(TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin()))+12; //distance between pads of cl1 and cl0 increased by 12 bins
2103 Double_t rmsPad1 = (cl1->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl1->GetSigmaY2())/padWidth);
2104 Double_t rmsPad0 = (cl0->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl0->GetSigmaY2())/padWidth);
265b5e37 2105
2106
a3986da2 2107
9a836cc2 2108 Double_t sumAmp1=0.;
2109 for (Int_t idelta =-2; idelta<=2;idelta++){
2110 sumAmp1+=TMath::Exp(-idelta*idelta/(2*rmsPad1));
2111 }
7bc5f28b 2112
9a836cc2 2113 Double_t sumAmp0=0.;
2114 for (Int_t idelta =-2; idelta<=2;idelta++){
2115 sumAmp0+=TMath::Exp(-idelta*idelta/(2*rmsPad0));
2116 }
7bc5f28b 2117
9a836cc2 2118 // Apply the correction --> cl1 corrects cl0 (loop over cl1's pads and find which pads of cl0 are going to be corrected)
2119 Int_t padScan=2; // +-2 pad-timebin window will be scanned
2120 for (Int_t ipad1=padcl1-padScan; ipad1<=padcl1+padScan; ipad1++) {
2121 //
2122 //
2123 Float_t deltaPad1 = TMath::Abs(cl1->GetPad()-(Float_t)ipad1);
2124 Double_t amp1 = (TMath::Exp(-(deltaPad1*deltaPad1)/(2*rmsPad1)))/sumAmp1; // normalized pad response function
2125 Float_t qTotPad1 = amp1*qTot1; // used as a factor to multipliy the response function
2126
2127 // find closest value of cl1 to COG (among the time response functions' amplitude array --> to select proper t.r.f.)
2128 Int_t ampIndex = 0;
2129 Float_t diffAmp = TMath::Abs(deltaPad1-indexAmpGraphs[0]);
2130 for (Int_t j=0;j<20;j++) {
2131 if (diffAmp > TMath::Abs(deltaPad1-indexAmpGraphs[j]) && indexAmpGraphs[j]!=0)
2132 {
2133 diffAmp = TMath::Abs(deltaPad1-indexAmpGraphs[j]);
2134 ampIndex = j;
2135 }
2136 }
2137 if (!graphRes[ampIndex]) continue;
b4742ddd 2138 if (deltaTimebin+2 >= graphRes[ampIndex]->GetN()) continue;
9a836cc2 2139 if (graphRes[ampIndex]->GetY()[deltaTimebin+2]>=0) continue;
2140
2141 for (Int_t ipad0=padcl0-padScan; ipad0<=padcl0+padScan; ipad0++) {
2142 //
2143 //
2144 if (ipad1!=ipad0) continue; // check if ipad1 channel sees ipad0 channel, if not no correction to be applied.
2145
2146 Float_t deltaPad0 = TMath::Abs(cl0->GetPad()-(Float_t)ipad0);
2147 Double_t amp0 = (TMath::Exp(-(deltaPad0*deltaPad0)/(2*rmsPad0)))/sumAmp0; // normalized pad resp function
2148 Float_t qMaxPad0 = amp0*qTot0;
2149
2150 // Add 5 timebin range contribution around the max peak (-+2 tb window)
2151 for (Int_t itb=deltaTimebin-2; itb<=deltaTimebin+2; itb++) {
2152
2153 if (itb<0) continue;
2154 if (itb>=graphRes[ampIndex]->GetN()) continue;
2155
2156 // calculate contribution to qTot
2157 Float_t tailCorr = TMath::Abs((qTotPad1*ampfactor)*(graphRes[ampIndex])->GetY()[itb]);
2158 if (ipad1!=padcl0) {
2159 ionTailTotal += TMath::Min(qMaxPad0,tailCorr); // for side pad
2160 } else {
2161 ionTailTotal += tailCorr; // for center pad
2162 }
2163 // calculate contribution to qMax
2164 if (itb == deltaTimebin && ipad1 == padcl0) ionTailMax += tailCorr;
2165
2166 } // end of tb correction loop which is applied over 5 tb range
2167
2168 } // end of cl0 loop
2169 } // end of cl1 loop
2170
2171}
7bc5f28b 2172
1c53abe2 2173//_____________________________________________________________________________
829455ad 2174Int_t AliTPCtracker::LoadOuterSectors() {
1c53abe2 2175 //-----------------------------------------------------------------
2176 // This function fills outer TPC sectors with clusters.
2177 //-----------------------------------------------------------------
91162307 2178 Int_t nrows = fOuterSec->GetNRows();
2179 UInt_t index=0;
2180 for (Int_t sec = 0;sec<fkNOS;sec++)
2181 for (Int_t row = 0;row<nrows;row++){
bd26fa83 2182 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
91162307 2183 Int_t sec2 = sec+2*fkNIS;
2184 //left
b9671574 2185 Int_t ncl = tpcrow->GetN1();
91162307 2186 while (ncl--) {
b9671574 2187 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
91162307 2188 index=(((sec2<<8)+row)<<16)+ncl;
2189 tpcrow->InsertCluster(c,index);
2190 }
2191 //right
b9671574 2192 ncl = tpcrow->GetN2();
91162307 2193 while (ncl--) {
b9671574 2194 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
91162307 2195 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
2196 tpcrow->InsertCluster(c,index);
2197 }
2198 //
2199 // write indexes for fast acces
2200 //
2201 for (Int_t i=0;i<510;i++)
b9671574 2202 tpcrow->SetFastCluster(i,-1);
91162307 2203 for (Int_t i=0;i<tpcrow->GetN();i++){
2204 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
b9671574 2205 tpcrow->SetFastCluster(zi,i); // write index
91162307 2206 }
2207 Int_t last = 0;
2208 for (Int_t i=0;i<510;i++){
b9671574 2209 if (tpcrow->GetFastCluster(i)<0)
2210 tpcrow->SetFastCluster(i,last);
91162307 2211 else
b9671574 2212 last = tpcrow->GetFastCluster(i);
91162307 2213 }
2214 }
1c53abe2 2215 fN=fkNOS;
2216 fSectors=fOuterSec;
91162307 2217 return 0;
1c53abe2 2218}
2219
2220
2221//_____________________________________________________________________________
829455ad 2222Int_t AliTPCtracker::LoadInnerSectors() {
1c53abe2 2223 //-----------------------------------------------------------------
2224 // This function fills inner TPC sectors with clusters.
2225 //-----------------------------------------------------------------
91162307 2226 Int_t nrows = fInnerSec->GetNRows();
2227 UInt_t index=0;
2228 for (Int_t sec = 0;sec<fkNIS;sec++)
2229 for (Int_t row = 0;row<nrows;row++){
bd26fa83 2230 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
91162307 2231 //
2232 //left
b9671574 2233 Int_t ncl = tpcrow->GetN1();
91162307 2234 while (ncl--) {
b9671574 2235 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
91162307 2236 index=(((sec<<8)+row)<<16)+ncl;
2237 tpcrow->InsertCluster(c,index);
2238 }
2239 //right
b9671574 2240 ncl = tpcrow->GetN2();
91162307 2241 while (ncl--) {
b9671574 2242 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
91162307 2243 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
2244 tpcrow->InsertCluster(c,index);
2245 }
2246 //
2247 // write indexes for fast acces
2248 //
2249 for (Int_t i=0;i<510;i++)
b9671574 2250 tpcrow->SetFastCluster(i,-1);
91162307 2251 for (Int_t i=0;i<tpcrow->GetN();i++){
2252 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
b9671574 2253 tpcrow->SetFastCluster(zi,i); // write index
91162307 2254 }
2255 Int_t last = 0;
2256 for (Int_t i=0;i<510;i++){
b9671574 2257 if (tpcrow->GetFastCluster(i)<0)
2258 tpcrow->SetFastCluster(i,last);
91162307 2259 else
b9671574 2260 last = tpcrow->GetFastCluster(i);
91162307 2261 }
1c53abe2 2262
91162307 2263 }
2264
1c53abe2 2265 fN=fkNIS;
2266 fSectors=fInnerSec;
91162307 2267 return 0;
2268}
2269
2270
2271
2272//_________________________________________________________________________
829455ad 2273AliTPCclusterMI *AliTPCtracker::GetClusterMI(Int_t index) const {
91162307 2274 //--------------------------------------------------------------------
2275 // Return pointer to a given cluster
2276 //--------------------------------------------------------------------
e7eb17e4 2277 if (index<0) return 0; // no cluster
91162307 2278 Int_t sec=(index&0xff000000)>>24;
2279 Int_t row=(index&0x00ff0000)>>16;
d26d9159 2280 Int_t ncl=(index&0x00007fff)>>00;
91162307 2281
bd26fa83 2282 const AliTPCtrackerRow * tpcrow=0;
bfa00fba 2283 TClonesArray * clrow =0;
3da363a1 2284
ad23441a 2285 if (sec<0 || sec>=fkNIS*4) {
2286 AliWarning(Form("Wrong sector %d",sec));
2287 return 0x0;
2288 }
2289
91162307 2290 if (sec<fkNIS*2){
462fb6e0 2291 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
2292 if (tracksec.GetNRows()<=row) return 0;
2293 tpcrow = &(tracksec[row]);
3da363a1 2294 if (tpcrow==0) return 0;
2295
2296 if (sec<fkNIS) {
b9671574 2297 if (tpcrow->GetN1()<=ncl) return 0;
2298 clrow = tpcrow->GetClusters1();
3da363a1 2299 }
2300 else {
b9671574 2301 if (tpcrow->GetN2()<=ncl) return 0;
2302 clrow = tpcrow->GetClusters2();
3da363a1 2303 }
91162307 2304 }
3da363a1 2305 else {
462fb6e0 2306 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
2307 if (tracksec.GetNRows()<=row) return 0;
2308 tpcrow = &(tracksec[row]);
3da363a1 2309 if (tpcrow==0) return 0;
2310
2311 if (sec-2*fkNIS<fkNOS) {
b9671574 2312 if (tpcrow->GetN1()<=ncl) return 0;
2313 clrow = tpcrow->GetClusters1();
3da363a1 2314 }
2315 else {
b9671574 2316 if (tpcrow->GetN2()<=ncl) return 0;
2317 clrow = tpcrow->GetClusters2();
3da363a1 2318 }
91162307 2319 }
3da363a1 2320
bfa00fba 2321 return (AliTPCclusterMI*)clrow->At(ncl);
91162307 2322
1c53abe2 2323}
2324
91162307 2325
2326
829455ad 2327Int_t AliTPCtracker::FollowToNext(AliTPCseed& t, Int_t nr) {
1c53abe2 2328 //-----------------------------------------------------------------
2329 // This function tries to find a track prolongation to next pad row
2330 //-----------------------------------------------------------------
1c53abe2 2331 //
91162307 2332 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
76d56fd6 2333 //
2334 //
4d158c36 2335 AliTPCclusterMI *cl=0;
2336 Int_t tpcindex= t.GetClusterIndex2(nr);
2337 //
2338 // update current shape info every 5 pad-row
2339 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
2340 GetShape(&t,nr);
2341 //}
2342 //
2343 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
2344 //
2345 if (tpcindex==-1) return 0; //track in dead zone
f124f8bf 2346 if (tpcindex >= 0){ //
b9671574 2347 cl = t.GetClusterPointer(nr);
97d77e7a 2348 //if (cl==0) cl = GetClusterMI(tpcindex);
2349 if (!cl) cl = GetClusterMI(tpcindex);
b9671574 2350 t.SetCurrentClusterIndex1(tpcindex);
4d158c36 2351 }
2352 if (cl){
2353 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
2354 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
2355 //
2356 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
2357 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
2358
2359 if (TMath::Abs(angle-t.GetAlpha())>0.001){
2360 Double_t rotation = angle-t.GetAlpha();
b9671574 2361 t.SetRelativeSector(relativesector);
1791d824 2362 if (!t.Rotate(rotation)) {
2363 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2364 return 0;
2365 }
2366 }
2367 if (!t.PropagateTo(x)) {
2368 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
2369 return 0;
4d158c36 2370 }
4d158c36 2371 //
b9671574 2372 t.SetCurrentCluster(cl);
2373 t.SetRow(nr);
00055a22 2374 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
4d158c36 2375 if ((tpcindex&0x8000)==0) accept =0;
2376 if (accept<3) {
2377 //if founded cluster is acceptible
2378 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
b9671574 2379 t.SetErrorY2(t.GetErrorY2()+0.03);
2380 t.SetErrorZ2(t.GetErrorZ2()+0.03);
2381 t.SetErrorY2(t.GetErrorY2()*3);
2382 t.SetErrorZ2(t.GetErrorZ2()*3);
4d158c36 2383 }
b9671574 2384 t.SetNFoundable(t.GetNFoundable()+1);
4d158c36 2385 UpdateTrack(&t,accept);
2386 return 1;
f124f8bf 2387 }
2388 else { // Remove old cluster from track
2389 t.SetClusterIndex(nr, -3);
2390 t.SetClusterPointer(nr, 0);
2391 }
4d158c36 2392 }
1627d1c4 2393 }
3f82c4f2 2394 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
76d56fd6 2395 if (fIteration>1 && IsFindable(t)){
3f82c4f2 2396 // not look for new cluster during refitting
b9671574 2397 t.SetNFoundable(t.GetNFoundable()+1);
3f82c4f2 2398 return 0;
2399 }
91162307 2400 //
4d158c36 2401 UInt_t index=0;
ca142b1f 2402 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
f124f8bf 2403 if (!t.PropagateTo(x)) {
2404 if (fIteration==0) t.SetRemoval(10);
2405 return 0;
2406 }
2407 Double_t y = t.GetY();
91162307 2408 if (TMath::Abs(y)>ymax){
2409 if (y > ymax) {
b9671574 2410 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
91162307 2411 if (!t.Rotate(fSectors->GetAlpha()))
2412 return 0;
2413 } else if (y <-ymax) {
b9671574 2414 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
91162307 2415 if (!t.Rotate(-fSectors->GetAlpha()))
2416 return 0;
2417 }
f124f8bf 2418 if (!t.PropagateTo(x)) {
2419 if (fIteration==0) t.SetRemoval(10);
2420 return 0;
2421 }
2422 y = t.GetY();
91162307 2423 }
2424 //
4d158c36 2425 Double_t z=t.GetZ();
2426 //
a3232aae 2427
b9671574 2428 if (!IsActive(t.GetRelativeSector(),nr)) {
2429 t.SetInDead(kTRUE);
a3232aae 2430 t.SetClusterIndex2(nr,-1);
2431 return 0;
2432 }
2433 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
b9671574 2434 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
2435 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
a3232aae 2436
2437 if (!isActive || !isActive2) return 0;
2438
bd26fa83 2439 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
91162307 2440 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
2441 Double_t roady =1.;
2442 Double_t roadz = 1.;
2443 //
b9671574 2444 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2445 t.SetInDead(kTRUE);
91162307 2446 t.SetClusterIndex2(nr,-1);
1c53abe2 2447 return 0;
2448 }
2449 else
2450 {
76d56fd6 2451 if (IsFindable(t))
47af7ca4 2452 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
b9671574 2453 t.SetNFoundable(t.GetNFoundable()+1);
1627d1c4 2454 else
2455 return 0;
1c53abe2 2456 }
2457 //calculate
91162307 2458 if (krow) {
2459 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
2460 cl = krow.FindNearest2(y,z,roady,roadz,index);
b9671574 2461 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
91162307 2462 }
91162307 2463 if (cl) {
b9671574 2464 t.SetCurrentCluster(cl);
2465 t.SetRow(nr);
4d158c36 2466 if (fIteration==2&&cl->IsUsed(10)) return 0;
00055a22 2467 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
4d158c36 2468 if (fIteration==2&&cl->IsUsed(11)) {
b9671574 2469 t.SetErrorY2(t.GetErrorY2()+0.03);
2470 t.SetErrorZ2(t.GetErrorZ2()+0.03);
2471 t.SetErrorY2(t.GetErrorY2()*3);
2472 t.SetErrorZ2(t.GetErrorZ2()*3);
4d158c36 2473 }
d26d9159 2474 /*
91162307 2475 if (t.fCurrentCluster->IsUsed(10)){
2476 //
2477 //
c9427e08 2478
91162307 2479 t.fNShared++;
2480 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
2481 t.fRemoval =10;
2482 return 0;
2483 }
2484 }
d26d9159 2485 */
91162307 2486 if (accept<3) UpdateTrack(&t,accept);
c9427e08 2487
91162307 2488 } else {
b9671574 2489 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
91162307 2490
2491 }
2492 return 1;
2493}
c9427e08 2494
1c53abe2 2495
91162307 2496
5d837844 2497//_________________________________________________________________________
829455ad 2498Bool_t AliTPCtracker::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
5d837844 2499{
2500 // Get track space point by index
2501 // return false in case the cluster doesn't exist
2502 AliTPCclusterMI *cl = GetClusterMI(index);
2503 if (!cl) return kFALSE;
2504 Int_t sector = (index&0xff000000)>>24;
0201b65c 2505 // Int_t row = (index&0x00ff0000)>>16;
5d837844 2506 Float_t xyz[3];
47af7ca4 2507 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
0201b65c 2508 xyz[0] = cl->GetX();
5d837844 2509 xyz[1] = cl->GetY();
2510 xyz[2] = cl->GetZ();
2511 Float_t sin,cos;
47af7ca4 2512 fkParam->AdjustCosSin(sector,cos,sin);
5d837844 2513 Float_t x = cos*xyz[0]-sin*xyz[1];
2514 Float_t y = cos*xyz[1]+sin*xyz[0];
2515 Float_t cov[6];
2516 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
47af7ca4 2517 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
5d837844 2518 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
47af7ca4 2519 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
5d837844 2520 cov[0] = sin*sin*sigmaY2;
2521 cov[1] = -sin*cos*sigmaY2;
2522 cov[2] = 0.;
2523 cov[3] = cos*cos*sigmaY2;
2524 cov[4] = 0.;
2525 cov[5] = sigmaZ2;
2526 p.SetXYZ(x,y,xyz[2],cov);
ae079791 2527 AliGeomManager::ELayerID iLayer;
5d837844 2528 Int_t idet;
47af7ca4 2529 if (sector < fkParam->GetNInnerSector()) {
ae079791 2530 iLayer = AliGeomManager::kTPC1;
5d837844 2531 idet = sector;
2532 }
2533 else {
ae079791 2534 iLayer = AliGeomManager::kTPC2;
47af7ca4 2535 idet = sector - fkParam->GetNInnerSector();
5d837844 2536 }
ae079791 2537 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
5d837844 2538 p.SetVolumeID(volid);
2539 return kTRUE;
2540}
2541
2542
2543
829455ad 2544Int_t AliTPCtracker::UpdateClusters(AliTPCseed& t, Int_t nr) {
1c53abe2 2545 //-----------------------------------------------------------------
2546 // This function tries to find a track prolongation to next pad row
2547 //-----------------------------------------------------------------
b9671574 2548 t.SetCurrentCluster(0);
f124f8bf 2549 t.SetCurrentClusterIndex1(-3);
1c53abe2 2550
2551 Double_t xt=t.GetX();
91162307 2552 Int_t row = GetRowNumber(xt)-1;
2553 Double_t ymax= GetMaxY(nr);
2554
1c53abe2 2555 if (row < nr) return 1; // don't prolongate if not information until now -
ca142b1f 2556// if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
2557// t.fRemoval =10;
2558// return 0; // not prolongate strongly inclined tracks
2559// }
2560// if (TMath::Abs(t.GetSnp())>0.95) {
2561// t.fRemoval =10;
2562// return 0; // not prolongate strongly inclined tracks
2563// }// patch 28 fev 06
91162307 2564
2565 Double_t x= GetXrow(nr);
2566 Double_t y,z;
2567 //t.PropagateTo(x+0.02);
2568 //t.PropagateTo(x+0.01);
1627d1c4 2569 if (!t.PropagateTo(x)){
1627d1c4 2570 return 0;
2571 }
1c53abe2 2572 //
91162307 2573 y=t.GetY();
2574 z=t.GetZ();
1c53abe2 2575
91162307 2576 if (TMath::Abs(y)>ymax){
2577 if (y > ymax) {
b9671574 2578 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
91162307 2579 if (!t.Rotate(fSectors->GetAlpha()))
2580 return 0;
2581 } else if (y <-ymax) {
b9671574 2582 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
91162307 2583 if (!t.Rotate(-fSectors->GetAlpha()))
2584 return 0;
2585 }
982aff31 2586 // if (!t.PropagateTo(x)){
2587 // return 0;
2588 //}
2589 return 1;
2590 //y = t.GetY();
1c53abe2 2591 }
91162307 2592 //
3f82c4f2 2593 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
a3232aae 2594
b9671574 2595 if (!IsActive(t.GetRelativeSector(),nr)) {
2596 t.SetInDead(kTRUE);
a3232aae 2597 t.SetClusterIndex2(nr,-1);
2598 return 0;
2599 }
2600 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
2601
bd26fa83 2602 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1c53abe2 2603
b9671574 2604 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
2605 t.SetInDead(kTRUE);
91162307 2606 t.SetClusterIndex2(nr,-1);
1c53abe2 2607 return 0;
2608 }
2609 else
2610 {
76d56fd6 2611
2612 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
2613 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1627d1c4 2614 else
2615 return 0;
1c53abe2 2616 }
1c53abe2 2617
91162307 2618 // update current
42aec1f1 2619 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
91162307 2620 // t.fCurrentSigmaY = GetSigmaY(&t);
2621 //t.fCurrentSigmaZ = GetSigmaZ(&t);
2622 GetShape(&t,nr);
2623 }
1c53abe2 2624
91162307 2625 AliTPCclusterMI *cl=0;
d184f295 2626 Int_t index=0;
91162307 2627 //
2628 Double_t roady = 1.;
2629 Double_t roadz = 1.;
2630 //
d26d9159 2631
2632 if (!cl){
2633 index = t.GetClusterIndex2(nr);
bad6eb00 2634 if ( (index >= 0) && (index&0x8000)==0){
b9671574 2635 cl = t.GetClusterPointer(nr);
bad6eb00 2636 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
b9671574 2637 t.SetCurrentClusterIndex1(index);
d26d9159 2638 if (cl) {
b9671574 2639 t.SetCurrentCluster(cl);
d26d9159 2640 return 1;
2641 }
2642 }
2643 }
2644
3d172d79 2645 // if (index<0) return 0;
2646 UInt_t uindex = TMath::Abs(index);
d184f295 2647
91162307 2648 if (krow) {
d184f295 2649 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
2650 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
91162307 2651 }
d26d9159 2652
b9671574 2653 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
2654 t.SetCurrentCluster(cl);
d26d9159 2655
91162307 2656 return 1;
2657}
1c53abe2 2658
1c53abe2 2659
829455ad 2660Int_t AliTPCtracker::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
91162307 2661 //-----------------------------------------------------------------
2662 // This function tries to find a track prolongation to next pad row
2663 //-----------------------------------------------------------------
1c53abe2 2664
91162307 2665 //update error according neighborhoud
1c53abe2 2666
b9671574 2667 if (t.GetCurrentCluster()) {
2668 t.SetRow(nr);
00055a22 2669 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
91162307 2670
b9671574 2671 if (t.GetCurrentCluster()->IsUsed(10)){
91162307 2672 //
2673 //
2674 // t.fErrorZ2*=2;
2675 // t.fErrorY2*=2;
b9671574 2676 t.SetNShared(t.GetNShared()+1);
2677 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
2678 t.SetRemoval(10);
91162307 2679 return 0;
2680 }
b364ca79 2681 }
d26d9159 2682 if (fIteration>0) accept = 0;
b364ca79 2683 if (accept<3) UpdateTrack(&t,accept);
2684
1c53abe2 2685 } else {
91162307 2686 if (fIteration==0){
f124f8bf 2687 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
2688 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
91162307 2689
b9671574 2690 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1c53abe2 2691 }
2692 }
2693 return 1;
2694}
2695
2696
2697
91162307 2698//_____________________________________________________________________________
829455ad 2699Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1c53abe2 2700 //-----------------------------------------------------------------
91162307 2701 // This function tries to find a track prolongation.
1c53abe2 2702 //-----------------------------------------------------------------
2703 Double_t xt=t.GetX();
2704 //
f124f8bf 2705 Double_t alpha=t.GetAlpha();
1c53abe2 2706 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2707 if (alpha < 0. ) alpha += 2.*TMath::Pi();
91162307 2708 //
a3f36f42 2709 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1c53abe2 2710
f124f8bf 2711 Int_t first = GetRowNumber(xt);
2712 if (!fromSeeds)
2713 first -= step;
a3f36f42 2714 if (first < 0)
2715 first = 0;
51ad6848 2716 for (Int_t nr= first; nr>=rf; nr-=step) {
2717 // update kink info
2718 if (t.GetKinkIndexes()[0]>0){
2719 for (Int_t i=0;i<3;i++){
2720 Int_t index = t.GetKinkIndexes()[i];
2721 if (index==0) break;
2722 if (index<0) continue;
2723 //
6c94f330 2724 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
51ad6848 2725 if (!kink){
2726 printf("PROBLEM\n");
2727 }
2728 else{
eea478d3 2729 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
51ad6848 2730 if (kinkrow==nr){
2731 AliExternalTrackParam paramd(t);
2732 kink->SetDaughter(paramd);
eea478d3 2733 kink->SetStatus(2,5);
51ad6848 2734 kink->Update();
51ad6848 2735 }
2736 }
2737 }
2738 }
2739
2740 if (nr==80) t.UpdateReference();
982aff31 2741 if (nr<fInnerSec->GetNRows())
2742 fSectors = fInnerSec;
2743 else
2744 fSectors = fOuterSec;
91162307 2745 if (FollowToNext(t,nr)==0)
4d158c36 2746 if (!t.IsActive())
2747 return 0;
91162307 2748
2749 }
2750 return 1;
2751}
2752
1c53abe2 2753
1c53abe2 2754
91162307 2755
2756
2757
829455ad 2758Int_t AliTPCtracker::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
1c53abe2 2759 //-----------------------------------------------------------------
2760 // This function tries to find a track prolongation.
2761 //-----------------------------------------------------------------
1c53abe2 2762 //
eea478d3 2763 Double_t xt=t.GetX();
f124f8bf 2764 Double_t alpha=t.GetAlpha();
1c53abe2 2765 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2766 if (alpha < 0. ) alpha += 2.*TMath::Pi();
bad6eb00 2767 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1c53abe2 2768
b9671574 2769 Int_t first = t.GetFirstPoint();
f124f8bf 2770 Int_t ri = GetRowNumber(xt);
2771 if (!fromSeeds)
2772 ri += 1;
2773
2774 if (first<ri) first = ri;
91162307 2775 //
2776 if (first<0) first=0;
4d158c36 2777 for (Int_t nr=first; nr<=rf; nr++) {
ca142b1f 2778 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
51ad6848 2779 if (t.GetKinkIndexes()[0]<0){
2780 for (Int_t i=0;i<3;i++){
2781 Int_t index = t.GetKinkIndexes()[i];
2782 if (index==0) break;
2783 if (index>0) continue;
2784 index = TMath::Abs(index);
6c94f330 2785 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
51ad6848 2786 if (!kink){
2787 printf("PROBLEM\n");
2788 }
2789 else{
eea478d3 2790 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
51ad6848 2791 if (kinkrow==nr){
2792 AliExternalTrackParam paramm(t);
2793 kink->SetMother(paramm);
eea478d3 2794 kink->SetStatus(2,1);
51ad6848 2795 kink->Update();
51ad6848 2796 }
2797 }
eea478d3 2798 }
51ad6848 2799 }
eea478d3 2800 //
d26d9159 2801 if (nr<fInnerSec->GetNRows())
2802 fSectors = fInnerSec;
2803 else
2804 fSectors = fOuterSec;
c274e255 2805
d26d9159 2806 FollowToNext(t,nr);
1c53abe2 2807 }
2808 return 1;
2809}
2810
2811
2812
2813
2814
829455ad 2815Float_t AliTPCtracker::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1c53abe2 2816{
544c295f 2817 // overlapping factor
1c53abe2 2818 //
2819 sum1=0;
2820 sum2=0;
2821 Int_t sum=0;
1c53abe2 2822 //
2823 Float_t dz2 =(s1->GetZ() - s2->GetZ());
c9427e08 2824 dz2*=dz2;
91162307 2825
c9427e08 2826 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1c53abe2 2827 dy2*=dy2;
2828 Float_t distance = TMath::Sqrt(dz2+dy2);
c9427e08 2829 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1c53abe2 2830
91162307 2831 // Int_t offset =0;
b9671574 2832 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2833 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
c9427e08 2834 if (lastpoint>160)
2835 lastpoint =160;
2836 if (firstpoint<0)
2837 firstpoint = 0;
91162307 2838 if (firstpoint>lastpoint) {
2839 firstpoint =lastpoint;
2840 // lastpoint =160;
c9427e08 2841 }
2842
2843
91162307 2844 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2845 if (s1->GetClusterIndex2(i)>0) sum1++;
2846 if (s2->GetClusterIndex2(i)>0) sum2++;
2847 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1c53abe2 2848 sum++;
2849 }
2850 }
91162307 2851 if (sum<5) return 0;
2852
1627d1c4 2853 Float_t summin = TMath::Min(sum1+1,sum2+1);
2854 Float_t ratio = (sum+1)/Float_t(summin);
1c53abe2 2855 return ratio;
2856}
2857
829455ad 2858void AliTPCtracker::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1c53abe2 2859{
544c295f 2860 // shared clusters
1c53abe2 2861 //
a0f4d6a6 2862 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2863 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2864 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2865 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
91162307 2866
1c53abe2 2867 //
91162307 2868 Int_t sumshared=0;
2869 //
a0f4d6a6 2870 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2871 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2872 Int_t firstpoint = 0;
2873 Int_t lastpoint = 160;
91162307 2874 //
a0f4d6a6 2875 // if (firstpoint>=lastpoint-5) return;;
1af5da7e 2876
91162307 2877 for (Int_t i=firstpoint;i<lastpoint;i++){
2878 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2879 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2880 sumshared++;
2881 }
2882 }
a0f4d6a6 2883 if (sumshared>cutN0){
91162307 2884 // sign clusters
2885 //
2886 for (Int_t i=firstpoint;i<lastpoint;i++){
2887 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2888 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2889 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2890 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2891 if (s1->IsActive()&&s2->IsActive()){
b9671574 2892 p1->SetShared(kTRUE);
2893 p2->SetShared(kTRUE);
91162307 2894 }
2895 }
2896 }
2897 }
2898 //
a0f4d6a6 2899 if (sumshared>cutN0){
91162307 2900 for (Int_t i=0;i<4;i++){
b9671574 2901 if (s1->GetOverlapLabel(3*i)==0){
2902 s1->SetOverlapLabel(3*i, s2->GetLabel());
2903 s1->SetOverlapLabel(3*i+1,sumshared);
2904 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
91162307 2905 break;
2906 }
2907 }
2908 for (Int_t i=0;i<4;i++){
b9671574 2909 if (s2->GetOverlapLabel(3*i)==0){
2910 s2->SetOverlapLabel(3*i, s1->GetLabel());
2911 s2->SetOverlapLabel(3*i+1,sumshared);
2912 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
91162307 2913 break;
2914 }
2915 }
2916 }
91162307 2917}
1c53abe2 2918
829455ad 2919void AliTPCtracker::SignShared(TObjArray * arr)
91162307 2920{
1c53abe2 2921 //
91162307 2922 //sort trackss according sectors
2923 //
c9427e08 2924 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2925 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 2926 if (!pt) continue;
2927 //if (pt) RotateToLocal(pt);
b9671574 2928 pt->SetSort(0);
c9427e08 2929 }
91162307 2930 arr->UnSort();
6d493ea0 2931 arr->Sort(); // sorting according relative sectors
1c53abe2 2932 arr->Expand(arr->GetEntries());
91162307 2933 //
2934 //
1c53abe2 2935 Int_t nseed=arr->GetEntriesFast();
1c53abe2 2936 for (Int_t i=0; i<nseed; i++) {
2937 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
91162307 2938 if (!pt) continue;
ec26e231 2939 for (Int_t j=0;j<12;j++){
b9671574 2940 pt->SetOverlapLabel(j,0);
1c53abe2 2941 }
91162307 2942 }
2943 for (Int_t i=0; i<nseed; i++) {
2944 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2945 if (!pt) continue;
b9671574 2946 if (pt->GetRemoval()>10) continue;
1c53abe2 2947 for (Int_t j=i+1; j<nseed; j++){
2948 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1af5da7e 2949 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
91162307 2950 // if (pt2){
b9671574 2951 if (pt2->GetRemoval()<=10) {
6d493ea0 2952 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
91162307 2953 SignShared(pt,pt2);
c9427e08 2954 }
91162307 2955 }
2956 }
2957}
2958
91162307 2959
829455ad 2960void AliTPCtracker::SortTracks(TObjArray * arr, Int_t mode) const
91162307 2961{
2962 //
2963 //sort tracks in array according mode criteria
2964 Int_t nseed = arr->GetEntriesFast();
2965 for (Int_t i=0; i<nseed; i++) {
2966 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2967 if (!pt) {
2968 continue;
2969 }
b9671574 2970 pt->SetSort(mode);
91162307 2971 }
2972 arr->UnSort();
2973 arr->Sort();
1c53abe2 2974}
c9427e08 2975
51ad6848 2976
829455ad 2977void AliTPCtracker::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
51ad6848 2978{
51ad6848 2979 //
6d493ea0 2980 // Loop over all tracks and remove overlaped tracks (with lower quality)
2981 // Algorithm:
2982 // 1. Unsign clusters
2983 // 2. Sort tracks according quality
2984 // Quality is defined by the number of cluster between first and last points
2985 //
2986 // 3. Loop over tracks - decreasing quality order
2987 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2988 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2989 // c.) if track accepted - sign clusters
51ad6848 2990 //
829455ad 2991 //Called in - AliTPCtracker::Clusters2Tracks()
2992 // - AliTPCtracker::PropagateBack()
2993 // - AliTPCtracker::RefitInward()
6d493ea0 2994 //
be34cb88 2995 // Arguments:
2996 // factor1 - factor for constrained
2997 // factor2 - for non constrained tracks
2998 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2999 //
51ad6848 3000 UnsignClusters();
3001 //
3002 Int_t nseed = arr->GetEntriesFast();
3003 Float_t * quality = new Float_t[nseed];
3004 Int_t * indexes = new Int_t[nseed];
3005 Int_t good =0;
3006 //
3007 //
3008 for (Int_t i=0; i<nseed; i++) {
3009 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
3010 if (!pt){
3011 quality[i]=-1;
3012 continue;
3013 }
3014 pt->UpdatePoints(); //select first last max dens points
3015 Float_t * points = pt->GetPoints();
3016 if (points[3]<0.8) quality[i] =-1;
51ad6848 3017 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
6d493ea0 3018 //prefer high momenta tracks if overlaps
3019 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
51ad6848 3020 }
3021 TMath::Sort(nseed,quality,indexes);
3022 //
3023 //
3024 for (Int_t itrack=0; itrack<nseed; itrack++) {
3025 Int_t trackindex = indexes[itrack];
6d493ea0 3026 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
3027 if (!pt) continue;
3028 //
51ad6848 3029 if (quality[trackindex]<0){
ddfbc51a 3030 MarkSeedFree( arr->RemoveAt(trackindex) );
f06a1ff6 3031 continue;
51ad6848 3032 }
3033 //
6d493ea0 3034 //
51ad6848 3035