]> 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
macros for space charge plots (MW)
[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
7f72a764 967 // 4 - preliminary eta
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
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;
7f72a764 986 if ( fCorrectionType == kPreliminaryEta ) xyz[2] = r/TMath::Tan(prelTheta)*sign;
600eaa0d 987 }
4a777885 988
d1cf83f5 989 if ( fCorrectionType == kIdeal ) xyz[2] = seedCluster[iseed]->GetZ();
48e13991 990
d1cf83f5 991 //!!! TODO: to be replaced with the proper correction
223d9e38 992 fTPCCorrection->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
d1cf83f5 993 }
994
600eaa0d 995 // after the correction set the time bin as z-Position in case of a T0 seed
996 if ( fCreateT0seed )
c484014f 997 xyz[2]=seedCluster[iseed]->GetTimeBin() + sign*( zBeforeCorr - xyz[2] )/GetVDrift();
48e13991 998 // xyz[2]=seedCluster[iseed]->GetTimeBin()*sign;
999
d1cf83f5 1000 seedPoint[iseed].SetXYZ(xyz);
1001 }
1002
1003 const Double_t kMaxSnp = 0.85;
1004 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1005
1006 AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
d1cf83f5 1007
48e13991 1008 if (fCreateT0seed&&!fLongT0seed){
1009 // only propagate to vertex if we don't create a long seed
d1cf83f5 1010 // if fTime0 < 0 we assume that we create a seed for the T0 estimate
1011 AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
1012 if (TMath::Abs(seed->GetX())>3) {
05da1b4e 1013// printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
d1cf83f5 1014 }
1015 }
48e13991 1016
1017 seed->ResetCovariance(10);
d1cf83f5 1018 return seed;
1019
1020}
1021
e83fd282 1022//____________________________________________________________________________________
1023AliExternalTrackParam* AliToyMCReconstruction::GetTrackRefit(const AliToyMCTrack * const tr, EDet det)
1024{
1025 //
1026 // Get the ITS or TRD track refitted from the toy track
1027 // type: 0=ITS; 1=TRD
1028 //
1029
1030 AliExternalTrackParam *track=GetSeedFromTrackIdeal(tr,det);
1031 if (!track) return 0x0;
1032
1033 const Double_t kMaxSnp = 0.85;
1034 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1035
1036 Int_t npoints=0;
1037 switch (det) {
1038 case kITS:
1039 npoints=tr->GetNumberOfITSPoints();
1040 break;
1041 case kTPC:
1042 npoints=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1043 break;
1044 case kTRD:
1045 npoints=tr->GetNumberOfTRDPoints();
1046 break;
1047 }
1048
1049 const AliCluster *cl=0x0;
1050
1051 for (Int_t ipoint=0; ipoint<npoints; ++ipoint) {
1052 AliTrackPoint pIn;
1053
1054 switch (det) {
1055 case kITS:
1056 pIn=(*tr->GetITSPoint(ipoint));
1057 break;
1058 case kTPC:
1059 cl=tr->GetSpacePoint(ipoint);
1060 if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
1061 AliTPCclusterMI::SetGlobalTrackPoint(*cl,pIn);
1062 break;
1063 case kTRD:
1064 pIn=(*tr->GetTRDPoint(ipoint));
1065 break;
1066 }
1067
1068
1069 const Double_t angle=pIn.GetAngle();
1070 track->Rotate(angle);
1071 AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1072
1073 if (!AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial)) {
1074 AliInfo(Form("Could not propagate track to x=%.2f (a=%.2f) for det %d",prot.GetX(),angle,det));
1075 }
1076 //
1077
1078 Double_t pointPos[2]={0,0};
1079 Double_t pointCov[3]={0,0,0};
1080 pointPos[0]=prot.GetY();//local y
1081 pointPos[1]=prot.GetZ();//local z
1082 pointCov[0]=prot.GetCov()[3];//simay^2
1083 pointCov[1]=prot.GetCov()[4];//sigmayz
1084 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1085
1086 if (!track->Update(pointPos,pointCov)) {
1087 AliInfo(Form("no update: det: %d",det));
1088 break;
1089 }
1090
1091 }
1092
1093 return track;
1094}
1095
79f54519 1096
d1cf83f5 1097//____________________________________________________________________________________
1098void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
1099{
1100 //
1101 // make AliTrackPoint out of AliTPCclusterMI
1102 //
1103
1104 if (!cl) return;
4a777885 1105 Float_t xyz[3]={0.,0.,0.};
d1cf83f5 1106 // ClusterToSpacePoint(cl,xyz);
1107 // cl->GetGlobalCov(cov);
1108 //TODO: what to do with the covariance matrix???
1109 //TODO: the problem is that it is used in GetAngle in AliTrackPoint
1110 //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
1111 //TODO: for the moment simply assign 1 permill squared
1112 // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
1113 // Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
1114 // xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
1115 // cl->GetGlobalXYZ(xyz);
1116 // cl->GetGlobalCov(cov);
1117 // voluem ID to add later ....
1118 // p.SetXYZ(xyz);
1119 // p.SetCov(cov);
0403120d 1120// AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
1121// p=*tp;
1122// delete tp;
e83fd282 1123// const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
1124 AliTPCclusterMI::SetGlobalTrackPoint(*cl,p);
d1cf83f5 1125 // cl->Print();
1126 // p.Print();
1127 p.SetVolumeID(cl->GetDetector());
600eaa0d 1128
1129
1130 if ( !fCreateT0seed && !fIdealTracking ) {
1131 p.GetXYZ(xyz);
1132 const Int_t sector=cl->GetDetector();
1133 const Int_t sign=1-2*((sector/18)%2);
1134 const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
05da1b4e 1135// printf(" z: %.2f %.2f\n",xyz[2],zT0);
600eaa0d 1136 xyz[2]=zT0;
1137 p.SetXYZ(xyz);
4a777885 1138 }
600eaa0d 1139
1140
d1cf83f5 1141 // p.Rotate(p.GetAngle()).Print();
1142}
1143
1144//____________________________________________________________________________________
1145void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
1146{
1147 //
1148 // convert the cluster to a space point in global coordinates
1149 //
1150 if (!cl) return;
1151 xyz[0]=cl->GetRow();
1152 xyz[1]=cl->GetPad();
1153 xyz[2]=cl->GetTimeBin(); // this will not be correct at all
1154 Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
1155 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
1156 fTPCParam->Transform8to4(xyz,i);
1157 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
1158 fTPCParam->Transform4to3(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->Transform2to1(xyz,i);
1161 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
1162}
1163
1164//____________________________________________________________________________________
1165AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
1166{
1167 //
1168 //
1169 //
1170
1171 // create track
1172 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1173
1174 Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1175
1176 const AliTPCROC * roc = AliTPCROC::Instance();
1177
1178 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
1179 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1180 const Double_t kMaxSnp = 0.85;
1181 const Double_t kMaxR = 500.;
1182 const Double_t kMaxZ = 500.;
1183
1184 // const Double_t kMaxZ0=220;
1185// const Double_t kZcut=3;
1186
1187 const Double_t refX = tr->GetX();
1188
1189 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1190
1191 // loop over all other points and add to the track
4a777885 1192 for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
d1cf83f5 1193 AliTrackPoint pIn;
1194 const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
1195 if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
8a2a61d3 1196 const Int_t globalRow = cl->GetRow() +(cl->GetDetector() >35)*63;
1197 if ( fCreateT0seed ){
1198 if ( globalRow<fSeedingRow || globalRow>fSeedingRow+2*fSeedingDist ) continue;
1199 }
1200
d1cf83f5 1201 SetTrackPointFromCluster(cl, pIn);
8a2a61d3 1202
1203 Float_t xyz[3]={0,0,0};
1204 pIn.GetXYZ(xyz);
48e13991 1205 Float_t zBeforeCorr = xyz[2];
1206
1207 const Int_t sector=cl->GetDetector();
1208 const Int_t sign=1-2*((sector/18)%2);
8a2a61d3 1209
d1cf83f5 1210 if (fCorrectionType != kNoCorrection){
7f72a764 1211
8a2a61d3 1212 const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
8a2a61d3 1213
1214 if ( fCreateT0seed ){
1215 if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
1216 //!!! TODO: is this the correct association?
1217 if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
1218 if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
1219 }
1220
223d9e38 1221 fTPCCorrection->CorrectPoint(xyz, cl->GetDetector());
d1cf83f5 1222 }
8a2a61d3 1223
1224 if ( fCreateT0seed )
c484014f 1225 xyz[2]=cl->GetTimeBin() + sign*( zBeforeCorr - xyz[2] )/GetVDrift();
48e13991 1226 // xyz[2]=cl->GetTimeBin();
8a2a61d3 1227 pIn.SetXYZ(xyz);
1228
d1cf83f5 1229 // rotate the cluster to the local detector frame
1230 track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
1231 AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1232 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1233 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1234 //
05da1b4e 1235 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
d1cf83f5 1236
1237 if (!res) break;
1238
1239 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1240 if (TMath::Abs(track->GetX())>kMaxR) break;
1241// if (TMath::Abs(track->GetZ())<kZcut)continue;
1242 //
1243 Double_t pointPos[2]={0,0};
1244 Double_t pointCov[3]={0,0,0};
1245 pointPos[0]=prot.GetY();//local y
1246 pointPos[1]=prot.GetZ();//local z
1247 pointCov[0]=prot.GetCov()[3];//simay^2
1248 pointCov[1]=prot.GetCov()[4];//sigmayz
1249 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1250
1251 if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
1252 }
1253
05da1b4e 1254 AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
d1cf83f5 1255
48e13991 1256 if (!fCreateT0seed){
1257 // rotate fittet track to the frame of the original track and propagate to same reference
1258 track->Rotate(tr->GetAlpha());
d1cf83f5 1259
48e13991 1260 AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1261 }
d1cf83f5 1262
1263 return track;
1264}
1265
510cfcff 1266
1267//____________________________________________________________________________________
c0ebc0e0 1268AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
510cfcff 1269{
1270 //
1271 // Tracking for given seed on an array of clusters
1272 //
1273
1274 // create track
1275 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
5993ed4f 1276
1277 const AliTPCROC * roc = AliTPCROC::Instance();
1278
1279 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
1280 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
c0ebc0e0 1281 const Int_t kNRowsTPC = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1282 const Int_t kIRowsTPC = roc->GetNRows(0) - 1;
5993ed4f 1283 const Double_t kMaxSnp = 0.85;
1284 const Double_t kMaxR = 500.;
1285 const Double_t kMaxZ = 500.;
c0ebc0e0 1286 const Double_t roady = 100.;
1287 const Double_t roadz = 100.;
5993ed4f 1288
1289 const Double_t refX = tr->GetX();
1290
1291 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
c0ebc0e0 1292
1293 Int_t secCur = -1;
1294 UInt_t indexCur = 0;
1295 Double_t xCur, yCur, zCur = 0.;
5993ed4f 1296
64ff24a8 1297 Float_t vDrift = GetVDrift();
1298
5993ed4f 1299 // first propagate seed to outermost row
1300 AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
510cfcff 1301
5993ed4f 1302 // Loop over rows and find the cluster candidates
1303 for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
c0ebc0e0 1304
5993ed4f 1305 // inner or outer sector
1306 Bool_t bInnerSector = kTRUE;
1307 if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1308
c0ebc0e0 1309 // nearest track point/cluster (to be found)
1310 AliTrackPoint nearestPoint;
1311 AliTPCclusterMI *nearestCluster = NULL;
1312
1313 // Inner Sector
5993ed4f 1314 if(bInnerSector){
c0ebc0e0 1315
1316 // Propagate to center of pad row and extract parameters
5993ed4f 1317 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
c0ebc0e0 1318 xCur = track->GetX();
1319 yCur = track->GetY();
1320 zCur = track->GetZ();
64ff24a8 1321 if ( !fIdealTracking ) {
1322 zCur = zCur/vDrift + fTime0; // Look at time, not at z!
1323 }
c0ebc0e0 1324 secCur = GetSector(track);
1325
1326 // Find the nearest cluster (TODO: correct road settings!)
7f1a1b08 1327 Printf("inner tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
c0ebc0e0 1328 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1329
1330 // Move to next row if now cluster found
1331 if(!nearestCluster) continue;
1332 //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1333
5993ed4f 1334 }
c0ebc0e0 1335
1336 // Outer sector
5993ed4f 1337 else{
c0ebc0e0 1338
1339 // Propagate to center of pad row and extract parameters
5993ed4f 1340 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
c0ebc0e0 1341 xCur = track->GetX();
1342 yCur = track->GetY();
1343 zCur = track->GetZ();
64ff24a8 1344 if ( !fIdealTracking ) {
1345 zCur = zCur/vDrift + fTime0; // Look at time, not at z!
1346 }
c0ebc0e0 1347 secCur = GetSector(track);
1348
1349 // Find the nearest cluster (TODO: correct road settings!)
7f1a1b08 1350 Printf("outer tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
c0ebc0e0 1351 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1352
1353 // Move to next row if now cluster found
1354 if(!nearestCluster) continue;
1355 //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
1356
5993ed4f 1357 }
c0ebc0e0 1358
1359 // create track point from cluster
1360 SetTrackPointFromCluster(nearestCluster,nearestPoint);
1361
1362 //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
1363
64ff24a8 1364 // correction
1365 // TODO: also correction when looking for the next cluster?
1366 if (fCorrectionType != kNoCorrection){
1367 Float_t xyz[3]={0,0,0};
1368 nearestPoint.GetXYZ(xyz);
223d9e38 1369 fTPCCorrection->CorrectPoint(xyz, nearestCluster->GetDetector());
64ff24a8 1370 nearestPoint.SetXYZ(xyz);
1371 }
1372
c0ebc0e0 1373 // rotate the cluster to the local detector frame
1374 track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1375 AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1376 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1377 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1378
64ff24a8 1379
c0ebc0e0 1380 //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
1381
1382 // update track with the nearest track point
223d9e38 1383 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
510cfcff 1384
5993ed4f 1385 if (!res) break;
1386
1387 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1388 if (TMath::Abs(track->GetX())>kMaxR) break;
1389 //if (TMath::Abs(track->GetZ())<kZcut)continue;
1390
c0ebc0e0 1391 Double_t pointPos[2]={0,0};
1392 Double_t pointCov[3]={0,0,0};
1393 pointPos[0]=prot.GetY();//local y
1394 pointPos[1]=prot.GetZ();//local z
1395 pointCov[0]=prot.GetCov()[3];//simay^2
1396 pointCov[1]=prot.GetCov()[4];//sigmayz
1397 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1398
1399 if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
5993ed4f 1400
c0ebc0e0 1401 ++nClus;
5993ed4f 1402 }
c0ebc0e0 1403
5993ed4f 1404
c0ebc0e0 1405 // propagation to refX
223d9e38 1406 AliTrackerBase::PropagateTrackTo(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
5993ed4f 1407
1408 // rotate fittet track to the frame of the original track and propagate to same reference
1409 track->Rotate(tr->GetAlpha());
1410
223d9e38 1411 AliTrackerBase::PropagateTrackTo(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
c0ebc0e0 1412
1413 Printf("We have %d clusters in this track!",nClus);
5993ed4f 1414
510cfcff 1415 return track;
1416}
1417
79f54519 1418//____________________________________________________________________________________
1419AliExternalTrackParam* AliToyMCReconstruction::ClusterToTrackAssociation(const AliTPCseed *seed, Int_t trackID, Int_t &nClus)
1420{
1421 //
1422 // Cluster to track association for given seed on an array of clusters
1423 //
1424
1425 // create track
1426 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1427
1428 const AliTPCROC * roc = AliTPCROC::Instance();
1429
1430 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
1431 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1432 const Int_t kNRowsTPC = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1433 const Int_t kIRowsTPC = roc->GetNRows(0) - 1;
1434 const Double_t kMaxSnp = 0.85;
1435 const Double_t kMaxR = 500.;
1436 const Double_t kMaxZ = 500.;
1437 const Double_t roady = 0.1;
1438 const Double_t roadz = 0.01;
1439
1440 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1441
6f2a2a54 1442 Int_t secCur, secOld = -1;
79f54519 1443 UInt_t indexCur = 0;
1444 Double_t xCur, yCur, zCur = 0.;
1445
cedbc66e 1446// Float_t vDrift = GetVDrift();
79f54519 1447
1448 // first propagate seed to outermost row
96495539 1449 Bool_t res0=AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
79f54519 1450
1451 // Loop over rows and find the cluster candidates
1452 for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
1453
1454 // inner or outer sector
1455 Bool_t bInnerSector = kTRUE;
1456 if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1457
1458 // nearest track point/cluster (to be found)
1459 AliTrackPoint nearestPoint;
1460 AliTPCclusterMI *nearestCluster = NULL;
1461
1462 // Inner Sector
1463 if(bInnerSector){
1464
1465 // Propagate to center of pad row and extract parameters
1466 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
1467 xCur = track->GetX();
1468 yCur = track->GetY();
1469 zCur = track->GetZ();
1470 secCur = GetSector(track);
1471
1472 // Find the nearest cluster (TODO: correct road settings!)
6f2a2a54 1473 //Printf("inner tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
79f54519 1474 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
6f2a2a54 1475
1476 // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1477 // Increase also the road in this case
96495539 1478 if(!nearestCluster && secCur != secOld && secOld > -1){
6f2a2a54 1479 //Printf("inner tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1480 nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1481 }
79f54519 1482
1483 // Move to next row if now cluster found
1484 if(!nearestCluster) continue;
1485 //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1486
1487 }
1488
1489 // Outer sector
1490 else{
1491
1492 // Propagate to center of pad row and extract parameters
1493 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
1494 xCur = track->GetX();
1495 yCur = track->GetY();
1496 zCur = track->GetZ();
1497 secCur = GetSector(track);
1498
1499 // Find the nearest cluster (TODO: correct road settings!)
96495539 1500 Printf("res0 = %d, outer tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",res0,xCur,yCur,zCur,iRow,secCur);
79f54519 1501 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1502
6f2a2a54 1503 // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1504 // Increase also the road in this case
96495539 1505 if(!nearestCluster && secCur != secOld && secOld > -1){
1506 Printf("outer tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
6f2a2a54 1507 nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1508 }
1509
1510
79f54519 1511 // Move to next row if now cluster found
1512 if(!nearestCluster) continue;
96495539 1513 Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
79f54519 1514
1515 }
1516
1517 // create track point from cluster
1518 SetTrackPointFromCluster(nearestCluster,nearestPoint);
1519
6f2a2a54 1520 //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
79f54519 1521
1522 // rotate the cluster to the local detector frame
1523 track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1524 AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
1525 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1526 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1527
1528
6f2a2a54 1529 //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
79f54519 1530
1531 // update track with the nearest track point
1532 Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1533
1534 if (!res) break;
1535
1536 if (TMath::Abs(track->GetZ())>kMaxZ) break;
1537 if (TMath::Abs(track->GetX())>kMaxR) break;
1538 //if (TMath::Abs(track->GetZ())<kZcut)continue;
1539
1540 Double_t pointPos[2]={0,0};
1541 Double_t pointCov[3]={0,0,0};
1542 pointPos[0]=prot.GetY();//local y
1543 pointPos[1]=prot.GetZ();//local z
1544 pointCov[0]=prot.GetCov()[3];//simay^2
1545 pointCov[1]=prot.GetCov()[4];//sigmayz
1546 pointCov[2]=prot.GetCov()[5];//sigmaz^2
1547
1548 if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
6f2a2a54 1549 secOld = secCur;
79f54519 1550
6f2a2a54 1551 //Printf("Cluster belongs to track = %d",nearestCluster->GetLabel(0));
79f54519 1552
1553 // only count as associate cluster if it belongs to correct track!
1554 if(nearestCluster->GetLabel(0) == trackID)
1555 ++nClus;
1556 }
1557
1558 Printf("We have %d clusters in this track!",nClus);
1559
1560 return track;
1561}
1562
9d548f3c 1563//____________________________________________________________________________________
1564void AliToyMCReconstruction::AssociateClusters(AliTPCseed &seed, Int_t firstRow, Int_t lastRow, Bool_t direction)
1565{
1566 //
1567 // do cluster to track association from first to last row
1568 // direction 0: outwards; 1: inwards
1569 //
1570
1571 Double_t roady=10.;
1572 Double_t roadz=10.;
1573
1574 AliRieman rieman1(160);
1575 AliRieman rieman2(160);
1576 SetRieman(seed,rieman1);
1577 CopyRieman(rieman1,rieman2);
1578
1579 Int_t sec=seed.GetSector();
1580 Int_t noLastPoint=0;
1581 //TODO: change to inward and outwar search?
1582 // -> better handling of non consecutive points
1583 if (direction){
1584 firstRow*=-1;
1585 lastRow*=-1;
1586 }
1587
1588 //always from inside out
1589 if (firstRow>lastRow){
1590 Int_t tmp=firstRow;
1591 firstRow=lastRow;
1592 lastRow=tmp;
1593 }
1594
1595 for (Int_t row=firstRow; row<=lastRow && noLastPoint<3;++row) {
1596 Int_t iRow=TMath::Abs(row);
1597 const AliTPCclusterMI *cl=seed.GetClusterPointer(iRow);
1598 if (cl) continue;
1599
1600 const Int_t secrow = iRow<63?iRow:iRow-63;
1601
1602 AliTPCtrackerSector *arrSec=(iRow<63)?fInnerSectorArray:fOuterSectorArray;
1603 const AliTPCtrackerRow& kr = arrSec[sec%36][secrow];
1604 const Double_t maxy=arrSec[sec%36].GetMaxY(secrow);
1605
1606 Double_t y=rieman1.GetYat(kr.GetX());
1607 Double_t z=rieman1.GetZat(kr.GetX());
1608
1609 if (TMath::Abs(y)>maxy) {
1610 AliError("Tracking over sector boundaries not implemented, yet");
1611 continue;
1612 }
1613
1614 AliTPCclusterMI *n=kr.FindNearest(y,z,roady,roadz);
1615 if (!n || n->IsUsed()) {
1616 ++noLastPoint;
1617 continue;
1618 }
1619 // check for quality of the cluster
1620 // TODO: better?
1621 rieman2.AddPoint(n->GetX(), n->GetY(), n->GetZ(),
1622 TMath::Sqrt(n->GetSigmaY2()), TMath::Sqrt(n->GetSigmaZ2()));
1623 rieman2.Update();
afc28cd1 1624// printf(" Riemann results: row=%d valid=%d, Chi2=%.2f (%.2f) %d (%d)",
1625// iRow, rieman2.IsValid(), rieman2.GetChi2(), rieman1.GetChi2(), n->GetLabel(0),seed.GetLabel());
9d548f3c 1626 Double_t limit=2*rieman1.GetChi2();
1627 if (fClusterType==0) limit=1000;
1628 if (rieman2.GetChi2()>limit) {
1629 CopyRieman(rieman1,rieman2);
1630 ++noLastPoint;
afc28cd1 1631// printf("\n");
9d548f3c 1632 continue;
1633 }
afc28cd1 1634// printf(" +++ \n");
9d548f3c 1635
1636 noLastPoint=0;
1637 //use point
1638 rieman1.AddPoint(n->GetX(), n->GetY(), n->GetZ(),
1639 TMath::Sqrt(n->GetSigmaY2()), TMath::Sqrt(n->GetSigmaZ2()));
1640 rieman1.Update();
1641
1642 seed.SetClusterPointer(iRow,n);
1643 // if (iRow<seed.GetSeed1()) seed.SetSeed1(iRow);
1644 // if (iRow>seed.GetSeed2()) seed.SetSeed2(iRow);
1645 n->Use();
1646
1647 }
1648}
1649
1650//____________________________________________________________________________________
1651void AliToyMCReconstruction::ClusterToTrackAssociation(AliTPCseed &seed)
1652{
1653 //
1654 //
1655 //
1656
afc28cd1 1657// printf("\n ============ \nnext Seed: %d\n",seed.GetLabel());
9d548f3c 1658 //assume seed is within one sector
1659 Int_t iMiddle=(seed.GetSeed1()+seed.GetSeed2())/2;
1660 //outward
1661 AssociateClusters(seed,iMiddle+1,158,kFALSE);
1662 //inward
1663 AssociateClusters(seed,0,iMiddle,kTRUE);
1664 seed.SetIsSeeding(kFALSE);
1665
1666 CookLabel(&seed,.6);
1667}
1668
510cfcff 1669
d1cf83f5 1670//____________________________________________________________________________________
1671void AliToyMCReconstruction::InitSpaceCharge()
1672{
1673 //
1674 // Init the space charge map
1675 //
1676
1677 TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
1678 if (fTree) {
1679 TList *l=fTree->GetUserInfo();
1680 for (Int_t i=0; i<l->GetEntries(); ++i) {
1681 TString s(l->At(i)->GetName());
1682 if (s.Contains("SC_")) {
1683 filename=s;
1684 break;
1685 }
1686 }
1687 }
1688
1689 printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
1690 TFile f(filename.Data());
223d9e38 1691 fTPCCorrection=(AliTPCCorrection*)f.Get("map");
d1cf83f5 1692
223d9e38 1693 // fTPCCorrection = new AliTPCSpaceCharge3D();
1694 // fTPCCorrection->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
1695 // fTPCCorrection->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1696 // // fTPCCorrection->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
1697 // fTPCCorrection->InitSpaceCharge3DDistortion();
d1cf83f5 1698
1699}
1700
1701//____________________________________________________________________________________
1702Double_t AliToyMCReconstruction::GetVDrift() const
1703{
1704 //
1705 //
1706 //
1707 return fTPCParam->GetDriftV();
1708}
1709
1710//____________________________________________________________________________________
1711Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
1712{
1713 //
1714 //
1715 //
1716 if (roc<0 || roc>71) return -1;
1717 return fTPCParam->GetZLength(roc);
1718}
1719
9e98dea8 1720//____________________________________________________________________________________
1721TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
1722 TString s=gSystem->GetFromPipe(Form("ls %s",files));
1723
1724 TTree *tFirst=0x0;
1725 TObjArray *arrFiles=s.Tokenize("\n");
1726
1727 for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
1728 TString name(arrFiles->At(ifile)->GetName());
1729
1730 TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
1731 TObjArray *arrMatch=0x0;
1732 arrMatch=reg.MatchS(name);
05da1b4e 1733 TString matchName;
1734 if (arrMatch && arrMatch->At(1)) matchName=arrMatch->At(1)->GetName();
1735 else matchName=Form("%02d",ifile);
1736 delete arrMatch;
9e98dea8 1737
1738 if (!tFirst) {
1739 TFile *f=TFile::Open(name.Data());
1740 if (!f) continue;
1741 TTree *t=(TTree*)f->Get("Tracks");
1742 if (!t) {
1743 delete f;
1744 continue;
1745 }
1746
05da1b4e 1747 t->SetName(matchName.Data());
9e98dea8 1748 tFirst=t;
1749 } else {
05da1b4e 1750 tFirst->AddFriend(Form("t%s=Tracks",matchName.Data()), name.Data());
9e98dea8 1751// tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
1752 }
1753 }
1754
c484014f 1755 if (tFirst->GetListOfFriends()) tFirst->GetListOfFriends()->Print();
9e98dea8 1756 return tFirst;
1757}
d1cf83f5 1758
0403120d 1759//____________________________________________________________________________________
5993ed4f 1760Int_t AliToyMCReconstruction::LoadOuterSectors() {
1761 //-----------------------------------------------------------------
1762 // This function fills outer TPC sectors with clusters.
1763 // Copy and paste from AliTPCtracker
1764 //-----------------------------------------------------------------
1765 Int_t nrows = fOuterSectorArray->GetNRows();
1766 UInt_t index=0;
1767 for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
1768 for (Int_t row = 0;row<nrows;row++){
1769 AliTPCtrackerRow* tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
1770 Int_t sec2 = sec+2*fkNSectorInner;
1771 //left
1772 Int_t ncl = tpcrow->GetN1();
1773 while (ncl--) {
1774 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1775 index=(((sec2<<8)+row)<<16)+ncl;
1776 tpcrow->InsertCluster(c,index);
1777 }
1778 //right
1779 ncl = tpcrow->GetN2();
1780 while (ncl--) {
1781 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1782 index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
1783 tpcrow->InsertCluster(c,index);
1784 }
1785 //
1786 // write indexes for fast acces
1787 //
1788 for (Int_t i=0;i<510;i++)
1789 tpcrow->SetFastCluster(i,-1);
1790 for (Int_t i=0;i<tpcrow->GetN();i++){
1791 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1792 tpcrow->SetFastCluster(zi,i); // write index
1793 }
1794 Int_t last = 0;
1795 for (Int_t i=0;i<510;i++){
1796 if (tpcrow->GetFastCluster(i)<0)
1797 tpcrow->SetFastCluster(i,last);
1798 else
1799 last = tpcrow->GetFastCluster(i);
1800 }
1801 }
1802 return 0;
1803}
1804
1805
0403120d 1806//____________________________________________________________________________________
5993ed4f 1807Int_t AliToyMCReconstruction::LoadInnerSectors() {
1808 //-----------------------------------------------------------------
1809 // This function fills inner TPC sectors with clusters.
1810 // Copy and paste from AliTPCtracker
1811 //-----------------------------------------------------------------
1812 Int_t nrows = fInnerSectorArray->GetNRows();
1813 UInt_t index=0;
1814 for (Int_t sec = 0;sec<fkNSectorInner;sec++)
1815 for (Int_t row = 0;row<nrows;row++){
1816 AliTPCtrackerRow* tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
1817 //
1818 //left
1819 Int_t ncl = tpcrow->GetN1();
1820 while (ncl--) {
1821 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1822 index=(((sec<<8)+row)<<16)+ncl;
1823 tpcrow->InsertCluster(c,index);
1824 }
1825 //right
1826 ncl = tpcrow->GetN2();
1827 while (ncl--) {
1828 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1829 index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
1830 tpcrow->InsertCluster(c,index);
1831 }
1832 //
1833 // write indexes for fast acces
1834 //
1835 for (Int_t i=0;i<510;i++)
1836 tpcrow->SetFastCluster(i,-1);
1837 for (Int_t i=0;i<tpcrow->GetN();i++){
1838 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1839 tpcrow->SetFastCluster(zi,i); // write index
1840 }
1841 Int_t last = 0;
1842 for (Int_t i=0;i<510;i++){
1843 if (tpcrow->GetFastCluster(i)<0)
1844 tpcrow->SetFastCluster(i,last);
1845 else
1846 last = tpcrow->GetFastCluster(i);
1847 }
1848
1849 }
1850 return 0;
1851}
c0ebc0e0 1852
0403120d 1853//____________________________________________________________________________________
c0ebc0e0 1854Int_t AliToyMCReconstruction::GetSector(AliExternalTrackParam *track) {
1855 //-----------------------------------------------------------------
1856 // This function returns the sector number for a given track
1857 //-----------------------------------------------------------------
1858
1859 Int_t sector = -1;
1860
1861 // get the sector number
1862 // rotate point to global system (track is already global!)
1863 Double_t xd[3];
1864 track->GetXYZ(xd);
1865 //track->Local2GlobalPosition(xd,track->GetAlpha());
1866
1867 // use TPCParams to get the sector number
1868 Float_t xyz[3] = {xd[0],xd[1],xd[2]};
1869 Int_t i[3] = {0,0,0};
1870 if(fTPCParam){
1871 sector = fTPCParam->Transform0to1(xyz,i);
1872 }
1873
1874 return sector;
1875}
1876
0403120d 1877//____________________________________________________________________________________
c0ebc0e0 1878void AliToyMCReconstruction::FillSectorStructure(Int_t maxev) {
1879 //-----------------------------------------------------------------
1880 // This function fills the sector structure of AliToyMCReconstruction
1881 //-----------------------------------------------------------------
1882
1883 // cluster array of all sectors
1884 fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
1885 fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
1886
1887 for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
1888 for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
1889
1890 Int_t count[72][96] = { {0} , {0} };
1891
1892 for (Int_t iev=0; iev<maxev; ++iev){
1893 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
1894 fTree->GetEvent(iev);
1895 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1896 printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
1897 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1898
1899 // number of clusters to loop over
1900 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1901
1902 for(Int_t icl=0; icl<ncls; ++icl){
1903
1904 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1905 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1906 if (!cl) continue;
1907
1908 Int_t sec = cl->GetDetector();
1909 Int_t row = cl->GetRow();
1910
0403120d 1911 // set Q of the cluster to 1, Q==0 does not work for the seeding
1912 cl->SetQ(1);
1913
64ff24a8 1914 // set cluster time to cluster Z (if not ideal tracking)
1915 if ( !fIdealTracking ) {
0403120d 1916 // a 'valid' position in z is needed for the seeding procedure
1917// cl->SetZ(cl->GetTimeBin()*GetVDrift());
1918 cl->SetZ(cl->GetTimeBin());
1919 }
c0ebc0e0 1920 //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1921
1922 // fill arrays for inner and outer sectors (A/C side handled internally)
1923 if (sec<fkNSectorInner*2){
7f1a1b08 1924 fInnerSectorArray[sec%fkNSectorInner].InsertCluster(cl, count[sec][row], fTPCParam);
c0ebc0e0 1925 }
1926 else{
7f1a1b08 1927 fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(cl, count[sec][row], fTPCParam);
c0ebc0e0 1928 }
1929
1930 ++count[sec][row];
1931 }
1932 }
1933 }
1934
1935 // fill the arrays completely
64ea683c 1936 // LoadOuterSectors();
1937 // LoadInnerSectors();
1938
1939 // // check the arrays
1940 // for (Int_t i=0; i<fkNSectorInner; ++i){
1941 // for (Int_t j=0; j<fInnerSectorArray[i].GetNRows(); ++j){
1942 // if(fInnerSectorArray[i][j].GetN()>0){
1943 // Printf("Inner: Sector %d Row %d : %d",i,j,fInnerSectorArray[i][j].GetN());
1944 // }
1945 // }
1946 // }
1947 // for (Int_t i=0; i<fkNSectorInner; ++i){
1948 // for (Int_t j=0; j<fOuterSectorArray[i].GetNRows(); ++j){
1949 // if(fOuterSectorArray[i][j].GetN()>0){
1950 // Printf("Outer: Sector %d Row %d : %d",i,j,fOuterSectorArray[i][j].GetN());
1951 // }
1952 // }
1953 // }
c0ebc0e0 1954}
0403120d 1955
38d9d609 1956//____________________________________________________________________________________
a06336b6 1957void AliToyMCReconstruction::FillSectorStructureAC() {
38d9d609 1958 //-----------------------------------------------------------------
1959 // This function fills the sector structure of AliToyMCReconstruction
1960 //-----------------------------------------------------------------
1961
1962 /*
1963 my god is the AliTPCtrackerSector stuff complicated!!!
1964 Ok, so here we will not fill the fClusters1 and fClusters2 of AliTPCtrackerRow,
1965 using InsertCluster of AliTPCtrackerSector, but only the fClusters via InsertCluster
1966 of AliTPCtrackerRow itself which then will not be owner, but we create an array in
1967 here (fAllClusters) which owns all clusters ...
1968 */
1969
1970 fIsAC=kTRUE;
1971 // cluster array of all sectors
1972 fInnerSectorArray = new AliTPCtrackerSector[2*fkNSectorInner];
1973 fOuterSectorArray = new AliTPCtrackerSector[2*fkNSectorOuter];
1974
1975 for (Int_t i=0; i<2*fkNSectorInner; ++i) {
1976 fInnerSectorArray[i].Setup(fTPCParam,0);
1977 }
1978
1979 for (Int_t i=0; i<2*fkNSectorOuter; ++i) {
1980 fOuterSectorArray[i].Setup(fTPCParam,1);
1981 }
1982
1983 Int_t count[72][96] = { {0} , {0} };
1984
a06336b6 1985 for (Int_t iev=0; iev<fNmaxEvents; ++iev){
38d9d609 1986 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
1987 fTree->GetEvent(iev);
1988 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
a06336b6 1989// printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
38d9d609 1990 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1991
1992 // number of clusters to loop over
1993 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1994
1995 // check if expansion of the cluster arrays is needed.
1996 if (fAllClusters.GetEntriesFast()+ncls>=fAllClusters.Capacity()) fAllClusters.Expand(2*fAllClusters.Capacity());
1997 for(Int_t icl=0; icl<ncls; ++icl){
1998
1999 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
2000 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
2001 if (!cl) continue;
2002
2003 // register copy to the cluster array
2004 cl = new(fAllClusters[fAllClusters.GetEntriesFast()]) AliTPCclusterMI(*cl);
2005
2006 Int_t sec = cl->GetDetector();
2007 Int_t row = cl->GetRow();
2008
2009 // set Q of the cluster to 1, Q==0 does not work for the seeding
2010 cl->SetQ(1);
2011
2012 // set cluster time to cluster Z (if not ideal tracking)
2013 if ( !fIdealTracking ) {
2014 // a 'valid' position in z is needed for the seeding procedure
9d548f3c 2015 Double_t sign=1;
2016 if (((sec/18)%2)==1) sign=-1;
48e13991 2017 cl->SetZ(cl->GetTimeBin()*GetVDrift());
9d548f3c 2018 //mark cluster to be time*vDrift by setting the type to 1
2019 cl->SetType(1);
38d9d609 2020// cl->SetZ(cl->GetTimeBin());
2021 }
2022 //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
2023
2024 // fill arrays for inner and outer sectors (A/C side handled internally)
2025 if (sec<fkNSectorInner*2){
2026 fInnerSectorArray[sec][row].InsertCluster(cl, count[sec][row]);
2027 }
2028 else{
2029 fOuterSectorArray[(sec-fkNSectorInner*2)][row].InsertCluster(cl, count[sec][row]);
2030 }
2031
2032 ++count[sec][row];
2033 }
2034 }
2035 }
2036
2037}
2038
0403120d 2039//____________________________________________________________________________________
2040AliToyMCTrack *AliToyMCReconstruction::ConvertTPCSeedToToyMCTrack(const AliTPCseed &seed)
2041{
2042 //
2043 //
2044 //
2045
2046
2047 AliToyMCTrack *tToy = new AliToyMCTrack(seed);
2048
2049 for (Int_t icl=0; icl<159; ++icl){
cedbc66e 2050 const AliTPCclusterMI *cl=seed.GetClusterFast(icl);
2051 if (!cl) continue;
2052 if (fClusterType==0){
2053 tToy->AddSpacePoint(*cl);
2054 } else {
2055 tToy->AddDistortedSpacePoint(*cl);
2056 }
0403120d 2057 }
2058
2059 return tToy;
2060}
2061
2062//____________________________________________________________________________________
2063AliExternalTrackParam* AliToyMCReconstruction::GetRefittedTrack(const AliTPCseed &seed)
2064{
2065 //
2066 //
2067 //
2068
cedbc66e 2069 AliExternalTrackParam *track=0x0;
0403120d 2070
38d9d609 2071 const Float_t vDrift=GetVDrift();
cedbc66e 2072 const Float_t zLength=GetZLength(0);
2073 const Double_t kMaxSnp = 0.85;
2074 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2075
cedbc66e 2076 AliToyMCTrack *toyTrack = ConvertTPCSeedToToyMCTrack(seed);
9d548f3c 2077
2078 fTime0 = 0;
2079
2080 //get t0 estimate
2081 fCreateT0seed = kTRUE;
2082 AliExternalTrackParam *t0seed = GetSeedFromTrack(toyTrack,kTRUE);
2083 if (!t0seed) return 0x0;
cedbc66e 2084
9d548f3c 2085 fTime0 = t0seed->GetZ()-zLength/vDrift;
2086 delete t0seed;
2087 t0seed=0x0;
2088
cedbc66e 2089 fCreateT0seed = kFALSE;
9d548f3c 2090 AliExternalTrackParam *dummy = GetSeedFromTrack(toyTrack,kTRUE);
2091 track = GetFittedTrackFromSeed(toyTrack, dummy);
2092 delete dummy;
2093 // propagate seed to 0
2094 AliTrackerBase::PropagateTrackTo(track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
0403120d 2095
2096 return track;
2097}
502eb90b 2098
2099//____________________________________________________________________________________
38d9d609 2100AliTPCclusterMI* AliToyMCReconstruction::FindMiddleCluster(const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
2101 Double_t roady, Double_t roadz,
2102 AliRieman &seedFit)
2103{
2104 //
2105 //
2106 //
2107
2108 const Int_t rocInner = clInner->GetDetector();
2109 const Int_t rocOuter = clOuter->GetDetector();
2110
2111 if ( (rocInner%18) != (rocOuter%18) ){
2112 AliError("Currently only same Sector implemented");
2113 return 0x0;
2114 }
2115
2116 const Int_t innerDet=clInner->GetDetector();
2117 const Int_t outerDet=clOuter->GetDetector();
2118 const Int_t globalRowInner = clInner->GetRow() +(innerDet >35)*63;
2119 const Int_t globalRowOuter = clOuter->GetRow() +(outerDet >35)*63;
2120
9d548f3c 2121 AliTPCclusterMI *n=0x0;
2122
2123 // allow flexibility of +- nRowsGrace Rows to find a middle cluster
2124 const Int_t nRowsGrace = 0;
2125 for (Int_t iter=0; iter<2*nRowsGrace+1; ++iter){
2126 Int_t iMiddle = (globalRowInner+globalRowOuter)/2;
2127 iMiddle+=((iter+1)/2)*(1-2*((iter+1)%2));
38d9d609 2128
9d548f3c 2129 Int_t localRow=iMiddle;
2130 Int_t roc = innerDet;
2131 if (iMiddle>62){
2132 localRow-=63;
2133 roc=outerDet;
2134 }
2135
2136 AliTPCtrackerSector *arrRow=(iMiddle<63)?fInnerSectorArray:fOuterSectorArray;
2137 const AliTPCtrackerRow& krMiddle = arrRow[roc%36][localRow]; // middle
2138 // initial guess use simple linear interpolation
2139 Double_t y=(clInner->GetY()+clOuter->GetY())/2;
2140 Double_t z=(clInner->GetZ()+clOuter->GetZ())/2;
2141 if (seedFit.IsValid()) {
2142 // update values once the fit was performed
2143 y=seedFit.GetYat(krMiddle.GetX());
2144 z=seedFit.GetZat(krMiddle.GetX());
2145 }
38d9d609 2146
9d548f3c 2147 n=krMiddle.FindNearest(y,z,roady,roadz);
2148 if (n) break;
2149 }
2150
38d9d609 2151// if (n)
2152// 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,
2153// n->GetX(), n->GetY(),n->GetZ(),n->GetLabel(0),
2154// clInner->GetX(), clInner->GetY(),clInner->GetZ(),clInner->GetLabel(0),
2155// clOuter->GetX(), clOuter->GetY(),clOuter->GetZ(),clOuter->GetLabel(0)
2156// );
2157 return n;
2158}
2159
2160//____________________________________________________________________________________
2161void AliToyMCReconstruction::AddMiddleClusters(AliTPCseed *seed,
2162 const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
2163 Double_t roady, Double_t roadz,
2164 Int_t &nTotalClusters, AliRieman &seedFit)
2165{
2166 //
2167 // Iteratively add the middle clusters
2168 //
2169
2170 // update rieman fit with every second point added
2171 AliTPCclusterMI *clMiddle=FindMiddleCluster(clInner,clOuter,roady,roadz,seedFit);
2172
2173 // break iterative process
9d548f3c 2174 if (!clMiddle || clMiddle->IsUsed()) return;
38d9d609 2175
2176 const Int_t globalRowInner = clInner->GetRow() +(clInner->GetDetector() >35)*63;
2177 const Int_t globalRowMiddle = clMiddle->GetRow()+(clMiddle->GetDetector()>35)*63;
2178 const Int_t globalRowOuter = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
2179
2180 seed->SetClusterPointer(globalRowMiddle,clMiddle);
2181 ++nTotalClusters;
2182// printf(" > Total clusters: %d\n",nTotalClusters);
2183 seedFit.AddPoint(clMiddle->GetX(), clMiddle->GetY(), clMiddle->GetZ(),
2184 TMath::Sqrt(clMiddle->GetSigmaY2()), TMath::Sqrt(clMiddle->GetSigmaZ2()));
2185
9d548f3c 2186 if (seedFit.GetN()>3) {
afc28cd1 2187// printf(" call update: %d (%d)\n",seedFit.GetN(),nTotalClusters);
2188// printf(" Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f -- %d\n",
2189// seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z(), clMiddle->GetLabel(0));
38d9d609 2190 seedFit.Update();
2191 }
2192 if ( seedFit.IsValid() && seedFit.GetChi2()>1000 ) return;
2193
2194 // Add middle clusters towards outer radius
2195 if (globalRowMiddle+1<globalRowOuter) AddMiddleClusters(seed,clMiddle,clOuter,roady,roadz,nTotalClusters,seedFit);
2196
2197 // Add middle clusters towards innter radius
2198 if (globalRowInner+1<globalRowMiddle) AddMiddleClusters(seed,clInner,clMiddle,roady,roadz,nTotalClusters,seedFit);
2199
2200 return;
2201}
2202
2203//____________________________________________________________________________________
9d548f3c 2204Int_t AliToyMCReconstruction::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t iRowInner, Int_t iRowOuter)
38d9d609 2205{
2206 //
2207 // find seeds in a sector, requires iRowInner < iRowOuter
2208 // iRowXXX runs from 0-159, so over IROC and OROC
2209 //
2210
2211 if (!fIsAC) {
2212 AliError("This function requires the sector arrays filled for AC tracking");
9d548f3c 2213 return 0;
38d9d609 2214 }
2215
2216 // swap rows in case they are in the wrong order
2217 if (iRowInner>iRowOuter) {
2218 Int_t tmp=iRowInner;
2219 iRowInner=iRowOuter;
2220 iRowOuter=tmp;
2221 }
2222
a06336b6 2223 if (iRowOuter>158) iRowOuter=158;
2224 if (iRowInner<0) iRowInner=0;
2225
38d9d609 2226 // Define the search roads:
2227 // timeRoadCombinatorics - the maximum time difference used for the
2228 // combinatorics. Since we don't have a 'z-Vertex' estimate this will
2229 // reduce the combinatorics significantly when iterating on one TF
2230 // use a little more than one full drift length of the TPC to allow for
2231 // distortions
2232 //
2233 // timeRoad - the time difference allowed when associating the cluster
2234 // in the middle of the other two used for the initial search
2235 //
2236 // padRoad - the local y difference allowed when associating the middle cluster
38d9d609 2237
2238// Double_t timeRoadCombinatorics = 270./vDrift;
2239// Double_t timeRoad = 20./vDrift;
2240 Double_t timeRoadCombinatorics = 270.;
9d548f3c 2241 Double_t timeRoad = 10.;
2242 Double_t padRoad = 5.;
38d9d609 2243
2244
2245 // fOuterSectorArray runs from 0-95, so from iRowXXX which runs from 0-159
2246 // the number of rows in the IROC has to be subtracted
2247 const Int_t innerRows=fInnerSectorArray->GetNRows();
9d548f3c 2248
2249 AliTPCtrackerSector *arrInnerRow=(iRowInner<63)?fInnerSectorArray:fOuterSectorArray;
2250 AliTPCtrackerSector *arrOuterRow=(iRowOuter<63)?fInnerSectorArray:fOuterSectorArray;
2251
2252 const AliTPCtrackerRow& krInner = arrInnerRow[sec][iRowInner - (iRowInner>62)*innerRows]; // up
2253 const AliTPCtrackerRow& krOuter = arrOuterRow[sec][iRowOuter - (iRowOuter>62)*innerRows]; // down
38d9d609 2254
2255 AliTPCseed *seed = new AliTPCseed;
2256
2257 const Int_t nMaxClusters=iRowOuter-iRowInner+1;
2258// Int_t nScannedClusters = 0;
9d548f3c 2259
2260 Int_t nseedAdded=0;
38d9d609 2261 // loop over all points in the firstand last search row
2262 for (Int_t iOuter=0; iOuter < krOuter; iOuter++) {
2263 const AliTPCclusterMI *clOuter = krOuter[iOuter];
a06336b6 2264 if (clOuter->IsUsed()) continue;
38d9d609 2265 for (Int_t iInner=0; iInner < krInner; iInner++) {
2266 const AliTPCclusterMI *clInner = krInner[iInner];
a06336b6 2267 if (clInner->IsUsed()) continue;
afc28cd1 2268// 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 2269 // check maximum distance for combinatorics
2270 if (TMath::Abs(clOuter->GetZ()-clInner->GetZ())>timeRoadCombinatorics) continue;
afc28cd1 2271// printf(" Is inside one drift\n");
38d9d609 2272
2273 // use rieman fit for seed description
2274 AliRieman seedFit(159);
2275 // add first two points
2276 seedFit.AddPoint(clInner->GetX(), clInner->GetY(), clInner->GetZ(),
2277 TMath::Sqrt(clInner->GetSigmaY2()), TMath::Sqrt(clInner->GetSigmaZ2()));
2278 seedFit.AddPoint(clOuter->GetX(), clOuter->GetY(), clOuter->GetZ(),
2279 TMath::Sqrt(clOuter->GetSigmaY2()), TMath::Sqrt(clOuter->GetSigmaZ2()));
2280
2281 // Iteratively add all clusters in the respective middle
2282 Int_t nFoundClusters=2;
2283 AddMiddleClusters(seed,clInner,clOuter,padRoad,timeRoad,nFoundClusters,seedFit);
afc28cd1 2284// printf(" Clusters attached: %d\n",nFoundClusters);
a06336b6 2285 if (nFoundClusters>2) seedFit.Update();
afc28cd1 2286// printf(" Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f\n",
2287// seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z());
38d9d609 2288
2289 // check for minimum number of assigned clusters and a decent chi2
2290 if ( nFoundClusters<0.5*nMaxClusters || seedFit.GetChi2()>1000 ){
2291 seed->Reset();
2292 continue;
2293 }
2294// printf(" >>> Seed accepted\n");
2295 // add original outer clusters
2296 const Int_t globalRowInner = clInner->GetRow() +(clInner->GetDetector() >35)*63;
2297 const Int_t globalRowOuter = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
2298
2299 seed->SetClusterPointer(globalRowInner,const_cast<AliTPCclusterMI*>(clInner));
2300 seed->SetClusterPointer(globalRowOuter,const_cast<AliTPCclusterMI*>(clOuter));
2301
2302 // set parameters retrieved from AliRieman
2303 Double_t params[5]={0};
2304 Double_t covar[15]={0};
2305 const Double_t alpha=TMath::DegToRad()*(clInner->GetDetector()%18*20.+10.);
2306 const Double_t x=clInner->GetX();
2307 seedFit.GetExternalParameters(x,params,covar);
2308
2309 seed->Set(x,alpha,params,covar);
2310
2311 // set label of the seed. At least 60% of the clusters need the correct label
a06336b6 2312 CookLabel(seed,.6);
2313// printf(" - Label: %d\n",seed->GetLabel());
2314 // mark clusters as being used
2315 MarkClustersUsed(seed);
9d548f3c 2316
2317 seed->SetSeed1(iRowInner);
2318 seed->SetSeed2(iRowOuter);
2319 seed->SetSector(sec);
2320 ++nseedAdded;
2321
2322 seed->SetUniqueID(arr->GetEntriesFast());
2323 seed->SetIsSeeding(kTRUE);
38d9d609 2324
2325 arr->Add(seed);
2326 seed=new AliTPCseed;
9d548f3c 2327
2328 break;
38d9d609 2329 }
2330 }
2331 //delete surplus seed
2332 delete seed;
2333 seed=0x0;
2334
9d548f3c 2335 return nseedAdded;
38d9d609 2336}
2337//____________________________________________________________________________________
2338void AliToyMCReconstruction::MakeSeeds(TObjArray * /*arr*/, Int_t sec, Int_t iRow1, Int_t iRow2)
502eb90b 2339{
2340 //
2341 // Create seeds between i1 and i2 (stored in arr) with TLinearFitter
2342 //
2343 // sec: sector number
2344 // iRow1: upper row
2345 // iRow2: lower row
2346 //
2347
2348 // Create Fitter
2349 static TLinearFitter fitter(3,"pol2");
2350
2351 // Get 3 padrows (iRow1,iMiddle=(iRow1+iRow2)/2,iRow2)
2352 Int_t iMiddle = (iRow1+iRow2)/2;
2353 const AliTPCtrackerRow& krd = fOuterSectorArray[sec][iRow2-fInnerSectorArray->GetNRows()]; // down
2354 const AliTPCtrackerRow& krm = fOuterSectorArray[sec][iMiddle-fInnerSectorArray->GetNRows()]; // middle
2355 const AliTPCtrackerRow& kru = fOuterSectorArray[sec][iRow1-fInnerSectorArray->GetNRows()]; // up
2356
2357 // Loop over 3 cluster possibilities
2358 for (Int_t iu=0; iu < kru; iu++) {
2359 for (Int_t im=0; im < krm; im++) {
2360 for (Int_t id=0; id < krd; id++) {
2361
2362 // clear all points
2363 fitter.ClearPoints();
2364
2365 // add all three points to fitter
2366 Double_t xy[2] = {kru[iu]->GetX(),kru[iu]->GetY()};
2367 Double_t z = kru[iu]->GetZ();
2368 fitter.AddPoint(xy,z);
2369
2370 xy[0] = krm[im]->GetX();
2371 xy[1] = krm[im]->GetY();
2372 z = krm[im]->GetZ();
2373 fitter.AddPoint(xy,z);
2374
2375 xy[0] = krd[id]->GetX();
2376 xy[1] = krd[id]->GetY();
2377 z = krd[id]->GetZ();
2378 fitter.AddPoint(xy,z);
2379
2380 // Evaluate and get parameters
2381 fitter.Eval();
2382
2383 // how to get the other clusters now?
2384 // ...
2385
2386 }
2387 }
2388 }
2389}
38d9d609 2390
2391//____________________________________________________________________________________
a06336b6 2392void AliToyMCReconstruction::InitStreamer(TString addName, Int_t level)
38d9d609 2393{
2394 //
2395 // initilise the debug streamer and set the logging level
2396 // use a default naming convention
2397 //
a06336b6 2398
2399 delete fStreamer;
2400 fStreamer=0x0;
2401
2402 if (addName.IsNull()) addName=".dummy";
38d9d609 2403
a06336b6 2404 if (!fTree) return;
2405
48e13991 2406 TString debugName=gSystem->BaseName(fInputFile->GetName());
a06336b6 2407 debugName.ReplaceAll(".root","");
2408 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
2409 fUseMaterial,fIdealTracking,fClusterType,
2410 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
2411 debugName.Append(addName);
2412 debugName.Append(".root");
2413
2414 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
2415 fStreamer=new TTreeSRedirector(debugName.Data());
2416 fStreamer->SetUniqueID(level);
2417
2418 gROOT->cd();
38d9d609 2419}
2420
2421//____________________________________________________________________________________
a06336b6 2422void AliToyMCReconstruction::ConnectInputFile(const char* file, Int_t nmaxEv)
38d9d609 2423{
2424 //
2425 // connect the tree and event pointer from the input file
2426 //
2427
2428 delete fInputFile;
2429 fInputFile=0x0;
2430 fEvent=0x0;
2431 fTree=0;
2432
2433 fInputFile=TFile::Open(file);
2434 if (!fInputFile || !fInputFile->IsOpen() || fInputFile->IsZombie()) {
2435 delete fInputFile;
2436 fInputFile=0x0;
afc28cd1 2437 AliError(Form("ERROR: couldn't open the file '%s'\n", file));
38d9d609 2438 return;
2439 }
2440
2441 fTree=(TTree*)fInputFile->Get("toyMCtree");
2442 if (!fTree) {
afc28cd1 2443 AliError(Form("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file));
38d9d609 2444 return;
2445 }
2446
2447 fEvent=0x0;
2448 fTree->SetBranchAddress("event",&fEvent);
2449
2450 gROOT->cd();
a06336b6 2451
2452 fNmaxEvents=fTree->GetEntries();
2453 if (nmaxEv>-1) fNmaxEvents=TMath::Min(nmaxEv,fNmaxEvents);
38d9d609 2454
a06336b6 2455 // setup space charge map from Userinfo of the tree
2456 InitSpaceCharge();
2457
2458 // setup the track maps
38d9d609 2459 SetupTrackMaps();
2460}
2461
2462//____________________________________________________________________________________
2463void AliToyMCReconstruction::Cleanup()
2464{
2465 //
2466 // Cleanup input data
2467 //
2468
a06336b6 2469 if (fStreamer) delete fStreamer;
38d9d609 2470 fStreamer=0x0;
2471
2472 delete fEvent;
2473 fEvent = 0x0;
2474
a06336b6 2475// delete fTree;
38d9d609 2476 fTree=0x0;
2477
2478 delete fInputFile;
2479 fInputFile=0x0;
2480
2481}
2482
2483//____________________________________________________________________________________
2484void AliToyMCReconstruction::SetupTrackMaps()
2485{
2486 //
2487 //
2488 //
2489
2490 fMapTrackEvent.Delete();
2491 fMapTrackTrackInEvent.Delete();
2492
2493 if (!fTree) {
2494 AliError("Tree not connected");
2495 return;
2496 }
2497
a06336b6 2498 Int_t nmaxEv=fTree->GetEntries();
2499 if (fNmaxEvents>-1) nmaxEv=fNmaxEvents;
2500
2501 for (Int_t iev=0; iev<nmaxEv; ++iev) {
38d9d609 2502 fTree->GetEvent(iev);
2503 if (!fEvent) continue;
2504
2505 const Int_t ntracks=fEvent->GetNumberOfTracks();
2506 if (fMapTrackEvent.GetSize()+ntracks>=fMapTrackEvent.Capacity()) fMapTrackEvent.Expand(2*fMapTrackEvent.Capacity());
2507 if (fMapTrackTrackInEvent.GetSize()+ntracks>=fMapTrackTrackInEvent.Capacity()) fMapTrackTrackInEvent.Expand(2*fMapTrackTrackInEvent.Capacity());
2508
2509 for (Int_t itrack=0; itrack<ntracks; ++itrack){
2510 Int_t label=fEvent->GetTrack(itrack)->GetUniqueID();
2511
2512 fMapTrackEvent.Add(label,iev);
2513 fMapTrackTrackInEvent.Add(label,itrack);
2514 }
2515 }
2516
2517}
a06336b6 2518
2519//____________________________________________________________________________________
2520void AliToyMCReconstruction::CookLabel(AliTPCseed *seed, Double_t fraction, Int_t info[5])
2521{
2522 //
2523 //
2524 //
2525
2526 Int_t labels[159]={0};
2527// Long64_t posMaxLabel=-1;
2528 Int_t nMaxLabel=0; // clusters from maximum label
2529 Int_t nMaxLabel2=0; // clusters from second maximum
2530 Int_t nlabels=0;
2531 Int_t maxLabel=-1; // label with most clusters
2532 Int_t maxLabel2=-1; // label with second most clusters
2533 Int_t nclusters=0;
2534 TExMap labelMap(159);
2535
2536 for (Int_t icl=0; icl<159; ++icl) {
2537 const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2538 if (!cl) continue;
2539 ++nclusters;
2540
2541 const Int_t clLabel=cl->GetLabel(0);
2542 // a not assinged value returns 0, so we need to add 1 and subtract it afterwards
2543 Long64_t labelPos=labelMap.GetValue(clLabel);
2544
2545 if (!labelPos) {
2546 labelPos=nlabels+1;
2547 labelMap.Add(clLabel,labelPos);
2548 ++nlabels;
2549 }
2550 --labelPos;
2551
2552 const Int_t nCurrentLabel=++labels[labelPos];
2553 if (nCurrentLabel>nMaxLabel) {
2554 nMaxLabel2=nMaxLabel;
2555 nMaxLabel=nCurrentLabel;
2556// posMaxLabel=labelPos;
2557 maxLabel2=maxLabel;
2558 maxLabel=clLabel;
2559 }
2560 }
2561
2562 if (Double_t(nMaxLabel)/nclusters<fraction) maxLabel=-maxLabel;
2563
2564 seed->SetLabel(maxLabel);
2565
2566 if (info) {
2567 info[0]=nMaxLabel;
2568 info[1]=nMaxLabel2;
2569 info[2]=maxLabel2;
2570 info[3]=nclusters;
2571 info[4]=nlabels;
2572 }
2573}
2574
2575
2576//____________________________________________________________________________________
9d548f3c 2577void AliToyMCReconstruction::DumpSeedInfo(TObjArray *arr)
a06336b6 2578{
2579
2580 // for debugging
2581 if (!fStreamer || !fTree) return;
2582 // swap rows in case they are in the wrong order
9d548f3c 2583 AliInfo(Form("Number of Seeds: %d (out of %d tracks)",arr->GetEntriesFast(),fMapTrackEvent.GetSize()));
a06336b6 2584
2585 //loop over all events and tracks and try to associate the seed to the track
2586 for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed){
2587 AliTPCseed *seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
a06336b6 2588
2589 // get original track
9d548f3c 2590 Int_t seedLabel=seed->GetLabel();
a06336b6 2591 Int_t iev=fMapTrackEvent.GetValue(TMath::Abs(seedLabel));
2592 Int_t itr=fMapTrackTrackInEvent.GetValue(TMath::Abs(seedLabel));
2593
2594 fTree->GetEvent(iev);
2595 const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
2596
9d548f3c 2597 DumpSeedInfo(toyTrack,seed);
2598 }
2599}
2600
2601//____________________________________________________________________________________
2602void AliToyMCReconstruction::DumpTrackInfo(TObjArray *arr)
2603{
2604
2605 // for debugging
2606 if (!fStreamer || !fTree) return;
2607 // swap rows in case they are in the wrong order
2608 AliInfo(Form("Number of Seeds: %d (out of %d tracks)",arr->GetEntriesFast(),fMapTrackEvent.GetSize()));
2609
2610 //loop over all events and tracks and try to associate the seed to the track
2611 AliTPCseed dummySeed;
2612 for (Int_t iev=0; iev<fNmaxEvents; ++iev) {
2613 fTree->GetEvent(iev);
2614 const Int_t ntracks=fEvent->GetNumberOfTracks();
2615 for (Int_t itr=0; itr<ntracks; ++itr) {
2616 const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
2617
2618 Bool_t foundSeed=kFALSE;
2619 for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed) {
2620 AliTPCseed *seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
2621 const UInt_t tmpLabel=TMath::Abs(seed->GetLabel());
2622 if (toyTrack->GetUniqueID()!=tmpLabel) continue;
2623
2624 // dump all seeds with the correct label
2625 DumpSeedInfo(toyTrack,seed);
2626 foundSeed=kTRUE;
2627 }
2628
2629 if (!foundSeed) DumpSeedInfo(toyTrack,&dummySeed);
a06336b6 2630
a06336b6 2631 }
9d548f3c 2632 }
2633}
a06336b6 2634
9d548f3c 2635//____________________________________________________________________________________
2636void AliToyMCReconstruction::DumpSeedInfo(const AliToyMCTrack *toyTrack, AliTPCseed *seed)
2637{
2638 //
2639 //
2640 //
2641
2642 const Double_t kMaxSnp = 0.85;
2643 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2644 Float_t vDrift=GetVDrift();
2645 Float_t zLength=GetZLength(0);
2646
2647 AliExternalTrackParam tOrig(*toyTrack);
2648 AliExternalTrackParam tOrigITS(*toyTrack);
2649
2650 // propagate original track to ITS last layer
e83fd282 2651// Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
2652 const Double_t iFCRadius = 83.5; //radius constants found in AliTPCCorrection.cxx
2653 Double_t lastLayerITS = iFCRadius; // its track propgated to inner TPC wall
9d548f3c 2654 AliTrackerBase::PropagateTrackTo(&tOrigITS,lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2655
2656 AliExternalTrackParam dummyParam;
2657 Bool_t isDummy=kFALSE;
2658 //create refitted track, this will also set the fTime0
2659 AliExternalTrackParam *track=GetRefittedTrack(*seed);
2660 if (!track) {
2661 track=&dummyParam;
2662 isDummy=kTRUE;
2663 }
2664 AliTrackerBase::PropagateTrackTo(track,0,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2665 track->Rotate(tOrig.GetAlpha());
2666 AliTrackerBase::PropagateTrackTo(track,0,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2667
2668 // rotate fitted track to the frame of the original track and propagate to same reference
2669 AliExternalTrackParam trackITS(*track);
2670 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2671 trackITS.Rotate(tOrigITS.GetAlpha());
2672 AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2673
2674 Int_t seedSec=seed->GetSector();
2675 Int_t seedID =seed->GetUniqueID();
2676 //
2677 //count findable and found clusters in the seed
2678 //
2679 Int_t iRowInner=seed->GetSeed1();
2680 Int_t iRowOuter=seed->GetSeed2();
2681
2682 Int_t nClustersSeedMax=iRowOuter-iRowInner+1;
2683 Int_t nClustersFindable=0;
2684 Int_t nClustersSeed=0;
2685
2686 Int_t ncls=(fClusterType==0)?toyTrack->GetNumberOfSpacePoints():toyTrack->GetNumberOfDistSpacePoints();
2687
2688 Int_t rowInner=iRowInner-(iRowInner>62)*63;
2689 Int_t rowOuter=iRowOuter-(iRowOuter>62)*63;
2690
2691 //findable in the current seed sector
2692 Int_t lastROC=-1;
2693 Int_t rocMaxCl=-1;
2694 Int_t nCrossedROCs=0;
2695 Int_t nMaxClROC=0;
2696 Int_t nclROC=0;
2697 Int_t row1=-1;
2698 Int_t row2=-1;
2699 Int_t row1Maxcl=-1;
2700 Int_t row2Maxcl=-1;
2701 for (Int_t icl=0; icl<ncls; ++icl) {
2702 const AliTPCclusterMI *cl=(fClusterType==0)?toyTrack->GetSpacePoint(icl):toyTrack->GetDistortedSpacePoint(icl);
2703 const Int_t roc=cl->GetDetector();
2704 const Int_t row=cl->GetRow();
2705 const Int_t rowGlobal=row+(roc>35)*63;
2706
2707 AliTPCclusterMI *seedCl = seed->GetClusterPointer(rowGlobal);
2708 if (seedCl) {
2709 AliTPCclusterMI *clc=const_cast<AliTPCclusterMI*>(cl);
2710 if (seedCl->GetDetector()==roc&&seedCl->IsUsed()) clc->Use();
2711 clc->SetLabel(seedID,1);
2712
a06336b6 2713 }
9d548f3c 2714
2715// if (row1<0) row1=rowGlobal;
2716
2717 if ( (roc%36) != (lastROC%36)) {
2718 ++nCrossedROCs;
2719 if (nclROC>nMaxClROC) {
2720 nMaxClROC=nclROC;
2721 rocMaxCl=lastROC;
2722 row1Maxcl=row1;
2723 row2Maxcl=row2;
2724 }
2725 lastROC=roc%36;
2726 nclROC=0;
2727 row1=rowGlobal;
2728 }
2729 ++nclROC;
2730 row2=rowGlobal;
a06336b6 2731
9d548f3c 2732 if ( (roc%36)!=(seedSec%36) ) continue;
2733// if ( (row<rowInner) || (row>rowOuter) ) continue;
2734 ++nClustersFindable;
2735
a06336b6 2736 }
9d548f3c 2737
2738 if (nclROC>nMaxClROC) {
2739 rocMaxCl=lastROC;
2740 nMaxClROC=nclROC;
2741 row1Maxcl=row1;
2742 row2Maxcl=row2;
2743 }
2744
2745 Int_t firstRow=160;
2746 Int_t lastRow=0;
2747 Int_t nClustersInTrack=0;
2748 //found in seed
2749 for (Int_t icl=0; icl<159; ++icl) {
2750 const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2751 if (!cl) continue;
2752 ++nClustersInTrack;
2753 const Int_t row=cl->GetRow();
2754 const Int_t rowGlobal=row+(cl->GetDetector()>35)*63;
2755 if (rowGlobal>lastRow) lastRow=rowGlobal;
2756 if (rowGlobal<firstRow) firstRow=rowGlobal;
2757 if ( (row<rowInner) || (row>rowOuter) ) continue;
2758 ++nClustersSeed;
2759 }
2760
2761 Float_t z0=fEvent->GetZ();
2762 Float_t t0=fEvent->GetT0();
2763
2764 Int_t ctype(fCorrectionType);
2765
2766 Int_t info[5]={0};
2767 CookLabel(seed,.6,info);
2768 Int_t seedLabel=seed->GetLabel();
2769
2770 Int_t labelOrig=toyTrack->GetUniqueID();
2771
2772 AliToyMCTrack *tr2 = const_cast<AliToyMCTrack*>(toyTrack);
2773
2774 (*fStreamer) << "Seeds" <<
2775 // "iev=" << iev <<
2776 // "iseed=" << iseed <<
2777 // "itr=" << itr <<
2778 "z0=" << z0 <<
2779 "t0=" << t0 <<
2780 "vDrift=" << vDrift <<
2781 "zLength=" << zLength <<
2782 "fTime0=" << fTime0 <<
2783 "clsType=" << fClusterType <<
2784 "corrType=" << ctype <<
2785
2786 "tOrig.=" << &tOrig <<
2787 "tOrigITS.=" << &tOrigITS <<
2788
2789 "to.nclFindable=" << nClustersFindable <<
2790 "to.nclTot=" << ncls <<
2791 "to.label=" << labelOrig <<
2792 "to.nCrossedROCs="<< nCrossedROCs <<
2793 "to.rocMax=" << rocMaxCl <<
2794 "to.rocMaxNcl=" << nMaxClROC <<
2795 "to.row1Max=" << row1Maxcl <<
2796 "to.row2Max=" << row2Maxcl <<
2797
2798 "track.=" << track <<
2799 "trackITS.=" << &trackITS <<
2800
2801 "s.RowInner=" << iRowInner <<
2802 "s.RowOuter=" << iRowOuter <<
2803 "s.nclMax=" << nClustersSeedMax <<
2804 "s.ncl=" << nClustersSeed <<
2805 "s.ID=" << seedID <<
2806
2807 "tr.firstClRow=" << firstRow <<
2808 "tr.lastClRow=" << lastRow <<
2809 "tr.ncl=" << nClustersInTrack <<
2810 "tr.label=" << seedLabel <<
2811
2812 "tr.LabelNcl=" << info[0] <<
2813 "tr.Label2Ncl=" << info[1] <<
2814 "tr.Label2=" << info[2] <<
2815 "tr.nclTot=" << info[3] <<
2816 "tr.Nlabels=" << info[4] <<
2817
2818 "tr.Sec=" << seedSec <<
2819
2820 "seed.=" << seed <<
2821 "toyTrack.=" << tr2 <<
2822 "\n";
2823
2824 if (!isDummy) delete track;
a06336b6 2825}
2826
2827//____________________________________________________________________________________
2828void AliToyMCReconstruction::MarkClustersUsed(AliTPCseed *seed)
2829{
2830 //
2831 //
2832 //
2833
2834 for (Int_t icl=0; icl<159; ++icl) {
2835 AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2836 if (cl) cl->Use();
2837 }
2838}
2839
9d548f3c 2840//____________________________________________________________________________________
2841void AliToyMCReconstruction::ResetClustersZtoTime(AliTPCseed *seed)
2842{
2843 //
2844 //
2845 //
2846
2847 for (Int_t icl=0; icl<159; ++icl) {
2848 AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2849 if (cl) cl->SetZ(cl->GetTimeBin());
2850 }
2851}
2852
a06336b6 2853//____________________________________________________________________________________
2854void AliToyMCReconstruction::DumpTracksToTree(const char* file)
2855{
2856 //
2857 //
2858 //
2859 ConnectInputFile(file);
2860 if (!fTree) return;
2861
2862 delete fStreamer;
2863 fStreamer=0x0;
2864
2865 TString debugName=fInputFile->GetName();
2866 debugName.ReplaceAll(".root",".AllTracks.root");
2867
2868 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
2869 fStreamer=new TTreeSRedirector(debugName.Data());
2870
2871 for (Int_t iev=0;iev<fNmaxEvents;++iev){
2872 fTree->GetEvent(iev);
2873 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks();++itr){
2874 AliToyMCTrack *toyTrack=const_cast<AliToyMCTrack*>(fEvent->GetTrack(itr));
2875 Int_t trackID=toyTrack->GetUniqueID();
2876
2877 (*fStreamer) << "Tracks" <<
2878 "iev=" << iev <<
2879 "itr=" << itr <<
2880 "trackID=" << trackID <<
2881 "track.=" << toyTrack <<
2882 "\n";
2883
2884 }
2885 }
2886
2887 Cleanup();
2888}
9d548f3c 2889
2890//____________________________________________________________________________________
2891void AliToyMCReconstruction::SetRieman(const AliTPCseed &seed, AliRieman &rieman)
2892{
2893 //
2894 //
2895 //
2896
2897 rieman.Reset();
2898 for (Int_t icl=0; icl<159; ++icl) {
2899 const AliTPCclusterMI *cl=seed.GetClusterPointer(icl);
2900 if (!cl) continue;
2901 rieman.AddPoint(cl->GetX(), cl->GetY(), cl->GetZ(),
2902 TMath::Sqrt(cl->GetSigmaY2()), TMath::Sqrt(cl->GetSigmaZ2()));
2903 }
2904 rieman.Update();
2905}
2906
2907//____________________________________________________________________________________
2908void AliToyMCReconstruction::CopyRieman(const AliRieman &from, AliRieman &to)
2909{
2910 //
2911 //
2912 //
2913
2914 if (to.GetCapacity()<from.GetCapacity()) return;
2915 to.Reset();
2916
2917 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]);
2918}
2919