]>
Commit | Line | Data |
---|---|---|
e8189707 | 1 | #include "iostream.h" |
2 | ||
3 | void SPDclusterTest (Int_t evNumber1=0,Int_t evNumber2=0) | |
4 | { | |
5 | ///////////////////////////////////////////////////////////////////////// | |
6 | // This macro is a small example of a ROOT macro | |
7 | // illustrating how to read the output of GALICE | |
8 | // and do some analysis. | |
9 | // | |
10 | ///////////////////////////////////////////////////////////////////////// | |
11 | ||
12 | // Dynamically link some shared libs | |
13 | ||
14 | if (gClassTable->GetID("AliRun") < 0) { | |
15 | gROOT->LoadMacro("loadlibs.C"); | |
16 | loadlibs(); | |
17 | } | |
18 | ||
19 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
20 | ||
21 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
22 | if (!file) file = new TFile("galice.root"); | |
23 | file->ls(); | |
24 | ||
25 | // Get AliRun object from file or create it if not on file | |
26 | ||
27 | if (!gAlice) { | |
28 | gAlice = (AliRun*)file->Get("gAlice"); | |
29 | if (gAlice) printf("AliRun object found on file\n"); | |
30 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
31 | } | |
32 | ||
33 | // | |
34 | // Loop over events | |
35 | // | |
36 | Int_t Nh=0; | |
37 | Int_t Nh1=0; | |
38 | for (int nev=0; nev<= evNumber2; nev++) { | |
39 | Int_t nparticles = gAlice->GetEvent(nev); | |
40 | cout << "nev " << nev <<endl; | |
41 | cout << "nparticles " << nparticles <<endl; | |
42 | if (nev < evNumber1) continue; | |
43 | if (nparticles <= 0) return; | |
44 | ||
45 | TTree *TH = gAlice->TreeH(); | |
46 | Int_t ntracks = TH->GetEntries(); | |
47 | cout<<"ntracks "<<ntracks<<endl; | |
48 | ||
49 | Int_t nbytes = 0; | |
50 | ||
51 | AliITSRawClusterSPD *ITSclust; | |
52 | ||
53 | // Get pointers to Alice detectors and Digits containers | |
54 | AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); | |
55 | TClonesArray *Particles = gAlice->Particles(); | |
56 | ||
57 | if (ITS) { | |
58 | // fill modules with sorted by module hits | |
59 | Int_t nmodules; | |
60 | ITS->InitModules(-1,nmodules); | |
61 | ITS->FillModules(nev,-1,evNumber2,nmodules," "," "); | |
62 | //get pointer to modules array | |
63 | TObjArray *ITSmodules = ITS->GetModules(); | |
64 | AliITShit *itsHit; | |
65 | ||
66 | // get the Tree for clusters | |
67 | ITS->GetTreeC(nev); | |
68 | TTree *TC=ITS->TreeC(); | |
69 | Int_t nent=TC->GetEntries(); | |
70 | printf("Found %d entries in the tree (must be one per module per event!)\n",nent); | |
71 | ||
72 | for (Int_t idettype=0;idettype<3;idettype++) { | |
73 | ||
74 | TClonesArray *ITSclusters = ITS->ClustersAddress(idettype); | |
75 | //printf ("ITSclusters %p \n",ITSclusters); | |
76 | ||
77 | if (idettype != 0) continue; | |
78 | ||
79 | ||
80 | // ------------ Cluster and point analysis histogramms ------------ | |
81 | ||
82 | TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.); | |
83 | TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.); | |
84 | TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.); | |
85 | TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.); | |
86 | TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.); | |
87 | TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.); | |
88 | TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.); | |
89 | TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.); | |
90 | ||
91 | TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.); | |
92 | TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.); | |
93 | TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.); | |
94 | TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.); | |
95 | ||
96 | ||
97 | // -------------- Create ntuples -------------------- | |
98 | ||
99 | // ntuple structures: | |
100 | ||
101 | struct { | |
102 | Int_t lay; | |
103 | Int_t nx; | |
104 | Int_t nz; | |
105 | Int_t hitprim; | |
106 | Int_t partcode; | |
107 | Float_t dx; | |
108 | Float_t dz; | |
109 | Float_t pmod; | |
110 | } ntuple_st; | |
111 | ||
112 | struct { | |
113 | Int_t lay; | |
114 | Int_t lad; | |
115 | Int_t det; | |
116 | Int_t nx; | |
117 | Int_t nz; | |
118 | Int_t noverlaps; | |
119 | Int_t noverprim; | |
120 | Float_t qcl; | |
121 | Float_t dx; | |
122 | Float_t dz; | |
123 | } ntuple1_st; | |
124 | ||
125 | ||
126 | struct { | |
127 | // Int_t lay; | |
128 | Int_t nx; | |
129 | Int_t nz; | |
130 | } ntuple2_st; | |
131 | ||
132 | ||
133 | ntuple = new TTree("ntuple","Demo ntuple"); | |
134 | ntuple->Branch("lay",&ntuple_st.lay,"lay/I"); | |
135 | ntuple->Branch("nx",&ntuple_st.nx,"nx/I"); | |
136 | ntuple->Branch("nz",&ntuple_st.nz,"nz/I"); | |
137 | ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I"); | |
138 | ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I"); | |
139 | ntuple->Branch("dx",&ntuple_st.dx,"dx/F"); | |
140 | ntuple->Branch("dz",&ntuple_st.dz,"dz/F"); | |
141 | ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F"); | |
142 | ||
143 | ntuple1 = new TTree("ntuple1","Demo ntuple1"); | |
144 | ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I"); | |
145 | ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I"); | |
146 | ntuple1->Branch("det",&ntuple1_st.det,"det/I"); | |
147 | ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I"); | |
148 | ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I"); | |
149 | ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F"); | |
150 | ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I"); | |
151 | ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I"); | |
152 | ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F"); | |
153 | ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F"); | |
154 | ||
155 | ||
156 | ntuple2 = new TTree("ntuple2","Demo ntuple2"); | |
157 | // ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I"); | |
158 | ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I"); | |
159 | ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I"); | |
160 | ||
161 | // ------------------------------------------------------------------------ | |
162 | ||
163 | // Module loop | |
164 | ||
165 | for (Int_t mod=0; mod<nent; mod++) { | |
166 | AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod); | |
167 | ||
168 | Int_t nhits = itsModule->GetNhits(); | |
169 | if(nhits) printf("module nhits %d %d\n",mod,nhits); | |
170 | if(!nhits) continue; | |
171 | ||
172 | ITS->ResetClusters(); | |
173 | TC->GetEvent(mod); | |
174 | Int_t nclust = ITSclusters->GetEntries(); | |
175 | if (nclust) printf("Found %d clust for module %d in det type %d \n",nclust,mod,idettype); | |
176 | if (!nclust) continue; | |
177 | ||
178 | // cluster/hit loops | |
179 | ||
180 | for (Int_t clu=0;clu<nclust;clu++) { | |
181 | itsclu = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu); | |
182 | printf("%d %d %f %f %f\n",itsclu->NclZ(),itsclu->NclX(),itsclu->Q(),itsclu->X(),itsclu->Z()); | |
183 | ||
184 | Int_t noverlaps = 0; | |
185 | Int_t noverprim = 0; | |
186 | ||
187 | Int_t clustersizex = itsclu->NclX(); | |
188 | Int_t clustersizez = itsclu->NclZ(); | |
189 | // Int_t xstart = itsclu->XStart(); | |
190 | // Int_t xstop = itsclu->XStop(); | |
191 | Int_t xstart = itsclu->XStartf(); | |
192 | Int_t xstop = itsclu->XStopf(); | |
193 | Float_t fxstart = xstart*50; | |
194 | Float_t fxstop = (xstop+1)*50; | |
195 | Float_t zstart = itsclu->ZStart(); | |
196 | Float_t zstop = itsclu->ZStop(); | |
197 | Int_t zend = itsclu->Zend(); | |
198 | Float_t clusterx = itsclu->X(); | |
199 | Float_t clusterz = itsclu->Z(); | |
200 | Float_t clusterQ = itsclu->Q(); | |
201 | ||
202 | ntuple2_st.nx = clustersizex; | |
203 | ntuple2_st.nz = clustersizez; | |
204 | ||
205 | ntuple2->Fill(); | |
206 | ||
207 | Int_t icl = 0; | |
208 | Float_t dxprimlast = 10.e+6; | |
209 | Float_t dzprimlast = 10.e+6; | |
210 | ||
211 | ||
212 | // if(module > 217 && module < 226) { | |
213 | cout<<"mod,nclust,clu,Nxpix,Nzpix ="<<mod<<","<<nclust<<","<<clu<<","<<clustersizex<<","<<clustersizez<<endl; | |
214 | cout<<"clusx,clusz ="<<clusterx<<","<<clusterz<<endl; | |
215 | cout<<"XStartf,XStopf,ZStart,ZStop ="<<fxstart<<","<<fxstop<<","<<zstart<<","<<zstop<<endl; | |
216 | // } | |
217 | ||
218 | ||
219 | Float_t SPDlength = 83600; | |
220 | Float_t SPDwidth = 12800; | |
221 | Float_t xhit0 = 1e+5; | |
222 | Float_t zhit0 = 1e+5; | |
223 | ||
224 | ||
225 | for (Int_t hit=0;hit<nhits;hit++) { | |
226 | ||
227 | // Find coordinate differences between the hit and cluster positions | |
228 | // for the resolution determination. | |
229 | ||
230 | itsHit = (AliITShit*)itsModule->GetHit(hit); | |
231 | ||
232 | Int_t hitlayer = itsHit->GetLayer(); | |
233 | Int_t hitladder= itsHit->GetLadder(); | |
234 | Int_t hitdet= itsHit->GetDetector(); | |
235 | ||
236 | Int_t clusterlayer = hitlayer; | |
237 | Int_t clusterladder= hitladder; | |
238 | Int_t clusterdetector = hitdet; | |
239 | ||
240 | Int_t track = itsHit->fTrack; | |
241 | Int_t dray = 0; | |
242 | Int_t hitstat = itsHit->GetTrackStatus(); | |
243 | ||
244 | ||
245 | Float_t zhit = 10000*itsHit->GetZL(); | |
246 | Float_t xhit = 10000*itsHit->GetXL(); | |
247 | ||
248 | if(abs(zhit) > SPDlength/2) { | |
249 | if(hitstat == 66) zhit0 = 1e+5; | |
250 | continue; | |
251 | } | |
252 | ||
253 | if(abs(xhit) > SPDwidth/2) { | |
254 | if(hitstat == 66) xhit0 = 1e+5; | |
255 | continue; | |
256 | } | |
257 | ||
258 | zhit += SPDlength/2; | |
259 | xhit += SPDwidth/2; | |
260 | Float_t yhit = 10000*itsHit->GetYL(); | |
261 | ||
262 | if(hitlayer == 1 && hitstat == 66 && yhit > 71) { | |
263 | xhit0 = xhit; | |
264 | zhit0 = zhit; | |
265 | } | |
266 | if(hitlayer == 2 && hitstat == 66 && yhit < -71) { | |
267 | xhit0 = xhit; | |
268 | zhit0 = zhit; | |
269 | } | |
270 | ||
271 | if(hitstat != 68) continue; // Take only the hit if the last | |
272 | // track point went out from | |
273 | // the detector. | |
274 | ||
275 | if(xhit0 > 9e+4 || zhit0 > 9e+4) continue; | |
276 | ||
277 | ||
278 | Float_t xmed = (xhit + xhit0)/2; | |
279 | Float_t zmed = (zhit + zhit0)/2; | |
280 | ||
281 | Float_t xdif = xmed - clusterx; | |
282 | Float_t zdif = zmed - clusterz; | |
283 | ||
284 | cout<<"clu,hit,xmed,fxstart,fxstop,zmed,zstart,zstop ="<<clu<<","<<hit<<","<<xmed<<","<<fxstart<<","<<fxstop<<","<<zmed<<","<<zstart<<","<<zstop<<endl; | |
285 | ||
286 | // Consider the hits inside of cluster region only | |
287 | ||
288 | if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) { | |
289 | ||
290 | icl = 1; | |
291 | ||
292 | // part = (TParticle *)particles.UncheckedAt(track); | |
293 | // Int_t partcode = part->GetPdgCode(); | |
294 | // Int_t primery = gAlice->GetPrimary(track); | |
295 | ||
296 | Int_t parent = itsHit->GetParticle()->GetFirstMother(); | |
297 | Int_t partcode = itsHit->GetParticle()->GetPdgCode(); | |
298 | ||
299 | // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+ | |
300 | // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda | |
301 | ||
302 | ||
303 | Float_t px = itsHit->GetPXL(); | |
304 | Float_t py = itsHit->GetPYL(); | |
305 | Float_t pz = itsHit->GetPZL(); | |
306 | Float_t pmod = 1000*sqrt(px*px+py*py+pz*pz); | |
307 | ||
308 | cout<<"track,partcode,pmod,parent ="<<track<<","<<partcode<<","<<pmod<<","<<parent<<endl; | |
309 | ||
310 | Int_t hitprim = 0; | |
311 | ||
312 | if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e- | |
313 | // at p < 6 MeV/c | |
314 | ||
315 | if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but | |
316 | // not for delta ray which | |
317 | // also went out from the | |
318 | // detector and returned | |
319 | // again | |
320 | ||
321 | if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery | |
322 | // particles | |
323 | ||
324 | if(hitprim > 0) noverprim = noverprim + 1; | |
325 | ||
326 | if(hitprim > 0) { | |
327 | dxprimlast = xdif; | |
328 | dzprimlast = zdif; | |
329 | } | |
330 | ||
331 | // fill ntuple | |
332 | ||
333 | ntuple_st.lay = hitlayer; | |
334 | ntuple_st.nx = clustersizex; | |
335 | ntuple_st.nz = clustersizez; | |
336 | ntuple_st.hitprim = hitprim; | |
337 | ntuple_st.partcode = partcode; | |
338 | ntuple_st.dx = xdif; | |
339 | ntuple_st.dz = zdif; | |
340 | ntuple_st.pmod = pmod; | |
341 | ||
342 | ntuple->Fill(); | |
343 | ||
344 | ||
345 | ||
346 | if(hitprim > 0) { // for primary particles | |
347 | if(hitlayer == 1) { | |
348 | cout<<"!!!!!! lay,hitprim,xdif,zdif ="<<hitlayer<<","<<hitprim<<","<<xdif<<","<<zdif<<endl; | |
349 | Xres1->Fill(xdif); | |
350 | Zres1->Fill(zdif); | |
351 | } | |
352 | if(hitlayer == 2) { | |
353 | cout<<"!!!!!! lay,hitprim,xdif,zdif ="<<hitlayer<<","<<hitprim<<","<<xdif<<","<<zdif<<endl; | |
354 | Xres2->Fill(xdif); | |
355 | Zres2->Fill(zdif); | |
356 | } | |
357 | } // primery particles | |
358 | ||
359 | } // end of cluster region | |
360 | } // end of hit loop | |
361 | ||
362 | ||
363 | if(icl == 1) { | |
364 | ||
365 | // fill ntuple1 | |
366 | ||
367 | // ntuple1->Fill(clusterlayer,clustersizex,clustersizez,noverlaps,\ | |
368 | noverprim,dx,dz); | |
369 | ||
370 | if(noverlaps == 0) noverlaps = 1; // cluster contains one or more | |
371 | // delta rays only | |
372 | ||
373 | ntuple1_st.lay = clusterlayer; | |
374 | ntuple1_st.lad = clusterladder; | |
375 | ntuple1_st.det = clusterdetector; | |
376 | ntuple1_st.nx = clustersizex; | |
377 | ntuple1_st.nz = clustersizez; | |
378 | ntuple1_st.qcl = clusterQ; | |
379 | ntuple1_st.noverlaps = noverlaps; | |
380 | ntuple1_st.noverprim = noverprim; | |
381 | ntuple1_st.dx = dxprimlast; | |
382 | ntuple1_st.dz = dzprimlast; | |
383 | ||
384 | ntuple1->Fill(); | |
385 | ||
386 | } // icl = 1 | |
387 | } // cluster loop | |
388 | } // module loop | |
389 | } // idettype loop | |
390 | } // end if ITS | |
391 | } // event loop | |
392 | ||
393 | ||
394 | // Write and Draw Histogramms and ntuples | |
395 | ||
396 | ||
397 | ||
398 | TFile fhistos("SPD_his.root","RECREATE"); | |
399 | ||
400 | ntuple->Write(); | |
401 | ntuple1->Write(); | |
402 | ntuple2->Write(); | |
403 | ||
404 | Nxpix1->Write(); | |
405 | Nzpix1->Write(); | |
406 | Nxpix2->Write(); | |
407 | Nzpix2->Write(); | |
408 | ||
409 | Xpix1->Write(); | |
410 | Zpix1->Write(); | |
411 | Xpix2->Write(); | |
412 | Zpix2->Write(); | |
413 | ||
414 | Xres1->Write(); | |
415 | Zres1->Write(); | |
416 | Xres2->Write(); | |
417 | Zres2->Write(); | |
418 | ||
419 | fhistos.Close(); | |
420 | cout<<"!!! Histogramms and ntuples were written"<<endl; | |
421 | ||
422 | TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700); | |
423 | c1->Divide(2,2); | |
424 | c1->cd(1); | |
425 | gPad->SetFillColor(33); | |
426 | Xres1->SetFillColor(42); | |
427 | Xres1->Draw(); | |
428 | c1->cd(2); | |
429 | gPad->SetFillColor(33); | |
430 | Zres1->SetFillColor(46); | |
431 | Zres1->Draw(); | |
432 | c1->cd(3); | |
433 | gPad->SetFillColor(33); | |
434 | Xres2->SetFillColor(42); | |
435 | Xres2->Draw(); | |
436 | c1->cd(4); | |
437 | gPad->SetFillColor(33); | |
438 | Zres2->SetFillColor(46); | |
439 | Zres2->Draw(); | |
440 | ||
441 | // cout<<"END test for clusters and hits "<<endl; | |
442 | ||
443 | // file->Close(); | |
444 | } | |
445 | ||
446 | ||
447 |