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