]>
Commit | Line | Data |
---|---|---|
1 | #ifndef __CINT__ | |
2 | #include <iostream.h> | |
3 | ||
4 | #include "AliTRDtracker.h" | |
5 | #include "AliTRDcluster.h" | |
6 | #include "AliTRDhit.h" | |
7 | #include "AliTRDv1.h" | |
8 | #include "AliTRDgeometry.h" | |
9 | #include "AliTRDparameter.h" | |
10 | ||
11 | #include "alles.h" | |
12 | #include "TFile.h" | |
13 | #include "TStopwatch.h" | |
14 | ||
15 | #endif | |
16 | ||
17 | void AliTRDclusterErrors() { | |
18 | ||
19 | Int_t No_of_tracks_to_analyze = 40; | |
20 | ||
21 | const Int_t tbpp = 15; | |
22 | const Int_t nPlanes = 6; | |
23 | const Int_t ntb = tbpp * nPlanes; | |
24 | ||
25 | TH1F *hy = new TH1F("delta R*phi","Cluster displacement in R*phi",200,-1.,1.); | |
26 | TH1F *hyp = new TH1F("delta R*phi pos","delta R*phi, positive",200,-1.,1.); | |
27 | TH1F *hym = new TH1F("delta R*phi neg","delta R*phi, negative",200,-1.,1.); | |
28 | ||
29 | TH1F *hyn = new TH1F("Norm., d(R*phi)","Norm. cluster displacement in R*phi",400,-8.,8.); | |
30 | TH1F *hz = new TH1F("delta Z","Cluster displacement in Z",300,-10.,50.); | |
31 | TH2F *hy2 = new TH2F("Amp vs delta R*phi","Amplitude versus delta R*phi",200,-5.,5.,200,0.,600.); | |
32 | TH2F *herr = new TH2F("sigmaY vs delta R*phi","sigmaY vs delta R*phi",200,-1,1,200,0.,0.1); | |
33 | TH2F *hy3 = new TH2F("Position within pad vs delta R*phi","Position within pad vs delta R*phi",200,-1.,1.,200,-0.5,1.5); | |
34 | TH2F *hy4 = new TH2F("local tb vs delta R*phi","local tb vs delta R*phi",200,-1.,1.,20,-2.5,17.5); | |
35 | ||
36 | hy->SetXTitle("Displacement, cm"); | |
37 | hyn->SetXTitle("Displacement, SigmaY"); | |
38 | hy2->SetXTitle("Displacement, cm"); | |
39 | hy2->SetYTitle("Amplitude"); | |
40 | hz->SetXTitle("Displacement, cm"); | |
41 | hy3->SetXTitle("Displacement, cm"); | |
42 | hy3->SetYTitle("Position, cm"); | |
43 | hy4->SetXTitle("Displacement, cm"); | |
44 | hy4->SetYTitle("local time bin"); | |
45 | ||
46 | /* | |
47 | // Dynamically link some shared libs | |
48 | if (gClassTable->GetID("AliRun") < 0) { | |
49 | gROOT->LoadMacro("loadlibs.C"); | |
50 | loadlibs(); | |
51 | cout << "Loaded shared libraries" << endl; | |
52 | } | |
53 | */ | |
54 | ||
55 | // Load clusters | |
56 | Char_t *clusterfile = "AliTRDclusters.root"; | |
57 | Int_t nEvent = 0; | |
58 | ||
59 | TObjArray carray(2000); | |
60 | TObjArray *ClustersArray = &carray; | |
61 | TFile *geofile =TFile::Open("AliTRDclusters.root"); | |
62 | ||
63 | AliTRDparameter *par = (AliTRDparameter*) geofile->Get("TRDparameter"); | |
64 | AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry"); | |
65 | ||
66 | AliTRDtracker *Tracker = new AliTRDtracker(geofile); | |
67 | Tracker->SetEventNumber(nEvent); | |
68 | Tracker->ReadClusters(ClustersArray,clusterfile); | |
69 | Int_t nClusters = carray.GetEntriesFast(); | |
70 | ||
71 | printf("Total number of clusters %d \n", nClusters); | |
72 | ||
73 | Char_t *alifile = "galice.root"; | |
74 | // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits | |
75 | TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); | |
76 | if (!gafl) { | |
77 | cout << "Open the ALIROOT-file " << alifile << endl; | |
78 | gafl = new TFile(alifile); | |
79 | } | |
80 | else { | |
81 | cout << alifile << " is already open" << endl; | |
82 | } | |
83 | ||
84 | // Get AliRun object from file or create it if not on file | |
85 | gAlice = (AliRun*) gafl->Get("gAlice"); | |
86 | if (gAlice) | |
87 | cout << "AliRun object found on file" << endl; | |
88 | else | |
89 | gAlice = new AliRun("gAlice","Alice test program"); | |
90 | ||
91 | AliTRDv1 *fTRD = (AliTRDv1*) gAlice->GetDetector("TRD"); | |
92 | ||
93 | // Import the Trees for the event nEvent in the file | |
94 | Int_t nparticles = gAlice->GetEvent(nEvent); | |
95 | ||
96 | TObjArray *particles=gAlice->Particles(); | |
97 | ||
98 | TTree *hitTree = gAlice->TreeH(); | |
99 | Int_t nBytes = 0; | |
100 | ||
101 | // Get the number of entries in the hit tree | |
102 | // (Number of primary particles creating a hit somewhere) | |
103 | Int_t nTrack = (Int_t) hitTree->GetEntries(); | |
104 | cout << " Found " << nTrack << " primary particles with hits" << endl; | |
105 | No_of_tracks_to_analyze = TMath::Min(No_of_tracks_to_analyze, nTrack); | |
106 | ||
107 | ||
108 | // Loop through particles and fill histoes | |
109 | ||
110 | Float_t hitY[ntb]; | |
111 | Float_t hitZ[ntb]; | |
112 | Float_t hitO[ntb]; | |
113 | Float_t clusterY[ntb]; | |
114 | Float_t clusterZ[ntb]; | |
115 | Float_t clusterQ[ntb]; | |
116 | Float_t clusterSigmaY[ntb]; | |
117 | Float_t pos[3]; | |
118 | Float_t rot[3]; | |
119 | Float_t global[3]; | |
120 | Int_t det = 0; | |
121 | Int_t track_index[3]; | |
122 | ||
123 | printf("\n"); | |
124 | ||
125 | for (Int_t ii=0; ii<No_of_tracks_to_analyze; ii++) { | |
126 | printf("track %d out of %d \n", ii+1 , No_of_tracks_to_analyze); | |
127 | ||
128 | for(Int_t plane = 0; plane < nPlanes; plane++) { | |
129 | for(Int_t ltb = 14; ltb > -1; ltb--) { | |
130 | ||
131 | if(ii >= nTrack) continue; | |
132 | ||
133 | TParticle *p = gAlice->Particle(ii); | |
134 | if (p->GetFirstMother()>=0) continue; | |
135 | TParticlePDG *pdg = p->GetPDG(); | |
136 | Float_t charge=pdg->Charge(); | |
137 | if(TMath::Abs(charge) < 0.5) continue; | |
138 | ||
139 | Int_t gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
140 | Double_t x = Tracker->GetX(0,plane,ltb); | |
141 | ||
142 | // loop through clusters | |
143 | ||
144 | Bool_t cluster_found = kFALSE; | |
145 | Int_t nw = 0; | |
146 | ||
147 | for (Int_t i = 0; i < nClusters; i++) { | |
148 | ||
149 | AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i); | |
150 | ||
151 | nw = cl->GetDetector(); | |
152 | ||
153 | nw = fGeo->GetPlane(nw); | |
154 | ||
155 | if(nw != plane) continue; | |
156 | for(Int_t j=0; j<3; j++) track_index[j] = cl->GetLabel(j); | |
157 | if((track_index[0] != ii) && | |
158 | (track_index[1] != ii) && | |
159 | (track_index[2] != ii)) continue; | |
160 | ||
161 | nw = cl->GetLocalTimeBin(); | |
162 | if(nw != ltb) continue; | |
163 | ||
164 | clusterY[gtb] = cl->GetY(); | |
165 | clusterZ[gtb] = cl->GetZ(); | |
166 | clusterQ[gtb] = cl->GetQ(); | |
167 | ||
168 | clusterSigmaY[gtb] = TMath::Sqrt(cl->GetSigmaY2()); | |
169 | cluster_found = kTRUE; | |
170 | break; | |
171 | } | |
172 | ||
173 | if(!cluster_found) continue; | |
174 | ||
175 | gAlice->ResetHits(); | |
176 | ||
177 | // nBytes += hitTree->GetEvent(nPrimaries - ii - 1); | |
178 | nBytes += hitTree->GetEvent(nTrack - ii - 1); | |
179 | ||
180 | // Loop through the TRD hits | |
181 | Bool_t found_hit = kFALSE; | |
182 | ||
183 | ||
184 | for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); | |
185 | hit; | |
186 | hit = (AliTRDhit *) fTRD->NextHit()) { | |
187 | nw = hit->Track(); | |
188 | if(nw != ii) continue; | |
189 | det = hit->GetDetector(); | |
190 | nw = fGeo->GetPlane(det); | |
191 | if(nw != plane) continue; | |
192 | ||
193 | ||
194 | pos[0]=hit->X(); | |
195 | pos[1]=hit->Y(); | |
196 | pos[2]=hit->Z(); | |
197 | fGeo->Rotate(det,pos,rot); | |
198 | ||
199 | if(TMath::Abs(rot[0]-x) > 0.01) continue; | |
200 | hitY[gtb] = rot[1]; | |
201 | hitZ[gtb] = rot[2]; | |
202 | ||
203 | Float_t col0 = par->GetCol0(plane); | |
204 | Float_t colPadSize = par->GetColPadSize(plane); | |
205 | Float_t colH = (Int_t ((rot[1] - col0)/colPadSize)) * colPadSize; | |
206 | hitO[gtb] = (rot[1] - col0) - colH; | |
207 | found_hit = kTRUE; | |
208 | break; | |
209 | } | |
210 | ||
211 | if(!found_hit) continue; | |
212 | ||
213 | /* | |
214 | printf("gtb: %d, x: %f, rot[0]: %f, Yhit: %f, Ycl: %f\n", | |
215 | gtb, x, rot[0], rot[1], clusterY[gtb]); | |
216 | printf("\n Zhit - Zcl = %f - %f = %f\n", | |
217 | rot[2], clusterZ[gtb], rot[2] - clusterZ[gtb]); | |
218 | ||
219 | ||
220 | ||
221 | printf("found hit within dx = %f - %f \n",rot[0],x); | |
222 | printf("pos: %f, %f, %f \n",pos[0],pos[1],pos[2]); | |
223 | printf("rot: %f, %f, %f \n",rot[0],rot[1],rot[2]); | |
224 | printf("cluster: %d, %f, %f \n",gtb,clusterY[gtb],clusterZ[gtb]); | |
225 | */ | |
226 | ||
227 | ||
228 | hy->Fill(hitY[gtb]-clusterY[gtb]); | |
229 | if(charge > 0) hyp->Fill(hitY[gtb]-clusterY[gtb]); | |
230 | else if(charge < 0) hym->Fill(hitY[gtb]-clusterY[gtb]); | |
231 | ||
232 | if((clusterQ[gtb]>10)&&(clusterSigmaY[gtb]>0)) | |
233 | hyn->Fill((hitY[gtb]-clusterY[gtb])/clusterSigmaY[gtb]); | |
234 | hz->Fill(hitZ[gtb]-clusterZ[gtb]); | |
235 | hy2->Fill(hitY[gtb]-clusterY[gtb],clusterQ[gtb]); | |
236 | hy3->Fill(hitY[gtb]-clusterY[gtb],hitO[gtb]); | |
237 | hy4->Fill(hitY[gtb]-clusterY[gtb],(Float_t)(tbpp - 1 - gtb%tbpp)); | |
238 | herr->Fill(hitY[gtb]-clusterY[gtb],clusterSigmaY[gtb]); | |
239 | } | |
240 | } | |
241 | } | |
242 | ||
243 | gStyle->SetOptStat(1); | |
244 | gStyle->SetOptFit(1); | |
245 | ||
246 | TCanvas* c = new TCanvas("c", "c", 110, 110, 810, 840); | |
247 | c->SetFillColor(10); | |
248 | c->Divide(2,2); | |
249 | c->cd(1); hy->SetLineWidth(2); hy->SetFillColor(29); hy->Draw(); | |
250 | c->cd(2); hz->SetLineWidth(2); hz->SetFillColor(29); hz->Draw(); | |
251 | c->cd(3); hyn->Draw(); | |
252 | c->cd(4); hy4->Draw(); | |
253 | ||
254 | TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940); | |
255 | c1->SetFillColor(10); | |
256 | c1->Divide(2,2); | |
257 | c1->cd(1); hyp->SetLineWidth(2); hyp->SetFillColor(29); hyp->Fit("gaus"); | |
258 | c1->cd(2); hym->SetLineWidth(2); hym->SetFillColor(29); hym->Fit("gaus"); | |
259 | c1->cd(3); hy3->Draw(); | |
260 | c1->cd(4); hy2->Draw(); | |
261 | ||
262 | } | |
263 | ||
264 | ||
265 | ||
266 | ||
267 | ||
268 |