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 ****************************************************************************/
18 #include "AliRunLoader.h"
19 #include "AliLoader.h"
20 #include "AliTPCLoader.h"
21 #include "AliITSLoader.h"
23 #include "AliTPCtracker.h"
25 #include "AliITSgeom.h"
26 #include "AliITSRecPoint.h"
27 #include "AliITSclusterV2.h"
28 #include "AliITSsimulationFastPoints.h"
29 #include "AliITStrackerV2.h"
34 Int_t TPCFindClusters( Int_t n);
35 Int_t TPCFindTracks(Int_t n);
36 Int_t TPCSortTracks(const Char_t *outname, Int_t n);
37 Int_t TPCPropagateBack();
39 Int_t ITSFindClusters(Int_t n);
40 Int_t ITSFindTracks(const Char_t *inname2, Int_t n);
41 Int_t ITSPropagateBack();
43 const char* TPCtrkNameS= "TPC.TracksSorted.root";
48 AliRunLoader *rl = 0x0;
49 AliTPCLoader *tpcl = 0x0;
50 AliTPCParam *param = 0x0;
51 Bool_t debug = kFALSE;
52 Int_t AliBarrelReconstruction(Int_t n=1)
57 delete gAlice->GetRunLoader();
58 delete gAlice;//if everything was OK here it is already NULL
62 rl = AliRunLoader::Open("galice.root");
65 cerr<<"Can not open session"<<endl;
71 cerr<<"Error occured while l"<<endl;
74 AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
76 tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
79 cerr<<"Can not get TPC Loader"<<endl;
83 // ********** Find TPC clusters *********** //
84 if ( TPCFindClusters(n) )
86 cerr<<"Failed to get TPC clusters !\n";
90 // ********** Find TPC tracks *********** //
93 cerr<<"Failed to get TPC tracks !\n";
97 // cout<<"Stopping tracking on TPC\n";
98 // ********** Sort and label TPC tracks *********** //
99 if (TPCSortTracks(TPCtrkNameS,n)) {
100 cerr<<"Failed to sort TPC tracks !\n";
104 // cout<<"Stopping tracking on TPC sorted T\n";
107 // ********** Find ITS clusters *********** //
108 if (ITSFindClusters(n))
110 cerr<<"Failed to get ITS clusters !\n";
114 // ********** Find ITS tracks *********** //
116 if (ITSFindTracks(TPCtrkNameS,n))
118 cerr<<"Failed to get ITS tracks !\n";
122 cout<<"Stopping on ITSFindTracks\n";
128 // ********** Back propagation of the ITS tracks *********** //
129 if ( ITSPropagateBack() ) {
130 cerr<<"Failed to propagate back the ITS tracks !\n";
135 // ********** Back propagation of the TPC tracks *********** //
137 if (TPCPropagateBack()) {
138 cerr<<"Failed to propagate back the TPC tracks !\n";
146 Int_t TPCFindClusters(Int_t n)
149 const Char_t *name="TPCFindClusters";
150 cerr<<'\n'<<name<<"...\n";
152 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
155 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
158 cerr<<"TPC parameters have not been found !\n";
165 gBenchmark->Start(name);
171 tpcl->LoadDigits("read");
172 tpcl->LoadRecPoints("recreate");
175 //tpc.Digits2Clusters(out); //MI change
176 for (Int_t i=0;i<n;i++)
178 printf("Processing event %d\n",i);
179 tpc.Digits2Clusters(i);
180 // AliTPCclusterer::Digits2Clusters(dig, out, i);
183 tpcl->UnloadDigits();
184 tpcl->UnloadRecPoints();
185 gBenchmark->Stop(name);
186 gBenchmark->Show(name);
191 Int_t TPCFindTracks(Int_t n) {
194 AliTPCtracker *tracker = 0x0;
195 const Char_t *name="TPCFindTracks";
198 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
202 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
205 cerr<<"TPC parameters have not been found !\n";
212 gBenchmark->Start(name);
214 tpcl->LoadRecPoints("read");
215 tpcl->LoadTracks("recreate");
217 //AliTPCtracker *tracker=new AliTPCtracker(param);
218 //rc=tracker->Clusters2Tracks(0,out);
221 for (Int_t i=0;i<n;i++){
223 printf("Processing event %d\n",i);
225 tracker = new AliTPCtracker(param, i, (AliConfig::GetDefaultEventFolderName()).Data());
227 tracker->Clusters2Tracks();
230 tpcl->UnloadRecPoints();
231 tpcl->UnloadTracks();
233 gBenchmark->Stop(name);
234 gBenchmark->Show(name);
240 Int_t TPCSortTracks(const Char_t * outname, Int_t eventn){
242 const Char_t *name="TPCSortTracks";
244 cerr<<'\n'<<name<<"...\n";
246 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
249 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
252 cerr<<"TPC parameters have not been found !\n";
258 gBenchmark->Start(name);
261 AliRunLoader* rl2 = AliRunLoader::Open("galice.root","tmp");
263 AliLoader* tpcl2 = (AliTPCLoader*)rl2->GetLoader("TPCLoader");
264 cout<<"tpcl2->SetTracksFileName("<<outname<<");\n";
265 tpcl2->SetTracksFileName(TString(outname));
266 tpcl2->LoadTracks("recreate");
269 for (Int_t event=0;event<eventn; event++)
273 rl2->GetEvent(event);
275 TObjArray tarray(10000);
276 AliTPCtrack *iotrack=0;
279 if (tpcl->TreeT() == 0x0) tpcl->LoadTracks("read");
280 TTree *tracktree=tpcl->TreeT();
281 if (tracktree == 0x0)
283 cerr<<"Can not get TreeT for event "<<event<<endl;
286 TBranch *tbranch=tracktree->GetBranch("tracks");
287 Int_t nentr=(Int_t)tracktree->GetEntries();
288 cout<<"IN Tracks nentr = "<<nentr<<endl;
289 for (i=0; i<nentr; i++)
291 iotrack=new AliTPCtrack;
292 tbranch->SetAddress(&iotrack);
293 tracktree->GetEvent(i);
294 tarray.AddLast(iotrack);
298 //assign thacks GEANT labels
299 cout<<"Running Tracker\n";
300 AliTPCtracker *tracker = new AliTPCtracker(param, event);
302 cout<<"Load Sectors\n";
304 tracker->LoadInnerSectors();
305 tracker->LoadOuterSectors();
307 cout<<"Cooking Labels\n";
309 for (i=0; i<nentr; i++)
311 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
312 tracker->CookLabel(iotrack,0.1);
314 cout<<"deleting tracker\n";
317 //end of GEANT label assignment
320 if (tpcl2->TreeT() == 0x0) tpcl2->MakeTree("T");
321 tracktree= tpcl2->TreeT();
322 if (tracktree == 0x0)
324 cerr<<"Can not get TreeT for Sorted tracks\n";
327 tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
328 for (i=0; i<nentr; i++) {
329 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
333 cout<<" Tracks File Name is "<<tpcl2->GetTracksFileName()<<endl;
334 tpcl2->WriteTracks("OVERWRITE");
340 tpcl->UnloadTracks();
341 tpcl->UnloadRecPoints();
343 gBenchmark->Stop(name);
344 gBenchmark->Show(name);
349 Int_t ITSFindClusters(Int_t n)
352 const Char_t *name="ITSFindClusters";
353 cerr<<'\n'<<name<<"...\n";
354 gBenchmark->Start(name);
362 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
365 cerr<<"Can not get ITS Loader from Run Loader\n";
369 gAlice=rl->GetAliRun();
371 cerr<<"Can't get gAlice !\n";
374 if (rl->TreeK() == 0x0) rl->LoadKinematics();
375 if (rl->TreeE() == 0x0) rl->LoadHeader();
377 itsl->LoadRawClusters("recreate");
379 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
380 if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
381 AliITSgeom *geom=ITS->GetITSgeom();
384 for (ev = 0; ev<n; ev++)
388 TBranch *branch = 0x0;
394 itsl->LoadRecPoints("read");
398 ::Error("AliBarrelReonstruction.C::ITSFindClusters",
399 "Can not get TreeR for event %d",ev);
405 branch=(pTree)?pTree->GetBranch("ITSRecPoints"):0x0;
408 //if not reconstructed ITS branch do reconstruction
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");
414 ITS->MakeBranch("RF");
415 ITS->SetTreeAddress();
418 //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
419 for (Int_t i=0;i<3;i++) {
420 ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
424 Int_t bgr_ev=Int_t(ev/nsignal);
426 ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
427 //////////////////////////////////////////////////////////////////////////
428 //MI comment - in HitsToFast... they reset treeR to 0
429 //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
432 ITS->SetTreeAddress();
435 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
437 if (itsl->TreeC() == 0x0) itsl->MakeTree("C");
438 TTree *cTree = itsl->TreeC();
439 TBranch* b = cTree->GetBranch("Clusters");
441 cTree->Branch("Clusters",&clusters);
442 else b->SetAddress(&clusters);
445 if (itsl->TreeR() == 0x0) itsl->LoadRecPoints();
448 {//it means that we produced FastRecPoints above
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;}
454 TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
455 branch->SetAddress(&points);
457 TClonesArray &cl=*clusters;
459 Int_t nentr=(Int_t)branch->GetEntries();
460 AliITSgeom *geom=ITS->GetITSgeom();
462 cout<<"\n\n Number of entries in TreeR = "<<nentr<<endl;
463 for (Int_t i=0; i<nentr; i++)
465 if (!branch->GetEvent(i))
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();
479 if (debug) ::Info("AliBarrelReconstruction.C",
480 "i=%d lay=%d lad=%d det=%d NRP=%d",
481 i, lay, lad, det, ncl);
483 for (Int_t j=0; j<ncl; j++)
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();
494 lab[0]=p->GetLabel(0);
495 lab[1]=p->GetLabel(1);
496 lab[2]=p->GetLabel(2);
500 if (label<0) continue;
501 TParticle *part=(TParticle*)rl->Stack()->Particle(label);
503 while (part->P() < 0.005) {
504 Int_t m=part->GetFirstMother();
505 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
507 part=(TParticle*)gAlice->Particle(label);
509 if (lab[1]<0) lab[1]=label;
510 else if (lab[2]<0) lab[2]=label;
511 else cerr<<"No empty labels !\n";
513 new(cl[j]) AliITSclusterV2(lab,lp);
515 cTree->Fill(); clusters->Delete();
518 itsl->WriteRawClusters("OVERWRITE");
520 cerr<<"Number of clusters: "<<nclusters<<endl;
521 delete clusters; delete points;
524 itsl->UnloadRecPoints();
525 itsl->UnloadRawClusters();
526 rl->UnloadKinematics();
528 gBenchmark->Stop(name);
529 gBenchmark->Show(name);
534 Int_t ITSFindTracks(const Char_t *inname2, Int_t n)
538 const Char_t *name="ITSFindTracks";
539 cerr<<'\n'<<name<<"...\n";
540 gBenchmark->Start(name);
544 AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
545 if ( (itsl == 0x0) || (tpcl == 0x0))
547 cerr<<"Can not get ITS Loader from Run Loader\n";
551 tpcl->SetTracksFileName(TString(inname2));
552 if(rl->GetAliRun() == 0x0) rl->LoadgAlice();
554 AliITS* ITS = rl->GetAliRun()->GetModule("ITS");
557 AliITSgeom *geom=ITS->GetITSgeom();
558 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
560 itsl->LoadTracks("recreate");
561 itsl->LoadRawClusters("read");
562 for (Int_t i=0;i<n;i++)
565 AliITStrackerV2 tracker(geom,i,AliConfig::GetDefaultEventFolderName());
566 rc=tracker.Clusters2Tracks();
570 gBenchmark->Stop(name);
571 gBenchmark->Show(name);
576 Int_t ITSPropagateBack()
581 const Char_t *name="ITSPropagateBack";
582 cerr<<'\n'<<name<<"...\n";
583 gBenchmark->Start(name);
585 AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
586 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
587 AliITStrackerV2 tracker(geom);
589 TFile *out=TFile::Open(outname,"update");
590 TFile *in =TFile::Open(inname);
591 rc=tracker.PropagateBack(in,out);
595 gBenchmark->Stop(name);
596 gBenchmark->Show(name);
601 Int_t TPCPropagateBack()
604 const Char_t *name="TPCPropagateBack";
605 cerr<<'\n'<<name<<"...\n";
606 gBenchmark->Start(name);
608 param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
611 param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
614 cerr<<"TPC parameters have not been found !\n";
621 AliTPCtracker *tracker = new AliTPCtracker(param, AliConfig::GetDefaultEventFolderName());
623 // TFile *out=TFile::Open(outname,"update");
624 // TFile *in =TFile::Open(inname);
626 rc=tracker->PropagateBack();
629 gBenchmark->Stop(name);
630 gBenchmark->Show(name);