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 ****************************************************************************/
17 #include "AliTPCtracker.h"
20 #include "AliITSgeom.h"
21 #include "AliITSRecPoint.h"
22 #include "AliITSclusterV2.h"
23 #include "AliITSsimulationFastPoints.h"
24 #include "AliITStrackerV2.h"
27 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
28 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n);
29 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
30 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname);
32 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
33 Int_t ITSFindTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
34 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname);
37 Int_t AliBarrelReconstruction(Int_t n=1) {
38 const Char_t *TPCdigName="rfio:galice.root";
39 const Char_t *TPCclsName="AliTPCclusters.root";
40 const Char_t *TPCtrkName="AliTPCtracks.root";
41 const Char_t *TPCtrkNameS="AliTPCtracksSorted.root";
44 const Char_t *ITSdigName="rfio:galice.root";
45 const Char_t *ITSclsName="AliITSclustersV2.root";
46 const Char_t *ITStrkName="AliITStracksV2.root";
48 AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
51 // ********** Find TPC clusters *********** //
52 if (TPCFindClusters(TPCdigName,TPCclsName, n)) {
53 cerr<<"Failed to get TPC clusters !\n";
57 // ********** Find TPC tracks *********** //
58 if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
59 cerr<<"Failed to get TPC tracks !\n";
63 // ********** Sort and label TPC tracks *********** //
64 if (TPCSortTracks(TPCclsName,TPCtrkName,TPCtrkNameS,n)) {
65 cerr<<"Failed to sort TPC tracks !\n";
69 // ********** Find ITS clusters *********** //
70 if (ITSFindClusters(ITSdigName,ITSclsName,n)) {
71 cerr<<"Failed to get ITS clusters !\n";
75 // ********** Find ITS tracks *********** //
77 //TFile *clsFile=TFile::Open(ITSclsName);
78 if (ITSFindTracks(TPCtrkNameS,ITSclsName,ITStrkName,n)) {
79 cerr<<"Failed to get ITS tracks !\n";
85 // ********** Back propagation of the ITS tracks *********** //
86 {TFile *clsFile=TFile::Open(ITSclsName);
87 if (ITSPropagateBack(ITStrkName,TPCtrkNameS)) {
88 cerr<<"Failed to propagate back the ITS tracks !\n";
94 // ********** Back propagation of the TPC tracks *********** //
95 {TFile *clsFile=TFile::Open(TPCclsName);
96 if (TPCPropagateBack(TPCtrkName,TPCtrkName)) {
97 cerr<<"Failed to propagate back the TPC tracks !\n";
106 Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
108 const Char_t *name="TPCFindClusters";
109 cerr<<'\n'<<name<<"...\n";
110 gBenchmark->Start(name);
112 TFile *out=TFile::Open(outname,"recreate");
113 TFile *in =TFile::Open(inname);
115 AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
116 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
120 //tpc.Digits2Clusters(out); //MI change
121 for (Int_t i=0;i<n;i++){
122 printf("Processing event %d\n",i);
123 tpc.Digits2Clusters(out,i);
124 // AliTPCclusterer::Digits2Clusters(dig, out, i);
128 gBenchmark->Stop(name);
129 gBenchmark->Show(name);
134 Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n) {
136 const Char_t *name="TPCFindTracks";
137 cerr<<'\n'<<name<<"...\n";
138 gBenchmark->Start(name);
139 TFile *out=TFile::Open(outname,"recreate");
140 TFile *in =TFile::Open(inname);
141 AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
142 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
144 //AliTPCtracker *tracker=new AliTPCtracker(param);
145 //rc=tracker->Clusters2Tracks(0,out);
148 for (Int_t i=0;i<n;i++){
149 printf("Processing event %d\n",i);
150 AliTPCtracker *tracker = new AliTPCtracker(param,i);
152 tracker->Clusters2Tracks(0,out);
158 gBenchmark->Stop(name);
159 gBenchmark->Show(name);
163 Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t * outname, Int_t eventn){
165 const Char_t *name="TPCSortTracks";
166 cerr<<'\n'<<name<<"...\n";
167 gBenchmark->Start(name);
169 TFile *out =TFile::Open(outname,"recreate");
170 TFile *in1 =TFile::Open(inname);
171 TFile *in2 =TFile::Open(inname2);
174 AliTPCParam *param=(AliTPCParam *)in1->Get("75x40_100x60");
175 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
180 for (Int_t event=0;event<eventn; event++){
181 TObjArray tarray(10000);
182 AliTPCtrack *iotrack=0;
188 sprintf(tname,"TreeT_TPC_%d",event);
190 TTree *tracktree=(TTree*)in2->Get(tname);
191 TBranch *tbranch=tracktree->GetBranch("tracks");
192 Int_t nentr=(Int_t)tracktree->GetEntries();
193 for (i=0; i<nentr; i++) {
194 iotrack=new AliTPCtrack;
195 tbranch->SetAddress(&iotrack);
196 tracktree->GetEvent(i);
197 tarray.AddLast(iotrack);
202 //assign thacks GEANT labels
204 AliTPCtracker *tracker = new AliTPCtracker(param,event);
205 tracker->LoadInnerSectors();
206 tracker->LoadOuterSectors();
207 for (i=0; i<nentr; i++) {
208 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
209 tracker->CookLabel(iotrack,0.1);
213 //end of GEANT label assignment
215 tracktree=new TTree(tname,"Tree with TPC tracks");
216 tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
217 for (i=0; i<nentr; i++) {
218 iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
229 gBenchmark->Stop(name);
230 gBenchmark->Show(name);
235 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
237 const Char_t *name="ITSFindClusters";
238 cerr<<'\n'<<name<<"...\n";
239 gBenchmark->Start(name);
240 TFile *out=TFile::Open(outname,"recreate");
241 TFile *in =TFile::Open(inname,"update");
243 if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
244 cerr<<"Can't get gAlice !\n";
248 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
249 if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
250 AliITSgeom *geom=ITS->GetITSgeom();
255 for (ev = 0; ev<n; ev++){
256 in->cd(); // !!!! MI directory must point to galice. - othervise problem with Tree -connection
257 gAlice->GetEvent(ev);
258 //gAlice->TreeR()->Reset(); //reset reconstructed tree
261 TTree *pTree=gAlice->TreeR();
263 gAlice->MakeTree("R");
264 pTree = gAlice->TreeR();
266 TBranch *branch=pTree->GetBranch("ITSRecPoints");
268 //if not reconstructed ITS branch do reconstruction
269 ITS->MakeBranch("R",0);
270 //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
271 for (Int_t i=0;i<3;i++) {
272 ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
276 Int_t bgr_ev=Int_t(ev/nsignal);
277 ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
278 //////////////////////////////////////////////////////////////////////////
279 gAlice->GetEvent(ev); //MI comment - in HitsToFast... they reset treeR to 0
280 //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
287 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
289 sprintf(cname,"TreeC_ITS_%d",ev);
291 TTree *cTree=new TTree(cname,"ITS clusters");
292 cTree->Branch("Clusters",&clusters);
294 pTree=gAlice->TreeR();
295 if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
296 branch=pTree->GetBranch("ITSRecPoints");
297 if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
298 TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
299 branch->SetAddress(&points);
301 TClonesArray &cl=*clusters;
303 Int_t nentr=(Int_t)pTree->GetEntries();
304 AliITSgeom *geom=ITS->GetITSgeom();
306 for (Int_t i=0; i<nentr; i++) {
307 if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
308 Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
309 Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
310 Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
311 Double_t yshift = x*rot[0] + y*rot[1];
312 Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
313 Int_t ncl=points->GetEntriesFast();
315 for (Int_t j=0; j<ncl; j++) {
316 AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
318 lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
319 lp[1]=p->GetZ()+zshift;
320 lp[2]=p->GetSigmaX2();
321 lp[3]=p->GetSigmaZ2();
324 lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
328 TParticle *part=(TParticle*)gAlice->Particle(label);
330 while (part->P() < 0.005) {
331 Int_t m=part->GetFirstMother();
332 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
334 part=(TParticle*)gAlice->Particle(label);
336 if (lab[1]<0) lab[1]=label;
337 else if (lab[2]<0) lab[2]=label;
338 else cerr<<"No empty labels !\n";
340 new(cl[j]) AliITSclusterV2(lab,lp);
342 cTree->Fill(); clusters->Delete();
346 cerr<<"Number of clusters: "<<nclusters<<endl;
347 delete cTree; delete clusters; delete points;
351 delete gAlice; gAlice=0;
354 gBenchmark->Stop(name);
355 gBenchmark->Show(name);
360 Int_t ITSFindTracks(const Char_t * inname, const Char_t *inname2, const Char_t *outname, Int_t n) {
362 const Char_t *name="ITSFindTracks";
363 cerr<<'\n'<<name<<"...\n";
364 gBenchmark->Start(name);
367 TFile *out=TFile::Open(outname,"recreate");
368 TFile *in =TFile::Open(inname);
369 TFile *in2 =TFile::Open(inname2);
371 AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
372 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
375 for (Int_t i=0;i<n;i++){
376 AliITStrackerV2 tracker(geom,i);
377 rc=tracker.Clusters2Tracks(in,out);
384 gBenchmark->Stop(name);
385 gBenchmark->Show(name);
390 Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname) {
394 const Char_t *name="ITSPropagateBack";
395 cerr<<'\n'<<name<<"...\n";
396 gBenchmark->Start(name);
398 AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
399 if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
400 AliITStrackerV2 tracker(geom);
402 TFile *out=TFile::Open(outname,"update");
403 TFile *in =TFile::Open(inname);
404 rc=tracker.PropagateBack(in,out);
408 gBenchmark->Stop(name);
409 gBenchmark->Show(name);
414 Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname) {
416 const Char_t *name="TPCPropagateBack";
417 cerr<<'\n'<<name<<"...\n";
418 gBenchmark->Start(name);
420 AliTPCParam *param=(AliTPCParam *)gFile->Get("75x40_100x60");
421 if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
422 AliTPCtracker *tracker=new AliTPCtracker(param);
424 TFile *out=TFile::Open(outname,"update");
425 TFile *in =TFile::Open(inname);
426 rc=tracker->PropagateBack(in,out);
431 gBenchmark->Stop(name);
432 gBenchmark->Show(name);