]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | void MUONtestnew (Int_t evNumber1=0,Int_t evNumber2=0) |
2 | { | |
3 | ///////////////////////////////////////////////////////////////////////// | |
4 | // This macro is a small example of a ROOT macro | |
5 | // illustrating how to read the output of GALICE | |
6 | // and do some analysis. | |
7 | // | |
8 | ///////////////////////////////////////////////////////////////////////// | |
9 | ||
10 | // Dynamically link some shared libs | |
11 | ||
12 | if (gClassTable->GetID("AliRun") < 0) { | |
13 | gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so"); | |
14 | gSystem->Load("$ALITOP/cern.so/lib/libPythia.so"); | |
15 | gSystem->Load("$ROOTSYS/lib/libEG.so"); | |
16 | gSystem->Load("$ROOTSYS/lib/libEGPythia.so"); | |
17 | gSystem->Load("libGeant3Dummy.so"); //a dummy version of Geant3 | |
18 | gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes | |
19 | gSystem->Load("libgalice.so"); // the standard Alice classes | |
20 | } | |
21 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
22 | ||
23 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
24 | ||
25 | ||
26 | if (!file) { | |
27 | printf("\n Creating galice.root \n"); | |
28 | file = new TFile("galice.root"); | |
29 | } else { | |
30 | printf("\n galice.root found in file list"); | |
31 | } | |
32 | file->ls(); | |
33 | ||
34 | // Get AliRun object from file or create it if not on file | |
35 | if (!gAlice) { | |
36 | gAlice = (AliRun*)(file->Get("gAlice")); | |
37 | if (gAlice) printf("AliRun object found on file\n"); | |
38 | if (!gAlice) { | |
39 | printf("\n create new gAlice object"); | |
40 | gAlice = new AliRun("gAlice","Alice test program"); | |
41 | } | |
42 | } | |
43 | ||
44 | // Create some histograms | |
45 | ||
46 | ||
47 | TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution" | |
48 | ,100,0.,500.); | |
49 | TH1F *pcharge1 = new TH1F("pcharge1","Pad Charge distribution" | |
50 | ,100,0.,200.); | |
51 | TH1F *xresid1 = new TH1F("xresid1","x-Residuals" | |
52 | ,100,-0.5,0.5); | |
53 | TH1F *yresid1 = new TH1F("yresid1","y-Residuals" | |
54 | ,100,-0.2,0.2); | |
55 | TH1F *npads1 = new TH1F("npads1" ,"Pads per Hit" | |
56 | ,20,-0.5,19.5); | |
57 | TH2F *xresys1 = new TH2F("xresys1","x-Residuals systematics" | |
58 | ,50,-2.0,2.0,100,-1.0,1.0); | |
59 | TH2F *yresys1 = new TH2F("yresys1","y-Residuals systematics" | |
60 | ,50,-0.4,0.4,100,-0.1,0.1); | |
61 | ||
62 | TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution" | |
63 | ,100,0.,500.); | |
64 | TH1F *pcharge2 = new TH1F("pcharge2","Pad Charge distribution" | |
65 | ,100,0.,200.); | |
66 | TH1F *xresid2 = new TH1F("xresid2","x-Residuals" | |
67 | ,100,-0.5,0.5); | |
68 | TH1F *yresid2 = new TH1F("yresid2","y-Residuals" | |
69 | ,100,-0.2,0.2); | |
70 | TH1F *npads2 = new TH1F("npads2" ,"Pads per Hit" | |
71 | ,20,-0.5,19.5); | |
72 | TH2F *xresys2 = new TH2F("xresys2","x-Residuals systematics" | |
73 | ,50,-2.0,2.0,100,-1.0,1.0); | |
74 | TH2F *yresys2 = new TH2F("yresys2","y-Residuals systematics" | |
75 | ,50,-0.4,0.4,100,-0.1,0.1); | |
76 | ||
77 | ||
78 | AliMUONchamber* iChamber; | |
79 | // | |
80 | // Loop over events | |
81 | // | |
82 | Int_t Nh=0; | |
83 | Int_t Nh1=0; | |
84 | for (Int_t nev=0; nev<= evNumber2; nev++) { | |
85 | cout << "nev " << nev <<endl; | |
86 | Int_t nparticles = gAlice->GetEvent(nev); | |
87 | cout << "nparticles " << nparticles <<endl; | |
88 | if (nev < evNumber1) continue; | |
89 | if (nparticles <= 0) return; | |
90 | ||
91 | AliMUON *MUON = (AliMUON*) gAlice->GetDetector("MUON"); | |
92 | printf("\n track %d %d \n ", nev, MUON); | |
93 | ||
94 | TTree *TH = gAlice->TreeH(); | |
95 | Int_t ntracks = TH->GetEntries(); | |
96 | Int_t Nc=0; | |
97 | Int_t npad[2]; | |
98 | Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2]; | |
99 | // | |
100 | // Loop over events | |
101 | // | |
102 | for (Int_t track=0; track<ntracks;track++) { | |
103 | ||
104 | gAlice->ResetHits(); | |
105 | Int_t nbytes += TH->GetEvent(track); | |
106 | if (MUON) { | |
107 | for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); | |
108 | mHit; | |
109 | mHit=(AliMUONhit*)MUON->NextHit()) | |
110 | { | |
111 | Int_t nch = mHit->fChamber; // chamber number | |
112 | Float_t x = mHit->fX; // x-pos of hit | |
113 | Float_t y = mHit->fY; // y-pos | |
114 | // | |
115 | // | |
116 | iChamber = & MUON->Chamber(nch-1); | |
117 | response=iChamber->GetResponseModel(); | |
118 | // | |
119 | // | |
120 | if (nch <= 1) { | |
121 | for (Int_t i=0; i<2; i++) { | |
122 | xmean[i]=0; | |
123 | ymean[i]=0; | |
124 | Q[i]=0; | |
125 | npad[i]=0; | |
126 | } | |
127 | for (AliMUONcluster *mPad=(AliMUONcluster*)MUON | |
128 | ->FirstPad(mHit); | |
129 | mPad; | |
130 | mPad=(AliMUONcluster*)MUON->NextPad() | |
131 | ) | |
132 | { | |
133 | Int_t nseg = mPad->fCathode; // segmentation | |
134 | Int_t nhit = mPad->fHitNumber; // hit number | |
135 | Int_t qtot = mPad->fQ; // charge | |
136 | Int_t ipx = mPad->fPadX; // pad number on X | |
137 | Int_t ipy = mPad->fPadY; // pad number on Y | |
138 | Int_t iqpad = mPad->fQpad; // charge per pad | |
139 | Int_t secpad = mPad->fRSec; // r-pos of pad | |
140 | // | |
141 | // | |
142 | segmentation=iChamber->GetSegmentationModel(nseg); | |
143 | Int_t ind=nseg-1; | |
144 | // printf("Pad of hit, padx, pady, iqpad, ncha %d %d %d %d %d %d\n", | |
145 | // nhit, ipx, ipy, iqpad, nseg, ind); | |
146 | ||
147 | ||
148 | if (nseg == 1) { | |
149 | pcharge1->Fill(iqpad,(float) 1); | |
150 | ccharge1->Fill(qtot ,(float) 1); | |
151 | } else { | |
152 | pcharge2->Fill(iqpad,(float) 1); | |
153 | ccharge2->Fill(qtot ,(float) 1); | |
154 | } | |
155 | // Calculate centre of gravity | |
156 | // | |
157 | if (iqpad == 0 && ind==0) { | |
158 | printf("\n attention iqpad 0"); | |
159 | } | |
160 | ||
161 | if (iqpad > 0 && secpad==1) { | |
162 | Float_t xc,yc; | |
163 | npad[ind]++; | |
164 | segmentation->GetPadCxy(ipx,ipy,xc,yc); | |
165 | // printf("\n pad %d %d ", ipx, ipy); | |
166 | // printf("\n pos %f %f ", xc, yc); | |
167 | xmean[ind]+=Float_t(iqpad*xc); | |
168 | ymean[ind]+=Float_t(iqpad*yc); | |
169 | Q[ind]+=iqpad; | |
170 | } | |
171 | ||
172 | } //pad hit loop | |
173 | for (Int_t i=0; i<2; i++) { | |
174 | segmentation = iChamber->GetSegmentationModel(i+1); | |
175 | ||
176 | // if (npad[i] ==0) { | |
177 | // printf("\n #pads=0 %f %f",x,y); | |
178 | // } | |
179 | ||
180 | if (Q[i] >0) { | |
181 | xmean[i] = xmean[i]/Q[i]; | |
182 | xres[i] = xmean[i]-x; | |
183 | ymean[i] = ymean[i]/Q[i]; | |
184 | yres[i] = ymean[i]-y; | |
185 | // printf("\n xmean, ymean %f %f",xmean[i],ymean[i]); | |
186 | // printf("\n xhit, yhit %f %f",x,y); | |
187 | ||
188 | // Systematic Error | |
189 | // | |
190 | Int_t icx, icy; | |
191 | segmentation-> | |
192 | GetPadIxy(x,y,icx,icy); | |
193 | segmentation-> | |
194 | GetPadCxy(icx,icy,xonpad[i],yonpad[i]); | |
195 | xonpad[i]-=x; | |
196 | yonpad[i]-=y; | |
197 | } | |
198 | } // plane loop | |
199 | xresid1->Fill(xres[0],(float) 1); | |
200 | yresid1->Fill(yres[0],(float) 1); | |
201 | npads1->Fill(npad[0],(float) 1); | |
202 | if (npad[0] >=2 && Q[0] > 6 ) { | |
203 | xresys1->Fill(xonpad[0],xres[0],(float) 1); | |
204 | yresys1->Fill(yonpad[0],yres[0],(float) 1); | |
205 | } | |
206 | ||
207 | xresid2->Fill(xres[1],(float) 1); | |
208 | yresid2->Fill(yres[1],(float) 1); | |
209 | npads2->Fill(npad[1],(float) 1); | |
210 | if (npad[1] >=2 && Q[1] > 6) { | |
211 | xresys2->Fill(xonpad[1],xres[1],(float) 1); | |
212 | yresys2->Fill(yonpad[1],yres[1],(float) 1); | |
213 | } | |
214 | } // chamber 1 | |
215 | } // hit loop | |
216 | } // if MUON | |
217 | } // track loop | |
218 | } // event loop | |
219 | ||
220 | //Create a canvas, set the view range, show histograms | |
221 | TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700); | |
222 | pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99); | |
223 | pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99); | |
224 | pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49); | |
225 | pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49); | |
226 | pad11->SetFillColor(11); | |
227 | pad12->SetFillColor(11); | |
228 | pad13->SetFillColor(11); | |
229 | pad14->SetFillColor(11); | |
230 | pad11->Draw(); | |
231 | pad12->Draw(); | |
232 | pad13->Draw(); | |
233 | pad14->Draw(); | |
234 | ||
235 | pad11->cd(); | |
236 | ccharge1->SetFillColor(42); | |
237 | ccharge1->SetXTitle("ADC units"); | |
238 | ccharge1->Draw(); | |
239 | ||
240 | pad12->cd(); | |
241 | pcharge1->SetFillColor(42); | |
242 | pcharge1->SetXTitle("ADC units"); | |
243 | pcharge1->Draw(); | |
244 | ||
245 | pad13->cd(); | |
246 | xresid1->SetFillColor(42); | |
247 | xresid1->Draw(); | |
248 | ||
249 | pad14->cd(); | |
250 | yresid1->SetFillColor(42); | |
251 | yresid1->Draw(); | |
252 | ||
253 | TCanvas *c2 = new TCanvas("c2","Cluster Size",400,10,600,700); | |
254 | pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99); | |
255 | pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99); | |
256 | pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49); | |
257 | pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49); | |
258 | pad21->SetFillColor(11); | |
259 | pad22->SetFillColor(11); | |
260 | pad23->SetFillColor(11); | |
261 | pad24->SetFillColor(11); | |
262 | pad21->Draw(); | |
263 | pad22->Draw(); | |
264 | pad23->Draw(); | |
265 | pad24->Draw(); | |
266 | ||
267 | pad21->cd(); | |
268 | npads1->SetFillColor(42); | |
269 | npads1->SetXTitle("Cluster Size"); | |
270 | npads1->Draw(); | |
271 | ||
272 | pad23->cd(); | |
273 | xresys1->SetXTitle("x on pad"); | |
274 | xresys1->SetYTitle("x-xcog"); | |
275 | xresys1->Draw(); | |
276 | ||
277 | pad24->cd(); | |
278 | yresys1->SetXTitle("y on pad"); | |
279 | yresys1->SetYTitle("y-ycog"); | |
280 | yresys1->Draw(); | |
281 | ||
282 | TCanvas *c3 = new TCanvas("c3","Charge and Residuals",400,10,600,700); | |
283 | pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99); | |
284 | pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99); | |
285 | pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49); | |
286 | pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49); | |
287 | pad31->SetFillColor(11); | |
288 | pad32->SetFillColor(11); | |
289 | pad33->SetFillColor(11); | |
290 | pad34->SetFillColor(11); | |
291 | pad31->Draw(); | |
292 | pad32->Draw(); | |
293 | pad33->Draw(); | |
294 | pad34->Draw(); | |
295 | ||
296 | pad31->cd(); | |
297 | ccharge2->SetFillColor(42); | |
298 | ccharge2->SetXTitle("ADC units"); | |
299 | ccharge2->Draw(); | |
300 | ||
301 | pad32->cd(); | |
302 | pcharge2->SetFillColor(42); | |
303 | pcharge2->SetXTitle("ADC units"); | |
304 | pcharge2->Draw(); | |
305 | ||
306 | pad33->cd(); | |
307 | xresid2->SetFillColor(42); | |
308 | xresid2->Draw(); | |
309 | ||
310 | pad34->cd(); | |
311 | yresid2->SetFillColor(42); | |
312 | yresid2->Draw(); | |
313 | ||
314 | TCanvas *c4 = new TCanvas("c4","Cluster Size",400,10,600,700); | |
315 | pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99); | |
316 | pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99); | |
317 | pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49); | |
318 | pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49); | |
319 | pad41->SetFillColor(11); | |
320 | pad42->SetFillColor(11); | |
321 | pad43->SetFillColor(11); | |
322 | pad44->SetFillColor(11); | |
323 | pad41->Draw(); | |
324 | pad42->Draw(); | |
325 | pad43->Draw(); | |
326 | pad44->Draw(); | |
327 | ||
328 | pad41->cd(); | |
329 | npads2->SetFillColor(42); | |
330 | npads2->SetXTitle("Cluster Size"); | |
331 | npads2->Draw(); | |
332 | ||
333 | pad43->cd(); | |
334 | xresys2->SetXTitle("x on pad"); | |
335 | xresys2->SetYTitle("x-xcog"); | |
336 | xresys2->Draw(); | |
337 | ||
338 | pad44->cd(); | |
339 | yresys2->SetXTitle("y on pad"); | |
340 | yresys2->SetYTitle("y-ycog"); | |
341 | yresys2->Draw(); | |
342 | ||
343 | } |