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