]> 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
larger canvases for better pngs
[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)
d1cf83f5 49, fStreamer(0x0)
38d9d609 50, fInputFile(0x0)
d1cf83f5 51, fTree(0x0)
4a777885 52, fEvent(0x0)
d1cf83f5 53, fTPCParam(0x0)
223d9e38 54, fTPCCorrection(0x0)
5993ed4f 55, fkNSectorInner(18) // hard-coded to avoid loading the parameters before
56, fInnerSectorArray(0x0)
57, fkNSectorOuter(18) // hard-coded to avoid loading the parameters before
58, fOuterSectorArray(0x0)
38d9d609 59, fAllClusters("AliTPCclusterMI",10000)
60, fMapTrackEvent(10000)
61, fMapTrackTrackInEvent(10000)
62, fIsAC(kFALSE)
d1cf83f5 63{
64 //
65 // ctor
66 //
67 fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
5993ed4f 68
d1cf83f5 69}
70
71//____________________________________________________________________________________
72AliToyMCReconstruction::~AliToyMCReconstruction()
73{
74 //
75 // dtor
76 //
77
38d9d609 78 Cleanup();
4a777885 79}
80
81//____________________________________________________________________________________
82void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
83{
84 //
510cfcff 85 // Recostruction from associated clusters
4a777885 86 //
87
a06336b6 88 ConnectInputFile(file, nmaxEv);
38d9d609 89 if (!fTree) return;
a06336b6 90
91 InitStreamer(".debug");
4a777885 92
93 gROOT->cd();
94
95 static AliExternalTrackParam dummySeedT0;
96 static AliExternalTrackParam dummySeed;
97 static AliExternalTrackParam dummyTrack;
98
99 AliExternalTrackParam t0seed;
100 AliExternalTrackParam seed;
101 AliExternalTrackParam track;
d5704a39 102 AliExternalTrackParam trackITS;
4a777885 103 AliExternalTrackParam tOrig;
d5704a39 104 AliExternalTrackParam tOrigITS;
4a777885 105
106 AliExternalTrackParam *dummy;
107
108 Int_t maxev=fTree->GetEntries();
109 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
110
111 for (Int_t iev=0; iev<maxev; ++iev){
112 printf("============== Processing Event %6d =================\n",iev);
113 fTree->GetEvent(iev);
114 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
0403120d 115// printf(" > ====== Processing Track %6d ======== \n",itr);
4a777885 116 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
117 tOrig = *tr;
118
119
120 // set dummy
121 t0seed = dummySeedT0;
122 seed = dummySeed;
123 track = dummyTrack;
124
125 Float_t z0=fEvent->GetZ();
126 Float_t t0=fEvent->GetT0();
127
38d9d609 128 Float_t vDrift=GetVDrift();
4a777885 129 Float_t zLength=GetZLength(0);
130
600eaa0d 131 // crate time0 seed, steered by fCreateT0seed
05da1b4e 132// printf("t0 seed\n");
4a777885 133 fTime0=-1.;
600eaa0d 134 fCreateT0seed=kTRUE;
4a777885 135 dummy = GetSeedFromTrack(tr);
136
137 if (dummy) {
138 t0seed = *dummy;
139 delete dummy;
140
141 // crate real seed using the time 0 from the first seed
600eaa0d 142 // set fCreateT0seed now to false to get the seed in z coordinates
38d9d609 143 fTime0 = t0seed.GetZ()-zLength/vDrift;
600eaa0d 144 fCreateT0seed = kFALSE;
05da1b4e 145// printf("seed (%.2g)\n",fTime0);
4a777885 146 dummy = GetSeedFromTrack(tr);
147 if (dummy) {
148 seed = *dummy;
149 delete dummy;
150
151 // create fitted track
152 if (fDoTrackFit){
05da1b4e 153// printf("track\n");
4a777885 154 dummy = GetFittedTrackFromSeed(tr, &seed);
155 track = *dummy;
156 delete dummy;
157 }
600eaa0d 158
d5704a39 159 // Copy original track and fitted track
160 // for extrapolation to ITS last layer
161 tOrigITS = *tr;
162 trackITS = track;
163
600eaa0d 164 // propagate seed to 0
165 const Double_t kMaxSnp = 0.85;
166 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
05da1b4e 167 AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
d5704a39 168
169 // propagate original track to ITS last layer
170 Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
171 AliTrackerBase::PropagateTrackTo(&tOrigITS,lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
172
173 // rotate fitted track to the frame of the original track and propagate to same reference
949d759d 174 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
d5704a39 175 trackITS.Rotate(tOrigITS.GetAlpha());
949d759d 176 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
177
d5704a39 178
4a777885 179 }
180 }
181
182 Int_t ctype(fCorrectionType);
183
184 if (fStreamer) {
185 (*fStreamer) << "Tracks" <<
186 "iev=" << iev <<
187 "z0=" << z0 <<
188 "t0=" << t0 <<
600eaa0d 189 "fTime0=" << fTime0 <<
4a777885 190 "itr=" << itr <<
191 "clsType=" << fClusterType <<
192 "corrType=" << ctype <<
193 "seedRow=" << fSeedingRow <<
194 "seedDist=" << fSeedingDist <<
38d9d609 195 "vDrift=" << vDrift <<
4a777885 196 "zLength=" << zLength <<
197 "t0seed.=" << &t0seed <<
198 "seed.=" << &seed <<
199 "track.=" << &track <<
200 "tOrig.=" << &tOrig <<
d5704a39 201 "trackITS.=" << &trackITS <<
202 "tOrigITS.=" << &tOrigITS <<
4a777885 203 "\n";
204 }
205
206
207 }
208 }
209
38d9d609 210 Cleanup();
d1cf83f5 211}
212
510cfcff 213
214//____________________________________________________________________________________
215void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
216{
217 //
5993ed4f 218 // Reconstruction for seed from associated clusters, but array of clusters:
219 // Step 1) Filling of cluster arrays
220 // Step 2) Seeding from clusters associated to tracks
221 // Step 3) Free track reconstruction using all clusters
510cfcff 222 //
223
224 TFile f(file);
225 if (!f.IsOpen() || f.IsZombie()) {
226 printf("ERROR: couldn't open the file '%s'\n", file);
227 return;
228 }
229
230 fTree=(TTree*)f.Get("toyMCtree");
231 if (!fTree) {
232 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
233 return;
234 }
235
236 fEvent=0x0;
237 fTree->SetBranchAddress("event",&fEvent);
238
239 // read spacecharge from the Userinfo ot the tree
240 InitSpaceCharge();
241
242 TString debugName=file;
243 debugName.ReplaceAll(".root","");
244 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
245 fUseMaterial,fIdealTracking,fClusterType,
246 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
247 debugName.Append(".allClusters.debug.root");
248
249 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
250 if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
251
252 gROOT->cd();
253
254 static AliExternalTrackParam dummySeedT0;
255 static AliExternalTrackParam dummySeed;
256 static AliExternalTrackParam dummyTrack;
257
258 AliExternalTrackParam t0seed;
259 AliExternalTrackParam seed;
5993ed4f 260 AliExternalTrackParam track;
510cfcff 261 AliExternalTrackParam tOrig;
262
510cfcff 263 AliExternalTrackParam *dummy;
264
265 Int_t maxev=fTree->GetEntries();
266 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
267
5993ed4f 268 // ===========================================================================================
269 // Loop 1: Fill AliTPCtrackerSector structure
270 // ===========================================================================================
c0ebc0e0 271 FillSectorStructure(maxev);
5993ed4f 272
c0ebc0e0 273 // settings (TODO: find the correct settings)
274 AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
275 tpcRecoParam->SetDoKinks(kFALSE);
276 AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
277 //tpcRecoParam->Print();
5993ed4f 278
c0ebc0e0 279 // need AliTPCReconstructor for parameter settings in AliTPCtracker
280 AliTPCReconstructor *tpcRec = new AliTPCReconstructor();
281 tpcRec->SetRecoParam(tpcRecoParam);
5993ed4f 282
283
284 // ===========================================================================================
285 // Loop 2: Seeding from clusters associated to tracks
286 // TODO: Implement tracking from given seed!
287 // ===========================================================================================
288 for (Int_t iev=0; iev<maxev; ++iev){
289 printf("============== Processing Event %6d =================\n",iev);
290 fTree->GetEvent(iev);
291 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
292 printf(" > ====== Processing Track %6d ======== \n",itr);
293 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
294 tOrig = *tr;
295
296
297 // set dummy
298 t0seed = dummySeedT0;
299 seed = dummySeed;
300 track = dummyTrack;
301
302 Float_t z0=fEvent->GetZ();
303 Float_t t0=fEvent->GetT0();
304
38d9d609 305 Float_t vDrift=GetVDrift();
5993ed4f 306 Float_t zLength=GetZLength(0);
307
c0ebc0e0 308 Int_t nClus = 0;
309
5993ed4f 310 // crate time0 seed, steered by fCreateT0seed
311 printf("t0 seed\n");
312 fTime0=-1.;
313 fCreateT0seed=kTRUE;
314 dummy = GetSeedFromTrack(tr);
315
316 if (dummy) {
317 t0seed = *dummy;
318 delete dummy;
319
320 // crate real seed using the time 0 from the first seed
321 // set fCreateT0seed now to false to get the seed in z coordinates
38d9d609 322 fTime0 = t0seed.GetZ()-zLength/vDrift;
5993ed4f 323 fCreateT0seed = kFALSE;
324 printf("seed (%.2g)\n",fTime0);
325 dummy = GetSeedFromTrack(tr);
326 if (dummy) {
327 seed = *dummy;
328 delete dummy;
329
330 // create fitted track
331 if (fDoTrackFit){
332 printf("track\n");
c0ebc0e0 333 dummy = GetFittedTrackFromSeedAllClusters(tr, &seed,nClus);
5993ed4f 334 track = *dummy;
335 delete dummy;
336 }
337
338 // propagate seed to 0
339 const Double_t kMaxSnp = 0.85;
340 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
341 AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
342
343 }
344 }
38d9d609 345
5993ed4f 346 Int_t ctype(fCorrectionType);
347
348 if (fStreamer) {
349 (*fStreamer) << "Tracks" <<
350 "iev=" << iev <<
351 "z0=" << z0 <<
352 "t0=" << t0 <<
353 "fTime0=" << fTime0 <<
354 "itr=" << itr <<
355 "clsType=" << fClusterType <<
356 "corrType=" << ctype <<
357 "seedRow=" << fSeedingRow <<
358 "seedDist=" << fSeedingDist <<
38d9d609 359 "vDrift=" << vDrift <<
5993ed4f 360 "zLength=" << zLength <<
c0ebc0e0 361 "nClus=" << nClus <<
5993ed4f 362 "t0seed.=" << &t0seed <<
363 "seed.=" << &seed <<
364 "track.=" << &track <<
365 "tOrig.=" << &tOrig <<
366 "\n";
367 }
368
369
370 }
371 }
372
373
374 delete fStreamer;
375 fStreamer=0x0;
376
377 delete fEvent;
378 fEvent = 0x0;
379
380 delete fTree;
381 fTree=0x0;
382 f.Close();
383}
384
385//____________________________________________________________________________________
386void AliToyMCReconstruction::RunRecoAllClustersStandardTracking(const char* file, Int_t nmaxEv)
387{
388 //
389 // Reconstruction for seed from associated clusters, but array of clusters
390 // Step 1) Filling of cluster arrays
391 // Step 2) Use the standard tracking: AliTPCtracker::Clusters2Tracks();
392 //
393
394 TFile f(file);
395 if (!f.IsOpen() || f.IsZombie()) {
396 printf("ERROR: couldn't open the file '%s'\n", file);
397 return;
398 }
399
400 fTree=(TTree*)f.Get("toyMCtree");
401 if (!fTree) {
402 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
403 return;
404 }
405
406 fEvent=0x0;
407 fTree->SetBranchAddress("event",&fEvent);
408
409 // read spacecharge from the Userinfo ot the tree
410 InitSpaceCharge();
411
412 TString debugName=file;
413 debugName.ReplaceAll(".root","");
414 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
415 fUseMaterial,fIdealTracking,fClusterType,
416 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
417 debugName.Append(".allClusters.debug.root");
418
419 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
420 if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
421
422 gROOT->cd();
79f54519 423
424 AliExternalTrackParam t0seed;
425 AliExternalTrackParam seed;
426 AliExternalTrackParam track;
64ea683c 427 AliExternalTrackParam tOrig;
502eb90b 428 AliToyMCTrack tOrigToy;
64ea683c 429
79f54519 430 AliExternalTrackParam *dummy;
431 AliTPCseed *seedBest;
432 AliTPCseed *seedTmp;
433 AliTPCclusterMI *cluster;
5993ed4f 434
435 Int_t maxev=fTree->GetEntries();
436 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
437
438
510cfcff 439 // ===========================================================================================
440 // Loop 1: Fill AliTPCtrackerSector structure
441 // ===========================================================================================
c0ebc0e0 442 FillSectorStructure(maxev);
510cfcff 443
510cfcff 444 // ===========================================================================================
502eb90b 445 // Loop 2: Use the TPC tracker for seeding (MakeSeeds3)
50f68b1a 446 // TODO: - check tracking configuration
447 // - add clusters and original tracks to output (how?)
510cfcff 448 // ===========================================================================================
510cfcff 449
50f68b1a 450 // settings (TODO: find the correct settings)
451 AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
452 tpcRecoParam->SetDoKinks(kFALSE);
453 AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
454 //tpcRecoParam->Print();
455
456 // need AliTPCReconstructor for parameter settings in AliTPCtracker
457 AliTPCReconstructor *tpcRec = new AliTPCReconstructor();
458 tpcRec->SetRecoParam(tpcRecoParam);
459
460 // AliTPCtracker
461 AliTPCtracker *tpcTracker = new AliTPCtracker(fTPCParam);
462 tpcTracker->SetDebug(10);
463
464 // set sector arrays
465 tpcTracker->SetTPCtrackerSectors(fInnerSectorArray,fOuterSectorArray);
466 tpcTracker->LoadInnerSectors();
467 tpcTracker->LoadOuterSectors();
468
7f1a1b08 469 // seeding
470 static TObjArray arrTracks;
79f54519 471 TObjArray * arr = &arrTracks;
472 TObjArray * seeds = new TObjArray;
7f1a1b08 473
64ea683c 474 // cuts for seeding
38d9d609 475// Float_t cuts[4];
476// cuts[0]=0.0070; // cuts[0] - fP4 cut
477// cuts[1] = 1.5; // cuts[1] - tan(phi) cut
478// cuts[2] = 3.; // cuts[2] - zvertex cut
479// cuts[3] = 3.; // cuts[3] - fP3 cut
64ea683c 480
481 // rows for seeding
482 Int_t lowerRow = fSeedingRow;
483 Int_t upperRow = fSeedingRow+2*fSeedingDist;
7f1a1b08 484 const AliTPCROC * roc = AliTPCROC::Instance();
64ea683c 485 const Int_t kNRowsInnerTPC = roc->GetNRows(0);
486 const Int_t kNRowsTPC = kNRowsInnerTPC + roc->GetNRows(36);
487 if(lowerRow < kNRowsInnerTPC){
488 Printf("Seeding row requested (%d) is lower than kNRowsInnerTPC --> use %d",lowerRow,kNRowsInnerTPC);
489 lowerRow = kNRowsInnerTPC;
490 upperRow = lowerRow + 20;
491 }
492 if(upperRow >= kNRowsTPC){
493 Printf("Seeding row requested (%d) is larger than kNRowsTPC --> use %d",upperRow,kNRowsTPC-1);
494 upperRow = kNRowsTPC-1;
495 lowerRow = upperRow-20;
496 }
497
498 // do the seeding
7f1a1b08 499 for (Int_t sec=0;sec<fkNSectorOuter;sec++){
500 //
502eb90b 501 //tpcTracker->MakeSeeds3(arr, sec,upperRow,lowerRow,cuts,-1,1);
d5704a39 502 MakeSeeds(arr, sec,upperRow,lowerRow); // own function (based on TLinearFitter)
502eb90b 503 //tpcTracker->SumTracks(seeds,arr);
504 //tpcTracker->SignClusters(seeds,3.0,3.0);
96495539 505
7f1a1b08 506 }
7f1a1b08 507
38d9d609 508 Printf("After seeding we have %d tracks",seeds->GetEntriesFast());
96495539 509
510 // Standard tracking
511 tpcTracker->SetSeeds(seeds);
512 tpcTracker->PropagateForward();
38d9d609 513 Printf("After trackinging we have %d tracks",seeds->GetEntriesFast());
514return;
50f68b1a 515
64ea683c 516 // Loop over all input tracks and connect to found seeds
517 for (Int_t iev=0; iev<maxev; ++iev){
518 printf("============== Fill Tracks: Processing Event %6d =================\n",iev);
519 fTree->GetEvent(iev);
520 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
521 printf(" > ====== Fill Tracks: Processing Track %6d ======== \n",itr);
522 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
523 tOrig = *tr;
502eb90b 524 tOrigToy = *tr;
64ea683c 525
526 Float_t z0=fEvent->GetZ();
527 Float_t t0=fEvent->GetT0();
38d9d609 528 Float_t vDrift=GetVDrift();
64ea683c 529 Float_t zLength=GetZLength(0);
7f1a1b08 530
79f54519 531 // find the corresponding seed (and track)
96495539 532 Int_t trackID = tr->GetUniqueID();
533 Int_t nClustersMC = tr->GetNumberOfSpacePoints(); // number of findable clusters (ideal)
534 if(fClusterType==1)
535 nClustersMC = tr->GetNumberOfDistSpacePoints(); // number of findable clusters (distorted)
38d9d609 536// Int_t idxSeed = 0; // index of best seed (best is with maximum number of clusters with correct ID)
96495539 537 Int_t nSeeds = 0; // number of seeds for MC track
538 Int_t nSeedClusters = 0; // number of clusters for best seed
539 Int_t nSeedClustersTmp = 0; // number of clusters for current seed
540 Int_t nSeedClustersID = 0; // number of clusters with correct ID for best seed
541 Int_t nSeedClustersIDTmp = 0; // number of clusters with correct ID for current seed
64ea683c 542 for(Int_t iSeeds = 0; iSeeds < seeds->GetEntriesFast(); ++iSeeds){
543
96495539 544 // set current seed and reset counters
64ea683c 545 seedTmp = (AliTPCseed*)(seeds->At(iSeeds));
96495539 546 nSeedClustersTmp = 0;
547 nSeedClustersIDTmp = 0;
548
502eb90b 549 if(!seedTmp) continue;
550
96495539 551 // loop over all rows
64ea683c 552 for(Int_t iRow = seedTmp->GetRow(); iRow < seedTmp->GetRow() + seedTmp->GetNumberOfClustersIndices(); iRow++ ){
7f1a1b08 553
96495539 554 // get cluster and increment counters
64ea683c 555 cluster = seedTmp->GetClusterFast(iRow);
556 if(cluster){
96495539 557 nSeedClustersTmp++;
64ea683c 558 if(cluster->GetLabel(0)==trackID){
96495539 559 nSeedClustersIDTmp++;
64ea683c 560 }
561 }
562 }
96495539 563
64ea683c 564 // if number of corresponding clusters > 0,
565 // increase nSeeds
566 if(nSeedClustersTmp > 0){
567 nSeeds++;
568 }
569
570 // if number of corresponding clusters bigger than current nSeedClusters,
96495539 571 // take this seed as the best one
572 if(nSeedClustersIDTmp > nSeedClustersID){
38d9d609 573// idxSeed = iSeeds;
79f54519 574 seedBest = seedTmp;
96495539 575 nSeedClusters = nSeedClustersTmp; // number of correctly assigned clusters
576 nSeedClustersID = nSeedClustersIDTmp; // number of all clusters
64ea683c 577 }
578
579 }
580
96495539 581 // cluster to track association (commented out, when used standard tracking)
79f54519 582 if (nSeeds>0&&nSeedClusters>0) {
96495539 583 t0seed = (AliExternalTrackParam)*seedBest;
38d9d609 584// fTime0 = t0seed.GetZ()-zLength/vDrift;
cedbc66e 585 // get the refitted track from the seed
586 // this will also set the fTime0 from the seed extrapolation
587 dummy=GetRefittedTrack(*seedBest);
588 track=*dummy;
589 delete dummy;
96495539 590 //printf("seed (%.2g): %d seeds with %d clusters\n",fTime0,nSeeds,nSeedClusters);
591
592 // // cluster to track association for all good seeds
593 // // set fCreateT0seed to true to get the seed in time coordinates
594 // fCreateT0seed = kTRUE;
595 // dummy = ClusterToTrackAssociation(seedBest,trackID,nClus);
596
597 //// propagate track to 0
598 //const Double_t kMaxSnp = 0.85;
599 //const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
600 //AliTrackerBase::PropagateTrackTo(&track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
79f54519 601
602 }
603
64ea683c 604 Int_t ctype(fCorrectionType);
605
606 if (fStreamer) {
607 (*fStreamer) << "Tracks" <<
96495539 608 "iev=" << iev <<
609 "z0=" << z0 <<
610 "t0=" << t0 <<
611 "fTime0=" << fTime0 <<
612 "itr=" << itr <<
613 "clsType=" << fClusterType <<
614 "corrType=" << ctype <<
615 "seedRow=" << fSeedingRow <<
616 "seedDist=" << fSeedingDist <<
38d9d609 617 "vDrift=" << vDrift <<
96495539 618 "zLength=" << zLength <<
619 "nClustersMC=" << nClustersMC <<
620 "nSeeds=" << nSeeds <<
621 "nSeedClusters=" << nSeedClusters <<
622 "nSeedClustersID=" << nSeedClustersID <<
623 "t0seed.=" << &t0seed <<
cedbc66e 624 "track.=" << &track <<
96495539 625 "tOrig.=" << &tOrig <<
502eb90b 626 "tOrigToy.=" << &tOrigToy <<
64ea683c 627 "\n";
628 }
629 }
630 }
631
50f68b1a 632
510cfcff 633 delete fStreamer;
634 fStreamer=0x0;
635
636 delete fEvent;
637 fEvent = 0x0;
638
639 delete fTree;
640 fTree=0x0;
641 f.Close();
642}
643
644
a06336b6 645//____________________________________________________________________________________
646void AliToyMCReconstruction::RunFullTracking(const char* file, Int_t nmaxEv)
647{
648 //
649 //
650 //
651
652 ConnectInputFile(file,nmaxEv);
653 if (!fTree) return;
654
655 InitStreamer(".fullTracking");
656
657 FillSectorStructureAC();
658
659 AliTPCReconstructor::SetStreamLevel(0);
660
661 TObjArray seeds;
662 seeds.SetOwner();
663 const Int_t lowerRow=fSeedingRow;
664 const Int_t upperRow=fSeedingRow+2*fSeedingDist;
665
666 // seeding, currently only for outer sectors
667 for (Int_t sec=0;sec<36;sec++){
668 printf("Seeding in sector: %d\n",sec);
669 MakeSeeds2(&seeds, sec,lowerRow,upperRow);
670 }
671
672 DumpSeedInfo(&seeds,lowerRow,upperRow);
673
674 Cleanup();
675}
676
d1cf83f5 677//____________________________________________________________________________________
678AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
679{
680 //
681 // if we don't have a valid time0 informaion (fTime0) available yet
682 // assume we create a seed for the time0 estimate
683 //
684
685 // seed point informaion
686 AliTrackPoint seedPoint[3];
687 const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
688
689 // number of clusters to loop over
690 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
691
692 UChar_t nextSeedRow=fSeedingRow;
693 Int_t nseeds=0;
694
695 //assumes sorted clusters
696 for (Int_t icl=0;icl<ncls;++icl) {
697 const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
698 if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
699 if (!cl) continue;
700 // use row in sector
701 const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
702 // skip clusters without proper pad row
703 if (row>200) continue;
704
705 //check seeding row
706 // if we are in the last row and still miss a seed we use the last row
707 // even if the row spacing will not be equal
708 if (row>=nextSeedRow || icl==ncls-1){
709 seedCluster[nseeds]=cl;
710 SetTrackPointFromCluster(cl, seedPoint[nseeds]);
711 ++nseeds;
712 nextSeedRow+=fSeedingDist;
713
714 if (nseeds==3) break;
715 }
716 }
717
718 // check we really have 3 seeds
719 if (nseeds!=3) {
720 AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
721 return 0x0;
722 }
723
724 // do cluster correction for fCorrectionType:
725 // 0 - no correction
726 // 1 - TPC center
727 // 2 - average eta
728 // 3 - ideal
729 // assign the cluster abs time as z component to all seeds
730 for (Int_t iseed=0; iseed<3; ++iseed) {
731 Float_t xyz[3]={0,0,0};
732 seedPoint[iseed].GetXYZ(xyz);
600eaa0d 733 const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
d1cf83f5 734
600eaa0d 735 const Int_t sector=seedCluster[iseed]->GetDetector();
736 const Int_t sign=1-2*((sector/18)%2);
d1cf83f5 737
738 if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
05da1b4e 739// printf("correction type: %d\n",(Int_t)fCorrectionType);
600eaa0d 740
741 // the settings below are for the T0 seed
742 // for known T0 the z position is already calculated in SetTrackPointFromCluster
743 if ( fCreateT0seed ){
744 if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
745 //!!! TODO: is this the correct association?
746 if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
747 }
4a777885 748
d1cf83f5 749 if ( fCorrectionType == kIdeal ) xyz[2] = seedCluster[iseed]->GetZ();
600eaa0d 750
d1cf83f5 751 //!!! TODO: to be replaced with the proper correction
223d9e38 752 fTPCCorrection->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
d1cf83f5 753 }
754
600eaa0d 755 // after the correction set the time bin as z-Position in case of a T0 seed
756 if ( fCreateT0seed )
757 xyz[2]=seedCluster[iseed]->GetTimeBin();
d1cf83f5 758
759 seedPoint[iseed].SetXYZ(xyz);
760 }
761
762 const Double_t kMaxSnp = 0.85;
763 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
764
765 AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
766 seed->ResetCovariance(10);
767
600eaa0d 768 if (fCreateT0seed){
d1cf83f5 769 // if fTime0 < 0 we assume that we create a seed for the T0 estimate
770 AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
771 if (TMath::Abs(seed->GetX())>3) {
05da1b4e 772// printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
d1cf83f5 773 }
774 }
775
776 return seed;
777
778}
779
79f54519 780
d1cf83f5 781//____________________________________________________________________________________
782void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
783{
784 //
785 // make AliTrackPoint out of AliTPCclusterMI
786 //
787
788 if (!cl) return;
4a777885 789 Float_t xyz[3]={0.,0.,0.};
d1cf83f5 790 // ClusterToSpacePoint(cl,xyz);
791 // cl->GetGlobalCov(cov);
792 //TODO: what to do with the covariance matrix???
793 //TODO: the problem is that it is used in GetAngle in AliTrackPoint
794 //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
795 //TODO: for the moment simply assign 1 permill squared
796 // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
797 // Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
798 // xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
799 // cl->GetGlobalXYZ(xyz);
800 // cl->GetGlobalCov(cov);
801 // voluem ID to add later ....
802 // p.SetXYZ(xyz);
803 // p.SetCov(cov);
0403120d 804// AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
805// p=*tp;
806// delete tp;
807 const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
d1cf83f5 808 // cl->Print();
809 // p.Print();
810 p.SetVolumeID(cl->GetDetector());
600eaa0d 811
812
813 if ( !fCreateT0seed && !fIdealTracking ) {
814 p.GetXYZ(xyz);
815 const Int_t sector=cl->GetDetector();
816 const Int_t sign=1-2*((sector/18)%2);
817 const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
05da1b4e 818// printf(" z: %.2f %.2f\n",xyz[2],zT0);
600eaa0d 819 xyz[2]=zT0;
820 p.SetXYZ(xyz);
4a777885 821 }
600eaa0d 822
823
d1cf83f5 824 // p.Rotate(p.GetAngle()).Print();
825}
826
827//____________________________________________________________________________________
828void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
829{
830 //
831 // convert the cluster to a space point in global coordinates
832 //
833 if (!cl) return;
834 xyz[0]=cl->GetRow();
835 xyz[1]=cl->GetPad();
836 xyz[2]=cl->GetTimeBin(); // this will not be correct at all
837 Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
838 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
839 fTPCParam->Transform8to4(xyz,i);
840 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
841 fTPCParam->Transform4to3(xyz,i);
842 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
843 fTPCParam->Transform2to1(xyz,i);
844 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
845}
846
847//____________________________________________________________________________________
848AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
849{
850 //
851 //
852 //
853
854 // create track
855 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
856
857 Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
858
859 const AliTPCROC * roc = AliTPCROC::Instance();
860
861 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
862 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
863 const Double_t kMaxSnp = 0.85;
864 const Double_t kMaxR = 500.;
865 const Double_t kMaxZ = 500.;
866
867 // const Double_t kMaxZ0=220;
868// const Double_t kZcut=3;
869
870 const Double_t refX = tr->GetX();
871
872 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
873
874 // loop over all other points and add to the track
4a777885 875 for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
d1cf83f5 876 AliTrackPoint pIn;
877 const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
878 if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
879 SetTrackPointFromCluster(cl, pIn);
880 if (fCorrectionType != kNoCorrection){
881 Float_t xyz[3]={0,0,0};
882 pIn.GetXYZ(xyz);
600eaa0d 883// if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
223d9e38 884 fTPCCorrection->CorrectPoint(xyz, cl->GetDetector());
d1cf83f5 885 pIn.SetXYZ(xyz);
886 }
887 // rotate the cluster to the local detector frame
888 track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
889 AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
890 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
891 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
892 //
05da1b4e 893 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
d1cf83f5 894
895 if (!res) break;
896
897 if (TMath::Abs(track->GetZ())>kMaxZ) break;
898 if (TMath::Abs(track->GetX())>kMaxR) break;
899// if (TMath::Abs(track->GetZ())<kZcut)continue;
900 //
901 Double_t pointPos[2]={0,0};
902 Double_t pointCov[3]={0,0,0};
903 pointPos[0]=prot.GetY();//local y
904 pointPos[1]=prot.GetZ();//local z
905 pointCov[0]=prot.GetCov()[3];//simay^2
906 pointCov[1]=prot.GetCov()[4];//sigmayz
907 pointCov[2]=prot.GetCov()[5];//sigmaz^2
908
909 if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
910 }
911
05da1b4e 912 AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
d1cf83f5 913
914 // rotate fittet track to the frame of the original track and propagate to same reference
915 track->Rotate(tr->GetAlpha());
916
05da1b4e 917 AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
d1cf83f5 918
919 return track;
920}
921
510cfcff 922
923//____________________________________________________________________________________
c0ebc0e0 924AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
510cfcff 925{
926 //
927 // Tracking for given seed on an array of clusters
928 //
929
930 // create track
931 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
5993ed4f 932
933 const AliTPCROC * roc = AliTPCROC::Instance();
934
935 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
936 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
c0ebc0e0 937 const Int_t kNRowsTPC = roc->GetNRows(0) + roc->GetNRows(36) - 1;
938 const Int_t kIRowsTPC = roc->GetNRows(0) - 1;
5993ed4f 939 const Double_t kMaxSnp = 0.85;
940 const Double_t kMaxR = 500.;
941 const Double_t kMaxZ = 500.;
c0ebc0e0 942 const Double_t roady = 100.;
943 const Double_t roadz = 100.;
5993ed4f 944
945 const Double_t refX = tr->GetX();
946
947 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
c0ebc0e0 948
949 Int_t secCur = -1;
950 UInt_t indexCur = 0;
951 Double_t xCur, yCur, zCur = 0.;
5993ed4f 952
64ff24a8 953 Float_t vDrift = GetVDrift();
954
5993ed4f 955 // first propagate seed to outermost row
956 AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
510cfcff 957
5993ed4f 958 // Loop over rows and find the cluster candidates
959 for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
c0ebc0e0 960
5993ed4f 961 // inner or outer sector
962 Bool_t bInnerSector = kTRUE;
963 if(iRow > kIRowsTPC) bInnerSector = kFALSE;
964
c0ebc0e0 965 // nearest track point/cluster (to be found)
966 AliTrackPoint nearestPoint;
967 AliTPCclusterMI *nearestCluster = NULL;
968
969 // Inner Sector
5993ed4f 970 if(bInnerSector){
c0ebc0e0 971
972 // Propagate to center of pad row and extract parameters
5993ed4f 973 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
c0ebc0e0 974 xCur = track->GetX();
975 yCur = track->GetY();
976 zCur = track->GetZ();
64ff24a8 977 if ( !fIdealTracking ) {
978 zCur = zCur/vDrift + fTime0; // Look at time, not at z!
979 }
c0ebc0e0 980 secCur = GetSector(track);
981
982 // Find the nearest cluster (TODO: correct road settings!)
7f1a1b08 983 Printf("inner tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
c0ebc0e0 984 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
985
986 // Move to next row if now cluster found
987 if(!nearestCluster) continue;
988 //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
989
5993ed4f 990 }
c0ebc0e0 991
992 // Outer sector
5993ed4f 993 else{
c0ebc0e0 994
995 // Propagate to center of pad row and extract parameters
5993ed4f 996 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
c0ebc0e0 997 xCur = track->GetX();
998 yCur = track->GetY();
999 zCur = track->GetZ();
64ff24a8 1000 if ( !fIdealTracking ) {
1001 zCur = zCur/vDrift + fTime0; // Look at time, not at z!
1002 }
c0ebc0e0 1003 secCur = GetSector(track);
1004
1005 // Find the nearest cluster (TODO: correct road settings!)
7f1a1b08 1006 Printf("outer tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
c0ebc0e0 1007 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1008
1009 // Move to next row if now cluster found
1010 if(!nearestCluster) continue;
1011 //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
1012
5993ed4f 1013 }
c0ebc0e0 1014
1015 // create track point from cluster
1016 SetTrackPointFromCluster(nearestCluster,nearestPoint);
1017
1018 //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
1019
64ff24a8 1020 // correction
1021 // TODO: also correction when looking for the next cluster?
1022 if (fCorrectionType != kNoCorrection){
1023 Float_t xyz[3]={0,0,0};
1024 nearestPoint.GetXYZ(xyz);
223d9e38 1025 fTPCCorrection->CorrectPoint(xyz, nearestCluster->GetDetector());
64ff24a8 1026 nearestPoint.SetXYZ(xyz);
1027 }
1028
c0ebc0e0 1029 // rotate the cluster to the local detector frame
1030 track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1031 AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1032 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1033 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1034
64ff24a8 1035
c0ebc0e0 1036 //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
1037
1038 // update track with the nearest track point
223d9e38 1039 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
510cfcff 1040
5993ed4f 1041 if (!res) break;
1042
1043 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1044 if (TMath::Abs(track->GetX())>kMaxR) break;
1045 //if (TMath::Abs(track->GetZ())<kZcut)continue;
1046
c0ebc0e0 1047 Double_t pointPos[2]={0,0};
1048 Double_t pointCov[3]={0,0,0};
1049 pointPos[0]=prot.GetY();//local y
1050 pointPos[1]=prot.GetZ();//local z
1051 pointCov[0]=prot.GetCov()[3];//simay^2
1052 pointCov[1]=prot.GetCov()[4];//sigmayz
1053 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1054
1055 if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
5993ed4f 1056
c0ebc0e0 1057 ++nClus;
5993ed4f 1058 }
c0ebc0e0 1059
5993ed4f 1060
c0ebc0e0 1061 // propagation to refX
223d9e38 1062 AliTrackerBase::PropagateTrackTo(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
5993ed4f 1063
1064 // rotate fittet track to the frame of the original track and propagate to same reference
1065 track->Rotate(tr->GetAlpha());
1066
223d9e38 1067 AliTrackerBase::PropagateTrackTo(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
c0ebc0e0 1068
1069 Printf("We have %d clusters in this track!",nClus);
5993ed4f 1070
510cfcff 1071 return track;
1072}
1073
79f54519 1074//____________________________________________________________________________________
1075AliExternalTrackParam* AliToyMCReconstruction::ClusterToTrackAssociation(const AliTPCseed *seed, Int_t trackID, Int_t &nClus)
1076{
1077 //
1078 // Cluster to track association for given seed on an array of clusters
1079 //
1080
1081 // create track
1082 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1083
1084 const AliTPCROC * roc = AliTPCROC::Instance();
1085
1086 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
1087 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1088 const Int_t kNRowsTPC = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1089 const Int_t kIRowsTPC = roc->GetNRows(0) - 1;
1090 const Double_t kMaxSnp = 0.85;
1091 const Double_t kMaxR = 500.;
1092 const Double_t kMaxZ = 500.;
1093 const Double_t roady = 0.1;
1094 const Double_t roadz = 0.01;
1095
1096 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1097
6f2a2a54 1098 Int_t secCur, secOld = -1;
79f54519 1099 UInt_t indexCur = 0;
1100 Double_t xCur, yCur, zCur = 0.;
1101
cedbc66e 1102// Float_t vDrift = GetVDrift();
79f54519 1103
1104 // first propagate seed to outermost row
96495539 1105 Bool_t res0=AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
79f54519 1106
1107 // Loop over rows and find the cluster candidates
1108 for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
1109
1110 // inner or outer sector
1111 Bool_t bInnerSector = kTRUE;
1112 if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1113
1114 // nearest track point/cluster (to be found)
1115 AliTrackPoint nearestPoint;
1116 AliTPCclusterMI *nearestCluster = NULL;
1117
1118 // Inner Sector
1119 if(bInnerSector){
1120
1121 // Propagate to center of pad row and extract parameters
1122 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
1123 xCur = track->GetX();
1124 yCur = track->GetY();
1125 zCur = track->GetZ();
1126 secCur = GetSector(track);
1127
1128 // Find the nearest cluster (TODO: correct road settings!)
6f2a2a54 1129 //Printf("inner tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
79f54519 1130 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
6f2a2a54 1131
1132 // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1133 // Increase also the road in this case
96495539 1134 if(!nearestCluster && secCur != secOld && secOld > -1){
6f2a2a54 1135 //Printf("inner tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1136 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1137 }
79f54519 1138
1139 // Move to next row if now cluster found
1140 if(!nearestCluster) continue;
1141 //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1142
1143 }
1144
1145 // Outer sector
1146 else{
1147
1148 // Propagate to center of pad row and extract parameters
1149 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
1150 xCur = track->GetX();
1151 yCur = track->GetY();
1152 zCur = track->GetZ();
1153 secCur = GetSector(track);
1154
1155 // Find the nearest cluster (TODO: correct road settings!)
96495539 1156 Printf("res0 = %d, outer tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",res0,xCur,yCur,zCur,iRow,secCur);
79f54519 1157 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1158
6f2a2a54 1159 // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1160 // Increase also the road in this case
96495539 1161 if(!nearestCluster && secCur != secOld && secOld > -1){
1162 Printf("outer tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
6f2a2a54 1163 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1164 }
1165
1166
79f54519 1167 // Move to next row if now cluster found
1168 if(!nearestCluster) continue;
96495539 1169 Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
79f54519 1170
1171 }
1172
1173 // create track point from cluster
1174 SetTrackPointFromCluster(nearestCluster,nearestPoint);
1175
6f2a2a54 1176 //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
79f54519 1177
1178 // rotate the cluster to the local detector frame
1179 track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1180 AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1181 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1182 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1183
1184
6f2a2a54 1185 //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
79f54519 1186
1187 // update track with the nearest track point
1188 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1189
1190 if (!res) break;
1191
1192 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1193 if (TMath::Abs(track->GetX())>kMaxR) break;
1194 //if (TMath::Abs(track->GetZ())<kZcut)continue;
1195
1196 Double_t pointPos[2]={0,0};
1197 Double_t pointCov[3]={0,0,0};
1198 pointPos[0]=prot.GetY();//local y
1199 pointPos[1]=prot.GetZ();//local z
1200 pointCov[0]=prot.GetCov()[3];//simay^2
1201 pointCov[1]=prot.GetCov()[4];//sigmayz
1202 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1203
1204 if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
6f2a2a54 1205 secOld = secCur;
79f54519 1206
6f2a2a54 1207 //Printf("Cluster belongs to track = %d",nearestCluster->GetLabel(0));
79f54519 1208
1209 // only count as associate cluster if it belongs to correct track!
1210 if(nearestCluster->GetLabel(0) == trackID)
1211 ++nClus;
1212 }
1213
1214 Printf("We have %d clusters in this track!",nClus);
1215
1216 return track;
1217}
1218
510cfcff 1219
d1cf83f5 1220//____________________________________________________________________________________
1221void AliToyMCReconstruction::InitSpaceCharge()
1222{
1223 //
1224 // Init the space charge map
1225 //
1226
1227 TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
1228 if (fTree) {
1229 TList *l=fTree->GetUserInfo();
1230 for (Int_t i=0; i<l->GetEntries(); ++i) {
1231 TString s(l->At(i)->GetName());
1232 if (s.Contains("SC_")) {
1233 filename=s;
1234 break;
1235 }
1236 }
1237 }
1238
1239 printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
1240 TFile f(filename.Data());
223d9e38 1241 fTPCCorrection=(AliTPCCorrection*)f.Get("map");
d1cf83f5 1242
223d9e38 1243 // fTPCCorrection = new AliTPCSpaceCharge3D();
1244 // fTPCCorrection->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
1245 // fTPCCorrection->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1246 // // fTPCCorrection->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
1247 // fTPCCorrection->InitSpaceCharge3DDistortion();
d1cf83f5 1248
1249}
1250
1251//____________________________________________________________________________________
1252Double_t AliToyMCReconstruction::GetVDrift() const
1253{
1254 //
1255 //
1256 //
1257 return fTPCParam->GetDriftV();
1258}
1259
1260//____________________________________________________________________________________
1261Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
1262{
1263 //
1264 //
1265 //
1266 if (roc<0 || roc>71) return -1;
1267 return fTPCParam->GetZLength(roc);
1268}
1269
9e98dea8 1270//____________________________________________________________________________________
1271TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
1272 TString s=gSystem->GetFromPipe(Form("ls %s",files));
1273
1274 TTree *tFirst=0x0;
1275 TObjArray *arrFiles=s.Tokenize("\n");
1276
1277 for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
1278 TString name(arrFiles->At(ifile)->GetName());
1279
1280 TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
1281 TObjArray *arrMatch=0x0;
1282 arrMatch=reg.MatchS(name);
05da1b4e 1283 TString matchName;
1284 if (arrMatch && arrMatch->At(1)) matchName=arrMatch->At(1)->GetName();
1285 else matchName=Form("%02d",ifile);
1286 delete arrMatch;
9e98dea8 1287
1288 if (!tFirst) {
1289 TFile *f=TFile::Open(name.Data());
1290 if (!f) continue;
1291 TTree *t=(TTree*)f->Get("Tracks");
1292 if (!t) {
1293 delete f;
1294 continue;
1295 }
1296
05da1b4e 1297 t->SetName(matchName.Data());
9e98dea8 1298 tFirst=t;
1299 } else {
05da1b4e 1300 tFirst->AddFriend(Form("t%s=Tracks",matchName.Data()), name.Data());
9e98dea8 1301// tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
1302 }
1303 }
1304
376089b6 1305 tFirst->GetListOfFriends()->Print();
9e98dea8 1306 return tFirst;
1307}
d1cf83f5 1308
0403120d 1309//____________________________________________________________________________________
5993ed4f 1310Int_t AliToyMCReconstruction::LoadOuterSectors() {
1311 //-----------------------------------------------------------------
1312 // This function fills outer TPC sectors with clusters.
1313 // Copy and paste from AliTPCtracker
1314 //-----------------------------------------------------------------
1315 Int_t nrows = fOuterSectorArray->GetNRows();
1316 UInt_t index=0;
1317 for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
1318 for (Int_t row = 0;row<nrows;row++){
1319 AliTPCtrackerRow* tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
1320 Int_t sec2 = sec+2*fkNSectorInner;
1321 //left
1322 Int_t ncl = tpcrow->GetN1();
1323 while (ncl--) {
1324 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1325 index=(((sec2<<8)+row)<<16)+ncl;
1326 tpcrow->InsertCluster(c,index);
1327 }
1328 //right
1329 ncl = tpcrow->GetN2();
1330 while (ncl--) {
1331 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1332 index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
1333 tpcrow->InsertCluster(c,index);
1334 }
1335 //
1336 // write indexes for fast acces
1337 //
1338 for (Int_t i=0;i<510;i++)
1339 tpcrow->SetFastCluster(i,-1);
1340 for (Int_t i=0;i<tpcrow->GetN();i++){
1341 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1342 tpcrow->SetFastCluster(zi,i); // write index
1343 }
1344 Int_t last = 0;
1345 for (Int_t i=0;i<510;i++){
1346 if (tpcrow->GetFastCluster(i)<0)
1347 tpcrow->SetFastCluster(i,last);
1348 else
1349 last = tpcrow->GetFastCluster(i);
1350 }
1351 }
1352 return 0;
1353}
1354
1355
0403120d 1356//____________________________________________________________________________________
5993ed4f 1357Int_t AliToyMCReconstruction::LoadInnerSectors() {
1358 //-----------------------------------------------------------------
1359 // This function fills inner TPC sectors with clusters.
1360 // Copy and paste from AliTPCtracker
1361 //-----------------------------------------------------------------
1362 Int_t nrows = fInnerSectorArray->GetNRows();
1363 UInt_t index=0;
1364 for (Int_t sec = 0;sec<fkNSectorInner;sec++)
1365 for (Int_t row = 0;row<nrows;row++){
1366 AliTPCtrackerRow* tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
1367 //
1368 //left
1369 Int_t ncl = tpcrow->GetN1();
1370 while (ncl--) {
1371 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1372 index=(((sec<<8)+row)<<16)+ncl;
1373 tpcrow->InsertCluster(c,index);
1374 }
1375 //right
1376 ncl = tpcrow->GetN2();
1377 while (ncl--) {
1378 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1379 index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
1380 tpcrow->InsertCluster(c,index);
1381 }
1382 //
1383 // write indexes for fast acces
1384 //
1385 for (Int_t i=0;i<510;i++)
1386 tpcrow->SetFastCluster(i,-1);
1387 for (Int_t i=0;i<tpcrow->GetN();i++){
1388 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1389 tpcrow->SetFastCluster(zi,i); // write index
1390 }
1391 Int_t last = 0;
1392 for (Int_t i=0;i<510;i++){
1393 if (tpcrow->GetFastCluster(i)<0)
1394 tpcrow->SetFastCluster(i,last);
1395 else
1396 last = tpcrow->GetFastCluster(i);
1397 }
1398
1399 }
1400 return 0;
1401}
c0ebc0e0 1402
0403120d 1403//____________________________________________________________________________________
c0ebc0e0 1404Int_t AliToyMCReconstruction::GetSector(AliExternalTrackParam *track) {
1405 //-----------------------------------------------------------------
1406 // This function returns the sector number for a given track
1407 //-----------------------------------------------------------------
1408
1409 Int_t sector = -1;
1410
1411 // get the sector number
1412 // rotate point to global system (track is already global!)
1413 Double_t xd[3];
1414 track->GetXYZ(xd);
1415 //track->Local2GlobalPosition(xd,track->GetAlpha());
1416
1417 // use TPCParams to get the sector number
1418 Float_t xyz[3] = {xd[0],xd[1],xd[2]};
1419 Int_t i[3] = {0,0,0};
1420 if(fTPCParam){
1421 sector = fTPCParam->Transform0to1(xyz,i);
1422 }
1423
1424 return sector;
1425}
1426
0403120d 1427//____________________________________________________________________________________
c0ebc0e0 1428void AliToyMCReconstruction::FillSectorStructure(Int_t maxev) {
1429 //-----------------------------------------------------------------
1430 // This function fills the sector structure of AliToyMCReconstruction
1431 //-----------------------------------------------------------------
1432
1433 // cluster array of all sectors
1434 fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
1435 fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
1436
1437 for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
1438 for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
1439
1440 Int_t count[72][96] = { {0} , {0} };
1441
1442 for (Int_t iev=0; iev<maxev; ++iev){
1443 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
1444 fTree->GetEvent(iev);
1445 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1446 printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
1447 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1448
1449 // number of clusters to loop over
1450 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1451
1452 for(Int_t icl=0; icl<ncls; ++icl){
1453
1454 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1455 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1456 if (!cl) continue;
1457
1458 Int_t sec = cl->GetDetector();
1459 Int_t row = cl->GetRow();
1460
0403120d 1461 // set Q of the cluster to 1, Q==0 does not work for the seeding
1462 cl->SetQ(1);
1463
64ff24a8 1464 // set cluster time to cluster Z (if not ideal tracking)
1465 if ( !fIdealTracking ) {
0403120d 1466 // a 'valid' position in z is needed for the seeding procedure
1467// cl->SetZ(cl->GetTimeBin()*GetVDrift());
1468 cl->SetZ(cl->GetTimeBin());
1469 }
c0ebc0e0 1470 //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1471
1472 // fill arrays for inner and outer sectors (A/C side handled internally)
1473 if (sec<fkNSectorInner*2){
7f1a1b08 1474 fInnerSectorArray[sec%fkNSectorInner].InsertCluster(cl, count[sec][row], fTPCParam);
c0ebc0e0 1475 }
1476 else{
7f1a1b08 1477 fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(cl, count[sec][row], fTPCParam);
c0ebc0e0 1478 }
1479
1480 ++count[sec][row];
1481 }
1482 }
1483 }
1484
1485 // fill the arrays completely
64ea683c 1486 // LoadOuterSectors();
1487 // LoadInnerSectors();
1488
1489 // // check the arrays
1490 // for (Int_t i=0; i<fkNSectorInner; ++i){
1491 // for (Int_t j=0; j<fInnerSectorArray[i].GetNRows(); ++j){
1492 // if(fInnerSectorArray[i][j].GetN()>0){
1493 // Printf("Inner: Sector %d Row %d : %d",i,j,fInnerSectorArray[i][j].GetN());
1494 // }
1495 // }
1496 // }
1497 // for (Int_t i=0; i<fkNSectorInner; ++i){
1498 // for (Int_t j=0; j<fOuterSectorArray[i].GetNRows(); ++j){
1499 // if(fOuterSectorArray[i][j].GetN()>0){
1500 // Printf("Outer: Sector %d Row %d : %d",i,j,fOuterSectorArray[i][j].GetN());
1501 // }
1502 // }
1503 // }
c0ebc0e0 1504}
0403120d 1505
38d9d609 1506//____________________________________________________________________________________
a06336b6 1507void AliToyMCReconstruction::FillSectorStructureAC() {
38d9d609 1508 //-----------------------------------------------------------------
1509 // This function fills the sector structure of AliToyMCReconstruction
1510 //-----------------------------------------------------------------
1511
1512 /*
1513 my god is the AliTPCtrackerSector stuff complicated!!!
1514 Ok, so here we will not fill the fClusters1 and fClusters2 of AliTPCtrackerRow,
1515 using InsertCluster of AliTPCtrackerSector, but only the fClusters via InsertCluster
1516 of AliTPCtrackerRow itself which then will not be owner, but we create an array in
1517 here (fAllClusters) which owns all clusters ...
1518 */
1519
1520 fIsAC=kTRUE;
1521 // cluster array of all sectors
1522 fInnerSectorArray = new AliTPCtrackerSector[2*fkNSectorInner];
1523 fOuterSectorArray = new AliTPCtrackerSector[2*fkNSectorOuter];
1524
1525 for (Int_t i=0; i<2*fkNSectorInner; ++i) {
1526 fInnerSectorArray[i].Setup(fTPCParam,0);
1527 }
1528
1529 for (Int_t i=0; i<2*fkNSectorOuter; ++i) {
1530 fOuterSectorArray[i].Setup(fTPCParam,1);
1531 }
1532
1533 Int_t count[72][96] = { {0} , {0} };
1534
a06336b6 1535 for (Int_t iev=0; iev<fNmaxEvents; ++iev){
38d9d609 1536 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
1537 fTree->GetEvent(iev);
1538 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
a06336b6 1539// printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
38d9d609 1540 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1541
1542 // number of clusters to loop over
1543 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1544
1545 // check if expansion of the cluster arrays is needed.
1546 if (fAllClusters.GetEntriesFast()+ncls>=fAllClusters.Capacity()) fAllClusters.Expand(2*fAllClusters.Capacity());
1547 for(Int_t icl=0; icl<ncls; ++icl){
1548
1549 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1550 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1551 if (!cl) continue;
1552
1553 // register copy to the cluster array
1554 cl = new(fAllClusters[fAllClusters.GetEntriesFast()]) AliTPCclusterMI(*cl);
1555
1556 Int_t sec = cl->GetDetector();
1557 Int_t row = cl->GetRow();
1558
1559 // set Q of the cluster to 1, Q==0 does not work for the seeding
1560 cl->SetQ(1);
1561
1562 // set cluster time to cluster Z (if not ideal tracking)
1563 if ( !fIdealTracking ) {
1564 // a 'valid' position in z is needed for the seeding procedure
1565 cl->SetZ(cl->GetTimeBin()*GetVDrift());
1566// cl->SetZ(cl->GetTimeBin());
1567 }
1568 //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1569
1570 // fill arrays for inner and outer sectors (A/C side handled internally)
1571 if (sec<fkNSectorInner*2){
1572 fInnerSectorArray[sec][row].InsertCluster(cl, count[sec][row]);
1573 }
1574 else{
1575 fOuterSectorArray[(sec-fkNSectorInner*2)][row].InsertCluster(cl, count[sec][row]);
1576 }
1577
1578 ++count[sec][row];
1579 }
1580 }
1581 }
1582
1583}
1584
0403120d 1585//____________________________________________________________________________________
1586AliToyMCTrack *AliToyMCReconstruction::ConvertTPCSeedToToyMCTrack(const AliTPCseed &seed)
1587{
1588 //
1589 //
1590 //
1591
1592
1593 AliToyMCTrack *tToy = new AliToyMCTrack(seed);
1594
1595 for (Int_t icl=0; icl<159; ++icl){
cedbc66e 1596 const AliTPCclusterMI *cl=seed.GetClusterFast(icl);
1597 if (!cl) continue;
1598 if (fClusterType==0){
1599 tToy->AddSpacePoint(*cl);
1600 } else {
1601 tToy->AddDistortedSpacePoint(*cl);
1602 }
0403120d 1603 }
1604
1605 return tToy;
1606}
1607
1608//____________________________________________________________________________________
1609AliExternalTrackParam* AliToyMCReconstruction::GetRefittedTrack(const AliTPCseed &seed)
1610{
1611 //
1612 //
1613 //
1614
cedbc66e 1615 AliExternalTrackParam *track=0x0;
0403120d 1616
38d9d609 1617 const Float_t vDrift=GetVDrift();
cedbc66e 1618 const Float_t zLength=GetZLength(0);
1619 const Double_t kMaxSnp = 0.85;
1620 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1621
1622 // for propagation use a copy
1623 AliExternalTrackParam t0seed(seed);
1624 AliTrackerBase::PropagateTrackTo(&t0seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
38d9d609 1625 fTime0 = t0seed.GetZ()-zLength/vDrift;
cedbc66e 1626 fCreateT0seed = kFALSE;
1627
1628 AliToyMCTrack *toyTrack = ConvertTPCSeedToToyMCTrack(seed);
1629
38d9d609 1630 fTime0 = t0seed.GetZ()-zLength/vDrift;
cedbc66e 1631 fCreateT0seed = kFALSE;
1632 // printf("seed (%.2g)\n",fTime0);
1633 AliExternalTrackParam *dummy = GetSeedFromTrack(toyTrack);
1634 if (dummy) {
1635 track = GetFittedTrackFromSeed(toyTrack, dummy);
1636 delete dummy;
1637 // propagate seed to 0
1638 AliTrackerBase::PropagateTrackTo(track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1639
1640 }
0403120d 1641
1642 return track;
1643}
502eb90b 1644
1645//____________________________________________________________________________________
38d9d609 1646AliTPCclusterMI* AliToyMCReconstruction::FindMiddleCluster(const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
1647 Double_t roady, Double_t roadz,
1648 AliRieman &seedFit)
1649{
1650 //
1651 //
1652 //
1653
1654 const Int_t rocInner = clInner->GetDetector();
1655 const Int_t rocOuter = clOuter->GetDetector();
1656
1657 if ( (rocInner%18) != (rocOuter%18) ){
1658 AliError("Currently only same Sector implemented");
1659 return 0x0;
1660 }
1661
1662 const Int_t innerDet=clInner->GetDetector();
1663 const Int_t outerDet=clOuter->GetDetector();
1664 const Int_t globalRowInner = clInner->GetRow() +(innerDet >35)*63;
1665 const Int_t globalRowOuter = clOuter->GetRow() +(outerDet >35)*63;
1666
1667 Int_t iMiddle = (globalRowInner+globalRowOuter)/2;
1668 Int_t roc = innerDet;
1669 if (iMiddle>62){
1670 iMiddle-=63;
1671 roc=outerDet;
1672 }
1673
1674 const AliTPCtrackerRow& krMiddle = fOuterSectorArray[roc%36][iMiddle]; // middle
1675 // initial guess use simple linear interpolation
1676 Double_t y=(clInner->GetY()+clOuter->GetY())/2;
1677 Double_t z=(clInner->GetZ()+clOuter->GetZ())/2;
1678 if (seedFit.IsValid()) {
1679 // update values once the fit was performed
1680 y=seedFit.GetYat(krMiddle.GetX());
1681 z=seedFit.GetZat(krMiddle.GetX());
1682 }
1683
1684 AliTPCclusterMI *n=krMiddle.FindNearest(y,z,roady,roadz);
1685// if (n)
1686// 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,
1687// n->GetX(), n->GetY(),n->GetZ(),n->GetLabel(0),
1688// clInner->GetX(), clInner->GetY(),clInner->GetZ(),clInner->GetLabel(0),
1689// clOuter->GetX(), clOuter->GetY(),clOuter->GetZ(),clOuter->GetLabel(0)
1690// );
1691 return n;
1692}
1693
1694//____________________________________________________________________________________
1695void AliToyMCReconstruction::AddMiddleClusters(AliTPCseed *seed,
1696 const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
1697 Double_t roady, Double_t roadz,
1698 Int_t &nTotalClusters, AliRieman &seedFit)
1699{
1700 //
1701 // Iteratively add the middle clusters
1702 //
1703
1704 // update rieman fit with every second point added
1705 AliTPCclusterMI *clMiddle=FindMiddleCluster(clInner,clOuter,roady,roadz,seedFit);
1706
1707 // break iterative process
1708 if (!clMiddle) return;
1709
1710 const Int_t globalRowInner = clInner->GetRow() +(clInner->GetDetector() >35)*63;
1711 const Int_t globalRowMiddle = clMiddle->GetRow()+(clMiddle->GetDetector()>35)*63;
1712 const Int_t globalRowOuter = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
1713
1714 seed->SetClusterPointer(globalRowMiddle,clMiddle);
1715 ++nTotalClusters;
1716// printf(" > Total clusters: %d\n",nTotalClusters);
1717 seedFit.AddPoint(clMiddle->GetX(), clMiddle->GetY(), clMiddle->GetZ(),
1718 TMath::Sqrt(clMiddle->GetSigmaY2()), TMath::Sqrt(clMiddle->GetSigmaZ2()));
1719
1720 if (!(seedFit.GetN()%3)) {
1721// printf(" call update: %d (%d)\n",seedFit.GetN(),nTotalClusters);
1722// printf(" Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f\n",
1723// seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z());
1724 seedFit.Update();
1725 }
1726 if ( seedFit.IsValid() && seedFit.GetChi2()>1000 ) return;
1727
1728 // Add middle clusters towards outer radius
1729 if (globalRowMiddle+1<globalRowOuter) AddMiddleClusters(seed,clMiddle,clOuter,roady,roadz,nTotalClusters,seedFit);
1730
1731 // Add middle clusters towards innter radius
1732 if (globalRowInner+1<globalRowMiddle) AddMiddleClusters(seed,clInner,clMiddle,roady,roadz,nTotalClusters,seedFit);
1733
1734 return;
1735}
1736
1737//____________________________________________________________________________________
1738void AliToyMCReconstruction::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t iRowInner, Int_t iRowOuter)
1739{
1740 //
1741 // find seeds in a sector, requires iRowInner < iRowOuter
1742 // iRowXXX runs from 0-159, so over IROC and OROC
1743 //
1744
1745 if (!fIsAC) {
1746 AliError("This function requires the sector arrays filled for AC tracking");
1747 return;
1748 }
1749
1750 // swap rows in case they are in the wrong order
1751 if (iRowInner>iRowOuter) {
1752 Int_t tmp=iRowInner;
1753 iRowInner=iRowOuter;
1754 iRowOuter=tmp;
1755 }
1756
a06336b6 1757 if (iRowOuter>158) iRowOuter=158;
1758 if (iRowInner<0) iRowInner=0;
1759
38d9d609 1760 // only for CookLabel
1761 AliTPCtracker tpcTracker(fTPCParam);
1762
1763 // Define the search roads:
1764 // timeRoadCombinatorics - the maximum time difference used for the
1765 // combinatorics. Since we don't have a 'z-Vertex' estimate this will
1766 // reduce the combinatorics significantly when iterating on one TF
1767 // use a little more than one full drift length of the TPC to allow for
1768 // distortions
1769 //
1770 // timeRoad - the time difference allowed when associating the cluster
1771 // in the middle of the other two used for the initial search
1772 //
1773 // padRoad - the local y difference allowed when associating the middle cluster
38d9d609 1774
1775// Double_t timeRoadCombinatorics = 270./vDrift;
1776// Double_t timeRoad = 20./vDrift;
1777 Double_t timeRoadCombinatorics = 270.;
1778 Double_t timeRoad = 20.;
1779 Double_t padRoad = 10.;
1780
1781
1782 // fOuterSectorArray runs from 0-95, so from iRowXXX which runs from 0-159
1783 // the number of rows in the IROC has to be subtracted
1784 const Int_t innerRows=fInnerSectorArray->GetNRows();
1785 const AliTPCtrackerRow& krOuter = fOuterSectorArray[sec][iRowOuter - innerRows]; // down
1786 const AliTPCtrackerRow& krInner = fOuterSectorArray[sec][iRowInner - innerRows]; // up
1787
1788 AliTPCseed *seed = new AliTPCseed;
1789
1790 const Int_t nMaxClusters=iRowOuter-iRowInner+1;
1791// Int_t nScannedClusters = 0;
1792
1793 // loop over all points in the firstand last search row
1794 for (Int_t iOuter=0; iOuter < krOuter; iOuter++) {
1795 const AliTPCclusterMI *clOuter = krOuter[iOuter];
a06336b6 1796 if (clOuter->IsUsed()) continue;
38d9d609 1797 for (Int_t iInner=0; iInner < krInner; iInner++) {
1798 const AliTPCclusterMI *clInner = krInner[iInner];
a06336b6 1799 if (clInner->IsUsed()) continue;
38d9d609 1800// printf("\n\n Check combination %d (%d), %d (%d)\n",iOuter, iInner, clOuter->GetLabel(0), clInner->GetLabel(0));
1801 // check maximum distance for combinatorics
1802 if (TMath::Abs(clOuter->GetZ()-clInner->GetZ())>timeRoadCombinatorics) continue;
1803// printf(" Is inside one drift\n");
1804
1805 // use rieman fit for seed description
1806 AliRieman seedFit(159);
1807 // add first two points
1808 seedFit.AddPoint(clInner->GetX(), clInner->GetY(), clInner->GetZ(),
1809 TMath::Sqrt(clInner->GetSigmaY2()), TMath::Sqrt(clInner->GetSigmaZ2()));
1810 seedFit.AddPoint(clOuter->GetX(), clOuter->GetY(), clOuter->GetZ(),
1811 TMath::Sqrt(clOuter->GetSigmaY2()), TMath::Sqrt(clOuter->GetSigmaZ2()));
1812
1813 // Iteratively add all clusters in the respective middle
1814 Int_t nFoundClusters=2;
1815 AddMiddleClusters(seed,clInner,clOuter,padRoad,timeRoad,nFoundClusters,seedFit);
1816// printf(" Clusters attached: %d\n",nFoundClusters);
a06336b6 1817 if (nFoundClusters>2) seedFit.Update();
38d9d609 1818// printf(" Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f\n",
1819// seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z());
1820
1821 // check for minimum number of assigned clusters and a decent chi2
1822 if ( nFoundClusters<0.5*nMaxClusters || seedFit.GetChi2()>1000 ){
1823 seed->Reset();
1824 continue;
1825 }
1826// printf(" >>> Seed accepted\n");
1827 // add original outer clusters
1828 const Int_t globalRowInner = clInner->GetRow() +(clInner->GetDetector() >35)*63;
1829 const Int_t globalRowOuter = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
1830
1831 seed->SetClusterPointer(globalRowInner,const_cast<AliTPCclusterMI*>(clInner));
1832 seed->SetClusterPointer(globalRowOuter,const_cast<AliTPCclusterMI*>(clOuter));
1833
1834 // set parameters retrieved from AliRieman
1835 Double_t params[5]={0};
1836 Double_t covar[15]={0};
1837 const Double_t alpha=TMath::DegToRad()*(clInner->GetDetector()%18*20.+10.);
1838 const Double_t x=clInner->GetX();
1839 seedFit.GetExternalParameters(x,params,covar);
1840
1841 seed->Set(x,alpha,params,covar);
1842
1843 // set label of the seed. At least 60% of the clusters need the correct label
a06336b6 1844 CookLabel(seed,.6);
1845// printf(" - Label: %d\n",seed->GetLabel());
1846 // mark clusters as being used
1847 MarkClustersUsed(seed);
38d9d609 1848
1849 arr->Add(seed);
1850 seed=new AliTPCseed;
1851 }
1852 }
1853 //delete surplus seed
1854 delete seed;
1855 seed=0x0;
1856
38d9d609 1857}
1858//____________________________________________________________________________________
1859void AliToyMCReconstruction::MakeSeeds(TObjArray * /*arr*/, Int_t sec, Int_t iRow1, Int_t iRow2)
502eb90b 1860{
1861 //
1862 // Create seeds between i1 and i2 (stored in arr) with TLinearFitter
1863 //
1864 // sec: sector number
1865 // iRow1: upper row
1866 // iRow2: lower row
1867 //
1868
1869 // Create Fitter
1870 static TLinearFitter fitter(3,"pol2");
1871
1872 // Get 3 padrows (iRow1,iMiddle=(iRow1+iRow2)/2,iRow2)
1873 Int_t iMiddle = (iRow1+iRow2)/2;
1874 const AliTPCtrackerRow& krd = fOuterSectorArray[sec][iRow2-fInnerSectorArray->GetNRows()]; // down
1875 const AliTPCtrackerRow& krm = fOuterSectorArray[sec][iMiddle-fInnerSectorArray->GetNRows()]; // middle
1876 const AliTPCtrackerRow& kru = fOuterSectorArray[sec][iRow1-fInnerSectorArray->GetNRows()]; // up
1877
1878 // Loop over 3 cluster possibilities
1879 for (Int_t iu=0; iu < kru; iu++) {
1880 for (Int_t im=0; im < krm; im++) {
1881 for (Int_t id=0; id < krd; id++) {
1882
1883 // clear all points
1884 fitter.ClearPoints();
1885
1886 // add all three points to fitter
1887 Double_t xy[2] = {kru[iu]->GetX(),kru[iu]->GetY()};
1888 Double_t z = kru[iu]->GetZ();
1889 fitter.AddPoint(xy,z);
1890
1891 xy[0] = krm[im]->GetX();
1892 xy[1] = krm[im]->GetY();
1893 z = krm[im]->GetZ();
1894 fitter.AddPoint(xy,z);
1895
1896 xy[0] = krd[id]->GetX();
1897 xy[1] = krd[id]->GetY();
1898 z = krd[id]->GetZ();
1899 fitter.AddPoint(xy,z);
1900
1901 // Evaluate and get parameters
1902 fitter.Eval();
1903
1904 // how to get the other clusters now?
1905 // ...
1906
1907 }
1908 }
1909 }
1910}
38d9d609 1911
1912//____________________________________________________________________________________
a06336b6 1913void AliToyMCReconstruction::InitStreamer(TString addName, Int_t level)
38d9d609 1914{
1915 //
1916 // initilise the debug streamer and set the logging level
1917 // use a default naming convention
1918 //
a06336b6 1919
1920 delete fStreamer;
1921 fStreamer=0x0;
1922
1923 if (addName.IsNull()) addName=".dummy";
38d9d609 1924
a06336b6 1925 if (!fTree) return;
1926
1927 TString debugName=fInputFile->GetName();
1928 debugName.ReplaceAll(".root","");
1929 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
1930 fUseMaterial,fIdealTracking,fClusterType,
1931 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
1932 debugName.Append(addName);
1933 debugName.Append(".root");
1934
1935 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
1936 fStreamer=new TTreeSRedirector(debugName.Data());
1937 fStreamer->SetUniqueID(level);
1938
1939 gROOT->cd();
38d9d609 1940}
1941
1942//____________________________________________________________________________________
a06336b6 1943void AliToyMCReconstruction::ConnectInputFile(const char* file, Int_t nmaxEv)
38d9d609 1944{
1945 //
1946 // connect the tree and event pointer from the input file
1947 //
1948
1949 delete fInputFile;
1950 fInputFile=0x0;
1951 fEvent=0x0;
1952 fTree=0;
1953
1954 fInputFile=TFile::Open(file);
1955 if (!fInputFile || !fInputFile->IsOpen() || fInputFile->IsZombie()) {
1956 delete fInputFile;
1957 fInputFile=0x0;
1958 printf("ERROR: couldn't open the file '%s'\n", file);
1959 return;
1960 }
1961
1962 fTree=(TTree*)fInputFile->Get("toyMCtree");
1963 if (!fTree) {
1964 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
1965 return;
1966 }
1967
1968 fEvent=0x0;
1969 fTree->SetBranchAddress("event",&fEvent);
1970
1971 gROOT->cd();
a06336b6 1972
1973 fNmaxEvents=fTree->GetEntries();
1974 if (nmaxEv>-1) fNmaxEvents=TMath::Min(nmaxEv,fNmaxEvents);
38d9d609 1975
a06336b6 1976 // setup space charge map from Userinfo of the tree
1977 InitSpaceCharge();
1978
1979 // setup the track maps
38d9d609 1980 SetupTrackMaps();
1981}
1982
1983//____________________________________________________________________________________
1984void AliToyMCReconstruction::Cleanup()
1985{
1986 //
1987 // Cleanup input data
1988 //
1989
a06336b6 1990 if (fStreamer) delete fStreamer;
38d9d609 1991 fStreamer=0x0;
1992
1993 delete fEvent;
1994 fEvent = 0x0;
1995
a06336b6 1996// delete fTree;
38d9d609 1997 fTree=0x0;
1998
1999 delete fInputFile;
2000 fInputFile=0x0;
2001
2002}
2003
2004//____________________________________________________________________________________
2005void AliToyMCReconstruction::SetupTrackMaps()
2006{
2007 //
2008 //
2009 //
2010
2011 fMapTrackEvent.Delete();
2012 fMapTrackTrackInEvent.Delete();
2013
2014 if (!fTree) {
2015 AliError("Tree not connected");
2016 return;
2017 }
2018
a06336b6 2019 Int_t nmaxEv=fTree->GetEntries();
2020 if (fNmaxEvents>-1) nmaxEv=fNmaxEvents;
2021
2022 for (Int_t iev=0; iev<nmaxEv; ++iev) {
38d9d609 2023 fTree->GetEvent(iev);
2024 if (!fEvent) continue;
2025
2026 const Int_t ntracks=fEvent->GetNumberOfTracks();
2027 if (fMapTrackEvent.GetSize()+ntracks>=fMapTrackEvent.Capacity()) fMapTrackEvent.Expand(2*fMapTrackEvent.Capacity());
2028 if (fMapTrackTrackInEvent.GetSize()+ntracks>=fMapTrackTrackInEvent.Capacity()) fMapTrackTrackInEvent.Expand(2*fMapTrackTrackInEvent.Capacity());
2029
2030 for (Int_t itrack=0; itrack<ntracks; ++itrack){
2031 Int_t label=fEvent->GetTrack(itrack)->GetUniqueID();
2032
2033 fMapTrackEvent.Add(label,iev);
2034 fMapTrackTrackInEvent.Add(label,itrack);
2035 }
2036 }
2037
2038}
a06336b6 2039
2040//____________________________________________________________________________________
2041void AliToyMCReconstruction::CookLabel(AliTPCseed *seed, Double_t fraction, Int_t info[5])
2042{
2043 //
2044 //
2045 //
2046
2047 Int_t labels[159]={0};
2048// Long64_t posMaxLabel=-1;
2049 Int_t nMaxLabel=0; // clusters from maximum label
2050 Int_t nMaxLabel2=0; // clusters from second maximum
2051 Int_t nlabels=0;
2052 Int_t maxLabel=-1; // label with most clusters
2053 Int_t maxLabel2=-1; // label with second most clusters
2054 Int_t nclusters=0;
2055 TExMap labelMap(159);
2056
2057 for (Int_t icl=0; icl<159; ++icl) {
2058 const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2059 if (!cl) continue;
2060 ++nclusters;
2061
2062 const Int_t clLabel=cl->GetLabel(0);
2063 // a not assinged value returns 0, so we need to add 1 and subtract it afterwards
2064 Long64_t labelPos=labelMap.GetValue(clLabel);
2065
2066 if (!labelPos) {
2067 labelPos=nlabels+1;
2068 labelMap.Add(clLabel,labelPos);
2069 ++nlabels;
2070 }
2071 --labelPos;
2072
2073 const Int_t nCurrentLabel=++labels[labelPos];
2074 if (nCurrentLabel>nMaxLabel) {
2075 nMaxLabel2=nMaxLabel;
2076 nMaxLabel=nCurrentLabel;
2077// posMaxLabel=labelPos;
2078 maxLabel2=maxLabel;
2079 maxLabel=clLabel;
2080 }
2081 }
2082
2083 if (Double_t(nMaxLabel)/nclusters<fraction) maxLabel=-maxLabel;
2084
2085 seed->SetLabel(maxLabel);
2086
2087 if (info) {
2088 info[0]=nMaxLabel;
2089 info[1]=nMaxLabel2;
2090 info[2]=maxLabel2;
2091 info[3]=nclusters;
2092 info[4]=nlabels;
2093 }
2094}
2095
2096
2097//____________________________________________________________________________________
2098void AliToyMCReconstruction::DumpSeedInfo(TObjArray *arr, Int_t iRowInner, Int_t iRowOuter)
2099{
2100
2101 // for debugging
2102 if (!fStreamer || !fTree) return;
2103 // swap rows in case they are in the wrong order
2104 if (iRowInner>iRowOuter) {
2105 Int_t tmp=iRowInner;
2106 iRowInner=iRowOuter;
2107 iRowOuter=tmp;
2108 }
2109
2110 AliInfo("");
2111
2112 const Double_t kMaxSnp = 0.85;
2113 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2114 Float_t vDrift=GetVDrift();
2115 Float_t zLength=GetZLength(0);
2116
2117 //loop over all events and tracks and try to associate the seed to the track
2118 for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed){
2119 AliTPCseed *seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
2120 Int_t seedLabel=seed->GetLabel();
2121
2122 // get original track
2123 Int_t iev=fMapTrackEvent.GetValue(TMath::Abs(seedLabel));
2124 Int_t itr=fMapTrackTrackInEvent.GetValue(TMath::Abs(seedLabel));
2125
2126 fTree->GetEvent(iev);
2127 const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
2128
2129 AliExternalTrackParam extTrack(*toyTrack);
2130
2131 //propagate to same reference frame
2132 AliExternalTrackParam extSeed(*seed);
2133 AliTrackerBase::PropagateTrackTo(&extSeed,0,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2134 extSeed.Rotate(extTrack.GetAlpha());
2135 AliTrackerBase::PropagateTrackTo(&extSeed,0,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2136 //create propagated track
2137 AliExternalTrackParam *extSeedRefit=GetRefittedTrack(*seed);
2138 AliTrackerBase::PropagateTrackTo(extSeedRefit,0,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2139 extSeedRefit->Rotate(extTrack.GetAlpha());
2140 AliTrackerBase::PropagateTrackTo(extSeedRefit,0,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2141
2142 Int_t roc=-1;
2143 //
2144 //count findable and found clusters in the seed
2145 //
2146 Int_t nClustersSeedMax=iRowOuter-iRowInner+1;
2147 Int_t nClustersFindable=0;
2148 Int_t nClustersSeed=0;
2149
2150 const Int_t ncls=(fClusterType==0)?toyTrack->GetNumberOfSpacePoints():toyTrack->GetNumberOfDistSpacePoints();
2151
2152 Int_t rowInner=iRowInner-(iRowInner>62)*63;
2153 Int_t rowOuter=iRowOuter-(iRowOuter>62)*63;
2154 //findable
2155 for (Int_t icl=0; icl<ncls; ++icl) {
2156 const AliTPCclusterMI *cl=(fClusterType==0)?toyTrack->GetSpacePoint(icl):toyTrack->GetDistortedSpacePoint(icl);
2157 roc=cl->GetDetector();
2158 Int_t row=cl->GetRow();
2159 // printf("row: %d, iRowInner: %d, iRowOuter: %d\n", row, iRowInner, iRowOuter);
2160 if ( (row<rowInner) || (row>rowOuter) ) continue;
2161 ++nClustersFindable;
2162 }
2163
2164 //found in seed
2165 for (Int_t icl=0; icl<159; ++icl) {
2166 const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2167 if (!cl) continue;
2168 const Int_t row=cl->GetRow();
2169 if ( (row<rowInner) || (row>rowOuter) ) continue;
2170 ++nClustersSeed;
2171 }
2172
2173 // convert back to time, since we made tracking in 'pseude z coordinates'
2174 fTime0=extSeed.GetZ()/vDrift;
2175
2176 Float_t z0=fEvent->GetZ();
2177 Float_t t0=fEvent->GetT0();
2178
2179 Int_t ctype(fCorrectionType);
2180
2181 Int_t info[5]={0};
2182 CookLabel(seed,.6,info);
2183
2184 (*fStreamer) << "Seeds" <<
2185 "iev=" << iev <<
2186 "iseed=" << iseed <<
2187 "z0=" << z0 <<
2188 "t0=" << t0 <<
2189 "fTime0=" << fTime0 <<
2190 "itr=" << itr <<
2191 "clsType=" << fClusterType <<
2192 "corrType=" << ctype <<
2193 "seedRowInner=" << iRowInner <<
2194 "seedRowOuter=" << iRowOuter <<
2195 "vDrift=" << vDrift <<
2196 "zLength=" << zLength <<
2197 "track.=" << &extSeed <<
2198 "track2.=" << extSeedRefit <<
2199 "tOrig.=" << &extTrack <<
2200 "seedLabel=" << seedLabel <<
2201 "nclMax=" << nClustersSeedMax <<
2202 "nclFindable=" << nClustersFindable <<
2203 "nclFound=" << nClustersSeed <<
2204 "nMaxLabel=" << info[0] <<
2205 "nMaxLabel2=" << info[1] <<
2206 "maxLabel2=" << info[2] <<
2207 "nclusters=" << info[3] <<
2208 "nlabels=" << info[4] <<
2209 "roc=" << roc <<
2210 "\n";
2211
2212 delete extSeedRefit;
2213 }
2214}
2215
2216//____________________________________________________________________________________
2217void AliToyMCReconstruction::MarkClustersUsed(AliTPCseed *seed)
2218{
2219 //
2220 //
2221 //
2222
2223 for (Int_t icl=0; icl<159; ++icl) {
2224 AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2225 if (cl) cl->Use();
2226 }
2227}
2228
2229//____________________________________________________________________________________
2230void AliToyMCReconstruction::DumpTracksToTree(const char* file)
2231{
2232 //
2233 //
2234 //
2235 ConnectInputFile(file);
2236 if (!fTree) return;
2237
2238 delete fStreamer;
2239 fStreamer=0x0;
2240
2241 TString debugName=fInputFile->GetName();
2242 debugName.ReplaceAll(".root",".AllTracks.root");
2243
2244 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
2245 fStreamer=new TTreeSRedirector(debugName.Data());
2246
2247 for (Int_t iev=0;iev<fNmaxEvents;++iev){
2248 fTree->GetEvent(iev);
2249 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks();++itr){
2250 AliToyMCTrack *toyTrack=const_cast<AliToyMCTrack*>(fEvent->GetTrack(itr));
2251 Int_t trackID=toyTrack->GetUniqueID();
2252
2253 (*fStreamer) << "Tracks" <<
2254 "iev=" << iev <<
2255 "itr=" << itr <<
2256 "trackID=" << trackID <<
2257 "track.=" << toyTrack <<
2258 "\n";
2259
2260 }
2261 }
2262
2263 Cleanup();
2264}