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