]> 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
First Steps for tracking from seed of associated clusters:
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCReconstruction.cxx
CommitLineData
d1cf83f5 1
2#include <TDatabasePDG.h>
3#include <TString.h>
4a777885 4#include <TSystem.h>
5#include <TROOT.h>
6#include <TFile.h>
9e98dea8 7#include <TPRegexp.h>
d1cf83f5 8
9#include <AliExternalTrackParam.h>
10#include <AliTPCcalibDB.h>
11#include <AliTPCclusterMI.h>
12#include <AliTPCSpaceCharge3D.h>
13#include <AliTrackerBase.h>
14#include <AliTrackPointArray.h>
15#include <AliLog.h>
16#include <AliTPCParam.h>
17#include <AliTPCROC.h>
18#include <TTreeStream.h>
50f68b1a 19#include <AliTPCReconstructor.h>
20#include <AliTPCTransform.h>
510cfcff 21#include <AliTPCtracker.h>
22#include <AliTPCtrackerSector.h>
d1cf83f5 23
24#include "AliToyMCTrack.h"
4a777885 25#include "AliToyMCEvent.h"
d1cf83f5 26
27#include "AliToyMCReconstruction.h"
28
29
30//____________________________________________________________________________________
31AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
32, fSeedingRow(140)
33, fSeedingDist(10)
34, fClusterType(0)
35, fCorrectionType(kNoCorrection)
36, fDoTrackFit(kTRUE)
37, fUseMaterial(kFALSE)
600eaa0d 38, fIdealTracking(kFALSE)
d1cf83f5 39, fTime0(-1)
600eaa0d 40, fCreateT0seed(kFALSE)
d1cf83f5 41, fStreamer(0x0)
42, fTree(0x0)
4a777885 43, fEvent(0x0)
d1cf83f5 44, fTPCParam(0x0)
45, fSpaceCharge(0x0)
5993ed4f 46, fkNSectorInner(18) // hard-coded to avoid loading the parameters before
47, fInnerSectorArray(0x0)
48, fkNSectorOuter(18) // hard-coded to avoid loading the parameters before
49, fOuterSectorArray(0x0)
d1cf83f5 50{
51 //
52 // ctor
53 //
54 fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
5993ed4f 55
d1cf83f5 56}
57
58//____________________________________________________________________________________
59AliToyMCReconstruction::~AliToyMCReconstruction()
60{
61 //
62 // dtor
63 //
64
65 delete fStreamer;
4a777885 66// delete fTree;
67}
68
69//____________________________________________________________________________________
70void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
71{
72 //
510cfcff 73 // Recostruction from associated clusters
4a777885 74 //
75
76 TFile f(file);
77 if (!f.IsOpen() || f.IsZombie()) {
78 printf("ERROR: couldn't open the file '%s'\n", file);
79 return;
80 }
81
82 fTree=(TTree*)f.Get("toyMCtree");
83 if (!fTree) {
84 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
85 return;
86 }
87
88 fEvent=0x0;
89 fTree->SetBranchAddress("event",&fEvent);
90
91 // read spacecharge from the Userinfo ot the tree
92 InitSpaceCharge();
93
94 TString debugName=file;
95 debugName.ReplaceAll(".root","");
9e98dea8 96 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
600eaa0d 97 fUseMaterial,fIdealTracking,fClusterType,
98 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
4a777885 99 debugName.Append(".debug.root");
100
101 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
102 if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
103
104 gROOT->cd();
105
106 static AliExternalTrackParam dummySeedT0;
107 static AliExternalTrackParam dummySeed;
108 static AliExternalTrackParam dummyTrack;
109
110 AliExternalTrackParam t0seed;
111 AliExternalTrackParam seed;
112 AliExternalTrackParam track;
113 AliExternalTrackParam tOrig;
114
115 AliExternalTrackParam *dummy;
116
117 Int_t maxev=fTree->GetEntries();
118 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
119
120 for (Int_t iev=0; iev<maxev; ++iev){
121 printf("============== Processing Event %6d =================\n",iev);
122 fTree->GetEvent(iev);
123 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
124 printf(" > ====== Processing Track %6d ======== \n",itr);
125 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
126 tOrig = *tr;
127
128
129 // set dummy
130 t0seed = dummySeedT0;
131 seed = dummySeed;
132 track = dummyTrack;
133
134 Float_t z0=fEvent->GetZ();
135 Float_t t0=fEvent->GetT0();
136
137 Float_t vdrift=GetVDrift();
138 Float_t zLength=GetZLength(0);
139
600eaa0d 140 // crate time0 seed, steered by fCreateT0seed
4a777885 141 printf("t0 seed\n");
142 fTime0=-1.;
600eaa0d 143 fCreateT0seed=kTRUE;
4a777885 144 dummy = GetSeedFromTrack(tr);
145
146 if (dummy) {
147 t0seed = *dummy;
148 delete dummy;
149
150 // crate real seed using the time 0 from the first seed
600eaa0d 151 // set fCreateT0seed now to false to get the seed in z coordinates
4a777885 152 fTime0 = t0seed.GetZ()-zLength/vdrift;
600eaa0d 153 fCreateT0seed = kFALSE;
154 printf("seed (%.2g)\n",fTime0);
4a777885 155 dummy = GetSeedFromTrack(tr);
156 if (dummy) {
157 seed = *dummy;
158 delete dummy;
159
160 // create fitted track
161 if (fDoTrackFit){
162 printf("track\n");
163 dummy = GetFittedTrackFromSeed(tr, &seed);
164 track = *dummy;
165 delete dummy;
166 }
600eaa0d 167
168 // propagate seed to 0
169 const Double_t kMaxSnp = 0.85;
170 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
171// AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
172
4a777885 173 }
174 }
175
176 Int_t ctype(fCorrectionType);
177
178 if (fStreamer) {
179 (*fStreamer) << "Tracks" <<
180 "iev=" << iev <<
181 "z0=" << z0 <<
182 "t0=" << t0 <<
600eaa0d 183 "fTime0=" << fTime0 <<
4a777885 184 "itr=" << itr <<
185 "clsType=" << fClusterType <<
186 "corrType=" << ctype <<
187 "seedRow=" << fSeedingRow <<
188 "seedDist=" << fSeedingDist <<
189 "vdrift=" << vdrift <<
190 "zLength=" << zLength <<
191 "t0seed.=" << &t0seed <<
192 "seed.=" << &seed <<
193 "track.=" << &track <<
194 "tOrig.=" << &tOrig <<
195 "\n";
196 }
197
198
199 }
200 }
201
202 delete fStreamer;
203 fStreamer=0x0;
600eaa0d 204
205 delete fEvent;
206 fEvent = 0x0;
4a777885 207
600eaa0d 208 delete fTree;
209 fTree=0x0;
210 f.Close();
d1cf83f5 211}
212
510cfcff 213
214//____________________________________________________________________________________
215void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
216{
217 //
5993ed4f 218 // Reconstruction for seed from associated clusters, but array of clusters:
219 // Step 1) Filling of cluster arrays
220 // Step 2) Seeding from clusters associated to tracks
221 // Step 3) Free track reconstruction using all clusters
510cfcff 222 //
223
224 TFile f(file);
225 if (!f.IsOpen() || f.IsZombie()) {
226 printf("ERROR: couldn't open the file '%s'\n", file);
227 return;
228 }
229
230 fTree=(TTree*)f.Get("toyMCtree");
231 if (!fTree) {
232 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
233 return;
234 }
235
236 fEvent=0x0;
237 fTree->SetBranchAddress("event",&fEvent);
238
239 // read spacecharge from the Userinfo ot the tree
240 InitSpaceCharge();
241
242 TString debugName=file;
243 debugName.ReplaceAll(".root","");
244 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
245 fUseMaterial,fIdealTracking,fClusterType,
246 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
247 debugName.Append(".allClusters.debug.root");
248
249 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
250 if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
251
252 gROOT->cd();
253
254 static AliExternalTrackParam dummySeedT0;
255 static AliExternalTrackParam dummySeed;
256 static AliExternalTrackParam dummyTrack;
257
258 AliExternalTrackParam t0seed;
259 AliExternalTrackParam seed;
5993ed4f 260 AliExternalTrackParam track;
510cfcff 261 AliExternalTrackParam tOrig;
262
263 // cluster array of all sectors
5993ed4f 264 fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
265 fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
510cfcff 266
267 for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
268 for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
269
270 Int_t count[72][96] = { {0} , {0} };
271
272
273
274 AliExternalTrackParam *dummy;
275
276 Int_t maxev=fTree->GetEntries();
277 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
278
279
5993ed4f 280 // ===========================================================================================
281 // Loop 1: Fill AliTPCtrackerSector structure
282 // ===========================================================================================
283 for (Int_t iev=0; iev<maxev; ++iev){
284 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
285 fTree->GetEvent(iev);
286 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
287 printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
288 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
289
290 // number of clusters to loop over
291 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
292
293 for(Int_t icl=0; icl<ncls; ++icl){
294
295 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
296 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
297 if (!cl) continue;
298
299 Int_t sec = cl->GetDetector();
300 Int_t row = cl->GetRow();
301
302 // set cluster time to cluster Z
303 cl->SetZ(cl->GetTimeBin());
304
305 // fill arrays for inner and outer sectors (A/C side handled internally)
306 if (sec<fkNSectorInner*2){
307 fInnerSectorArray[sec%fkNSectorInner].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
308 }
309 else{
310 fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
311 }
312
313 ++count[sec][row];
314 }
315 }
316 }
317
318 // fill the arrays completely
319 LoadOuterSectors();
320 LoadInnerSectors();
321
322
323 // ===========================================================================================
324 // Loop 2: Seeding from clusters associated to tracks
325 // TODO: Implement tracking from given seed!
326 // ===========================================================================================
327 for (Int_t iev=0; iev<maxev; ++iev){
328 printf("============== Processing Event %6d =================\n",iev);
329 fTree->GetEvent(iev);
330 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
331 printf(" > ====== Processing Track %6d ======== \n",itr);
332 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
333 tOrig = *tr;
334
335
336 // set dummy
337 t0seed = dummySeedT0;
338 seed = dummySeed;
339 track = dummyTrack;
340
341 Float_t z0=fEvent->GetZ();
342 Float_t t0=fEvent->GetT0();
343
344 Float_t vdrift=GetVDrift();
345 Float_t zLength=GetZLength(0);
346
347 // crate time0 seed, steered by fCreateT0seed
348 printf("t0 seed\n");
349 fTime0=-1.;
350 fCreateT0seed=kTRUE;
351 dummy = GetSeedFromTrack(tr);
352
353 if (dummy) {
354 t0seed = *dummy;
355 delete dummy;
356
357 // crate real seed using the time 0 from the first seed
358 // set fCreateT0seed now to false to get the seed in z coordinates
359 fTime0 = t0seed.GetZ()-zLength/vdrift;
360 fCreateT0seed = kFALSE;
361 printf("seed (%.2g)\n",fTime0);
362 dummy = GetSeedFromTrack(tr);
363 if (dummy) {
364 seed = *dummy;
365 delete dummy;
366
367 // create fitted track
368 if (fDoTrackFit){
369 printf("track\n");
370 dummy = GetFittedTrackFromSeedAllClusters(tr, &seed);
371 track = *dummy;
372 delete dummy;
373 }
374
375 // propagate seed to 0
376 const Double_t kMaxSnp = 0.85;
377 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
378 AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
379
380 }
381 }
382
383 Int_t ctype(fCorrectionType);
384
385 if (fStreamer) {
386 (*fStreamer) << "Tracks" <<
387 "iev=" << iev <<
388 "z0=" << z0 <<
389 "t0=" << t0 <<
390 "fTime0=" << fTime0 <<
391 "itr=" << itr <<
392 "clsType=" << fClusterType <<
393 "corrType=" << ctype <<
394 "seedRow=" << fSeedingRow <<
395 "seedDist=" << fSeedingDist <<
396 "vdrift=" << vdrift <<
397 "zLength=" << zLength <<
398 "t0seed.=" << &t0seed <<
399 "seed.=" << &seed <<
400 "track.=" << &track <<
401 "tOrig.=" << &tOrig <<
402 "\n";
403 }
404
405
406 }
407 }
408
409
410 delete fStreamer;
411 fStreamer=0x0;
412
413 delete fEvent;
414 fEvent = 0x0;
415
416 delete fTree;
417 fTree=0x0;
418 f.Close();
419}
420
421//____________________________________________________________________________________
422void AliToyMCReconstruction::RunRecoAllClustersStandardTracking(const char* file, Int_t nmaxEv)
423{
424 //
425 // Reconstruction for seed from associated clusters, but array of clusters
426 // Step 1) Filling of cluster arrays
427 // Step 2) Use the standard tracking: AliTPCtracker::Clusters2Tracks();
428 //
429
430 TFile f(file);
431 if (!f.IsOpen() || f.IsZombie()) {
432 printf("ERROR: couldn't open the file '%s'\n", file);
433 return;
434 }
435
436 fTree=(TTree*)f.Get("toyMCtree");
437 if (!fTree) {
438 printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
439 return;
440 }
441
442 fEvent=0x0;
443 fTree->SetBranchAddress("event",&fEvent);
444
445 // read spacecharge from the Userinfo ot the tree
446 InitSpaceCharge();
447
448 TString debugName=file;
449 debugName.ReplaceAll(".root","");
450 debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
451 fUseMaterial,fIdealTracking,fClusterType,
452 Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
453 debugName.Append(".allClusters.debug.root");
454
455 gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
456 if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
457
458 gROOT->cd();
459
460 // cluster array of all sectors
461 fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
462 fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
463
464 for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
465 for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
466
467 Int_t count[72][96] = { {0} , {0} };
468
469 AliExternalTrackParam *dummy;
470 AliExternalTrackParam *track;
471
472 Int_t maxev=fTree->GetEntries();
473 if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
474
475
510cfcff 476 // ===========================================================================================
477 // Loop 1: Fill AliTPCtrackerSector structure
478 // ===========================================================================================
479 for (Int_t iev=0; iev<maxev; ++iev){
480 printf("============== Fill Clusters: Processing Event %6d =================\n",iev);
481 fTree->GetEvent(iev);
482 for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
483 printf(" > ====== Fill Clusters: Processing Track %6d ======== \n",itr);
484 const AliToyMCTrack *tr=fEvent->GetTrack(itr);
485
486 // number of clusters to loop over
487 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
488
489 for(Int_t icl=0; icl<ncls; ++icl){
490
50f68b1a 491 AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
492 if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
510cfcff 493 if (!cl) continue;
494
495 Int_t sec = cl->GetDetector();
496 Int_t row = cl->GetRow();
497
50f68b1a 498 // set cluster time to cluster Z
499 cl->SetZ(cl->GetTimeBin());
500
510cfcff 501 // fill arrays for inner and outer sectors (A/C side handled internally)
502 if (sec<fkNSectorInner*2){
503 fInnerSectorArray[sec%fkNSectorInner].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
504 }
505 else{
506 fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
507 }
508
509 ++count[sec][row];
510 }
511 }
512 }
513
510cfcff 514 // ===========================================================================================
5993ed4f 515 // Loop 2: Use the full TPC tracker
50f68b1a 516 // TODO: - check tracking configuration
517 // - add clusters and original tracks to output (how?)
510cfcff 518 // ===========================================================================================
510cfcff 519
50f68b1a 520 // settings (TODO: find the correct settings)
521 AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
522 tpcRecoParam->SetDoKinks(kFALSE);
523 AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
524 //tpcRecoParam->Print();
525
526 // need AliTPCReconstructor for parameter settings in AliTPCtracker
527 AliTPCReconstructor *tpcRec = new AliTPCReconstructor();
528 tpcRec->SetRecoParam(tpcRecoParam);
529
530 // AliTPCtracker
531 AliTPCtracker *tpcTracker = new AliTPCtracker(fTPCParam);
532 tpcTracker->SetDebug(10);
533
534 // set sector arrays
535 tpcTracker->SetTPCtrackerSectors(fInnerSectorArray,fOuterSectorArray);
536 tpcTracker->LoadInnerSectors();
537 tpcTracker->LoadOuterSectors();
538
539 // tracking
540 tpcTracker->Clusters2Tracks();
5993ed4f 541 //tpcTracker->PropagateForward();
50f68b1a 542 TObjArray *trackArray = tpcTracker->GetSeeds();
543
544 for(Int_t iTracks = 0; iTracks < trackArray->GetEntriesFast(); ++iTracks){
545 printf(" > ====== Fill Track %6d ======== \n",iTracks);
546
547 track = (AliExternalTrackParam*)(trackArray->At(iTracks));
548
549 if (fStreamer) {
550 (*fStreamer) << "Tracks" <<
551 "track.=" << track <<
552 "\n";
553 }
554 }
555
510cfcff 556 delete fStreamer;
557 fStreamer=0x0;
558
559 delete fEvent;
560 fEvent = 0x0;
561
562 delete fTree;
563 fTree=0x0;
564 f.Close();
565}
566
567
d1cf83f5 568//____________________________________________________________________________________
569AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
570{
571 //
572 // if we don't have a valid time0 informaion (fTime0) available yet
573 // assume we create a seed for the time0 estimate
574 //
575
576 // seed point informaion
577 AliTrackPoint seedPoint[3];
578 const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
579
580 // number of clusters to loop over
581 const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
582
583 UChar_t nextSeedRow=fSeedingRow;
584 Int_t nseeds=0;
585
586 //assumes sorted clusters
587 for (Int_t icl=0;icl<ncls;++icl) {
588 const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
589 if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
590 if (!cl) continue;
591 // use row in sector
592 const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
593 // skip clusters without proper pad row
594 if (row>200) continue;
595
596 //check seeding row
597 // if we are in the last row and still miss a seed we use the last row
598 // even if the row spacing will not be equal
599 if (row>=nextSeedRow || icl==ncls-1){
600 seedCluster[nseeds]=cl;
601 SetTrackPointFromCluster(cl, seedPoint[nseeds]);
602 ++nseeds;
603 nextSeedRow+=fSeedingDist;
604
605 if (nseeds==3) break;
606 }
607 }
608
609 // check we really have 3 seeds
610 if (nseeds!=3) {
611 AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
612 return 0x0;
613 }
614
615 // do cluster correction for fCorrectionType:
616 // 0 - no correction
617 // 1 - TPC center
618 // 2 - average eta
619 // 3 - ideal
620 // assign the cluster abs time as z component to all seeds
621 for (Int_t iseed=0; iseed<3; ++iseed) {
622 Float_t xyz[3]={0,0,0};
623 seedPoint[iseed].GetXYZ(xyz);
600eaa0d 624 const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
d1cf83f5 625
600eaa0d 626 const Int_t sector=seedCluster[iseed]->GetDetector();
627 const Int_t sign=1-2*((sector/18)%2);
d1cf83f5 628
629 if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
4a777885 630 printf("correction type: %d\n",(Int_t)fCorrectionType);
600eaa0d 631
632 // the settings below are for the T0 seed
633 // for known T0 the z position is already calculated in SetTrackPointFromCluster
634 if ( fCreateT0seed ){
635 if ( fCorrectionType == kTPCCenter ) xyz[2] = 125.*sign;
636 //!!! TODO: is this the correct association?
637 if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
638 }
4a777885 639
d1cf83f5 640 if ( fCorrectionType == kIdeal ) xyz[2] = seedCluster[iseed]->GetZ();
600eaa0d 641
d1cf83f5 642 //!!! TODO: to be replaced with the proper correction
643 fSpaceCharge->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
644 }
645
600eaa0d 646 // after the correction set the time bin as z-Position in case of a T0 seed
647 if ( fCreateT0seed )
648 xyz[2]=seedCluster[iseed]->GetTimeBin();
d1cf83f5 649
650 seedPoint[iseed].SetXYZ(xyz);
651 }
652
653 const Double_t kMaxSnp = 0.85;
654 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
655
656 AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
657 seed->ResetCovariance(10);
658
600eaa0d 659 if (fCreateT0seed){
d1cf83f5 660 // if fTime0 < 0 we assume that we create a seed for the T0 estimate
661 AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
662 if (TMath::Abs(seed->GetX())>3) {
663 printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
664 }
665 }
666
667 return seed;
668
669}
670
671//____________________________________________________________________________________
672void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
673{
674 //
675 // make AliTrackPoint out of AliTPCclusterMI
676 //
677
678 if (!cl) return;
4a777885 679 Float_t xyz[3]={0.,0.,0.};
d1cf83f5 680 // ClusterToSpacePoint(cl,xyz);
681 // cl->GetGlobalCov(cov);
682 //TODO: what to do with the covariance matrix???
683 //TODO: the problem is that it is used in GetAngle in AliTrackPoint
684 //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
685 //TODO: for the moment simply assign 1 permill squared
686 // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
687 // Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
688 // xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
689 // cl->GetGlobalXYZ(xyz);
690 // cl->GetGlobalCov(cov);
691 // voluem ID to add later ....
692 // p.SetXYZ(xyz);
693 // p.SetCov(cov);
694 AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint();
695 p=*tp;
696 delete tp;
697 // cl->Print();
698 // p.Print();
699 p.SetVolumeID(cl->GetDetector());
600eaa0d 700
701
702 if ( !fCreateT0seed && !fIdealTracking ) {
703 p.GetXYZ(xyz);
704 const Int_t sector=cl->GetDetector();
705 const Int_t sign=1-2*((sector/18)%2);
706 const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
707 printf(" z: %.2f %.2f\n",xyz[2],zT0);
708 xyz[2]=zT0;
709 p.SetXYZ(xyz);
4a777885 710 }
600eaa0d 711
712
d1cf83f5 713 // p.Rotate(p.GetAngle()).Print();
714}
715
716//____________________________________________________________________________________
717void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
718{
719 //
720 // convert the cluster to a space point in global coordinates
721 //
722 if (!cl) return;
723 xyz[0]=cl->GetRow();
724 xyz[1]=cl->GetPad();
725 xyz[2]=cl->GetTimeBin(); // this will not be correct at all
726 Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
727 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
728 fTPCParam->Transform8to4(xyz,i);
729 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
730 fTPCParam->Transform4to3(xyz,i);
731 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
732 fTPCParam->Transform2to1(xyz,i);
733 // printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
734}
735
736//____________________________________________________________________________________
737AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
738{
739 //
740 //
741 //
742
743 // create track
744 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
745
746 Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
747
748 const AliTPCROC * roc = AliTPCROC::Instance();
749
750 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
751 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
752 const Double_t kMaxSnp = 0.85;
753 const Double_t kMaxR = 500.;
754 const Double_t kMaxZ = 500.;
755
756 // const Double_t kMaxZ0=220;
757// const Double_t kZcut=3;
758
759 const Double_t refX = tr->GetX();
760
761 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
762
763 // loop over all other points and add to the track
4a777885 764 for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
d1cf83f5 765 AliTrackPoint pIn;
766 const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
767 if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
768 SetTrackPointFromCluster(cl, pIn);
769 if (fCorrectionType != kNoCorrection){
770 Float_t xyz[3]={0,0,0};
771 pIn.GetXYZ(xyz);
600eaa0d 772// if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
d1cf83f5 773 fSpaceCharge->CorrectPoint(xyz, cl->GetDetector());
774 pIn.SetXYZ(xyz);
775 }
776 // rotate the cluster to the local detector frame
777 track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
778 AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
779 if (TMath::Abs(prot.GetX())<kRTPC0) continue;
780 if (TMath::Abs(prot.GetX())>kRTPC1) continue;
781 //
782 Bool_t res=kTRUE;
783 if (fUseMaterial) res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp);
784 else res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
785
786 if (!res) break;
787
788 if (TMath::Abs(track->GetZ())>kMaxZ) break;
789 if (TMath::Abs(track->GetX())>kMaxR) break;
790// if (TMath::Abs(track->GetZ())<kZcut)continue;
791 //
792 Double_t pointPos[2]={0,0};
793 Double_t pointCov[3]={0,0,0};
794 pointPos[0]=prot.GetY();//local y
795 pointPos[1]=prot.GetZ();//local z
796 pointCov[0]=prot.GetCov()[3];//simay^2
797 pointCov[1]=prot.GetCov()[4];//sigmayz
798 pointCov[2]=prot.GetCov()[5];//sigmaz^2
799
800 if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
801 }
802
803 if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp);
804 else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
805
806 // rotate fittet track to the frame of the original track and propagate to same reference
807 track->Rotate(tr->GetAlpha());
808
809 if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp);
810 else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
811
812 return track;
813}
814
510cfcff 815
816//____________________________________________________________________________________
5993ed4f 817AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
510cfcff 818{
819 //
820 // Tracking for given seed on an array of clusters
821 //
822
823 // create track
824 AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
5993ed4f 825
826 const AliTPCROC * roc = AliTPCROC::Instance();
827
828 const Double_t kRTPC0 = roc->GetPadRowRadii(0,0);
829 const Double_t kRTPC1 = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
830 const Double_t kNRowsTPC = roc->GetNRows(0) + roc->GetNRows(36) - 1;
831 const Double_t kIRowsTPC = roc->GetNRows(0) - 1;
832 const Double_t kMaxSnp = 0.85;
833 const Double_t kMaxR = 500.;
834 const Double_t kMaxZ = 500.;
835
836 const Double_t refX = tr->GetX();
837
838 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
839
840 // first propagate seed to outermost row
841 AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
510cfcff 842
5993ed4f 843 // Loop over rows and find the cluster candidates
844 for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
845
846 // inner or outer sector
847 Bool_t bInnerSector = kTRUE;
848 if(iRow > kIRowsTPC) bInnerSector = kFALSE;
849
850 //AliTrackPoint prot = pIn.Rotate(track->GetAlpha()); // rotate to the local frame - non distoted point
851
852 if(bInnerSector){
853 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
854 Double_t x = track->GetX();
855 Int_t sec = (Int_t)((track->Phi() * TMath::RadToDeg() - 10.) / 20. );
856 Printf("inner tracking here: %.2f %d %d from %.2f",x,iRow,sec,track->Phi() * TMath::RadToDeg());
857 //Printf("N1 = %d",fInnerSectorArray[sec%fkNSectorInner][iRow].GetN1());
858 }
859 else{
860 AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
861 Double_t x = track->GetX();
862 Printf("outer tracking here: %.2f %d",x,iRow);
863 }
864
865 Bool_t res=kTRUE;
866 //if (fUseMaterial) res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp);
867 //else res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
510cfcff 868
5993ed4f 869 if (!res) break;
870
871 if (TMath::Abs(track->GetZ())>kMaxZ) break;
872 if (TMath::Abs(track->GetX())>kMaxR) break;
873 //if (TMath::Abs(track->GetZ())<kZcut)continue;
874
875
876 }
877
878 //
879 // Double_t pointPos[2]={0,0};
880 // Double_t pointCov[3]={0,0,0};
881 // pointPos[0]=prot.GetY();//local y
882 // pointPos[1]=prot.GetZ();//local z
883 // pointCov[0]=prot.GetCov()[3];//simay^2
884 // pointCov[1]=prot.GetCov()[4];//sigmayz
885 // pointCov[2]=prot.GetCov()[5];//sigmaz^2
886
887 // if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
888 // }
889
890 if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp);
891 else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
892
893 // rotate fittet track to the frame of the original track and propagate to same reference
894 track->Rotate(tr->GetAlpha());
895
896 if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp);
897 else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
898
510cfcff 899 return track;
900}
901
902
d1cf83f5 903//____________________________________________________________________________________
904void AliToyMCReconstruction::InitSpaceCharge()
905{
906 //
907 // Init the space charge map
908 //
909
910 TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
911 if (fTree) {
912 TList *l=fTree->GetUserInfo();
913 for (Int_t i=0; i<l->GetEntries(); ++i) {
914 TString s(l->At(i)->GetName());
915 if (s.Contains("SC_")) {
916 filename=s;
917 break;
918 }
919 }
920 }
921
922 printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
923 TFile f(filename.Data());
924 fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
925
926 // fSpaceCharge = new AliTPCSpaceCharge3D();
927 // fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
928 // fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
929 // // fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
930 // fSpaceCharge->InitSpaceCharge3DDistortion();
931
932}
933
934//____________________________________________________________________________________
935Double_t AliToyMCReconstruction::GetVDrift() const
936{
937 //
938 //
939 //
940 return fTPCParam->GetDriftV();
941}
942
943//____________________________________________________________________________________
944Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
945{
946 //
947 //
948 //
949 if (roc<0 || roc>71) return -1;
950 return fTPCParam->GetZLength(roc);
951}
952
9e98dea8 953//____________________________________________________________________________________
954TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
955 TString s=gSystem->GetFromPipe(Form("ls %s",files));
956
957 TTree *tFirst=0x0;
958 TObjArray *arrFiles=s.Tokenize("\n");
959
960 for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
961 TString name(arrFiles->At(ifile)->GetName());
962
963 TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
964 TObjArray *arrMatch=0x0;
965 arrMatch=reg.MatchS(name);
966
967 if (!tFirst) {
968 TFile *f=TFile::Open(name.Data());
969 if (!f) continue;
970 TTree *t=(TTree*)f->Get("Tracks");
971 if (!t) {
972 delete f;
973 continue;
974 }
975
976 t->SetName(arrMatch->At(1)->GetName());
977 tFirst=t;
978 } else {
979 tFirst->AddFriend(Form("t%s=Tracks",arrMatch->At(1)->GetName()), name.Data());
980// tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
981 }
982 }
983
376089b6 984 tFirst->GetListOfFriends()->Print();
9e98dea8 985 return tFirst;
986}
d1cf83f5 987
5993ed4f 988//_____________________________________________________________________________
989Int_t AliToyMCReconstruction::LoadOuterSectors() {
990 //-----------------------------------------------------------------
991 // This function fills outer TPC sectors with clusters.
992 // Copy and paste from AliTPCtracker
993 //-----------------------------------------------------------------
994 Int_t nrows = fOuterSectorArray->GetNRows();
995 UInt_t index=0;
996 for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
997 for (Int_t row = 0;row<nrows;row++){
998 AliTPCtrackerRow* tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
999 Int_t sec2 = sec+2*fkNSectorInner;
1000 //left
1001 Int_t ncl = tpcrow->GetN1();
1002 while (ncl--) {
1003 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1004 index=(((sec2<<8)+row)<<16)+ncl;
1005 tpcrow->InsertCluster(c,index);
1006 }
1007 //right
1008 ncl = tpcrow->GetN2();
1009 while (ncl--) {
1010 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1011 index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
1012 tpcrow->InsertCluster(c,index);
1013 }
1014 //
1015 // write indexes for fast acces
1016 //
1017 for (Int_t i=0;i<510;i++)
1018 tpcrow->SetFastCluster(i,-1);
1019 for (Int_t i=0;i<tpcrow->GetN();i++){
1020 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1021 tpcrow->SetFastCluster(zi,i); // write index
1022 }
1023 Int_t last = 0;
1024 for (Int_t i=0;i<510;i++){
1025 if (tpcrow->GetFastCluster(i)<0)
1026 tpcrow->SetFastCluster(i,last);
1027 else
1028 last = tpcrow->GetFastCluster(i);
1029 }
1030 }
1031 return 0;
1032}
1033
1034
1035//_____________________________________________________________________________
1036Int_t AliToyMCReconstruction::LoadInnerSectors() {
1037 //-----------------------------------------------------------------
1038 // This function fills inner TPC sectors with clusters.
1039 // Copy and paste from AliTPCtracker
1040 //-----------------------------------------------------------------
1041 Int_t nrows = fInnerSectorArray->GetNRows();
1042 UInt_t index=0;
1043 for (Int_t sec = 0;sec<fkNSectorInner;sec++)
1044 for (Int_t row = 0;row<nrows;row++){
1045 AliTPCtrackerRow* tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
1046 //
1047 //left
1048 Int_t ncl = tpcrow->GetN1();
1049 while (ncl--) {
1050 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1051 index=(((sec<<8)+row)<<16)+ncl;
1052 tpcrow->InsertCluster(c,index);
1053 }
1054 //right
1055 ncl = tpcrow->GetN2();
1056 while (ncl--) {
1057 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1058 index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
1059 tpcrow->InsertCluster(c,index);
1060 }
1061 //
1062 // write indexes for fast acces
1063 //
1064 for (Int_t i=0;i<510;i++)
1065 tpcrow->SetFastCluster(i,-1);
1066 for (Int_t i=0;i<tpcrow->GetN();i++){
1067 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1068 tpcrow->SetFastCluster(zi,i); // write index
1069 }
1070 Int_t last = 0;
1071 for (Int_t i=0;i<510;i++){
1072 if (tpcrow->GetFastCluster(i)<0)
1073 tpcrow->SetFastCluster(i,last);
1074 else
1075 last = tpcrow->GetFastCluster(i);
1076 }
1077
1078 }
1079 return 0;
1080}