]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/TPCupgrade/AliToyMCReconstruction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TPC module
[u/mrichter/AliRoot.git] / TPC / TPCupgrade / AliToyMCReconstruction.cxx
CommitLineData
d1cf83f5 1
2#include <TDatabasePDG.h>
3#include <TString.h>
4a777885 4#include <TSystem.h>
5#include <TROOT.h>
6#include <TFile.h>
9e98dea8 7#include <TPRegexp.h>
36a79b0d 8#include <TVectorF.h>
d1cf83f5 9
10#include <AliExternalTrackParam.h>
11#include <AliTPCcalibDB.h>
12#include <AliTPCclusterMI.h>
223d9e38 13#include <AliTPCCorrection.h>
d1cf83f5 14#include <AliTrackerBase.h>
15#include <AliTrackPointArray.h>
16#include <AliLog.h>
17#include <AliTPCParam.h>
18#include <AliTPCROC.h>
19#include <TTreeStream.h>
50f68b1a 20#include <AliTPCReconstructor.h>
21#include <AliTPCTransform.h>
7f1a1b08 22#include <AliTPCseed.h>
510cfcff 23#include <AliTPCtracker.h>
24#include <AliTPCtrackerSector.h>
38d9d609 25#include <AliRieman.h>
d1cf83f5 26
27#include "AliToyMCTrack.h"
4a777885 28#include "AliToyMCEvent.h"
d1cf83f5 29
30#include "AliToyMCReconstruction.h"
31
05da1b4e 32/*
33
34
35
36*/
d1cf83f5 37
38//____________________________________________________________________________________
39AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
40, fSeedingRow(140)
41, fSeedingDist(10)
42, fClusterType(0)
43, fCorrectionType(kNoCorrection)
44, fDoTrackFit(kTRUE)
45, fUseMaterial(kFALSE)
600eaa0d 46, fIdealTracking(kFALSE)
a06336b6 47, fNmaxEvents(-1)
d1cf83f5 48, fTime0(-1)
600eaa0d 49, fCreateT0seed(kFALSE)
48e13991 50, fLongT0seed(kTRUE)
4388aa4a 51, fFillClusterRes(kFALSE)
a1d7f156 52, fUseT0list(kFALSE)
53, fUseZ0list(kFALSE)
36a79b0d 54, fForceAlpha(kFALSE)
70cc5429 55, fRecoInfo(-1)
d1cf83f5 56, fStreamer(0x0)
38d9d609 57, fInputFile(0x0)
d1cf83f5 58, fTree(0x0)
4a777885 59, fEvent(0x0)
d1cf83f5 60, fTPCParam(0x0)
223d9e38 61, fTPCCorrection(0x0)
5993ed4f 62, fkNSectorInner(18) // hard-coded to avoid loading the parameters before
63, fInnerSectorArray(0x0)
64, fkNSectorOuter(18) // hard-coded to avoid loading the parameters before
65, fOuterSectorArray(0x0)
38d9d609 66, fAllClusters("AliTPCclusterMI",10000)
67, fMapTrackEvent(10000)
68, fMapTrackTrackInEvent(10000)
5414766f 69, fHnDelta(0x0)
38d9d609 70, fIsAC(kFALSE)
d1cf83f5 71{
72 //
73 // ctor
74 //
75 fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
5993ed4f 76
d1cf83f5 77}
78
79//____________________________________________________________________________________
80AliToyMCReconstruction::~AliToyMCReconstruction()
81{
82 //
83 // dtor
84 //
85
38d9d609 86 Cleanup();
4a777885 87}
88
89//____________________________________________________________________________________
90void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
91{
92 //
510cfcff 93 // Recostruction from associated clusters
4a777885 94 //
95
a06336b6 96 ConnectInputFile(file, nmaxEv);
38d9d609 97 if (!fTree) return;
a06336b6 98
36a79b0d 99 Int_t maxev=fTree->GetEntries();
100 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
101
a06336b6 102 InitStreamer(".debug");
4a777885 103
104 gROOT->cd();
105
e83fd282 106 static AliExternalTrackParam resetParam;
4a777885 107
108 AliExternalTrackParam t0seed;
109 AliExternalTrackParam seed;
110 AliExternalTrackParam track;
111 AliExternalTrackParam tOrig;
112
e83fd282 113 // at ITS
114 AliExternalTrackParam tOrigITS; // ideal track
115 AliExternalTrackParam tRealITS; // ITS track with realistic space point resolution
116 AliExternalTrackParam trackITS; // TPC refitted track
117
118 //between TPC inner wall and ITS
119 AliExternalTrackParam tOrigITS1;
120 AliExternalTrackParam tRealITS1;
121 AliExternalTrackParam trackITS1;
122
123 //at TPC inner wall
124 AliExternalTrackParam tOrigITS2;
125 AliExternalTrackParam tRealITS2;
126 AliExternalTrackParam trackITS2;
127
4a777885 128 AliExternalTrackParam *dummy;
4388aa4a 129
36a79b0d 130 // prepare list of T0s
f45d00f9 131 TVectorF t0list(maxev);
132 TVectorF z0list(maxev);
133 if (fUseT0list || fUseZ0list) {
36a79b0d 134 for (Int_t iev=0; iev<maxev; ++iev){
135 fTree->GetEvent(iev);
136 const Float_t t0=fEvent->GetT0();
f45d00f9 137 const Float_t z0=fEvent->GetZ();
36a79b0d 138 t0list[iev]=t0;
f45d00f9 139 z0list[iev]=z0;
36a79b0d 140 }
141 }
142
143 // array with cluster residuals
4388aa4a 144 TClonesArray *arrClustRes=0x0;
145 if (fFillClusterRes){
146 arrClustRes=new TClonesArray("AliTPCclusterMI",160);
147 }
4a777885 148
e83fd282 149 const Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
150 const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
151 const Double_t betweeTPCITS = (lastLayerITS+iFCRadius)/2.; // its track propgated to inner TPC wall
152
153 const Double_t kMaxSnp = 0.85;
154 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
155
36a79b0d 156 Double_t lastT0=0;
5414766f 157
158 // residuals
159 // binning r, phi, z, delta
160 const Int_t nbins=4;
161 Int_t bins[nbins] = {16, 18*5, 50, 80};
162 Double_t xmin[nbins] = {86. , 0., -250., -2.};
163 Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250., 2.};
164 fHnDelta = new THnF("hn", "hn", nbins, bins, xmin, xmax);
3abb3e1e 165
166 // fill streamer?
167 Bool_t fillStreamer=(fStreamer!=0x0);
168 if (fRecoInfo>-1 && ((fRecoInfo&kFillNoTrackInfo)==kFillNoTrackInfo)) fillStreamer=kFALSE;
4a777885 169
170 for (Int_t iev=0; iev<maxev; ++iev){
171 printf("============== Processing Event %6d =================\n",iev);
172 fTree->GetEvent(iev);
36a79b0d 173
174 Float_t z0=fEvent->GetZ();
175 Float_t t0=fEvent->GetT0();
862220e2 176
177 // set SC scaling factor
178 fTPCCorrection->SetCorrScaleFactor(fEvent->GetSCscale());
36a79b0d 179
4a777885 180 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
0403120d 181// printf(" > ====== Processing Track %6d ======== \n",itr);
4a777885 182 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
183 tOrig = *tr;
70cc5429 184
e83fd282 185 // propagate original track to ITS comparison points
70cc5429 186 if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ) {
187 tOrigITS = *tr;
188 AliTrackerBase::PropagateTrackTo(&tOrigITS, lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
189 tRealITS = resetParam;
190 }
191 if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ) {
192 tOrigITS1 = *tr;
193 AliTrackerBase::PropagateTrackTo(&tOrigITS1,betweeTPCITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
194 tRealITS1 = resetParam;
195 }
196 if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ) {
197 tOrigITS2 = *tr;
198 AliTrackerBase::PropagateTrackTo(&tOrigITS2,iFCRadius, kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
199 tRealITS2 = resetParam;
200 }
e83fd282 201
202 // realistic ITS track propagated to reference points
e83fd282 203 dummy = GetTrackRefit(tr,kITS);
204 if (dummy){
e83fd282 205 // propagate realistic track to ITS comparison points
70cc5429 206 if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ) {
207 tRealITS = *dummy;
208 AliTrackerBase::PropagateTrackTo(&tRealITS, lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
209 }
210 if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ) {
211 tRealITS1 = *dummy;
212 AliTrackerBase::PropagateTrackTo(&tRealITS1,betweeTPCITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
213 }
214 if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ) {
215 tRealITS2 = *dummy;
216 AliTrackerBase::PropagateTrackTo(&tRealITS2,iFCRadius, kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
217 }
e83fd282 218 //
219 delete dummy;
220 dummy=0x0;
221 }
222
4a777885 223
224
e83fd282 225 // resetParam
226 t0seed = resetParam;
227 seed = resetParam;
228 track = resetParam;
70cc5429 229
230 if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ) trackITS = resetParam;
231 if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ) trackITS1 = resetParam;
232 if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ) trackITS2 = resetParam;
4a777885 233
38d9d609 234 Float_t vDrift=GetVDrift();
4a777885 235 Float_t zLength=GetZLength(0);
236
600eaa0d 237 // crate time0 seed, steered by fCreateT0seed
05da1b4e 238// printf("t0 seed\n");
4a777885 239 fTime0=-1.;
600eaa0d 240 fCreateT0seed=kTRUE;
4a777885 241 dummy = GetSeedFromTrack(tr);
242
243 if (dummy) {
244 t0seed = *dummy;
245 delete dummy;
246
8a2a61d3 247 // Long seed
48e13991 248 if (fLongT0seed){
249 dummy = GetFittedTrackFromSeed(tr,&t0seed);
250 t0seed = *dummy;
251 delete dummy;
252 }
36a79b0d 253
254 // set the T0 from the seed
255 // in case the match with the real T0 infor is requested, find the
256 // closes T0 from the list of T0s
38d9d609 257 fTime0 = t0seed.GetZ()-zLength/vDrift;
f45d00f9 258 if (fUseT0list || fUseZ0list) {
259 fTime0 = FindClosestT0(t0list, z0list, t0seed);
36a79b0d 260 }
261 // create real seed using the time 0 from the first seed
262 // set fCreateT0seed now to false to get the seed in z coordinates
600eaa0d 263 fCreateT0seed = kFALSE;
05da1b4e 264// printf("seed (%.2g)\n",fTime0);
4a777885 265 dummy = GetSeedFromTrack(tr);
266 if (dummy) {
267 seed = *dummy;
268 delete dummy;
70cc5429 269 dummy=0x0;
4a777885 270
271 // create fitted track
272 if (fDoTrackFit){
05da1b4e 273// printf("track\n");
4388aa4a 274 dummy = GetFittedTrackFromSeed(tr, &seed, arrClustRes);
4a777885 275 track = *dummy;
70cc5429 276 dummy=0x0;
4a777885 277 delete dummy;
278 }
e83fd282 279
600eaa0d 280 // propagate seed to 0
05da1b4e 281 AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
d5704a39 282
e83fd282 283 //
284 // ITS comparison
285 //
286
e83fd282 287 // rotate fitted track to the frame of the original track and propagate to same reference
70cc5429 288 if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ){
289 trackITS = track;
290 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
291 trackITS.Rotate(tOrigITS.GetAlpha());
292 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
293 }
e83fd282 294
295 // rotate fitted track to the frame of the original track and propagate to same reference
70cc5429 296 if (fRecoInfo<0 || (fRecoInfo&kFillITS1)==kFillITS1 ){
297 trackITS1 = track;
298 AliTrackerBase::PropagateTrackTo(&trackITS1,betweeTPCITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
299 trackITS1.Rotate(tOrigITS1.GetAlpha());
300 AliTrackerBase::PropagateTrackTo(&trackITS1,betweeTPCITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
301 }
8a2a61d3 302
303 // rotate fitted track to the frame of the original track and propagate to same reference
70cc5429 304 if (fRecoInfo<0 || (fRecoInfo&kFillITS2)==kFillITS2 ){
305 trackITS2 = track;
306 AliTrackerBase::PropagateTrackTo(&trackITS2,iFCRadius,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
307 trackITS2.Rotate(tOrigITS2.GetAlpha());
308 AliTrackerBase::PropagateTrackTo(&trackITS2,iFCRadius,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
309 }
4a777885 310 }
311 }
312
313 Int_t ctype(fCorrectionType);
70cc5429 314
3abb3e1e 315 if (fillStreamer){
4a777885 316 (*fStreamer) << "Tracks" <<
317 "iev=" << iev <<
318 "z0=" << z0 <<
319 "t0=" << t0 <<
36a79b0d 320 "lastt0=" << lastT0 <<
600eaa0d 321 "fTime0=" << fTime0 <<
4a777885 322 "itr=" << itr <<
323 "clsType=" << fClusterType <<
324 "corrType=" << ctype <<
325 "seedRow=" << fSeedingRow <<
326 "seedDist=" << fSeedingDist <<
38d9d609 327 "vDrift=" << vDrift <<
4a777885 328 "zLength=" << zLength <<
329 "t0seed.=" << &t0seed <<
330 "seed.=" << &seed <<
e83fd282 331
4a777885 332 "tOrig.=" << &tOrig <<
70cc5429 333 "track.=" << &track;
334
838d33a8 335
70cc5429 336 // ITS match
337 if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ){
338 (*fStreamer) << "Tracks" <<
339 "tOrigITS.=" << &tOrigITS <<
340 "tRealITS.=" << &tRealITS <<
341 "trackITS.=" << &trackITS;
342 }
838d33a8 343
70cc5429 344 if (fRecoInfo<0 || (fRecoInfo&kFillITS1) ==kFillITS1 ){
345 (*fStreamer) << "Tracks" <<
346 "tOrigITS1.=" << &tOrigITS1 <<
347 "tRealITS1.=" << &tRealITS1 <<
348 "trackITS1.=" << &trackITS1;
349 }
350
351 if (fRecoInfo<0 || (fRecoInfo&kFillITS) ==kFillITS ){
352 (*fStreamer) << "Tracks" <<
353 "tOrigITS2.=" << &tOrigITS2 <<
354 "tRealITS2.=" << &tRealITS2 <<
355 "trackITS2.=" << &trackITS2;
356 }
4388aa4a 357
358 if (arrClustRes) {
359 const Int_t nCl=arrClustRes->GetEntriesFast();
360 // fracktion of outliers from track extrapolation
12c02f0f 361 // for 3, 3.5, 4, 4.5 and 5 sigma of the cluster resolution (~1mm)
4388aa4a 362 Float_t fracY[5]={0.};
363 Float_t fracZ[5]={0.};
364
365 for (Int_t icl=0; icl<nCl; ++icl) {
366 AliTPCclusterMI *cl=static_cast<AliTPCclusterMI*>(arrClustRes->At(icl));
b7106990 367// const Float_t sigmaY=TMath::Sqrt(cl->GetSigmaY2());
368// const Float_t sigmaZ=TMath::Sqrt(cl->GetSigmaZ2());
4388aa4a 369 for (Int_t inSig=0; inSig<5; ++inSig) {
b7106990 370 fracY[inSig] += cl->GetY()>(3+inSig*.5)/**sigmaY*/;
371 fracZ[inSig] += cl->GetZ()>(3+inSig*.5)/**sigmaZ*/;
4388aa4a 372 }
373 }
374
375 if (nCl>0) {
376 for (Int_t inSig=0; inSig<5; ++inSig) {
377 fracY[inSig]/=nCl;
378 fracZ[inSig]/=nCl;
379 }
380 }
381
382 (*fStreamer) << "Tracks" <<
383 "clustRes.=" << arrClustRes;
384 for (Int_t inSig=0; inSig<5; ++inSig) {
12c02f0f 385 const char* fracYname=Form("clFracY%02d=", 30+inSig*5);
386 const char* fracZname=Form("clFracZ%02d=", 30+inSig*5);
4388aa4a 387 (*fStreamer) << "Tracks" <<
388 fracYname << fracY[inSig] <<
389 fracZname << fracZ[inSig];
390 }
391 }
392
393 (*fStreamer) << "Tracks" <<
4a777885 394 "\n";
395 }
396
397
398 }
36a79b0d 399 lastT0=t0;
4a777885 400 }
401
5414766f 402 fStreamer->GetFile()->cd();
403 fHnDelta->Write();
404
4388aa4a 405 delete arrClustRes;
38d9d609 406 Cleanup();
d1cf83f5 407}
408
510cfcff 409
410//____________________________________________________________________________________
411void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
412{
413 //
5993ed4f 414 // Reconstruction for seed from associated clusters, but array of clusters:
415 // Step 1) Filling of cluster arrays
416 // Step 2) Seeding from clusters associated to tracks
417 // Step 3) Free track reconstruction using all clusters
510cfcff 418 //
419
420 TFile f(file);
421 if (!f.IsOpen() || f.IsZombie()) {
422 printf("ERROR: couldn't open the file '%s'\n", file);
423 return;
424 }
425
426 fTree=(TTree*)f.Get("toyMCtree");
427 if (!fTree) {
428 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
429 return;
430 }
431
432 fEvent=0x0;
433 fTree->SetBranchAddress("event",&fEvent);
434
435 // read spacecharge from the Userinfo ot the tree
436 InitSpaceCharge();
437
438 TString debugName=file;
439 debugName.ReplaceAll(".root","");
440 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
441 fUseMaterial,fIdealTracking,fClusterType,
442 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
443 debugName.Append(".allClusters.debug.root");
444
445 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
446 if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
447
448 gROOT->cd();
449
450 static AliExternalTrackParam dummySeedT0;
451 static AliExternalTrackParam dummySeed;
452 static AliExternalTrackParam dummyTrack;
453
454 AliExternalTrackParam t0seed;
455 AliExternalTrackParam seed;
5993ed4f 456 AliExternalTrackParam track;
510cfcff 457 AliExternalTrackParam tOrig;
458
510cfcff 459 AliExternalTrackParam *dummy;
460
461 Int_t maxev=fTree->GetEntries();
462 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
463
5993ed4f 464 // ===========================================================================================
465 // Loop 1: Fill AliTPCtrackerSector structure
466 // ===========================================================================================
c0ebc0e0 467 FillSectorStructure(maxev);
5993ed4f 468
c0ebc0e0 469 // settings (TODO: find the correct settings)
470 AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
471 tpcRecoParam->SetDoKinks(kFALSE);
472 AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
473 //tpcRecoParam->Print();
5993ed4f 474
c0ebc0e0 475 // need AliTPCReconstructor for parameter settings in AliTPCtracker
476 AliTPCReconstructor *tpcRec = new AliTPCReconstructor();
477 tpcRec->SetRecoParam(tpcRecoParam);
5993ed4f 478
479
480 // ===========================================================================================
481 // Loop 2: Seeding from clusters associated to tracks
482 // TODO: Implement tracking from given seed!
483 // ===========================================================================================
484 for (Int_t iev=0; iev<maxev; ++iev){
485 printf("============== Processing Event %6d =================\n",iev);
486 fTree->GetEvent(iev);
487 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
488 printf(" > ====== Processing Track %6d ======== \n",itr);
489 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
490 tOrig = *tr;
491
492
493 // set dummy
494 t0seed = dummySeedT0;
495 seed = dummySeed;
496 track = dummyTrack;
497
498 Float_t z0=fEvent->GetZ();
499 Float_t t0=fEvent->GetT0();
500
38d9d609 501 Float_t vDrift=GetVDrift();
5993ed4f 502 Float_t zLength=GetZLength(0);
503
c0ebc0e0 504 Int_t nClus = 0;
505
5993ed4f 506 // crate time0 seed, steered by fCreateT0seed
507 printf("t0 seed\n");
508 fTime0=-1.;
509 fCreateT0seed=kTRUE;
510 dummy = GetSeedFromTrack(tr);
511
512 if (dummy) {
513 t0seed = *dummy;
514 delete dummy;
515
516 // crate real seed using the time 0 from the first seed
517 // set fCreateT0seed now to false to get the seed in z coordinates
38d9d609 518 fTime0 = t0seed.GetZ()-zLength/vDrift;
5993ed4f 519 fCreateT0seed = kFALSE;
520 printf("seed (%.2g)\n",fTime0);
521 dummy = GetSeedFromTrack(tr);
522 if (dummy) {
523 seed = *dummy;
524 delete dummy;
525
526 // create fitted track
527 if (fDoTrackFit){
528 printf("track\n");
c0ebc0e0 529 dummy = GetFittedTrackFromSeedAllClusters(tr, &seed,nClus);
5993ed4f 530 track = *dummy;
531 delete dummy;
532 }
533
534 // propagate seed to 0
535 const Double_t kMaxSnp = 0.85;
536 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
537 AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
538
539 }
540 }
38d9d609 541
5993ed4f 542 Int_t ctype(fCorrectionType);
543
544 if (fStreamer) {
545 (*fStreamer) << "Tracks" <<
546 "iev=" << iev <<
547 "z0=" << z0 <<
548 "t0=" << t0 <<
549 "fTime0=" << fTime0 <<
550 "itr=" << itr <<
551 "clsType=" << fClusterType <<
552 "corrType=" << ctype <<
553 "seedRow=" << fSeedingRow <<
554 "seedDist=" << fSeedingDist <<
38d9d609 555 "vDrift=" << vDrift <<
5993ed4f 556 "zLength=" << zLength <<
c0ebc0e0 557 "nClus=" << nClus <<
5993ed4f 558 "t0seed.=" << &t0seed <<
559 "seed.=" << &seed <<
560 "track.=" << &track <<
561 "tOrig.=" << &tOrig <<
562 "\n";
563 }
564
565
566 }
567 }
568
569
570 delete fStreamer;
571 fStreamer=0x0;
572
573 delete fEvent;
574 fEvent = 0x0;
575
576 delete fTree;
577 fTree=0x0;
578 f.Close();
579}
580
581//____________________________________________________________________________________
582void AliToyMCReconstruction::RunRecoAllClustersStandardTracking(const char* file, Int_t nmaxEv)
583{
584 //
585 // Reconstruction for seed from associated clusters, but array of clusters
586 // Step 1) Filling of cluster arrays
587 // Step 2) Use the standard tracking: AliTPCtracker::Clusters2Tracks();
588 //
589
590 TFile f(file);
591 if (!f.IsOpen() || f.IsZombie()) {
592 printf("ERROR: couldn't open the file '%s'\n", file);
593 return;
594 }
595
596 fTree=(TTree*)f.Get("toyMCtree");
597 if (!fTree) {
598 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
599 return;
600 }
601
602 fEvent=0x0;
603 fTree->SetBranchAddress("event",&fEvent);
604
605 // read spacecharge from the Userinfo ot the tree
606 InitSpaceCharge();
607
608 TString debugName=file;
609 debugName.ReplaceAll(".root","");
610 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
611 fUseMaterial,fIdealTracking,fClusterType,
612 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
613 debugName.Append(".allClusters.debug.root");
614
615 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
616 if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
617
618 gROOT->cd();
79f54519 619
620 AliExternalTrackParam t0seed;
621 AliExternalTrackParam seed;
622 AliExternalTrackParam track;
64ea683c 623 AliExternalTrackParam tOrig;
502eb90b 624 AliToyMCTrack tOrigToy;
64ea683c 625
79f54519 626 AliExternalTrackParam *dummy;
627 AliTPCseed *seedBest;
628 AliTPCseed *seedTmp;
629 AliTPCclusterMI *cluster;
5993ed4f 630
631 Int_t maxev=fTree->GetEntries();
632 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
633
634
510cfcff 635 // ===========================================================================================
636 // Loop 1: Fill AliTPCtrackerSector structure
637 // ===========================================================================================
c0ebc0e0 638 FillSectorStructure(maxev);
510cfcff 639
510cfcff 640 // ===========================================================================================
502eb90b 641 // Loop 2: Use the TPC tracker for seeding (MakeSeeds3)
50f68b1a 642 // TODO: - check tracking configuration
643 // - add clusters and original tracks to output (how?)
510cfcff 644 // ===========================================================================================
510cfcff 645
50f68b1a 646 // settings (TODO: find the correct settings)
647 AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
648 tpcRecoParam->SetDoKinks(kFALSE);
649 AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
650 //tpcRecoParam->Print();
651
652 // need AliTPCReconstructor for parameter settings in AliTPCtracker
653 AliTPCReconstructor *tpcRec = new AliTPCReconstructor();
654 tpcRec->SetRecoParam(tpcRecoParam);
655
656 // AliTPCtracker
657 AliTPCtracker *tpcTracker = new AliTPCtracker(fTPCParam);
658 tpcTracker->SetDebug(10);
659
660 // set sector arrays
661 tpcTracker->SetTPCtrackerSectors(fInnerSectorArray,fOuterSectorArray);
662 tpcTracker->LoadInnerSectors();
663 tpcTracker->LoadOuterSectors();
664
7f1a1b08 665 // seeding
666 static TObjArray arrTracks;
79f54519 667 TObjArray * arr = &arrTracks;
668 TObjArray * seeds = new TObjArray;
7f1a1b08 669
64ea683c 670 // cuts for seeding
38d9d609 671// Float_t cuts[4];
672// cuts[0]=0.0070; // cuts[0] - fP4 cut
673// cuts[1] = 1.5; // cuts[1] - tan(phi) cut
674// cuts[2] = 3.; // cuts[2] - zvertex cut
675// cuts[3] = 3.; // cuts[3] - fP3 cut
64ea683c 676
677 // rows for seeding
678 Int_t lowerRow = fSeedingRow;
679 Int_t upperRow = fSeedingRow+2*fSeedingDist;
7f1a1b08 680 const AliTPCROC * roc = AliTPCROC::Instance();
64ea683c 681 const Int_t kNRowsInnerTPC = roc->GetNRows(0);
682 const Int_t kNRowsTPC = kNRowsInnerTPC + roc->GetNRows(36);
683 if(lowerRow < kNRowsInnerTPC){
684 Printf("Seeding row requested (%d) is lower than kNRowsInnerTPC --> use %d",lowerRow,kNRowsInnerTPC);
685 lowerRow = kNRowsInnerTPC;
686 upperRow = lowerRow + 20;
687 }
688 if(upperRow >= kNRowsTPC){
689 Printf("Seeding row requested (%d) is larger than kNRowsTPC --> use %d",upperRow,kNRowsTPC-1);
690 upperRow = kNRowsTPC-1;
691 lowerRow = upperRow-20;
692 }
693
694 // do the seeding
7f1a1b08 695 for (Int_t sec=0;sec<fkNSectorOuter;sec++){
696 //
502eb90b 697 //tpcTracker->MakeSeeds3(arr, sec,upperRow,lowerRow,cuts,-1,1);
d5704a39 698 MakeSeeds(arr, sec,upperRow,lowerRow); // own function (based on TLinearFitter)
502eb90b 699 //tpcTracker->SumTracks(seeds,arr);
700 //tpcTracker->SignClusters(seeds,3.0,3.0);
96495539 701
7f1a1b08 702 }
7f1a1b08 703
38d9d609 704 Printf("After seeding we have %d tracks",seeds->GetEntriesFast());
96495539 705
706 // Standard tracking
707 tpcTracker->SetSeeds(seeds);
708 tpcTracker->PropagateForward();
38d9d609 709 Printf("After trackinging we have %d tracks",seeds->GetEntriesFast());
710return;
50f68b1a 711
64ea683c 712 // Loop over all input tracks and connect to found seeds
713 for (Int_t iev=0; iev<maxev; ++iev){
714 printf("============== Fill Tracks: Processing Event %6d =================\n",iev);
715 fTree->GetEvent(iev);
716 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
717 printf(" > ====== Fill Tracks: Processing Track %6d ======== \n",itr);
718 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
719 tOrig = *tr;
502eb90b 720 tOrigToy = *tr;
64ea683c 721
722 Float_t z0=fEvent->GetZ();
723 Float_t t0=fEvent->GetT0();
38d9d609 724 Float_t vDrift=GetVDrift();
64ea683c 725 Float_t zLength=GetZLength(0);
7f1a1b08 726
79f54519 727 // find the corresponding seed (and track)
96495539 728 Int_t trackID = tr->GetUniqueID();
729 Int_t nClustersMC = tr->GetNumberOfSpacePoints(); // number of findable clusters (ideal)
730 if(fClusterType==1)
731 nClustersMC = tr->GetNumberOfDistSpacePoints(); // number of findable clusters (distorted)
38d9d609 732// Int_t idxSeed = 0; // index of best seed (best is with maximum number of clusters with correct ID)
96495539 733 Int_t nSeeds = 0; // number of seeds for MC track
734 Int_t nSeedClusters = 0; // number of clusters for best seed
735 Int_t nSeedClustersTmp = 0; // number of clusters for current seed
736 Int_t nSeedClustersID = 0; // number of clusters with correct ID for best seed
737 Int_t nSeedClustersIDTmp = 0; // number of clusters with correct ID for current seed
64ea683c 738 for(Int_t iSeeds = 0; iSeeds < seeds->GetEntriesFast(); ++iSeeds){
739
96495539 740 // set current seed and reset counters
64ea683c 741 seedTmp = (AliTPCseed*)(seeds->At(iSeeds));
96495539 742 nSeedClustersTmp = 0;
743 nSeedClustersIDTmp = 0;
744
502eb90b 745 if(!seedTmp) continue;
746
96495539 747 // loop over all rows
64ea683c 748 for(Int_t iRow = seedTmp->GetRow(); iRow < seedTmp->GetRow() + seedTmp->GetNumberOfClustersIndices(); iRow++ ){
7f1a1b08 749
96495539 750 // get cluster and increment counters
64ea683c 751 cluster = seedTmp->GetClusterFast(iRow);
752 if(cluster){
96495539 753 nSeedClustersTmp++;
64ea683c 754 if(cluster->GetLabel(0)==trackID){
96495539 755 nSeedClustersIDTmp++;
64ea683c 756 }
757 }
758 }
96495539 759
64ea683c 760 // if number of corresponding clusters > 0,
761 // increase nSeeds
762 if(nSeedClustersTmp > 0){
763 nSeeds++;
764 }
765
766 // if number of corresponding clusters bigger than current nSeedClusters,
96495539 767 // take this seed as the best one
768 if(nSeedClustersIDTmp > nSeedClustersID){
38d9d609 769// idxSeed = iSeeds;
79f54519 770 seedBest = seedTmp;
96495539 771 nSeedClusters = nSeedClustersTmp; // number of correctly assigned clusters
772 nSeedClustersID = nSeedClustersIDTmp; // number of all clusters
64ea683c 773 }
774
775 }
776
96495539 777 // cluster to track association (commented out, when used standard tracking)
79f54519 778 if (nSeeds>0&&nSeedClusters>0) {
96495539 779 t0seed = (AliExternalTrackParam)*seedBest;
38d9d609 780// fTime0 = t0seed.GetZ()-zLength/vDrift;
cedbc66e 781 // get the refitted track from the seed
782 // this will also set the fTime0 from the seed extrapolation
783 dummy=GetRefittedTrack(*seedBest);
784 track=*dummy;
785 delete dummy;
96495539 786 //printf("seed (%.2g): %d seeds with %d clusters\n",fTime0,nSeeds,nSeedClusters);
787
788 // // cluster to track association for all good seeds
789 // // set fCreateT0seed to true to get the seed in time coordinates
790 // fCreateT0seed = kTRUE;
791 // dummy = ClusterToTrackAssociation(seedBest,trackID,nClus);
792
793 //// propagate track to 0
794 //const Double_t kMaxSnp = 0.85;
795 //const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
796 //AliTrackerBase::PropagateTrackTo(&track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
79f54519 797
798 }
799
64ea683c 800 Int_t ctype(fCorrectionType);
801
802 if (fStreamer) {
803 (*fStreamer) << "Tracks" <<
96495539 804 "iev=" << iev <<
805 "z0=" << z0 <<
806 "t0=" << t0 <<
807 "fTime0=" << fTime0 <<
808 "itr=" << itr <<
809 "clsType=" << fClusterType <<
810 "corrType=" << ctype <<
811 "seedRow=" << fSeedingRow <<
812 "seedDist=" << fSeedingDist <<
38d9d609 813 "vDrift=" << vDrift <<
96495539 814 "zLength=" << zLength <<
815 "nClustersMC=" << nClustersMC <<
816 "nSeeds=" << nSeeds <<
817 "nSeedClusters=" << nSeedClusters <<
818 "nSeedClustersID=" << nSeedClustersID <<
819 "t0seed.=" << &t0seed <<
cedbc66e 820 "track.=" << &track <<
96495539 821 "tOrig.=" << &tOrig <<
502eb90b 822 "tOrigToy.=" << &tOrigToy <<
64ea683c 823 "\n";
824 }
825 }
826 }
827
50f68b1a 828
510cfcff 829 delete fStreamer;
830 fStreamer=0x0;
831
832 delete fEvent;
833 fEvent = 0x0;
834
835 delete fTree;
836 fTree=0x0;
837 f.Close();
838}
839
840
a06336b6 841//____________________________________________________________________________________
842void AliToyMCReconstruction::RunFullTracking(const char* file, Int_t nmaxEv)
843{
844 //
845 //
846 //
847
848 ConnectInputFile(file,nmaxEv);
849 if (!fTree) return;
850
851 InitStreamer(".fullTracking");
852
853 FillSectorStructureAC();
854
855 AliTPCReconstructor::SetStreamLevel(0);
856
857 TObjArray seeds;
858 seeds.SetOwner();
9d548f3c 859 Int_t lowerRow=130;
860 Int_t upperRow=150;
861 if (lowerRow>upperRow){
862 Int_t tmp=lowerRow;
863 lowerRow=upperRow;
864 upperRow=tmp;
865 }
a06336b6 866
9d548f3c 867 // seeding.
868 // NOTE: the z position is set to GetTimeBin*vDrift
869 // therefore it is not possible to simply propagate
870 // the track using AliTrackerBase::Propagate, since a
871 // wrong B-Field will be assinged...
872 printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
a06336b6 873 for (Int_t sec=0;sec<36;sec++){
9d548f3c 874 printf(" in sector: %d\n",sec);
875 Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
876 printf(" -> Added Seeds: %d\n",nAdded);
877 nAdded=MakeSeeds2(&seeds, sec,lowerRow-2,upperRow-2);
878 printf(" -> Added Seeds: %d\n",nAdded);
879 nAdded=MakeSeeds2(&seeds, sec,lowerRow-4,upperRow-4);
880 printf(" -> Added Seeds: %d\n",nAdded);
a06336b6 881 }
882
9d548f3c 883 printf("Run Tracking in %3d - %3d\n",lowerRow,upperRow);
884 Int_t firstSeed=0;
885 for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
886 //first seed is used to not run the tracking twice on a seed
887 firstSeed=seeds.GetEntriesFast();
888// DumpTrackInfo(&seeds);
a06336b6 889
9d548f3c 890 lowerRow=110;
891 upperRow=130;
892
893 printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
894 for (Int_t sec=0;sec<36;sec++){
895 printf(" in sector: %d\n",sec);
896 Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
897 printf(" -> Added Seeds: %d\n",nAdded);
898 }
899 printf("Run Tracking in %3d - %3d\n",lowerRow,upperRow);
900 for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
901 firstSeed=seeds.GetEntriesFast();
902
903 //now seeding also at more central rows with shorter seeds
904 lowerRow=70;
905 upperRow=90;
906
907 printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
908 for (Int_t sec=0;sec<36;sec++){
909 printf(" in sector: %d\n",sec);
910 Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
911 printf(" -> Added Seeds: %d\n",nAdded);
912 }
913 printf("Run Tracking in %3d - %3d\n",lowerRow,upperRow);
914 for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
915 firstSeed=seeds.GetEntriesFast();
916
917 //shorter seeds
afc28cd1 918 Int_t startUpper=upperRow-10;
919 Int_t startLower=lowerRow-5;
9d548f3c 920 for (Int_t sec=0;sec<36;sec++){
afc28cd1 921 upperRow=startUpper;
922 lowerRow=startLower;
9d548f3c 923 printf(" in sector: %d\n",sec);
924 while (lowerRow>0){
925 printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
926 Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
927 printf(" -> Added Seeds: %d\n",nAdded);
928 for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
929 firstSeed=seeds.GetEntriesFast();
930 lowerRow-=5;
931 upperRow-=5;
932 }
933 }
934
935 //track remaining
936
937 DumpTrackInfo(&seeds);
938
939// TObjArray seedsCentral2;
940// lowerRow=45;
941// upperRow=62;
942//
943// for (Int_t sec=0;sec<36;sec++){
944// Int_t nAdded=MakeSeeds2(&seedsCentral2, sec,lowerRow,upperRow);
945// printf(" -> Added Seeds: %d\n",nAdded);
946// nAdded=MakeSeeds2(&seedsCentral2, sec,lowerRow-2,upperRow-2);
947// printf(" -> Added Seeds: %d\n",nAdded);
948// nAdded=MakeSeeds2(&seedsCentral2, sec,lowerRow-4,upperRow-4);
949// printf(" -> Added Seeds: %d\n",nAdded);
950// }
951// for (Int_t iseed=0; iseed<seedsCentral2.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seedsCentral2.UncheckedAt(iseed));
952// DumpTrackInfo(&seedsCentral2);
953
954 //dump clusters
e83fd282 955// (*fStreamer) << "clusters" <<
956// "cl.=" << &fAllClusters << "\n";
9d548f3c 957
a06336b6 958 Cleanup();
959}
960
e83fd282 961//____________________________________________________________________________________
962AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrackIdeal(const AliToyMCTrack * const tr, EDet det )
963{
964 //
965 // crate a seed from the track points of the respective detector
966 //
967 AliTrackPoint seedPoint[3];
968
969 Int_t npoints=0;
970 switch (det) {
971 case kITS:
972 npoints=tr->GetNumberOfITSPoints();
973 break;
974 case kTPC:
975 npoints=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
976 break;
977 case kTRD:
978 npoints=tr->GetNumberOfTRDPoints();
979 break;
980 }
981
982 if (npoints<3) return 0x0;
983
984 Int_t pos[3]={0,npoints/2,npoints-1};
985 const AliCluster *cl=0x0;
986
987 for (Int_t ipoint=0;ipoint<3;++ipoint){
988 Int_t cluster=pos[ipoint];
989 switch (det) {
990 case kITS:
991 seedPoint[ipoint]=(*tr->GetITSPoint(cluster));
992 break;
993 case kTPC:
994 cl=tr->GetSpacePoint(cluster);
995 if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(cluster);
996 AliTPCclusterMI::SetGlobalTrackPoint(*cl,seedPoint[ipoint]);
997 break;
998 case kTRD:
999 seedPoint[ipoint]=(*tr->GetTRDPoint(cluster));
1000 break;
1001 }
1002 }
1003
1004 AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
1005 seed->ResetCovariance(10);
1006
1007 return seed;
1008}
1009
d1cf83f5 1010//____________________________________________________________________________________
9d548f3c 1011AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr, Bool_t forceSeed)
d1cf83f5 1012{
1013 //
1014 // if we don't have a valid time0 informaion (fTime0) available yet
1015 // assume we create a seed for the time0 estimate
1016 //
1017
1018 // seed point informaion
1019 AliTrackPoint seedPoint[3];
1020 const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
1021
1022 // number of clusters to loop over
1023 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
9d548f3c 1024 if (ncls<3){
1025 AliError(Form("Not enough points to create a seed: %d",ncls));
1026 return 0x0;
1027 }
d1cf83f5 1028 UChar_t nextSeedRow=fSeedingRow;
1029 Int_t nseeds=0;
1030
1031 //assumes sorted clusters
9d548f3c 1032 if (forceSeed){
1033 // force the seed creation, using the first, middle and last cluster
1034 Int_t npoints[3]={0,ncls/2,ncls-1};
1035 for (Int_t icl=0;icl<3;++icl){
1036 const AliTPCclusterMI *cl=tr->GetSpacePoint(npoints[icl]);
1037 if (fClusterType==1) cl=tr->GetDistortedSpacePoint(npoints[icl]);
d1cf83f5 1038 seedCluster[nseeds]=cl;
1039 SetTrackPointFromCluster(cl, seedPoint[nseeds]);
1040 ++nseeds;
9d548f3c 1041 }
1042 }else{
1043 // create seeds according to the reco settings
1044 for (Int_t icl=0;icl<ncls;++icl) {
1045 const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
1046 if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
1047 if (!cl) continue;
1048 // use row in sector
1049 const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
1050 // skip clusters without proper pad row
1051 if (row>200) continue;
1052
1053 //check seeding row
1054 // if we are in the last row and still miss a seed we use the last row
1055 // even if the row spacing will not be equal
1056 if (row>=nextSeedRow || icl==ncls-1){
1057 seedCluster[nseeds]=cl;
1058 SetTrackPointFromCluster(cl, seedPoint[nseeds]);
1059 ++nseeds;
d1cf83f5 1060 nextSeedRow+=fSeedingDist;
1061
1062 if (nseeds==3) break;
9d548f3c 1063 }
d1cf83f5 1064 }
1065 }
1066
1067 // check we really have 3 seeds
1068 if (nseeds!=3) {
1069 AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
1070 return 0x0;
1071 }
c66b6495 1072
7f72a764 1073 // determine preliminary theta
1074 Float_t xyz1[3]={0,0,0};
1075 Float_t xyz2[3]={0,0,0};
1076 seedPoint[0].GetXYZ(xyz1);
1077 seedPoint[2].GetXYZ(xyz2);
1078 Float_t prelDeltaR = TMath::Sqrt(xyz2[0]*xyz2[0]+xyz2[1]*xyz2[1]) - TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]) ;
1079 Float_t prelDeltaZ = ( seedCluster[0]->GetTimeBin() - seedCluster[2]->GetTimeBin() ) * GetVDrift();
1080 Float_t prelTheta = TMath::ATan(prelDeltaR/prelDeltaZ);
1081 if(prelTheta > TMath::Pi()/2) prelTheta = TMath::Pi() - prelTheta;
d1cf83f5 1082
1083 // do cluster correction for fCorrectionType:
1084 // 0 - no correction
1085 // 1 - TPC center
1086 // 2 - average eta
1087 // 3 - ideal
88386c5c 1088 // 4 - preliminary eta (needs fixing!!! Not yet in full code!!!)
d1cf83f5 1089 // assign the cluster abs time as z component to all seeds
1090 for (Int_t iseed=0; iseed<3; ++iseed) {
1091 Float_t xyz[3]={0,0,0};
1092 seedPoint[iseed].GetXYZ(xyz);
600eaa0d 1093 const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
d1cf83f5 1094
600eaa0d 1095 const Int_t sector=seedCluster[iseed]->GetDetector();
1096 const Int_t sign=1-2*((sector/18)%2);
c66b6495 1097
88386c5c 1098 Float_t zBeforeCorr = xyz[2];
d1cf83f5 1099
1100 if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
600eaa0d 1101 // the settings below are for the T0 seed
1102 // for known T0 the z position is already calculated in SetTrackPointFromCluster
1103 if ( fCreateT0seed ){
1104 if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
1105 //!!! TODO: is this the correct association?
1106 if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
88386c5c 1107 if ( fCorrectionType == kPreliminaryEta ) xyz[2] = r/TMath::Tan(prelTheta)*sign;//(needs fixing!!! Not yet in full code!!!)
600eaa0d 1108 }
4a777885 1109
d1cf83f5 1110 if ( fCorrectionType == kIdeal ) xyz[2] = seedCluster[iseed]->GetZ();
48e13991 1111
88386c5c 1112 // Store xyz only here!!! To get the Delta z from the correction...
1113 zBeforeCorr = xyz[2];
1114
d1cf83f5 1115 //!!! TODO: to be replaced with the proper correction
223d9e38 1116 fTPCCorrection->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
d1cf83f5 1117 }
1118
600eaa0d 1119 // after the correction set the time bin as z-Position in case of a T0 seed
1120 if ( fCreateT0seed )
c484014f 1121 xyz[2]=seedCluster[iseed]->GetTimeBin() + sign*( zBeforeCorr - xyz[2] )/GetVDrift();
48e13991 1122
d1cf83f5 1123 seedPoint[iseed].SetXYZ(xyz);
1124 }
1125
1126 const Double_t kMaxSnp = 0.85;
1127 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1128
1129 AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
d1cf83f5 1130
48e13991 1131 if (fCreateT0seed&&!fLongT0seed){
1132 // only propagate to vertex if we don't create a long seed
d1cf83f5 1133 // if fTime0 < 0 we assume that we create a seed for the T0 estimate
36a79b0d 1134 Int_t ret=AliTrackerBase::PropagateTrackTo2(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
d1cf83f5 1135 if (TMath::Abs(seed->GetX())>3) {
36a79b0d 1136// printf("Could not propagate track to 0, x:%.2f, a:%.2f (%.2f), snp:%.2f (%.2f), pt:%.2f (%.2f), %d\n",seed->GetX(),seed->GetAlpha(),tr->GetAlpha(), seed->GetSnp(), tr->GetSnp(), seed->Pt(), tr->Pt(), ret);
1137 }
1138 if (fForceAlpha) {
1139 seed->Rotate(tr->GetAlpha());
1140 AliTrackerBase::PropagateTrackTo2(seed,0,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
d1cf83f5 1141 }
1142 }
48e13991 1143
1144 seed->ResetCovariance(10);
d1cf83f5 1145 return seed;
1146
1147}
1148
e83fd282 1149//____________________________________________________________________________________
1150AliExternalTrackParam* AliToyMCReconstruction::GetTrackRefit(const AliToyMCTrack * const tr, EDet det)
1151{
1152 //
1153 // Get the ITS or TRD track refitted from the toy track
1154 // type: 0=ITS; 1=TRD
1155 //
1156
1157 AliExternalTrackParam *track=GetSeedFromTrackIdeal(tr,det);
1158 if (!track) return 0x0;
1159
1160 const Double_t kMaxSnp = 0.85;
1161 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1162
1163 Int_t npoints=0;
1164 switch (det) {
1165 case kITS:
1166 npoints=tr->GetNumberOfITSPoints();
1167 break;
1168 case kTPC:
1169 npoints=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1170 break;
1171 case kTRD:
1172 npoints=tr->GetNumberOfTRDPoints();
1173 break;
1174 }
1175
1176 const AliCluster *cl=0x0;
1177
1178 for (Int_t ipoint=0; ipoint<npoints; ++ipoint) {
1179 AliTrackPoint pIn;
1180
1181 switch (det) {
1182 case kITS:
1183 pIn=(*tr->GetITSPoint(ipoint));
1184 break;
1185 case kTPC:
1186 cl=tr->GetSpacePoint(ipoint);
1187 if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
1188 AliTPCclusterMI::SetGlobalTrackPoint(*cl,pIn);
1189 break;
1190 case kTRD:
1191 pIn=(*tr->GetTRDPoint(ipoint));
1192 break;
1193 }
1194
1195
1196 const Double_t angle=pIn.GetAngle();
1197 track->Rotate(angle);
1198 AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1199
1200 if (!AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial)) {
1201 AliInfo(Form("Could not propagate track to x=%.2f (a=%.2f) for det %d",prot.GetX(),angle,det));
1202 }
1203 //
1204
1205 Double_t pointPos[2]={0,0};
1206 Double_t pointCov[3]={0,0,0};
1207 pointPos[0]=prot.GetY();//local y
1208 pointPos[1]=prot.GetZ();//local z
1209 pointCov[0]=prot.GetCov()[3];//simay^2
1210 pointCov[1]=prot.GetCov()[4];//sigmayz
1211 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1212
1213 if (!track->Update(pointPos,pointCov)) {
1214 AliInfo(Form("no update: det: %d",det));
1215 break;
1216 }
1217
1218 }
1219
1220 return track;
1221}
1222
79f54519 1223
d1cf83f5 1224//____________________________________________________________________________________
1225void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
1226{
1227 //
1228 // make AliTrackPoint out of AliTPCclusterMI
1229 //
1230
1231 if (!cl) return;
4a777885 1232 Float_t xyz[3]={0.,0.,0.};
d1cf83f5 1233 // ClusterToSpacePoint(cl,xyz);
1234 // cl->GetGlobalCov(cov);
1235 //TODO: what to do with the covariance matrix???
1236 //TODO: the problem is that it is used in GetAngle in AliTrackPoint
1237 //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
1238 //TODO: for the moment simply assign 1 permill squared
1239 // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
1240 // Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
1241 // xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
1242 // cl->GetGlobalXYZ(xyz);
1243 // cl->GetGlobalCov(cov);
1244 // voluem ID to add later ....
1245 // p.SetXYZ(xyz);
1246 // p.SetCov(cov);
0403120d 1247// AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
1248// p=*tp;
1249// delete tp;
e83fd282 1250// const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
1251 AliTPCclusterMI::SetGlobalTrackPoint(*cl,p);
d1cf83f5 1252 // cl->Print();
1253 // p.Print();
1254 p.SetVolumeID(cl->GetDetector());
600eaa0d 1255
1256
1257 if ( !fCreateT0seed && !fIdealTracking ) {
1258 p.GetXYZ(xyz);
1259 const Int_t sector=cl->GetDetector();
1260 const Int_t sign=1-2*((sector/18)%2);
1261 const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
05da1b4e 1262// printf(" z: %.2f %.2f\n",xyz[2],zT0);
600eaa0d 1263 xyz[2]=zT0;
1264 p.SetXYZ(xyz);
4a777885 1265 }
600eaa0d 1266
1267
d1cf83f5 1268 // p.Rotate(p.GetAngle()).Print();
1269}
1270
1271//____________________________________________________________________________________
1272void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
1273{
1274 //
1275 // convert the cluster to a space point in global coordinates
1276 //
1277 if (!cl) return;
1278 xyz[0]=cl->GetRow();
1279 xyz[1]=cl->GetPad();
1280 xyz[2]=cl->GetTimeBin(); // this will not be correct at all
1281 Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
1282 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
1283 fTPCParam->Transform8to4(xyz,i);
1284 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
1285 fTPCParam->Transform4to3(xyz,i);
1286 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
1287 fTPCParam->Transform2to1(xyz,i);
1288 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
1289}
1290
1291//____________________________________________________________________________________
4388aa4a 1292AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, TClonesArray *arrClustRes)
d1cf83f5 1293{
1294 //
1295 //
1296 //
1297
4388aa4a 1298 if (arrClustRes) {
1299 arrClustRes->Clear();
1300 }
1301
d1cf83f5 1302 // create track
1303 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
5414766f 1304 // track copy for propagation
1305 AliExternalTrackParam trCopy(*tr);
d1cf83f5 1306
1307 Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1308
1309 const AliTPCROC * roc = AliTPCROC::Instance();
1310
1311 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
1312 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1313 const Double_t kMaxSnp = 0.85;
1314 const Double_t kMaxR = 500.;
1315 const Double_t kMaxZ = 500.;
1316
1317 // const Double_t kMaxZ0=220;
1318// const Double_t kZcut=3;
1319
1320 const Double_t refX = tr->GetX();
1321
1322 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
5414766f 1323
1324 // parametrised track resolution
1325 Double_t trackRes=gRandom->Gaus();
d1cf83f5 1326
1327 // loop over all other points and add to the track
4a777885 1328 for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
d1cf83f5 1329 AliTrackPoint pIn;
1330 const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
1331 if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
8a2a61d3 1332 const Int_t globalRow = cl->GetRow() +(cl->GetDetector() >35)*63;
1333 if ( fCreateT0seed ){
1334 if ( globalRow<fSeedingRow || globalRow>fSeedingRow+2*fSeedingDist ) continue;
1335 }
1336
d1cf83f5 1337 SetTrackPointFromCluster(cl, pIn);
8a2a61d3 1338
1339 Float_t xyz[3]={0,0,0};
1340 pIn.GetXYZ(xyz);
48e13991 1341 Float_t zBeforeCorr = xyz[2];
1342
1343 const Int_t sector=cl->GetDetector();
1344 const Int_t sign=1-2*((sector/18)%2);
8a2a61d3 1345
d1cf83f5 1346 if (fCorrectionType != kNoCorrection){
7f72a764 1347
8a2a61d3 1348 const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
8a2a61d3 1349
1350 if ( fCreateT0seed ){
1351 if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
1352 //!!! TODO: is this the correct association?
1353 if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
1354 if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
1355 }
88386c5c 1356
1357 // Store xyz only here!!! To get the Delta z from the correction...
1358 zBeforeCorr = xyz[2];
8a2a61d3 1359
223d9e38 1360 fTPCCorrection->CorrectPoint(xyz, cl->GetDetector());
d1cf83f5 1361 }
8a2a61d3 1362
1363 if ( fCreateT0seed )
c484014f 1364 xyz[2]=cl->GetTimeBin() + sign*( zBeforeCorr - xyz[2] )/GetVDrift();
48e13991 1365 // xyz[2]=cl->GetTimeBin();
8a2a61d3 1366 pIn.SetXYZ(xyz);
1367
d1cf83f5 1368 // rotate the cluster to the local detector frame
1369 track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
1370 AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1371 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1372 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1373 //
05da1b4e 1374 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
d1cf83f5 1375
1376 if (!res) break;
1377
1378 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1379 if (TMath::Abs(track->GetX())>kMaxR) break;
1380// if (TMath::Abs(track->GetZ())<kZcut)continue;
1381 //
4388aa4a 1382
1383 // add residuals
1384 if (arrClustRes) {
1385 TClonesArray &arrDummy=*arrClustRes;
1386 AliTPCclusterMI *clRes = new(arrDummy[arrDummy.GetEntriesFast()]) AliTPCclusterMI(*cl);
1387 clRes->SetX(prot.GetX());
b7106990 1388 // residuals in terms of sigma cl and track
1389 clRes->SetY((track->GetY()-prot.GetY())/( sqrt ( prot.GetCov()[3] + track->GetSigmaY2()) ) );
1390 clRes->SetZ((track->GetZ()-prot.GetZ())/( sqrt ( prot.GetCov()[5] + track->GetSigmaZ2()) ) );
4388aa4a 1391 }
5414766f 1392
1393 // fill cluster residuals to ideal track for calibration studies
1394 // ideal cluster position
70cc5429 1395 // require at least 2 TRD points
842d682d 1396 if (fRecoInfo<0 || (fRecoInfo&kFillDeltas) ==kFillDeltas ) {
70cc5429 1397 trCopy.Rotate(track->GetAlpha());
1398 AliTrackerBase::PropagateTrackTo(&trCopy,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1399 // binning r, phi, z, delta (0=rphi, 1=z)
1400 // resolution parametrisation
842d682d 1401 Float_t soneOverPt= trCopy.GetSigned1Pt();
1402 Float_t oneOverPt = TMath::Abs(soneOverPt);
1403 Float_t radius = trCopy.GetX();
1404 Float_t trackY = trCopy.GetY();
1405 Float_t trackZ = trCopy.GetZ();
1406 Float_t trackPhi = trCopy.Phi();
9c4e14d7 1407 Float_t alpha = trCopy.GetAlpha();
70cc5429 1408
842d682d 1409 Float_t pointY = prot.GetY();
1410 Float_t pointZ = prot.GetZ();
1411
1412 Float_t resRphi = 0.004390 + oneOverPt*(-0.136403) + oneOverPt*radius*(0.002266) + oneOverPt*radius*radius*(-0.000006);
1413
1414 Float_t resRphiRandom = resRphi*trackRes;
1415 Float_t deviation = trackY+resRphiRandom-pointY;
1416 Short_t npTRD = tr->GetNumberOfTRDPoints();
70cc5429 1417
1418 // rphi residuals
c1e1ba89 1419 Double_t xx[4]={radius, trackPhi, trackZ ,deviation};
842d682d 1420 if (npTRD>=2){
1421 fHnDelta->Fill(xx);
1422 }
1423
1424 Short_t event=fTree->GetReadEntry();
1425
1426 if (fStreamer) {
1427 (*fStreamer) << "delta" <<
1428 "soneOverPt=" << soneOverPt <<
1429 "r=" << radius <<
1430 "trackPhi=" << trackPhi <<
1431 "trackY=" << trackY <<
1432 "trackZ=" << trackZ <<
9c4e14d7 1433 "alpha=" << alpha <<
842d682d 1434 "resRphi=" << resRphi <<
1435 "trackRes=" << trackRes <<
1436 "pointY=" << pointY <<
1437 "pointZ=" << pointZ <<
1438 "npTRD=" << npTRD <<
1439 "event=" << event <<
1440 "\n";
1441// "point.=" << &prot <<
1442// "track.=" << track <<
1443 }
70cc5429 1444 }
4388aa4a 1445
d1cf83f5 1446 Double_t pointPos[2]={0,0};
1447 Double_t pointCov[3]={0,0,0};
1448 pointPos[0]=prot.GetY();//local y
1449 pointPos[1]=prot.GetZ();//local z
1450 pointCov[0]=prot.GetCov()[3];//simay^2
1451 pointCov[1]=prot.GetCov()[4];//sigmayz
1452 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1453
1454 if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
1455 }
1456
05da1b4e 1457 AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
d1cf83f5 1458
36a79b0d 1459 if (!fCreateT0seed || fForceAlpha){
48e13991 1460 // rotate fittet track to the frame of the original track and propagate to same reference
1461 track->Rotate(tr->GetAlpha());
d1cf83f5 1462
48e13991 1463 AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1464 }
d1cf83f5 1465
1466 return track;
1467}
1468
510cfcff 1469//____________________________________________________________________________________
c0ebc0e0 1470AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
510cfcff 1471{
1472 //
1473 // Tracking for given seed on an array of clusters
1474 //
1475
1476 // create track
1477 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
5993ed4f 1478
1479 const AliTPCROC * roc = AliTPCROC::Instance();
1480
1481 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
1482 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
c0ebc0e0 1483 const Int_t kNRowsTPC = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1484 const Int_t kIRowsTPC = roc->GetNRows(0) - 1;
5993ed4f 1485 const Double_t kMaxSnp = 0.85;
1486 const Double_t kMaxR = 500.;
1487 const Double_t kMaxZ = 500.;
c0ebc0e0 1488 const Double_t roady = 100.;
1489 const Double_t roadz = 100.;
5993ed4f 1490
1491 const Double_t refX = tr->GetX();
1492
1493 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
c0ebc0e0 1494
1495 Int_t secCur = -1;
1496 UInt_t indexCur = 0;
1497 Double_t xCur, yCur, zCur = 0.;
5993ed4f 1498
64ff24a8 1499 Float_t vDrift = GetVDrift();
1500
5993ed4f 1501 // first propagate seed to outermost row
1502 AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
510cfcff 1503
5993ed4f 1504 // Loop over rows and find the cluster candidates
1505 for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
c0ebc0e0 1506
5993ed4f 1507 // inner or outer sector
1508 Bool_t bInnerSector = kTRUE;
1509 if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1510
c0ebc0e0 1511 // nearest track point/cluster (to be found)
1512 AliTrackPoint nearestPoint;
1513 AliTPCclusterMI *nearestCluster = NULL;
1514
1515 // Inner Sector
5993ed4f 1516 if(bInnerSector){
c0ebc0e0 1517
1518 // Propagate to center of pad row and extract parameters
5993ed4f 1519 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
c0ebc0e0 1520 xCur = track->GetX();
1521 yCur = track->GetY();
1522 zCur = track->GetZ();
64ff24a8 1523 if ( !fIdealTracking ) {
1524 zCur = zCur/vDrift + fTime0; // Look at time, not at z!
1525 }
c0ebc0e0 1526 secCur = GetSector(track);
1527
1528 // Find the nearest cluster (TODO: correct road settings!)
7f1a1b08 1529 Printf("inner tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
c0ebc0e0 1530 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1531
1532 // Move to next row if now cluster found
1533 if(!nearestCluster) continue;
1534 //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1535
5993ed4f 1536 }
c0ebc0e0 1537
1538 // Outer sector
5993ed4f 1539 else{
c0ebc0e0 1540
1541 // Propagate to center of pad row and extract parameters
5993ed4f 1542 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
c0ebc0e0 1543 xCur = track->GetX();
1544 yCur = track->GetY();
1545 zCur = track->GetZ();
64ff24a8 1546 if ( !fIdealTracking ) {
1547 zCur = zCur/vDrift + fTime0; // Look at time, not at z!
1548 }
c0ebc0e0 1549 secCur = GetSector(track);
1550
1551 // Find the nearest cluster (TODO: correct road settings!)
7f1a1b08 1552 Printf("outer tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
c0ebc0e0 1553 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1554
1555 // Move to next row if now cluster found
1556 if(!nearestCluster) continue;
1557 //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
1558
5993ed4f 1559 }
c0ebc0e0 1560
1561 // create track point from cluster
1562 SetTrackPointFromCluster(nearestCluster,nearestPoint);
1563
1564 //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
1565
64ff24a8 1566 // correction
1567 // TODO: also correction when looking for the next cluster?
1568 if (fCorrectionType != kNoCorrection){
1569 Float_t xyz[3]={0,0,0};
1570 nearestPoint.GetXYZ(xyz);
223d9e38 1571 fTPCCorrection->CorrectPoint(xyz, nearestCluster->GetDetector());
64ff24a8 1572 nearestPoint.SetXYZ(xyz);
1573 }
1574
c0ebc0e0 1575 // rotate the cluster to the local detector frame
1576 track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1577 AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1578 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1579 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1580
64ff24a8 1581
c0ebc0e0 1582 //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
1583
1584 // update track with the nearest track point
223d9e38 1585 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
510cfcff 1586
5993ed4f 1587 if (!res) break;
1588
1589 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1590 if (TMath::Abs(track->GetX())>kMaxR) break;
1591 //if (TMath::Abs(track->GetZ())<kZcut)continue;
1592
c0ebc0e0 1593 Double_t pointPos[2]={0,0};
1594 Double_t pointCov[3]={0,0,0};
1595 pointPos[0]=prot.GetY();//local y
1596 pointPos[1]=prot.GetZ();//local z
1597 pointCov[0]=prot.GetCov()[3];//simay^2
1598 pointCov[1]=prot.GetCov()[4];//sigmayz
1599 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1600
1601 if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
5993ed4f 1602
c0ebc0e0 1603 ++nClus;
5993ed4f 1604 }
c0ebc0e0 1605
5993ed4f 1606
c0ebc0e0 1607 // propagation to refX
223d9e38 1608 AliTrackerBase::PropagateTrackTo(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
5993ed4f 1609
1610 // rotate fittet track to the frame of the original track and propagate to same reference
1611 track->Rotate(tr->GetAlpha());
1612
223d9e38 1613 AliTrackerBase::PropagateTrackTo(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
c0ebc0e0 1614
1615 Printf("We have %d clusters in this track!",nClus);
5993ed4f 1616
510cfcff 1617 return track;
1618}
1619
79f54519 1620//____________________________________________________________________________________
1621AliExternalTrackParam* AliToyMCReconstruction::ClusterToTrackAssociation(const AliTPCseed *seed, Int_t trackID, Int_t &nClus)
1622{
1623 //
1624 // Cluster to track association for given seed on an array of clusters
1625 //
1626
1627 // create track
1628 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1629
1630 const AliTPCROC * roc = AliTPCROC::Instance();
1631
1632 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
1633 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1634 const Int_t kNRowsTPC = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1635 const Int_t kIRowsTPC = roc->GetNRows(0) - 1;
1636 const Double_t kMaxSnp = 0.85;
1637 const Double_t kMaxR = 500.;
1638 const Double_t kMaxZ = 500.;
1639 const Double_t roady = 0.1;
1640 const Double_t roadz = 0.01;
1641
1642 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1643
6f2a2a54 1644 Int_t secCur, secOld = -1;
79f54519 1645 UInt_t indexCur = 0;
1646 Double_t xCur, yCur, zCur = 0.;
1647
cedbc66e 1648// Float_t vDrift = GetVDrift();
79f54519 1649
1650 // first propagate seed to outermost row
96495539 1651 Bool_t res0=AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
79f54519 1652
1653 // Loop over rows and find the cluster candidates
1654 for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
1655
1656 // inner or outer sector
1657 Bool_t bInnerSector = kTRUE;
1658 if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1659
1660 // nearest track point/cluster (to be found)
1661 AliTrackPoint nearestPoint;
1662 AliTPCclusterMI *nearestCluster = NULL;
1663
1664 // Inner Sector
1665 if(bInnerSector){
1666
1667 // Propagate to center of pad row and extract parameters
1668 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
1669 xCur = track->GetX();
1670 yCur = track->GetY();
1671 zCur = track->GetZ();
1672 secCur = GetSector(track);
1673
1674 // Find the nearest cluster (TODO: correct road settings!)
6f2a2a54 1675 //Printf("inner tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
79f54519 1676 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
6f2a2a54 1677
1678 // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1679 // Increase also the road in this case
96495539 1680 if(!nearestCluster && secCur != secOld && secOld > -1){
6f2a2a54 1681 //Printf("inner tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1682 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1683 }
79f54519 1684
1685 // Move to next row if now cluster found
1686 if(!nearestCluster) continue;
1687 //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1688
1689 }
1690
1691 // Outer sector
1692 else{
1693
1694 // Propagate to center of pad row and extract parameters
1695 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
1696 xCur = track->GetX();
1697 yCur = track->GetY();
1698 zCur = track->GetZ();
1699 secCur = GetSector(track);
1700
1701 // Find the nearest cluster (TODO: correct road settings!)
96495539 1702 Printf("res0 = %d, outer tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",res0,xCur,yCur,zCur,iRow,secCur);
79f54519 1703 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1704
6f2a2a54 1705 // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1706 // Increase also the road in this case
96495539 1707 if(!nearestCluster && secCur != secOld && secOld > -1){
1708 Printf("outer tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
6f2a2a54 1709 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1710 }
1711
1712
79f54519 1713 // Move to next row if now cluster found
1714 if(!nearestCluster) continue;
96495539 1715 Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
79f54519 1716
1717 }
1718
1719 // create track point from cluster
1720 SetTrackPointFromCluster(nearestCluster,nearestPoint);
1721
6f2a2a54 1722 //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
79f54519 1723
1724 // rotate the cluster to the local detector frame
1725 track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1726 AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1727 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1728 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1729
1730
6f2a2a54 1731 //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
79f54519 1732
1733 // update track with the nearest track point
1734 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1735
1736 if (!res) break;
1737
1738 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1739 if (TMath::Abs(track->GetX())>kMaxR) break;
1740 //if (TMath::Abs(track->GetZ())<kZcut)continue;
1741
1742 Double_t pointPos[2]={0,0};
1743 Double_t pointCov[3]={0,0,0};
1744 pointPos[0]=prot.GetY();//local y
1745 pointPos[1]=prot.GetZ();//local z
1746 pointCov[0]=prot.GetCov()[3];//simay^2
1747 pointCov[1]=prot.GetCov()[4];//sigmayz
1748 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1749
1750 if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
6f2a2a54 1751 secOld = secCur;
79f54519 1752
6f2a2a54 1753 //Printf("Cluster belongs to track = %d",nearestCluster->GetLabel(0));
79f54519 1754
1755 // only count as associate cluster if it belongs to correct track!
1756 if(nearestCluster->GetLabel(0) == trackID)
1757 ++nClus;
1758 }
1759
1760 Printf("We have %d clusters in this track!",nClus);
1761
1762 return track;
1763}
1764
9d548f3c 1765//____________________________________________________________________________________
1766void AliToyMCReconstruction::AssociateClusters(AliTPCseed &seed, Int_t firstRow, Int_t lastRow, Bool_t direction)
1767{
1768 //
1769 // do cluster to track association from first to last row
1770 // direction 0: outwards; 1: inwards
1771 //
1772
1773 Double_t roady=10.;
1774 Double_t roadz=10.;
1775
1776 AliRieman rieman1(160);
1777 AliRieman rieman2(160);
1778 SetRieman(seed,rieman1);
1779 CopyRieman(rieman1,rieman2);
1780
1781 Int_t sec=seed.GetSector();
1782 Int_t noLastPoint=0;
1783 //TODO: change to inward and outwar search?
1784 // -> better handling of non consecutive points
1785 if (direction){
1786 firstRow*=-1;
1787 lastRow*=-1;
1788 }
1789
1790 //always from inside out
1791 if (firstRow>lastRow){
1792 Int_t tmp=firstRow;
1793 firstRow=lastRow;
1794 lastRow=tmp;
1795 }
1796
1797 for (Int_t row=firstRow; row<=lastRow && noLastPoint<3;++row) {
1798 Int_t iRow=TMath::Abs(row);
1799 const AliTPCclusterMI *cl=seed.GetClusterPointer(iRow);
1800 if (cl) continue;
1801
1802 const Int_t secrow = iRow<63?iRow:iRow-63;
1803
1804 AliTPCtrackerSector *arrSec=(iRow<63)?fInnerSectorArray:fOuterSectorArray;
1805 const AliTPCtrackerRow& kr = arrSec[sec%36][secrow];
1806 const Double_t maxy=arrSec[sec%36].GetMaxY(secrow);
1807
1808 Double_t y=rieman1.GetYat(kr.GetX());
1809 Double_t z=rieman1.GetZat(kr.GetX());
1810
1811 if (TMath::Abs(y)>maxy) {
1812 AliError("Tracking over sector boundaries not implemented, yet");
1813 continue;
1814 }
1815
1816 AliTPCclusterMI *n=kr.FindNearest(y,z,roady,roadz);
1817 if (!n || n->IsUsed()) {
1818 ++noLastPoint;
1819 continue;
1820 }
1821 // check for quality of the cluster
1822 // TODO: better?
1823 rieman2.AddPoint(n->GetX(), n->GetY(), n->GetZ(),
1824 TMath::Sqrt(n->GetSigmaY2()), TMath::Sqrt(n->GetSigmaZ2()));
1825 rieman2.Update();
afc28cd1 1826// printf(" Riemann results: row=%d valid=%d, Chi2=%.2f (%.2f) %d (%d)",
1827// iRow, rieman2.IsValid(), rieman2.GetChi2(), rieman1.GetChi2(), n->GetLabel(0),seed.GetLabel());
9d548f3c 1828 Double_t limit=2*rieman1.GetChi2();
1829 if (fClusterType==0) limit=1000;
1830 if (rieman2.GetChi2()>limit) {
1831 CopyRieman(rieman1,rieman2);
1832 ++noLastPoint;
afc28cd1 1833// printf("\n");
9d548f3c 1834 continue;
1835 }
afc28cd1 1836// printf(" +++ \n");
9d548f3c 1837
1838 noLastPoint=0;
1839 //use point
1840 rieman1.AddPoint(n->GetX(), n->GetY(), n->GetZ(),
1841 TMath::Sqrt(n->GetSigmaY2()), TMath::Sqrt(n->GetSigmaZ2()));
1842 rieman1.Update();
1843
1844 seed.SetClusterPointer(iRow,n);
1845 // if (iRow<seed.GetSeed1()) seed.SetSeed1(iRow);
1846 // if (iRow>seed.GetSeed2()) seed.SetSeed2(iRow);
1847 n->Use();
1848
1849 }
1850}
1851
1852//____________________________________________________________________________________
1853void AliToyMCReconstruction::ClusterToTrackAssociation(AliTPCseed &seed)
1854{
1855 //
1856 //
1857 //
1858
afc28cd1 1859// printf("\n ============ \nnext Seed: %d\n",seed.GetLabel());
9d548f3c 1860 //assume seed is within one sector
1861 Int_t iMiddle=(seed.GetSeed1()+seed.GetSeed2())/2;
1862 //outward
1863 AssociateClusters(seed,iMiddle+1,158,kFALSE);
1864 //inward
1865 AssociateClusters(seed,0,iMiddle,kTRUE);
1866 seed.SetIsSeeding(kFALSE);
1867
1868 CookLabel(&seed,.6);
1869}
1870
510cfcff 1871
d1cf83f5 1872//____________________________________________________________________________________
1873void AliToyMCReconstruction::InitSpaceCharge()
1874{
1875 //
1876 // Init the space charge map
1877 //
1878
0291b879 1879// TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
1880 TString filename;
d1cf83f5 1881 if (fTree) {
1882 TList *l=fTree->GetUserInfo();
1883 for (Int_t i=0; i<l->GetEntries(); ++i) {
2dd02504 1884 TObject *o=l->At(i);
1885 if (o->IsA() == TObjString::Class()){
1886 TString s(o->GetName());
0291b879 1887 if (s.Contains("lookup.root")) {
2dd02504 1888 filename=s;
1889 break;
1890 }
d1cf83f5 1891 }
1892 }
1893 }
1894
0291b879 1895 if (filename.IsNull()) {
1896 AliFatal("No SC map provided in the Userinfo of the simulation tree");
1897 return;
1898 }
1899
1900 AliInfo(Form("Initialising the space charge map using the file: '%s'\n",filename.Data()));
d1cf83f5 1901 TFile f(filename.Data());
223d9e38 1902 fTPCCorrection=(AliTPCCorrection*)f.Get("map");
d1cf83f5 1903
223d9e38 1904 // fTPCCorrection = new AliTPCSpaceCharge3D();
1905 // fTPCCorrection->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
1906 // fTPCCorrection->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1907 // // fTPCCorrection->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
1908 // fTPCCorrection->InitSpaceCharge3DDistortion();
d1cf83f5 1909
1910}
1911
1912//____________________________________________________________________________________
1913Double_t AliToyMCReconstruction::GetVDrift() const
1914{
1915 //
1916 //
1917 //
1918 return fTPCParam->GetDriftV();
1919}
1920
1921//____________________________________________________________________________________
1922Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
1923{
1924 //
1925 //
1926 //
1927 if (roc<0 || roc>71) return -1;
1928 return fTPCParam->GetZLength(roc);
1929}
1930
9e98dea8 1931//____________________________________________________________________________________
1932TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
1933 TString s=gSystem->GetFromPipe(Form("ls %s",files));
1934
1935 TTree *tFirst=0x0;
1936 TObjArray *arrFiles=s.Tokenize("\n");
1937
1938 for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
1939 TString name(arrFiles->At(ifile)->GetName());
1940
1941 TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
1942 TObjArray *arrMatch=0x0;
1943 arrMatch=reg.MatchS(name);
05da1b4e 1944 TString matchName;
1945 if (arrMatch && arrMatch->At(1)) matchName=arrMatch->At(1)->GetName();
1946 else matchName=Form("%02d",ifile);
1947 delete arrMatch;
9e98dea8 1948
1949 if (!tFirst) {
1950 TFile *f=TFile::Open(name.Data());
1951 if (!f) continue;
1952 TTree *t=(TTree*)f->Get("Tracks");
1953 if (!t) {
1954 delete f;
1955 continue;
1956 }
1957
05da1b4e 1958 t->SetName(matchName.Data());
9e98dea8 1959 tFirst=t;
1960 } else {
05da1b4e 1961 tFirst->AddFriend(Form("t%s=Tracks",matchName.Data()), name.Data());
9e98dea8 1962// tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
1963 }
1964 }
1965
c484014f 1966 if (tFirst->GetListOfFriends()) tFirst->GetListOfFriends()->Print();
9e98dea8 1967 return tFirst;
1968}
d1cf83f5 1969
0403120d 1970//____________________________________________________________________________________
5993ed4f 1971Int_t AliToyMCReconstruction::LoadOuterSectors() {
1972 //-----------------------------------------------------------------
1973 // This function fills outer TPC sectors with clusters.
1974 // Copy and paste from AliTPCtracker
1975 //-----------------------------------------------------------------
1976 Int_t nrows = fOuterSectorArray->GetNRows();
1977 UInt_t index=0;
1978 for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
1979 for (Int_t row = 0;row<nrows;row++){
1980 AliTPCtrackerRow* tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
1981 Int_t sec2 = sec+2*fkNSectorInner;
1982 //left
1983 Int_t ncl = tpcrow->GetN1();
1984 while (ncl--) {
1985 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1986 index=(((sec2<<8)+row)<<16)+ncl;
1987 tpcrow->InsertCluster(c,index);
1988 }
1989 //right
1990 ncl = tpcrow->GetN2();
1991 while (ncl--) {
1992 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1993 index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
1994 tpcrow->InsertCluster(c,index);
1995 }
1996 //
1997 // write indexes for fast acces
1998 //
1999 for (Int_t i=0;i<510;i++)
2000 tpcrow->SetFastCluster(i,-1);
2001 for (Int_t i=0;i<tpcrow->GetN();i++){
2002 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
2003 tpcrow->SetFastCluster(zi,i); // write index
2004 }
2005 Int_t last = 0;
2006 for (Int_t i=0;i<510;i++){
2007 if (tpcrow->GetFastCluster(i)<0)
2008 tpcrow->SetFastCluster(i,last);
2009 else
2010 last = tpcrow->GetFastCluster(i);
2011 }
2012 }
2013 return 0;
2014}
2015
2016
0403120d 2017//____________________________________________________________________________________
5993ed4f 2018Int_t AliToyMCReconstruction::LoadInnerSectors() {
2019 //-----------------------------------------------------------------
2020 // This function fills inner TPC sectors with clusters.
2021 // Copy and paste from AliTPCtracker
2022 //-----------------------------------------------------------------
2023 Int_t nrows = fInnerSectorArray->GetNRows();
2024 UInt_t index=0;
2025 for (Int_t sec = 0;sec<fkNSectorInner;sec++)
2026 for (Int_t row = 0;row<nrows;row++){
2027 AliTPCtrackerRow* tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
2028 //
2029 //left
2030 Int_t ncl = tpcrow->GetN1();
2031 while (ncl--) {
2032 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
2033 index=(((sec<<8)+row)<<16)+ncl;
2034 tpcrow->InsertCluster(c,index);
2035 }
2036 //right
2037 ncl = tpcrow->GetN2();
2038 while (ncl--) {
2039 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
2040 index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
2041 tpcrow->InsertCluster(c,index);
2042 }
2043 //
2044 // write indexes for fast acces
2045 //
2046 for (Int_t i=0;i<510;i++)
2047 tpcrow->SetFastCluster(i,-1);
2048 for (Int_t i=0;i<tpcrow->GetN();i++){
2049 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
2050 tpcrow->SetFastCluster(zi,i); // write index
2051 }
2052 Int_t last = 0;
2053 for (Int_t i=0;i<510;i++){
2054 if (tpcrow->GetFastCluster(i)<0)
2055 tpcrow->SetFastCluster(i,last);
2056 else
2057 last = tpcrow->GetFastCluster(i);
2058 }
2059
2060 }
2061 return 0;
2062}
c0ebc0e0 2063
0403120d 2064//____________________________________________________________________________________
c0ebc0e0 2065Int_t AliToyMCReconstruction::GetSector(AliExternalTrackParam *track) {
2066 //-----------------------------------------------------------------
2067 // This function returns the sector number for a given track
2068 //-----------------------------------------------------------------
2069
2070 Int_t sector = -1;
2071
2072 // get the sector number
2073 // rotate point to global system (track is already global!)
2074 Double_t xd[3];
2075 track->GetXYZ(xd);
2076 //track->Local2GlobalPosition(xd,track->GetAlpha());
2077
2078 // use TPCParams to get the sector number
2942f542 2079 Float_t xyz[3] = {static_cast<Float_t>(xd[0]),static_cast<Float_t>(xd[1]),static_cast<Float_t>(xd[2])};
c0ebc0e0 2080 Int_t i[3] = {0,0,0};
2081 if(fTPCParam){
2082 sector = fTPCParam->Transform0to1(xyz,i);
2083 }
2084
2085 return sector;
2086}
2087
0403120d 2088//____________________________________________________________________________________
c0ebc0e0 2089void AliToyMCReconstruction::FillSectorStructure(Int_t maxev) {
2090 //-----------------------------------------------------------------
2091 // This function fills the sector structure of AliToyMCReconstruction
2092 //-----------------------------------------------------------------
2093
2094 // cluster array of all sectors
2095 fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
2096 fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
2097
2098 for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
2099 for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
2100
2101 Int_t count[72][96] = { {0} , {0} };
2102
2103 for (Int_t iev=0; iev<maxev; ++iev){
2104 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
2105 fTree->GetEvent(iev);
2106 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
2107 printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
2108 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
2109
2110 // number of clusters to loop over
2111 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
2112
2113 for(Int_t icl=0; icl<ncls; ++icl){
2114
2115 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
2116 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
2117 if (!cl) continue;
2118
2119 Int_t sec = cl->GetDetector();
2120 Int_t row = cl->GetRow();
2121
0403120d 2122 // set Q of the cluster to 1, Q==0 does not work for the seeding
2123 cl->SetQ(1);
2124
64ff24a8 2125 // set cluster time to cluster Z (if not ideal tracking)
2126 if ( !fIdealTracking ) {
0403120d 2127 // a 'valid' position in z is needed for the seeding procedure
2128// cl->SetZ(cl->GetTimeBin()*GetVDrift());
2129 cl->SetZ(cl->GetTimeBin());
2130 }
c0ebc0e0 2131 //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
2132
2133 // fill arrays for inner and outer sectors (A/C side handled internally)
2134 if (sec<fkNSectorInner*2){
7f1a1b08 2135 fInnerSectorArray[sec%fkNSectorInner].InsertCluster(cl, count[sec][row], fTPCParam);
c0ebc0e0 2136 }
2137 else{
7f1a1b08 2138 fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(cl, count[sec][row], fTPCParam);
c0ebc0e0 2139 }
2140
2141 ++count[sec][row];
2142 }
2143 }
2144 }
2145
2146 // fill the arrays completely
64ea683c 2147 // LoadOuterSectors();
2148 // LoadInnerSectors();
2149
2150 // // check the arrays
2151 // for (Int_t i=0; i<fkNSectorInner; ++i){
2152 // for (Int_t j=0; j<fInnerSectorArray[i].GetNRows(); ++j){
2153 // if(fInnerSectorArray[i][j].GetN()>0){
2154 // Printf("Inner: Sector %d Row %d : %d",i,j,fInnerSectorArray[i][j].GetN());
2155 // }
2156 // }
2157 // }
2158 // for (Int_t i=0; i<fkNSectorInner; ++i){
2159 // for (Int_t j=0; j<fOuterSectorArray[i].GetNRows(); ++j){
2160 // if(fOuterSectorArray[i][j].GetN()>0){
2161 // Printf("Outer: Sector %d Row %d : %d",i,j,fOuterSectorArray[i][j].GetN());
2162 // }
2163 // }
2164 // }
c0ebc0e0 2165}
0403120d 2166
38d9d609 2167//____________________________________________________________________________________
a06336b6 2168void AliToyMCReconstruction::FillSectorStructureAC() {
38d9d609 2169 //-----------------------------------------------------------------
2170 // This function fills the sector structure of AliToyMCReconstruction
2171 //-----------------------------------------------------------------
2172
2173 /*
2174 my god is the AliTPCtrackerSector stuff complicated!!!
2175 Ok, so here we will not fill the fClusters1 and fClusters2 of AliTPCtrackerRow,
2176 using InsertCluster of AliTPCtrackerSector, but only the fClusters via InsertCluster
2177 of AliTPCtrackerRow itself which then will not be owner, but we create an array in
2178 here (fAllClusters) which owns all clusters ...
2179 */
2180
2181 fIsAC=kTRUE;
2182 // cluster array of all sectors
2183 fInnerSectorArray = new AliTPCtrackerSector[2*fkNSectorInner];
2184 fOuterSectorArray = new AliTPCtrackerSector[2*fkNSectorOuter];
2185
2186 for (Int_t i=0; i<2*fkNSectorInner; ++i) {
2187 fInnerSectorArray[i].Setup(fTPCParam,0);
2188 }
2189
2190 for (Int_t i=0; i<2*fkNSectorOuter; ++i) {
2191 fOuterSectorArray[i].Setup(fTPCParam,1);
2192 }
2193
2194 Int_t count[72][96] = { {0} , {0} };
2195
a06336b6 2196 for (Int_t iev=0; iev<fNmaxEvents; ++iev){
38d9d609 2197 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
2198 fTree->GetEvent(iev);
2199 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
a06336b6 2200// printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
38d9d609 2201 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
2202
2203 // number of clusters to loop over
2204 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
2205
2206 // check if expansion of the cluster arrays is needed.
2207 if (fAllClusters.GetEntriesFast()+ncls>=fAllClusters.Capacity()) fAllClusters.Expand(2*fAllClusters.Capacity());
2208 for(Int_t icl=0; icl<ncls; ++icl){
2209
2210 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
2211 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
2212 if (!cl) continue;
2213
2214 // register copy to the cluster array
2215 cl = new(fAllClusters[fAllClusters.GetEntriesFast()]) AliTPCclusterMI(*cl);
2216
2217 Int_t sec = cl->GetDetector();
2218 Int_t row = cl->GetRow();
2219
2220 // set Q of the cluster to 1, Q==0 does not work for the seeding
2221 cl->SetQ(1);
2222
2223 // set cluster time to cluster Z (if not ideal tracking)
2224 if ( !fIdealTracking ) {
2225 // a 'valid' position in z is needed for the seeding procedure
9d548f3c 2226 Double_t sign=1;
2227 if (((sec/18)%2)==1) sign=-1;
48e13991 2228 cl->SetZ(cl->GetTimeBin()*GetVDrift());
9d548f3c 2229 //mark cluster to be time*vDrift by setting the type to 1
2230 cl->SetType(1);
38d9d609 2231// cl->SetZ(cl->GetTimeBin());
2232 }
2233 //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
2234
2235 // fill arrays for inner and outer sectors (A/C side handled internally)
2236 if (sec<fkNSectorInner*2){
2237 fInnerSectorArray[sec][row].InsertCluster(cl, count[sec][row]);
2238 }
2239 else{
2240 fOuterSectorArray[(sec-fkNSectorInner*2)][row].InsertCluster(cl, count[sec][row]);
2241 }
2242
2243 ++count[sec][row];
2244 }
2245 }
2246 }
2247
2248}
2249
0403120d 2250//____________________________________________________________________________________
2251AliToyMCTrack *AliToyMCReconstruction::ConvertTPCSeedToToyMCTrack(const AliTPCseed &seed)
2252{
2253 //
2254 //
2255 //
2256
2257
2258 AliToyMCTrack *tToy = new AliToyMCTrack(seed);
2259
2260 for (Int_t icl=0; icl<159; ++icl){
cedbc66e 2261 const AliTPCclusterMI *cl=seed.GetClusterFast(icl);
2262 if (!cl) continue;
2263 if (fClusterType==0){
2264 tToy->AddSpacePoint(*cl);
2265 } else {
2266 tToy->AddDistortedSpacePoint(*cl);
2267 }
0403120d 2268 }
2269
2270 return tToy;
2271}
2272
2273//____________________________________________________________________________________
2274AliExternalTrackParam* AliToyMCReconstruction::GetRefittedTrack(const AliTPCseed &seed)
2275{
2276 //
2277 //
2278 //
2279
cedbc66e 2280 AliExternalTrackParam *track=0x0;
0403120d 2281
38d9d609 2282 const Float_t vDrift=GetVDrift();
cedbc66e 2283 const Float_t zLength=GetZLength(0);
2284 const Double_t kMaxSnp = 0.85;
2285 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2286
cedbc66e 2287 AliToyMCTrack *toyTrack = ConvertTPCSeedToToyMCTrack(seed);
9d548f3c 2288
2289 fTime0 = 0;
2290
2291 //get t0 estimate
2292 fCreateT0seed = kTRUE;
2293 AliExternalTrackParam *t0seed = GetSeedFromTrack(toyTrack,kTRUE);
2294 if (!t0seed) return 0x0;
cedbc66e 2295
9d548f3c 2296 fTime0 = t0seed->GetZ()-zLength/vDrift;
2297 delete t0seed;
2298 t0seed=0x0;
2299
cedbc66e 2300 fCreateT0seed = kFALSE;
9d548f3c 2301 AliExternalTrackParam *dummy = GetSeedFromTrack(toyTrack,kTRUE);
2302 track = GetFittedTrackFromSeed(toyTrack, dummy);
2303 delete dummy;
2304 // propagate seed to 0
2305 AliTrackerBase::PropagateTrackTo(track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
0403120d 2306
2307 return track;
2308}
502eb90b 2309
2310//____________________________________________________________________________________
38d9d609 2311AliTPCclusterMI* AliToyMCReconstruction::FindMiddleCluster(const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
2312 Double_t roady, Double_t roadz,
2313 AliRieman &seedFit)
2314{
2315 //
2316 //
2317 //
2318
2319 const Int_t rocInner = clInner->GetDetector();
2320 const Int_t rocOuter = clOuter->GetDetector();
2321
2322 if ( (rocInner%18) != (rocOuter%18) ){
2323 AliError("Currently only same Sector implemented");
2324 return 0x0;
2325 }
2326
2327 const Int_t innerDet=clInner->GetDetector();
2328 const Int_t outerDet=clOuter->GetDetector();
2329 const Int_t globalRowInner = clInner->GetRow() +(innerDet >35)*63;
2330 const Int_t globalRowOuter = clOuter->GetRow() +(outerDet >35)*63;
2331
9d548f3c 2332 AliTPCclusterMI *n=0x0;
2333
2334 // allow flexibility of +- nRowsGrace Rows to find a middle cluster
2335 const Int_t nRowsGrace = 0;
2336 for (Int_t iter=0; iter<2*nRowsGrace+1; ++iter){
2337 Int_t iMiddle = (globalRowInner+globalRowOuter)/2;
2338 iMiddle+=((iter+1)/2)*(1-2*((iter+1)%2));
38d9d609 2339
9d548f3c 2340 Int_t localRow=iMiddle;
2341 Int_t roc = innerDet;
2342 if (iMiddle>62){
2343 localRow-=63;
2344 roc=outerDet;
2345 }
2346
2347 AliTPCtrackerSector *arrRow=(iMiddle<63)?fInnerSectorArray:fOuterSectorArray;
2348 const AliTPCtrackerRow& krMiddle = arrRow[roc%36][localRow]; // middle
2349 // initial guess use simple linear interpolation
2350 Double_t y=(clInner->GetY()+clOuter->GetY())/2;
2351 Double_t z=(clInner->GetZ()+clOuter->GetZ())/2;
2352 if (seedFit.IsValid()) {
2353 // update values once the fit was performed
2354 y=seedFit.GetYat(krMiddle.GetX());
2355 z=seedFit.GetZat(krMiddle.GetX());
2356 }
38d9d609 2357
9d548f3c 2358 n=krMiddle.FindNearest(y,z,roady,roadz);
2359 if (n) break;
2360 }
2361
38d9d609 2362// if (n)
2363// printf(" Nearest cluster (%.2f, %.2f, %.2f) = m(%.2f, %.2f, %.2f : %d) i(%.2f, %.2f , %.2f : %d) o(%.2f, %.2f, %.2f : %d)\n",krMiddle.GetX(),y,z,
2364// n->GetX(), n->GetY(),n->GetZ(),n->GetLabel(0),
2365// clInner->GetX(), clInner->GetY(),clInner->GetZ(),clInner->GetLabel(0),
2366// clOuter->GetX(), clOuter->GetY(),clOuter->GetZ(),clOuter->GetLabel(0)
2367// );
2368 return n;
2369}
2370
2371//____________________________________________________________________________________
2372void AliToyMCReconstruction::AddMiddleClusters(AliTPCseed *seed,
2373 const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
2374 Double_t roady, Double_t roadz,
2375 Int_t &nTotalClusters, AliRieman &seedFit)
2376{
2377 //
2378 // Iteratively add the middle clusters
2379 //
2380
2381 // update rieman fit with every second point added
2382 AliTPCclusterMI *clMiddle=FindMiddleCluster(clInner,clOuter,roady,roadz,seedFit);
2383
2384 // break iterative process
9d548f3c 2385 if (!clMiddle || clMiddle->IsUsed()) return;
38d9d609 2386
2387 const Int_t globalRowInner = clInner->GetRow() +(clInner->GetDetector() >35)*63;
2388 const Int_t globalRowMiddle = clMiddle->GetRow()+(clMiddle->GetDetector()>35)*63;
2389 const Int_t globalRowOuter = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
2390
2391 seed->SetClusterPointer(globalRowMiddle,clMiddle);
2392 ++nTotalClusters;
2393// printf(" > Total clusters: %d\n",nTotalClusters);
2394 seedFit.AddPoint(clMiddle->GetX(), clMiddle->GetY(), clMiddle->GetZ(),
2395 TMath::Sqrt(clMiddle->GetSigmaY2()), TMath::Sqrt(clMiddle->GetSigmaZ2()));
2396
9d548f3c 2397 if (seedFit.GetN()>3) {
afc28cd1 2398// printf(" call update: %d (%d)\n",seedFit.GetN(),nTotalClusters);
2399// printf(" Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f -- %d\n",
2400// seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z(), clMiddle->GetLabel(0));
38d9d609 2401 seedFit.Update();
2402 }
2403 if ( seedFit.IsValid() && seedFit.GetChi2()>1000 ) return;
2404
2405 // Add middle clusters towards outer radius
2406 if (globalRowMiddle+1<globalRowOuter) AddMiddleClusters(seed,clMiddle,clOuter,roady,roadz,nTotalClusters,seedFit);
2407
2408 // Add middle clusters towards innter radius
2409 if (globalRowInner+1<globalRowMiddle) AddMiddleClusters(seed,clInner,clMiddle,roady,roadz,nTotalClusters,seedFit);
2410
2411 return;
2412}
2413
2414//____________________________________________________________________________________
9d548f3c 2415Int_t AliToyMCReconstruction::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t iRowInner, Int_t iRowOuter)
38d9d609 2416{
2417 //
2418 // find seeds in a sector, requires iRowInner < iRowOuter
2419 // iRowXXX runs from 0-159, so over IROC and OROC
2420 //
2421
2422 if (!fIsAC) {
2423 AliError("This function requires the sector arrays filled for AC tracking");
9d548f3c 2424 return 0;
38d9d609 2425 }
2426
2427 // swap rows in case they are in the wrong order
2428 if (iRowInner>iRowOuter) {
2429 Int_t tmp=iRowInner;
2430 iRowInner=iRowOuter;
2431 iRowOuter=tmp;
2432 }
2433
a06336b6 2434 if (iRowOuter>158) iRowOuter=158;
2435 if (iRowInner<0) iRowInner=0;
2436
38d9d609 2437 // Define the search roads:
2438 // timeRoadCombinatorics - the maximum time difference used for the
2439 // combinatorics. Since we don't have a 'z-Vertex' estimate this will
2440 // reduce the combinatorics significantly when iterating on one TF
2441 // use a little more than one full drift length of the TPC to allow for
2442 // distortions
2443 //
2444 // timeRoad - the time difference allowed when associating the cluster
2445 // in the middle of the other two used for the initial search
2446 //
2447 // padRoad - the local y difference allowed when associating the middle cluster
38d9d609 2448
2449// Double_t timeRoadCombinatorics = 270./vDrift;
2450// Double_t timeRoad = 20./vDrift;
2451 Double_t timeRoadCombinatorics = 270.;
9d548f3c 2452 Double_t timeRoad = 10.;
2453 Double_t padRoad = 5.;
38d9d609 2454
2455
2456 // fOuterSectorArray runs from 0-95, so from iRowXXX which runs from 0-159
2457 // the number of rows in the IROC has to be subtracted
2458 const Int_t innerRows=fInnerSectorArray->GetNRows();
9d548f3c 2459
2460 AliTPCtrackerSector *arrInnerRow=(iRowInner<63)?fInnerSectorArray:fOuterSectorArray;
2461 AliTPCtrackerSector *arrOuterRow=(iRowOuter<63)?fInnerSectorArray:fOuterSectorArray;
2462
2463 const AliTPCtrackerRow& krInner = arrInnerRow[sec][iRowInner - (iRowInner>62)*innerRows]; // up
2464 const AliTPCtrackerRow& krOuter = arrOuterRow[sec][iRowOuter - (iRowOuter>62)*innerRows]; // down
38d9d609 2465
2466 AliTPCseed *seed = new AliTPCseed;
2467
2468 const Int_t nMaxClusters=iRowOuter-iRowInner+1;
2469// Int_t nScannedClusters = 0;
9d548f3c 2470
2471 Int_t nseedAdded=0;
38d9d609 2472 // loop over all points in the firstand last search row
2473 for (Int_t iOuter=0; iOuter < krOuter; iOuter++) {
2474 const AliTPCclusterMI *clOuter = krOuter[iOuter];
a06336b6 2475 if (clOuter->IsUsed()) continue;
38d9d609 2476 for (Int_t iInner=0; iInner < krInner; iInner++) {
2477 const AliTPCclusterMI *clInner = krInner[iInner];
a06336b6 2478 if (clInner->IsUsed()) continue;
afc28cd1 2479// printf("\n\n Check combination %d (%d), %d (%d) -- %d (%d) -- %d\n",iOuter, iInner, clOuter->GetLabel(0), clInner->GetLabel(0),iRowOuter,iRowInner,sec);
38d9d609 2480 // check maximum distance for combinatorics
2481 if (TMath::Abs(clOuter->GetZ()-clInner->GetZ())>timeRoadCombinatorics) continue;
afc28cd1 2482// printf(" Is inside one drift\n");
38d9d609 2483
2484 // use rieman fit for seed description
2485 AliRieman seedFit(159);
2486 // add first two points
2487 seedFit.AddPoint(clInner->GetX(), clInner->GetY(), clInner->GetZ(),
2488 TMath::Sqrt(clInner->GetSigmaY2()), TMath::Sqrt(clInner->GetSigmaZ2()));
2489 seedFit.AddPoint(clOuter->GetX(), clOuter->GetY(), clOuter->GetZ(),
2490 TMath::Sqrt(clOuter->GetSigmaY2()), TMath::Sqrt(clOuter->GetSigmaZ2()));
2491
2492 // Iteratively add all clusters in the respective middle
2493 Int_t nFoundClusters=2;
2494 AddMiddleClusters(seed,clInner,clOuter,padRoad,timeRoad,nFoundClusters,seedFit);
afc28cd1 2495// printf(" Clusters attached: %d\n",nFoundClusters);
a06336b6 2496 if (nFoundClusters>2) seedFit.Update();
afc28cd1 2497// printf(" Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f\n",
2498// seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z());
38d9d609 2499
2500 // check for minimum number of assigned clusters and a decent chi2
2501 if ( nFoundClusters<0.5*nMaxClusters || seedFit.GetChi2()>1000 ){
2502 seed->Reset();
2503 continue;
2504 }
2505// printf(" >>> Seed accepted\n");
2506 // add original outer clusters
2507 const Int_t globalRowInner = clInner->GetRow() +(clInner->GetDetector() >35)*63;
2508 const Int_t globalRowOuter = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
2509
2510 seed->SetClusterPointer(globalRowInner,const_cast<AliTPCclusterMI*>(clInner));
2511 seed->SetClusterPointer(globalRowOuter,const_cast<AliTPCclusterMI*>(clOuter));
2512
2513 // set parameters retrieved from AliRieman
2514 Double_t params[5]={0};
2515 Double_t covar[15]={0};
2516 const Double_t alpha=TMath::DegToRad()*(clInner->GetDetector()%18*20.+10.);
2517 const Double_t x=clInner->GetX();
2518 seedFit.GetExternalParameters(x,params,covar);
2519
2520 seed->Set(x,alpha,params,covar);
2521
2522 // set label of the seed. At least 60% of the clusters need the correct label
a06336b6 2523 CookLabel(seed,.6);
2524// printf(" - Label: %d\n",seed->GetLabel());
2525 // mark clusters as being used
2526 MarkClustersUsed(seed);
9d548f3c 2527
2528 seed->SetSeed1(iRowInner);
2529 seed->SetSeed2(iRowOuter);
2530 seed->SetSector(sec);
2531 ++nseedAdded;
2532
2533 seed->SetUniqueID(arr->GetEntriesFast());
2534 seed->SetIsSeeding(kTRUE);
38d9d609 2535
2536 arr->Add(seed);
2537 seed=new AliTPCseed;
9d548f3c 2538
2539 break;
38d9d609 2540 }
2541 }
2542 //delete surplus seed
2543 delete seed;
2544 seed=0x0;
2545
9d548f3c 2546 return nseedAdded;
38d9d609 2547}
2548//____________________________________________________________________________________
2549void AliToyMCReconstruction::MakeSeeds(TObjArray * /*arr*/, Int_t sec, Int_t iRow1, Int_t iRow2)
502eb90b 2550{
2551 //
2552 // Create seeds between i1 and i2 (stored in arr) with TLinearFitter
2553 //
2554 // sec: sector number
2555 // iRow1: upper row
2556 // iRow2: lower row
2557 //
2558
2559 // Create Fitter
2560 static TLinearFitter fitter(3,"pol2");
2561
2562 // Get 3 padrows (iRow1,iMiddle=(iRow1+iRow2)/2,iRow2)
2563 Int_t iMiddle = (iRow1+iRow2)/2;
2564 const AliTPCtrackerRow& krd = fOuterSectorArray[sec][iRow2-fInnerSectorArray->GetNRows()]; // down
2565 const AliTPCtrackerRow& krm = fOuterSectorArray[sec][iMiddle-fInnerSectorArray->GetNRows()]; // middle
2566 const AliTPCtrackerRow& kru = fOuterSectorArray[sec][iRow1-fInnerSectorArray->GetNRows()]; // up
2567
2568 // Loop over 3 cluster possibilities
2569 for (Int_t iu=0; iu < kru; iu++) {
2570 for (Int_t im=0; im < krm; im++) {
2571 for (Int_t id=0; id < krd; id++) {
2572
2573 // clear all points
2574 fitter.ClearPoints();
2575
2576 // add all three points to fitter
2577 Double_t xy[2] = {kru[iu]->GetX(),kru[iu]->GetY()};
2578 Double_t z = kru[iu]->GetZ();
2579 fitter.AddPoint(xy,z);
2580
2581 xy[0] = krm[im]->GetX();
2582 xy[1] = krm[im]->GetY();
2583 z = krm[im]->GetZ();
2584 fitter.AddPoint(xy,z);
2585
2586 xy[0] = krd[id]->GetX();
2587 xy[1] = krd[id]->GetY();
2588 z = krd[id]->GetZ();
2589 fitter.AddPoint(xy,z);
2590
2591 // Evaluate and get parameters
2592 fitter.Eval();
2593
2594 // how to get the other clusters now?
2595 // ...
2596
2597 }
2598 }
2599 }
2600}
38d9d609 2601
2602//____________________________________________________________________________________
a06336b6 2603void AliToyMCReconstruction::InitStreamer(TString addName, Int_t level)
38d9d609 2604{
2605 //
2606 // initilise the debug streamer and set the logging level
2607 // use a default naming convention
2608 //
a06336b6 2609
2610 delete fStreamer;
2611 fStreamer=0x0;
2612
2613 if (addName.IsNull()) addName=".dummy";
38d9d609 2614
a06336b6 2615 if (!fTree) return;
2616
48e13991 2617 TString debugName=gSystem->BaseName(fInputFile->GetName());
a06336b6 2618 debugName.ReplaceAll(".root","");
2619 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
2620 fUseMaterial,fIdealTracking,fClusterType,
2621 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
2622 debugName.Append(addName);
2623 debugName.Append(".root");
2624
2625 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
2626 fStreamer=new TTreeSRedirector(debugName.Data());
2627 fStreamer->SetUniqueID(level);
2628
2629 gROOT->cd();
38d9d609 2630}
2631
2632//____________________________________________________________________________________
a06336b6 2633void AliToyMCReconstruction::ConnectInputFile(const char* file, Int_t nmaxEv)
38d9d609 2634{
2635 //
2636 // connect the tree and event pointer from the input file
2637 //
2638
2639 delete fInputFile;
2640 fInputFile=0x0;
2641 fEvent=0x0;
2642 fTree=0;
2643
2644 fInputFile=TFile::Open(file);
2645 if (!fInputFile || !fInputFile->IsOpen() || fInputFile->IsZombie()) {
2646 delete fInputFile;
2647 fInputFile=0x0;
afc28cd1 2648 AliError(Form("ERROR: couldn't open the file '%s'\n", file));
38d9d609 2649 return;
2650 }
2651
2652 fTree=(TTree*)fInputFile->Get("toyMCtree");
2653 if (!fTree) {
afc28cd1 2654 AliError(Form("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file));
38d9d609 2655 return;
2656 }
2657
2658 fEvent=0x0;
2659 fTree->SetBranchAddress("event",&fEvent);
2660
2661 gROOT->cd();
a06336b6 2662
2663 fNmaxEvents=fTree->GetEntries();
2664 if (nmaxEv>-1) fNmaxEvents=TMath::Min(nmaxEv,fNmaxEvents);
38d9d609 2665
a06336b6 2666 // setup space charge map from Userinfo of the tree
2667 InitSpaceCharge();
2668
2669 // setup the track maps
38d9d609 2670 SetupTrackMaps();
2671}
2672
2673//____________________________________________________________________________________
2674void AliToyMCReconstruction::Cleanup()
2675{
2676 //
2677 // Cleanup input data
2678 //
2679
a06336b6 2680 if (fStreamer) delete fStreamer;
38d9d609 2681 fStreamer=0x0;
2682
2683 delete fEvent;
2684 fEvent = 0x0;
5414766f 2685
2686 delete fHnDelta;
2687 fHnDelta=0x0;
a06336b6 2688// delete fTree;
38d9d609 2689 fTree=0x0;
2690
2691 delete fInputFile;
2692 fInputFile=0x0;
2693
2694}
2695
2696//____________________________________________________________________________________
2697void AliToyMCReconstruction::SetupTrackMaps()
2698{
2699 //
2700 //
2701 //
2702
2703 fMapTrackEvent.Delete();
2704 fMapTrackTrackInEvent.Delete();
2705
2706 if (!fTree) {
2707 AliError("Tree not connected");
2708 return;
2709 }
2710
a06336b6 2711 Int_t nmaxEv=fTree->GetEntries();
2712 if (fNmaxEvents>-1) nmaxEv=fNmaxEvents;
2713
2714 for (Int_t iev=0; iev<nmaxEv; ++iev) {
38d9d609 2715 fTree->GetEvent(iev);
2716 if (!fEvent) continue;
2717
2718 const Int_t ntracks=fEvent->GetNumberOfTracks();
2719 if (fMapTrackEvent.GetSize()+ntracks>=fMapTrackEvent.Capacity()) fMapTrackEvent.Expand(2*fMapTrackEvent.Capacity());
2720 if (fMapTrackTrackInEvent.GetSize()+ntracks>=fMapTrackTrackInEvent.Capacity()) fMapTrackTrackInEvent.Expand(2*fMapTrackTrackInEvent.Capacity());
2721
2722 for (Int_t itrack=0; itrack<ntracks; ++itrack){
2723 Int_t label=fEvent->GetTrack(itrack)->GetUniqueID();
2724
2725 fMapTrackEvent.Add(label,iev);
2726 fMapTrackTrackInEvent.Add(label,itrack);
2727 }
2728 }
2729
2730}
a06336b6 2731
2732//____________________________________________________________________________________
2733void AliToyMCReconstruction::CookLabel(AliTPCseed *seed, Double_t fraction, Int_t info[5])
2734{
2735 //
2736 //
2737 //
2738
2739 Int_t labels[159]={0};
2740// Long64_t posMaxLabel=-1;
2741 Int_t nMaxLabel=0; // clusters from maximum label
2742 Int_t nMaxLabel2=0; // clusters from second maximum
2743 Int_t nlabels=0;
2744 Int_t maxLabel=-1; // label with most clusters
2745 Int_t maxLabel2=-1; // label with second most clusters
2746 Int_t nclusters=0;
2747 TExMap labelMap(159);
2748
2749 for (Int_t icl=0; icl<159; ++icl) {
2750 const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2751 if (!cl) continue;
2752 ++nclusters;
2753
2754 const Int_t clLabel=cl->GetLabel(0);
2755 // a not assinged value returns 0, so we need to add 1 and subtract it afterwards
2756 Long64_t labelPos=labelMap.GetValue(clLabel);
2757
2758 if (!labelPos) {
2759 labelPos=nlabels+1;
2760 labelMap.Add(clLabel,labelPos);
2761 ++nlabels;
2762 }
2763 --labelPos;
2764
2765 const Int_t nCurrentLabel=++labels[labelPos];
2766 if (nCurrentLabel>nMaxLabel) {
2767 nMaxLabel2=nMaxLabel;
2768 nMaxLabel=nCurrentLabel;
2769// posMaxLabel=labelPos;
2770 maxLabel2=maxLabel;
2771 maxLabel=clLabel;
2772 }
2773 }
2774
2775 if (Double_t(nMaxLabel)/nclusters<fraction) maxLabel=-maxLabel;
2776
2777 seed->SetLabel(maxLabel);
2778
2779 if (info) {
2780 info[0]=nMaxLabel;
2781 info[1]=nMaxLabel2;
2782 info[2]=maxLabel2;
2783 info[3]=nclusters;
2784 info[4]=nlabels;
2785 }
2786}
2787
2788
2789//____________________________________________________________________________________
9d548f3c 2790void AliToyMCReconstruction::DumpSeedInfo(TObjArray *arr)
a06336b6 2791{
2792
2793 // for debugging
2794 if (!fStreamer || !fTree) return;
2795 // swap rows in case they are in the wrong order
9d548f3c 2796 AliInfo(Form("Number of Seeds: %d (out of %d tracks)",arr->GetEntriesFast(),fMapTrackEvent.GetSize()));
a06336b6 2797
2798 //loop over all events and tracks and try to associate the seed to the track
2799 for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed){
2800 AliTPCseed *seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
a06336b6 2801
2802 // get original track
9d548f3c 2803 Int_t seedLabel=seed->GetLabel();
a06336b6 2804 Int_t iev=fMapTrackEvent.GetValue(TMath::Abs(seedLabel));
2805 Int_t itr=fMapTrackTrackInEvent.GetValue(TMath::Abs(seedLabel));
2806
2807 fTree->GetEvent(iev);
2808 const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
2809
9d548f3c 2810 DumpSeedInfo(toyTrack,seed);
2811 }
2812}
2813
2814//____________________________________________________________________________________
2815void AliToyMCReconstruction::DumpTrackInfo(TObjArray *arr)
2816{
2817
2818 // for debugging
2819 if (!fStreamer || !fTree) return;
2820 // swap rows in case they are in the wrong order
2821 AliInfo(Form("Number of Seeds: %d (out of %d tracks)",arr->GetEntriesFast(),fMapTrackEvent.GetSize()));
2822
2823 //loop over all events and tracks and try to associate the seed to the track
2824 AliTPCseed dummySeed;
2825 for (Int_t iev=0; iev<fNmaxEvents; ++iev) {
2826 fTree->GetEvent(iev);
2827 const Int_t ntracks=fEvent->GetNumberOfTracks();
2828 for (Int_t itr=0; itr<ntracks; ++itr) {
2829 const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
2830
2831 Bool_t foundSeed=kFALSE;
2832 for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed) {
2833 AliTPCseed *seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
2834 const UInt_t tmpLabel=TMath::Abs(seed->GetLabel());
2835 if (toyTrack->GetUniqueID()!=tmpLabel) continue;
2836
2837 // dump all seeds with the correct label
2838 DumpSeedInfo(toyTrack,seed);
2839 foundSeed=kTRUE;
2840 }
2841
2842 if (!foundSeed) DumpSeedInfo(toyTrack,&dummySeed);
a06336b6 2843
a06336b6 2844 }
9d548f3c 2845 }
2846}
a06336b6 2847
9d548f3c 2848//____________________________________________________________________________________
2849void AliToyMCReconstruction::DumpSeedInfo(const AliToyMCTrack *toyTrack, AliTPCseed *seed)
2850{
2851 //
2852 //
2853 //
2854
2855 const Double_t kMaxSnp = 0.85;
2856 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2857 Float_t vDrift=GetVDrift();
2858 Float_t zLength=GetZLength(0);
2859
2860 AliExternalTrackParam tOrig(*toyTrack);
2861 AliExternalTrackParam tOrigITS(*toyTrack);
2862
2863 // propagate original track to ITS last layer
e83fd282 2864// Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
2865 const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
2866 Double_t lastLayerITS = iFCRadius; // its track propgated to inner TPC wall
9d548f3c 2867 AliTrackerBase::PropagateTrackTo(&tOrigITS,lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2868
2869 AliExternalTrackParam dummyParam;
2870 Bool_t isDummy=kFALSE;
2871 //create refitted track, this will also set the fTime0
2872 AliExternalTrackParam *track=GetRefittedTrack(*seed);
2873 if (!track) {
2874 track=&dummyParam;
2875 isDummy=kTRUE;
2876 }
2877 AliTrackerBase::PropagateTrackTo(track,0,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2878 track->Rotate(tOrig.GetAlpha());
2879 AliTrackerBase::PropagateTrackTo(track,0,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2880
2881 // rotate fitted track to the frame of the original track and propagate to same reference
2882 AliExternalTrackParam trackITS(*track);
2883 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2884 trackITS.Rotate(tOrigITS.GetAlpha());
2885 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2886
2887 Int_t seedSec=seed->GetSector();
2888 Int_t seedID =seed->GetUniqueID();
2889 //
2890 //count findable and found clusters in the seed
2891 //
2892 Int_t iRowInner=seed->GetSeed1();
2893 Int_t iRowOuter=seed->GetSeed2();
2894
2895 Int_t nClustersSeedMax=iRowOuter-iRowInner+1;
2896 Int_t nClustersFindable=0;
2897 Int_t nClustersSeed=0;
2898
2899 Int_t ncls=(fClusterType==0)?toyTrack->GetNumberOfSpacePoints():toyTrack->GetNumberOfDistSpacePoints();
2900
2901 Int_t rowInner=iRowInner-(iRowInner>62)*63;
2902 Int_t rowOuter=iRowOuter-(iRowOuter>62)*63;
2903
2904 //findable in the current seed sector
2905 Int_t lastROC=-1;
2906 Int_t rocMaxCl=-1;
2907 Int_t nCrossedROCs=0;
2908 Int_t nMaxClROC=0;
2909 Int_t nclROC=0;
2910 Int_t row1=-1;
2911 Int_t row2=-1;
2912 Int_t row1Maxcl=-1;
2913 Int_t row2Maxcl=-1;
2914 for (Int_t icl=0; icl<ncls; ++icl) {
2915 const AliTPCclusterMI *cl=(fClusterType==0)?toyTrack->GetSpacePoint(icl):toyTrack->GetDistortedSpacePoint(icl);
2916 const Int_t roc=cl->GetDetector();
2917 const Int_t row=cl->GetRow();
2918 const Int_t rowGlobal=row+(roc>35)*63;
2919
2920 AliTPCclusterMI *seedCl = seed->GetClusterPointer(rowGlobal);
2921 if (seedCl) {
2922 AliTPCclusterMI *clc=const_cast<AliTPCclusterMI*>(cl);
2923 if (seedCl->GetDetector()==roc&&seedCl->IsUsed()) clc->Use();
2924 clc->SetLabel(seedID,1);
2925
a06336b6 2926 }
9d548f3c 2927
2928// if (row1<0) row1=rowGlobal;
2929
2930 if ( (roc%36) != (lastROC%36)) {
2931 ++nCrossedROCs;
2932 if (nclROC>nMaxClROC) {
2933 nMaxClROC=nclROC;
2934 rocMaxCl=lastROC;
2935 row1Maxcl=row1;
2936 row2Maxcl=row2;
2937 }
2938 lastROC=roc%36;
2939 nclROC=0;
2940 row1=rowGlobal;
2941 }
2942 ++nclROC;
2943 row2=rowGlobal;
a06336b6 2944
9d548f3c 2945 if ( (roc%36)!=(seedSec%36) ) continue;
2946// if ( (row<rowInner) || (row>rowOuter) ) continue;
2947 ++nClustersFindable;
2948
a06336b6 2949 }
9d548f3c 2950
2951 if (nclROC>nMaxClROC) {
2952 rocMaxCl=lastROC;
2953 nMaxClROC=nclROC;
2954 row1Maxcl=row1;
2955 row2Maxcl=row2;
2956 }
2957
2958 Int_t firstRow=160;
2959 Int_t lastRow=0;
2960 Int_t nClustersInTrack=0;
2961 //found in seed
2962 for (Int_t icl=0; icl<159; ++icl) {
2963 const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2964 if (!cl) continue;
2965 ++nClustersInTrack;
2966 const Int_t row=cl->GetRow();
2967 const Int_t rowGlobal=row+(cl->GetDetector()>35)*63;
2968 if (rowGlobal>lastRow) lastRow=rowGlobal;
2969 if (rowGlobal<firstRow) firstRow=rowGlobal;
2970 if ( (row<rowInner) || (row>rowOuter) ) continue;
2971 ++nClustersSeed;
2972 }
2973
2974 Float_t z0=fEvent->GetZ();
2975 Float_t t0=fEvent->GetT0();
2976
2977 Int_t ctype(fCorrectionType);
2978
2979 Int_t info[5]={0};
2980 CookLabel(seed,.6,info);
2981 Int_t seedLabel=seed->GetLabel();
2982
2983 Int_t labelOrig=toyTrack->GetUniqueID();
2984
2985 AliToyMCTrack *tr2 = const_cast<AliToyMCTrack*>(toyTrack);
2986
2987 (*fStreamer) << "Seeds" <<
2988 // "iev=" << iev <<
2989 // "iseed=" << iseed <<
2990 // "itr=" << itr <<
2991 "z0=" << z0 <<
2992 "t0=" << t0 <<
2993 "vDrift=" << vDrift <<
2994 "zLength=" << zLength <<
2995 "fTime0=" << fTime0 <<
2996 "clsType=" << fClusterType <<
2997 "corrType=" << ctype <<
2998
2999 "tOrig.=" << &tOrig <<
3000 "tOrigITS.=" << &tOrigITS <<
3001
3002 "to.nclFindable=" << nClustersFindable <<
3003 "to.nclTot=" << ncls <<
3004 "to.label=" << labelOrig <<
3005 "to.nCrossedROCs="<< nCrossedROCs <<
3006 "to.rocMax=" << rocMaxCl <<
3007 "to.rocMaxNcl=" << nMaxClROC <<
3008 "to.row1Max=" << row1Maxcl <<
3009 "to.row2Max=" << row2Maxcl <<
3010
3011 "track.=" << track <<
3012 "trackITS.=" << &trackITS <<
3013
3014 "s.RowInner=" << iRowInner <<
3015 "s.RowOuter=" << iRowOuter <<
3016 "s.nclMax=" << nClustersSeedMax <<
3017 "s.ncl=" << nClustersSeed <<
3018 "s.ID=" << seedID <<
3019
3020 "tr.firstClRow=" << firstRow <<
3021 "tr.lastClRow=" << lastRow <<
3022 "tr.ncl=" << nClustersInTrack <<
3023 "tr.label=" << seedLabel <<
3024
3025 "tr.LabelNcl=" << info[0] <<
3026 "tr.Label2Ncl=" << info[1] <<
3027 "tr.Label2=" << info[2] <<
3028 "tr.nclTot=" << info[3] <<
3029 "tr.Nlabels=" << info[4] <<
3030
3031 "tr.Sec=" << seedSec <<
3032
3033 "seed.=" << seed <<
3034 "toyTrack.=" << tr2 <<
3035 "\n";
3036
3037 if (!isDummy) delete track;
a06336b6 3038}
3039
3040//____________________________________________________________________________________
3041void AliToyMCReconstruction::MarkClustersUsed(AliTPCseed *seed)
3042{
3043 //
3044 //
3045 //
3046
3047 for (Int_t icl=0; icl<159; ++icl) {
3048 AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
3049 if (cl) cl->Use();
3050 }
3051}
3052
9d548f3c 3053//____________________________________________________________________________________
3054void AliToyMCReconstruction::ResetClustersZtoTime(AliTPCseed *seed)
3055{
3056 //
3057 //
3058 //
3059
3060 for (Int_t icl=0; icl<159; ++icl) {
3061 AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
3062 if (cl) cl->SetZ(cl->GetTimeBin());
3063 }
3064}
3065
a06336b6 3066//____________________________________________________________________________________
3067void AliToyMCReconstruction::DumpTracksToTree(const char* file)
3068{
3069 //
3070 //
3071 //
3072 ConnectInputFile(file);
3073 if (!fTree) return;
3074
3075 delete fStreamer;
3076 fStreamer=0x0;
3077
3078 TString debugName=fInputFile->GetName();
3079 debugName.ReplaceAll(".root",".AllTracks.root");
3080
3081 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
3082 fStreamer=new TTreeSRedirector(debugName.Data());
3083
3084 for (Int_t iev=0;iev<fNmaxEvents;++iev){
3085 fTree->GetEvent(iev);
3086 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks();++itr){
3087 AliToyMCTrack *toyTrack=const_cast<AliToyMCTrack*>(fEvent->GetTrack(itr));
3088 Int_t trackID=toyTrack->GetUniqueID();
3089
3090 (*fStreamer) << "Tracks" <<
3091 "iev=" << iev <<
3092 "itr=" << itr <<
3093 "trackID=" << trackID <<
3094 "track.=" << toyTrack <<
3095 "\n";
3096
3097 }
3098 }
3099
3100 Cleanup();
3101}
9d548f3c 3102
3103//____________________________________________________________________________________
3104void AliToyMCReconstruction::SetRieman(const AliTPCseed &seed, AliRieman &rieman)
3105{
3106 //
3107 //
3108 //
3109
3110 rieman.Reset();
3111 for (Int_t icl=0; icl<159; ++icl) {
3112 const AliTPCclusterMI *cl=seed.GetClusterPointer(icl);
3113 if (!cl) continue;
3114 rieman.AddPoint(cl->GetX(), cl->GetY(), cl->GetZ(),
3115 TMath::Sqrt(cl->GetSigmaY2()), TMath::Sqrt(cl->GetSigmaZ2()));
3116 }
3117 rieman.Update();
3118}
3119
3120//____________________________________________________________________________________
3121void AliToyMCReconstruction::CopyRieman(const AliRieman &from, AliRieman &to)
3122{
3123 //
3124 //
3125 //
3126
3127 if (to.GetCapacity()<from.GetCapacity()) return;
3128 to.Reset();
3129
3130 for (Int_t i=0;i<from.GetN();++i) to.AddPoint(from.GetX()[i],from.GetY()[i],from.GetZ()[i],from.GetSy()[i],from.GetSz()[i]);
3131}
3132
36a79b0d 3133//____________________________________________________________________________________
f45d00f9 3134Float_t AliToyMCReconstruction::FindClosestT0(const TVectorF &t0list, const TVectorF &z0list, AliExternalTrackParam &t0seed)
36a79b0d 3135{
3136 //
3137 // find closes T0 in a list of T0s
3138 //
3139
3140 Long64_t size=t0list.GetNrows();
3141 const Float_t *array=t0list.GetMatrixArray();
3142
f45d00f9 3143 Float_t vDrift=GetVDrift();
3144 Float_t zLength=GetZLength(0);
3145
3146 Float_t sign=(1-2*(t0seed.GetTgl()>0));
3147
3148 Float_t vtxCorr=0.;
3149 Float_t t0=t0seed.GetZ()-zLength/vDrift;
3150
36a79b0d 3151 Int_t index=0;
3152 Float_t minDist=1e20;
3153 for (Int_t it0=0; it0<size; ++it0) {
f45d00f9 3154 if (fUseZ0list) vtxCorr=sign*z0list[it0]/vDrift;
61d165ed 3155// printf("vtxcorr %d: %.2g, %.2g\n",it0, vtxCorr, z0list[it0]);
f45d00f9 3156 const Float_t dist=TMath::Abs(t0list[it0]-t0-vtxCorr);
36a79b0d 3157 if (dist<minDist) {
3158 index=it0;
3159 minDist=dist;
3160 }
3161 }
3162
3163 return array[index];
3164}