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