Using getter instead of global constant
[u/mrichter/AliRoot.git] / TPC / AliBarrelReconstruction.C
CommitLineData
b9de75e1 1/****************************************************************************
2 * This macro is supposed to do reconstruction in the barrel ALICE trackers *
3 * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
4 * It does the following steps (April 12, 2001): *
5 * 1) TPC cluster finding *
6 * 2) TPC track finding *
7 * 3) ITS cluster finding V2 (via fast points !) *
8 * 4) ITS track finding V2 *
9 * 5) ITS back track propagation V2 *
10 * 6) TPC back track propagation *
11 * (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch) *
12 ****************************************************************************/
13
14#ifndef __CINT__
15 #include "alles.h"
88cb7938 16 #include "iostream.h"
17 #include "AliRun.h"
18 #include "AliRunLoader.h"
19 #include "AliLoader.h"
20 #include "AliTPCLoader.h"
21 #include "AliITSLoader.h"
b9de75e1 22 #include "AliMagF.h"
23 #include "AliTPCtracker.h"
b9de75e1 24 #include "AliITS.h"
25 #include "AliITSgeom.h"
26 #include "AliITSRecPoint.h"
27 #include "AliITSclusterV2.h"
28 #include "AliITSsimulationFastPoints.h"
29 #include "AliITStrackerV2.h"
30#endif
31
b9de75e1 32
b9de75e1 33
88cb7938 34Int_t TPCFindClusters( Int_t n);
35Int_t TPCFindTracks(Int_t n);
36Int_t TPCSortTracks(const Char_t *outname, Int_t n);
37Int_t TPCPropagateBack();
38
39Int_t ITSFindClusters(Int_t n);
40Int_t ITSFindTracks(const Char_t *inname2, Int_t n);
41Int_t ITSPropagateBack();
42
43const char* TPCtrkNameS= "TPC.TracksSorted.root";
44
45class AliRunLoader;
46class AliTPCLoader;
47class AliTPCParam;
48AliRunLoader *rl = 0x0;
49AliTPCLoader *tpcl = 0x0;
50AliTPCParam *param = 0x0;
51Bool_t debug = kFALSE;
52Int_t AliBarrelReconstruction(Int_t n=1)
53 {
54
55 if (gAlice)
56 {
57 delete gAlice->GetRunLoader();
58 delete gAlice;//if everything was OK here it is already NULL
59 gAlice = 0x0;
60 }
743a19f2 61
88cb7938 62 rl = AliRunLoader::Open("galice.root");
63 if (rl == 0x0)
64 {
65 cerr<<"Can not open session"<<endl;
66 return 1;
67 }
68
69 if (rl->LoadgAlice())
70 {
71 cerr<<"Error occured while l"<<endl;
72 return 1;
73 }
74 AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
75
76 tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
77 if (tpcl == 0x0)
78 {
79 cerr<<"Can not get TPC Loader"<<endl;
80 return 1;
81 }
82
b9de75e1 83// ********** Find TPC clusters *********** //
88cb7938 84 if ( TPCFindClusters(n) )
85 {
b9de75e1 86 cerr<<"Failed to get TPC clusters !\n";
87 return 1;
88cb7938 88 }
89
b9de75e1 90// ********** Find TPC tracks *********** //
88cb7938 91 if (TPCFindTracks(n))
92 {
b9de75e1 93 cerr<<"Failed to get TPC tracks !\n";
94 return 1;
88cb7938 95 }
96
97// cout<<"Stopping tracking on TPC\n";
b9de75e1 98// ********** Sort and label TPC tracks *********** //
88cb7938 99 if (TPCSortTracks(TPCtrkNameS,n)) {
b9de75e1 100 cerr<<"Failed to sort TPC tracks !\n";
101 return 1;
88cb7938 102 }
103
104// cout<<"Stopping tracking on TPC sorted T\n";
105
743a19f2 106
b9de75e1 107// ********** Find ITS clusters *********** //
88cb7938 108 if (ITSFindClusters(n))
109 {
b9de75e1 110 cerr<<"Failed to get ITS clusters !\n";
111 return 1;
88cb7938 112 }
743a19f2 113
b9de75e1 114// ********** Find ITS tracks *********** //
88cb7938 115
116 if (ITSFindTracks(TPCtrkNameS,n))
117 {
b9de75e1 118 cerr<<"Failed to get ITS tracks !\n";
119 return 1;
88cb7938 120 }
743a19f2 121 //clsFile->Close();
88cb7938 122 cout<<"Stopping on ITSFindTracks\n";
123 delete rl;
124 rl = 0x0;
125 return 0;
743a19f2 126 return 1;
88cb7938 127
b9de75e1 128// ********** Back propagation of the ITS tracks *********** //
88cb7938 129 if ( ITSPropagateBack() ) {
b9de75e1 130 cerr<<"Failed to propagate back the ITS tracks !\n";
131 return 1;
132 }
b9de75e1 133
134
135// ********** Back propagation of the TPC tracks *********** //
88cb7938 136
137 if (TPCPropagateBack()) {
b9de75e1 138 cerr<<"Failed to propagate back the TPC tracks !\n";
139 return 1;
140 }
b9de75e1 141
142 return 0;
143}
144
145
88cb7938 146Int_t TPCFindClusters(Int_t n)
147 {
b9de75e1 148 Int_t rc=0;
149 const Char_t *name="TPCFindClusters";
150 cerr<<'\n'<<name<<"...\n";
88cb7938 151 rl->CdGAFile();
152 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
153 if (!param)
154 {
155 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
156 if (!param)
157 {
158 cerr<<"TPC parameters have not been found !\n";
159 return 1;
160 }
161 }
162// param->Dump();
163 param->Update();
b9de75e1 164
88cb7938 165 gBenchmark->Start(name);
b9de75e1 166
b9de75e1 167 AliTPCv2 tpc;
168 tpc.SetParam(param);
88cb7938 169 tpc.SetLoader(tpcl);
170
171 tpcl->LoadDigits("read");
172 tpcl->LoadRecPoints("recreate");
173
b9de75e1 174
743a19f2 175 //tpc.Digits2Clusters(out); //MI change
88cb7938 176 for (Int_t i=0;i<n;i++)
177 {
178 printf("Processing event %d\n",i);
179 tpc.Digits2Clusters(i);
180 // AliTPCclusterer::Digits2Clusters(dig, out, i);
181 }
182
183 tpcl->UnloadDigits();
184 tpcl->UnloadRecPoints();
b9de75e1 185 gBenchmark->Stop(name);
186 gBenchmark->Show(name);
187
188 return rc;
189}
190
88cb7938 191Int_t TPCFindTracks(Int_t n) {
192
b9de75e1 193 Int_t rc=0;
88cb7938 194 AliTPCtracker *tracker = 0x0;
b9de75e1 195 const Char_t *name="TPCFindTracks";
88cb7938 196 rl->GetEvent(0);
197 rl->CdGAFile();
198 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
199
200 if (!param)
201 {
202 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
203 if (!param)
204 {
205 cerr<<"TPC parameters have not been found !\n";
206 return 1;
207 }
208 }
209// param->Dump();
210 param->Update();
211
b9de75e1 212 gBenchmark->Start(name);
88cb7938 213
214 tpcl->LoadRecPoints("read");
215 tpcl->LoadTracks("recreate");
b9de75e1 216
743a19f2 217 //AliTPCtracker *tracker=new AliTPCtracker(param);
218 //rc=tracker->Clusters2Tracks(0,out);
219 //delete tracker;
88cb7938 220
743a19f2 221 for (Int_t i=0;i<n;i++){
88cb7938 222 rl->GetEvent(i);
743a19f2 223 printf("Processing event %d\n",i);
88cb7938 224
f5a857b2 225 tracker = new AliTPCtracker(param, i, (AliConfig::GetDefaultEventFolderName()).Data());
743a19f2 226 //Int_t rc=
88cb7938 227 tracker->Clusters2Tracks();
743a19f2 228 delete tracker;
229 }
88cb7938 230 tpcl->UnloadRecPoints();
231 tpcl->UnloadTracks();
743a19f2 232
b9de75e1 233 gBenchmark->Stop(name);
234 gBenchmark->Show(name);
235
236 return rc;
237}
743a19f2 238
239
88cb7938 240Int_t TPCSortTracks(const Char_t * outname, Int_t eventn){
241 Int_t rc=0;
242 const Char_t *name="TPCSortTracks";
b9de75e1 243
88cb7938 244 cerr<<'\n'<<name<<"...\n";
245 rl->CdGAFile();
246 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
247 if (!param)
248 {
249 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
250 if (!param)
251 {
252 cerr<<"TPC parameters have not been found !\n";
253 return 1;
254 }
255 }
256 param->Update();
743a19f2 257
88cb7938 258 gBenchmark->Start(name);
259
260
261 AliRunLoader* rl2 = AliRunLoader::Open("galice.root","tmp");
262
263 AliLoader* tpcl2 = (AliTPCLoader*)rl2->GetLoader("TPCLoader");
264 cout<<"tpcl2->SetTracksFileName("<<outname<<");\n";
265 tpcl2->SetTracksFileName(TString(outname));
266 tpcl2->LoadTracks("recreate");
267
743a19f2 268 // loop over events
88cb7938 269 for (Int_t event=0;event<eventn; event++)
270 {
271
272 rl->GetEvent(event);
273 rl2->GetEvent(event);
274
743a19f2 275 TObjArray tarray(10000);
276 AliTPCtrack *iotrack=0;
277 Int_t i;
88cb7938 278
279 if (tpcl->TreeT() == 0x0) tpcl->LoadTracks("read");
280 TTree *tracktree=tpcl->TreeT();
281 if (tracktree == 0x0)
282 {
283 cerr<<"Can not get TreeT for event "<<event<<endl;
284 continue;
285 }
743a19f2 286 TBranch *tbranch=tracktree->GetBranch("tracks");
287 Int_t nentr=(Int_t)tracktree->GetEntries();
88cb7938 288 cout<<"IN Tracks nentr = "<<nentr<<endl;
289 for (i=0; i<nentr; i++)
290 {
b9de75e1 291 iotrack=new AliTPCtrack;
292 tbranch->SetAddress(&iotrack);
293 tracktree->GetEvent(i);
294 tarray.AddLast(iotrack);
88cb7938 295 }
743a19f2 296 tarray.Sort();
743a19f2 297
298 //assign thacks GEANT labels
88cb7938 299 cout<<"Running Tracker\n";
300 AliTPCtracker *tracker = new AliTPCtracker(param, event);
301
302 cout<<"Load Sectors\n";
303
743a19f2 304 tracker->LoadInnerSectors();
305 tracker->LoadOuterSectors();
88cb7938 306
307 cout<<"Cooking Labels\n";
308
309 for (i=0; i<nentr; i++)
310 {
b9de75e1 311 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
312 tracker->CookLabel(iotrack,0.1);
88cb7938 313 }
314 cout<<"deleting tracker\n";
743a19f2 315 delete tracker;
316 //in->Close();
317 //end of GEANT label assignment
318
88cb7938 319
320 if (tpcl2->TreeT() == 0x0) tpcl2->MakeTree("T");
321 tracktree= tpcl2->TreeT();
322 if (tracktree == 0x0)
323 {
324 cerr<<"Can not get TreeT for Sorted tracks\n";
325 return 1;
326 }
743a19f2 327 tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
328 for (i=0; i<nentr; i++) {
b9de75e1 329 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
330 tracktree->Fill();
743a19f2 331 }
88cb7938 332
333 cout<<" Tracks File Name is "<<tpcl2->GetTracksFileName()<<endl;
334 tpcl2->WriteTracks("OVERWRITE");
b9de75e1 335 }
88cb7938 336
743a19f2 337
88cb7938 338 delete rl2;
b9de75e1 339
88cb7938 340 tpcl->UnloadTracks();
341 tpcl->UnloadRecPoints();
342
b9de75e1 343 gBenchmark->Stop(name);
344 gBenchmark->Show(name);
345
346 return rc;
347}
348
88cb7938 349Int_t ITSFindClusters(Int_t n)
350 {
b9de75e1 351 Int_t rc=0;
352 const Char_t *name="ITSFindClusters";
353 cerr<<'\n'<<name<<"...\n";
354 gBenchmark->Start(name);
b9de75e1 355
88cb7938 356 if (rl->GetEvent(0))
357 {
358 cerr<<"Problems\n";
359 return 1;
360 }
361
362 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
363 if (itsl == 0x0)
364 {
365 cerr<<"Can not get ITS Loader from Run Loader\n";
366 return 1;
367 }
368
369 gAlice=rl->GetAliRun();
370 if (!gAlice) {
b9de75e1 371 cerr<<"Can't get gAlice !\n";
372 return 1;
373 }
88cb7938 374 if (rl->TreeK() == 0x0) rl->LoadKinematics();
375 if (rl->TreeE() == 0x0) rl->LoadHeader();
376
377 itsl->LoadRawClusters("recreate");
378
b9de75e1 379 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
380 if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
b9de75e1 381 AliITSgeom *geom=ITS->GetITSgeom();
88cb7938 382
743a19f2 383 Int_t ev=0;
88cb7938 384 for (ev = 0; ev<n; ev++)
385 {
743a19f2 386
88cb7938 387 rl->GetEvent(ev);
388 TBranch *branch = 0x0;
389 TTree *pTree = 0x0;
743a19f2 390
88cb7938 391 pTree=itsl->TreeR();
392 if (pTree == 0x0)
393 {
394 itsl->LoadRecPoints("read");
395 pTree=itsl->TreeR();
396 if (pTree == 0x0)
397 {
398 ::Error("AliBarrelReonstruction.C::ITSFindClusters",
399 "Can not get TreeR for event %d",ev);
400 return 1;
401 }
402 }
403
404 pTree=itsl->TreeR();
405 branch=(pTree)?pTree->GetBranch("ITSRecPoints"):0x0;
406
407 if (branch== 0x0) {
743a19f2 408 //if not reconstructed ITS branch do reconstruction
88cb7938 409 if (debug) ::Info("AliBarrelReconstruction.C","Did not get ITSRecPoints from TreeR.");
410 if (debug) ::Info("AliBarrelReconstruction.C","Making branch and Producing RecPoints");
411 itsl->SetRecPointsFileOption("recreate");
412 if (itsl->TreeR()==0x0) itsl->MakeTree("R");
413
414 ITS->MakeBranch("RF");
415 ITS->SetTreeAddress();
416 itsl->LoadHits();
417
743a19f2 418 //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
8d9c76b5 419 for (Int_t i=0;i<3;i++) {
420 ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
421 }
743a19f2 422 Int_t nsignal=25;
423 Int_t size=-1;
424 Int_t bgr_ev=Int_t(ev/nsignal);
88cb7938 425
743a19f2 426 ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
427 //////////////////////////////////////////////////////////////////////////
88cb7938 428 //MI comment - in HitsToFast... they reset treeR to 0
743a19f2 429 //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
430 //time
88cb7938 431 itsl->UnloadHits();
432 ITS->SetTreeAddress();
743a19f2 433 }
434
743a19f2 435 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
88cb7938 436
437 if (itsl->TreeC() == 0x0) itsl->MakeTree("C");
438 TTree *cTree = itsl->TreeC();
439 TBranch* b = cTree->GetBranch("Clusters");
440 if (b == 0x0)
441 cTree->Branch("Clusters",&clusters);
442 else b->SetAddress(&clusters);
743a19f2 443
88cb7938 444
445 if (itsl->TreeR() == 0x0) itsl->LoadRecPoints();
446
447 if (branch == 0x0)
448 {//it means that we produced FastRecPoints above
449 pTree=itsl->TreeR();
450 if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
451 branch=pTree->GetBranch("ITSRecPointsF");
452 if (!branch) { cerr<<"Can't get Fast RecPoints branch named ITSRecPointsF !\n"; return 1;}
453 }
743a19f2 454 TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
455 branch->SetAddress(&points);
456
457 TClonesArray &cl=*clusters;
458 Int_t nclusters=0;
88cb7938 459 Int_t nentr=(Int_t)branch->GetEntries();
743a19f2 460 AliITSgeom *geom=ITS->GetITSgeom();
88cb7938 461
462 cout<<"\n\n Number of entries in TreeR = "<<nentr<<endl;
463 for (Int_t i=0; i<nentr; i++)
464 {
465 if (!branch->GetEvent(i))
466 {
467 cTree->Fill();
468 continue;
469 }
470
b9de75e1 471 Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
472 Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
473 Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
474 Double_t yshift = x*rot[0] + y*rot[1];
475 Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
476 Int_t ncl=points->GetEntriesFast();
477 nclusters+=ncl;
88cb7938 478
479 if (debug) ::Info("AliBarrelReconstruction.C",
480 "i=%d lay=%d lad=%d det=%d NRP=%d",
481 i, lay, lad, det, ncl);
482
483 for (Int_t j=0; j<ncl; j++)
484 {
485 Float_t lp[5];
486 Int_t lab[4];
487
488 AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
489 lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
490 lp[1]=p->GetZ()+zshift;
491 lp[2]=p->GetSigmaX2();
492 lp[3]=p->GetSigmaZ2();
493 lp[4]=p->GetQ();
494 lab[0]=p->GetLabel(0);
495 lab[1]=p->GetLabel(1);
496 lab[2]=p->GetLabel(2);
497 lab[3]=ndet;
498
499 Int_t label=lab[0];
500 if (label<0) continue;
501 TParticle *part=(TParticle*)rl->Stack()->Particle(label);
502 label=-3;
503 while (part->P() < 0.005) {
504 Int_t m=part->GetFirstMother();
505 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
506 label=m;
507 part=(TParticle*)gAlice->Particle(label);
508 }
509 if (lab[1]<0) lab[1]=label;
510 else if (lab[2]<0) lab[2]=label;
511 else cerr<<"No empty labels !\n";
512
513 new(cl[j]) AliITSclusterV2(lab,lp);
514 }
b9de75e1 515 cTree->Fill(); clusters->Delete();
516 points->Delete();
743a19f2 517 }
88cb7938 518 itsl->WriteRawClusters("OVERWRITE");
519
743a19f2 520 cerr<<"Number of clusters: "<<nclusters<<endl;
88cb7938 521 delete clusters; delete points;
b9de75e1 522
743a19f2 523 }
88cb7938 524 itsl->UnloadRecPoints();
525 itsl->UnloadRawClusters();
526 rl->UnloadKinematics();
b9de75e1 527
b9de75e1 528 gBenchmark->Stop(name);
529 gBenchmark->Show(name);
530
531 return rc;
532}
533
88cb7938 534Int_t ITSFindTracks(const Char_t *inname2, Int_t n)
535{
536
b9de75e1 537 Int_t rc=0;
538 const Char_t *name="ITSFindTracks";
539 cerr<<'\n'<<name<<"...\n";
540 gBenchmark->Start(name);
541
88cb7938 542 rl->GetEvent(0);
b9de75e1 543
88cb7938 544 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
545 if ( (itsl == 0x0) || (tpcl == 0x0))
546 {
547 cerr<<"Can not get ITS Loader from Run Loader\n";
548 return 1;
549 }
550
551 tpcl->SetTracksFileName(TString(inname2));
552 if(rl->GetAliRun() == 0x0) rl->LoadgAlice();
553
554 AliITS* ITS = rl->GetAliRun()->GetModule("ITS");
555
556 rl->CdGAFile();
557 AliITSgeom *geom=ITS->GetITSgeom();
743a19f2 558 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
559
88cb7938 560 itsl->LoadTracks("recreate");
561 itsl->LoadRawClusters("read");
562 for (Int_t i=0;i<n;i++)
563 {
564 rl->GetEvent(i);
f5a857b2 565 AliITStrackerV2 tracker(geom,i,AliConfig::GetDefaultEventFolderName());
88cb7938 566 rc=tracker.Clusters2Tracks();
567 }
743a19f2 568
b9de75e1 569
570 gBenchmark->Stop(name);
571 gBenchmark->Show(name);
572
573 return rc;
574}
575
88cb7938 576Int_t ITSPropagateBack()
577{
743a19f2 578
579 Int_t rc=0;
580 /*
b9de75e1 581 const Char_t *name="ITSPropagateBack";
582 cerr<<'\n'<<name<<"...\n";
583 gBenchmark->Start(name);
584
585 AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
586 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
587 AliITStrackerV2 tracker(geom);
588
589 TFile *out=TFile::Open(outname,"update");
590 TFile *in =TFile::Open(inname);
591 rc=tracker.PropagateBack(in,out);
592 in->Close();
593 out->Close();
594
595 gBenchmark->Stop(name);
596 gBenchmark->Show(name);
743a19f2 597 */
b9de75e1 598 return rc;
599}
600
88cb7938 601Int_t TPCPropagateBack()
602{
b9de75e1 603 Int_t rc=0;
604 const Char_t *name="TPCPropagateBack";
605 cerr<<'\n'<<name<<"...\n";
606 gBenchmark->Start(name);
607
88cb7938 608 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
609 if (!param)
610 {
611 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
612 if (!param)
613 {
614 cerr<<"TPC parameters have not been found !\n";
615 return 1;
616 }
617 }
618// param->Dump();
619 param->Update();
b9de75e1 620
f5a857b2 621 AliTPCtracker *tracker = new AliTPCtracker(param, AliConfig::GetDefaultEventFolderName());
88cb7938 622
623// TFile *out=TFile::Open(outname,"update");
624// TFile *in =TFile::Open(inname);
625
626 rc=tracker->PropagateBack();
b9de75e1 627 delete tracker;
b9de75e1 628
629 gBenchmark->Stop(name);
630 gBenchmark->Show(name);
631
632 return rc;
633}