Update to track display by Chuncheng
[u/mrichter/AliRoot.git] / TRD / AliTRDdisplayTracks.C
1 void AliTRDdisplayTracks(Int_t track=157) 
2 {
3
4   c1 = new TCanvas( "RecPoints", "RecPoints Display", 10, 10, 710, 740);
5   c1->SetFillColor(1);
6   TView *v=new TView(1);
7   
8       v->SetRange(-350.,-350.,-400.,350.,350.,400.); // full
9   //  v->SetRange(0.,0.,0.,350.,350.,400.);  // top right
10   //  v->SetRange(-350.,0.,0.,0.,350.,400.); // top left
11   //  v->SetRange(0.,-350.,0.,350.,0.,400.); // bottom right 
12   //  v->SetRange(-350.,-350.,0.,0.,0.,400.); // bottom left
13   //  v->Side();
14   v->Top();
15   cerr<<"psi = "<<v->GetPsi()<<endl;
16
17   // Dynamically link some shared libs
18   if (gClassTable->GetID("AliRun") < 0) {
19     gROOT->LoadMacro("loadlibs.C");
20     loadlibs();
21     cout << "Loaded shared libraries" << endl;
22   }       
23
24
25 // Load tracks
26
27    TFile *tf=TFile::Open("AliTRDtracks.root");
28    if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return 3;}
29    TObjArray tarray(2000);
30    //  tf->ls();
31    TTree *tracktree=(TTree*)tf->Get("TreeT0_TRD");
32    //   tracktree->ls();
33    TBranch *tbranch=tracktree->GetBranch("tracks");
34    Int_t nentr=tracktree->GetEntries();
35    cerr<<"Found "<<nentr<<" entries in the track tree"<<endl;
36    for (Int_t i=0; i<nentr; i++) {
37        AliTRDtrack *iotrack=new AliTRDtrack;
38        tbranch->SetAddress(&iotrack);
39        tracktree->GetEvent(i);
40        tarray.AddLast(iotrack);
41        cerr<<"got track "<<i<<": index ="<<iotrack->GetLabel()<<endl;
42    }
43    tf->Close();                 
44
45
46   // Load clusters
47   Char_t *alifile = "TRDclusters.root";
48   Int_t   nEvent  = 0;
49   TObjArray rparray(2000);
50   TObjArray *RecPointsArray = &rparray;
51   const TFile *geofile =TFile::Open(alifile);
52   //  AliTRDtracker *Tracker = new AliTRDtracker("dummy","dummy");
53   AliTRDtracker *Tracker = new AliTRDtracker(geofile);
54    Tracker->ReadClusters(RecPointsArray,alifile); 
55   Int_t nRecPoints = RecPointsArray->GetEntriesFast();
56   cerr<<"Found "<<nRecPoints<<" rec. points"<<endl;
57
58
59   // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
60       alifile = "galice.root";
61
62    TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
63   if (!gafl) {
64     cout << "Open the ALIROOT-file " << alifile << endl;
65     gafl = new TFile(alifile);
66   }
67   else {
68     cout << alifile << " is already open" << endl;
69   }
70   
71
72   // Get AliRun object from file or create it if not on file
73   gAlice = (AliRun*) gafl->Get("gAlice");
74   if (gAlice)
75     cout << "AliRun object found on file" << endl;
76   else
77     gAlice = new AliRun("gAlice","Alice test program");
78         
79
80        AliTRDparameter *par =  ( AliTRDparameter *par) gafl->Get("TRDparameter");
81        AliTRDv1       *TRD           = (AliTRDv1*) gAlice->GetDetector("TRD");
82        AliTRDgeometry *fGeom   = TRD->GetGeometry();   
83        
84   Int_t i,j,index,det,sector, ti[3];
85   Double_t x,y,z, cs,sn,tmp;
86   Float_t global[3], local[3];
87
88   // Display all TRD RecPoints
89   TPolyMarker3D *pm = new TPolyMarker3D(nRecPoints);
90
91   for (Int_t i = 0; i < nRecPoints; i++) {
92     printf("\r point %d out of %d",i,nRecPoints);
93     AliTRDcluster *rp = (AliTRDcluster *) RecPointsArray->UncheckedAt(i);    
94   
95    Int_t  idet=rp->GetDetector();
96    Int_t iplan = fGeom->GetPlane(idet);    
97    Int_t itt=rp->GetLocalTimeBin();
98    Float_t  timeSlice = itt+0.5; 
99    Float_t  time0     = par->GetTime0(iplan);
100
101   // calculate (x,y,z) position in rotated chamber
102    local[0] = time0 - (timeSlice - par->GetTimeBefore()) 
103          * par->GetTimeBinSize();
104    local[1]=rp->GetY();
105    local[2]=rp->GetZ(); 
106
107     for (j = 0; j < 3; j++) { ti[j] = rp->GetTrackIndex(j); }
108
109     if((track < 0) ||
110        ((ti[0]==track)||(ti[1]==track)||(ti[2]==track))) {
111
112         if(fGeom->RotateBack(idet,local,global)) {     
113         x=global[0]; y=global[1]; z=global[2];
114         pm->SetPoint(i,x,y,z);
115          }
116
117     }
118   }
119
120   pm->SetMarkerSize(1); pm->SetMarkerColor(10); pm->SetMarkerStyle(1);
121   pm->Draw();
122   
123    AliTRDparameter *par =  ( AliTRDparameter *par) gafl->Get("TRDparameter");
124
125   Int_t ntracks = tarray.GetEntriesFast();
126
127   for (i = 0; i < ntracks; i++) {
128     AliTRDtrack *pt = (AliTRDtrack *) tarray.UncheckedAt(i), &t=*pt;
129     Int_t nclusters = t.GetNumberOfClusters();
130     cerr<<"in track "<<i<<" found "<<nclusters<<" clusters"<<endl;
131    
132     TPolyMarker3D *pm = new TPolyMarker3D(nclusters);
133     for(j = 0; j < nclusters; j++) {
134       index = t.GetClusterIndex(j);
135       AliTRDcluster *rp = (AliTRDcluster *) RecPointsArray->UncheckedAt(index); 
136
137    Int_t  idet=rp->GetDetector();
138    Int_t iplan = fGeom->GetPlane(idet);    
139    Int_t itt=rp->GetLocalTimeBin();
140    Float_t  timeSlice = itt+0.5; 
141    Float_t  time0     = par->GetTime0(iplan);
142
143
144   // calculate (x,y,z) position in rotated chamber
145        local[0] = time0 - (timeSlice - par->GetTimeBefore()) 
146          * par->GetTimeBinSize();
147        local[1]=rp->GetY(); 
148        local[2]=rp->GetZ();
149
150             if(fGeom->RotateBack(idet,local,global)) {     
151         x=global[0]; y=global[1]; z=global[2];
152         pm->SetPoint(j,x,y,z);
153        }
154     } 
155     pm->SetMarkerSize(1); pm->SetMarkerColor(i%6+3); pm->SetMarkerStyle(1);
156     //    pm->Draw();
157   }
158   
159   TGeometry *geom=(TGeometry*)gafl->Get("AliceGeom");
160   geom->Draw("same");
161
162   c1->Modified(); 
163   
164   c1->Update();
165   
166   gafl->Close();
167
168 }