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