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